Gerya T. Introduction to Numerical Geodynamic Modelling
Подождите немного. Документ загружается.
1.4 Derivation of the Eulerian continuity equation 17
The balance between the old (m
0
) and new (m
1
) mass results in the following
relation
m
1
= m
0
+ m
in
− m
out
,
m
in
= m
A
+ m
C
+ m
E
,
m
out
= m
B
+ m
D
+ m
F
,
m
A
= ρ
A
v
xA
yzt,
m
B
= ρ
B
v
xB
yzt, (1.7)
m
C
= ρ
C
v
yC
xzt,
m
D
= ρ
D
v
yD
xzt,
m
E
= ρ
E
v
zE
xyt,
m
F
= ρ
F
v
zF
xyt,
where m
in
and m
out
are the incoming and outgoing mass, respectively; m
A
– m
F
are
masses that passed through the respective boundaries during the time t; ρ
A
– ρ
F
is the density at the respective boundaries; v
x A
– v
zD
are the velocity components
responsible for material fluxes through the boundaries (Fig. 1.6(a)). If t is small,
we can now write an approximate expression for the Eulerian time derivative of
the average density in the volume as:
∂ρ
∂t
≈
ρ
t
=
ρ
1
− ρ
0
t
=
m
1
− m
0
x y z t
. (1.8)
By using Equation (1.7) the following expression can be further obtained (verify
as an exercise)
ρ
t
=−
ρ
B
v
xB
− ρ
A
v
xA
x
−
ρ
D
v
yD
− ρ
C
v
yC
y
−
ρ
F
v
zF
− ρ
E
v
zE
z
, (1.8a)
or
ρ
t
=−
(ρv
x
)
x
−
(ρv
y
)
y
−
(ρv
z
)
z
, (1.8b)
or
ρ
t
+
(ρv
x
)
x
+
(ρv
y
)
y
+
(ρv
z
)
z
= 0,
(ρv
x
) = ρ
B
v
xB
− ρ
A
v
xA
,
(1.8c)
(ρv
y
) = ρ
D
v
yD
− ρ
C
v
yC
,
(ρv
z
) = ρ
F
v
zF
− ρ
E
v
zE
,
where (ρv
x
), (ρv
y
)and(ρv
z
) are differences in the mass fluxes in x, y and z
directions respectively (i.e. through the respective pairs of boundaries, Fig. 1.6(a)).
18 The continuity equation
Obviously, in such cases when t, x, y and z all tend to zero, the differences
can be replaced by derivatives and we obtain the Eulerian continuity equation
∂ρ
∂t
+
∂(ρv
x
)
∂x
+
∂(ρv
y
)
∂y
+
∂(ρv
z
)
∂z
= 0, (1.9a)
or
∂ρ
∂t
+ div(ρ v) = 0. (1.9b)
1.5 Derivation of the Lagrangian continuity equation
Similarly, the Lagrangian continuity equation can be derived by analysing the
motion of a small, mobile, rectangular Lagrangian (material) volume of initial
dimensions x
0
, y
0
and z
0
(Fig. 1.6(b)). In contrast to the steady Eulerian
volume, the amount of mass m in the moving Lagrangian volume remains constant
(since this volume is always related to the same material points), but the dimensions
of the volume may change due to internal expansion/contraction processes. The
initial average fluid density (ρ
0
), within the volume is given by
ρ
0
=
m
x
0
y
0
z
0
. (1.10)
Internal expansion or contraction affects the dimensions of the Lagrangian volume,
and after a small period of time t, these dimensions become x
1
, y
1
and z
1
and the average fluid density (ρ
1
) changes to
ρ
1
=
m
x
1
y
1
z
1
. (1.11)
The relationship between the old (x
0
, y
0
, z
0
)andthenew(x
1
, y
1
, z
1
)
dimensions of the Lagrangian volume can be established on the basis of the relative
movements of the boundaries of the volume (A, B, C, D, E, F, Fig. 1.6(b))which
leads to the following relations
x
1
= x
0
+ t v
x
, (1.12)
v
x
= v
xB
− v
xA
y
1
= y
0
+ t v
y
, (1.13)
v
y
= v
yD
− v
yC
z
1
= z
0
+ t v
z
, (1.14)
v
z
= v
zF
− v
zE
where v
x
, v
y
and v
z
are the differences in the velocity components that corre-
spond to the relative movements of their respective pairs of boundaries (Fig. 1.6(b)).
Taking t to be small, we can now write an approximate expression for the
1.5 Derivation of the Lagrangian continuity equation 19
Lagrangian time derivative of the average density in the volume as
Dρ
Dt
≈
ρ
t
=
ρ
1
− ρ
0
t
=
m
x
1
y
1
z
1
t
−
m
x
0
y
0
z
0
t
. (1.15)
By using the equations derived for x
1
, y
1
and z
1
(Eqs. (1.12)–(1.14))the
following expression can be obtained (verify as an exercise)
ρ
t
= ρ
0
1 −
1 + t
v
x
x
0
1 + t
v
y
y
0
1 + t
v
z
z
0
t
1 + t
v
x
x
0
1 + t
v
y
y
0
1 + t
v
z
z
0
, (1.16a)
or
ρ
t
= ρ
0
−
v
x
x
0
−
v
y
y
0
−
v
z
z
0
− t
v
x
x
0
v
y
y
0
+
v
x
x
0
v
z
z
0
+
v
y
y
0
v
z
z
0
+ t
v
x
x
0
v
y
y
0
v
z
z
0
1 +t
v
x
x
0
1 +t
v
y
y
0
1 + t
v
z
z
0
,
(1.16b)
or
ρ
t
+ ρ
0
v
x
x
0
+
v
y
y
0
+
v
z
z
0
+ K
1
K
2
= 0,
K
1
= t
v
x
x
0
v
y
y
0
+
v
x
x
0
v
z
z
0
+
v
y
y
0
v
z
z
0
+ t
v
x
x
0
v
y
y
0
v
z
z
0
, (1.16c)
K
2
=
1 + t
v
x
x
0
1 + t
v
y
y
0
1 + t
v
z
z
0
,
where K
1
and K
2
are coefficients which respectively tend to zero and unity when t
tends to zero. Obviously in the case when t, x
0
, y
0
and z
0
all tend towards
zero, the differences in Eq. (1.16c) can be replaced by derivatives and taking K
1
=0
and K
2
=1 we obtain the Lagrangian continuity equation
Dρ
Dt
+ ρ
∂v
x
∂x
+ ρ
∂v
y
∂y
+ ρ
∂v
z
∂z
= 0, (1.17a)
or
Dρ
Dt
+ ρdiv(v) = 0. (1.17b)
20 The continuity equation
1.6 Comparing Eulerian and Lagrangian continuity equations.
Advective transport term
Let’s now perform transformations of the Eulerian continuity equation (Eq. (1.1))
in order to decipher its structure and to establish a relationship with the Lagrangian
continuity equation (Eq. (1.3)). It is convenient to decompose div(ρ v)usingthe
standard product rule (also called Leibniz’s law)
(
u · v
)
= u
· v +v
· u,orinan
alternative notation
∂
∂x
(
uv
)
=
∂u
∂x
· v +
∂v
∂x
· u.
div(ρ v) = ρdiv(v) +v grad(ρ), (1.18a)
or in a different symbolic notation
∇·(ρ v) = ρ∇·v +v∇ρ, (1.18b)
or ‘deciphering’ what we actually are doing in three dimensions
∂
∂x
(
ρv
x
)
+
∂
∂y
(ρv
y
) +
∂
∂z
(
ρv
z
)
=
ρ
∂v
x
∂x
+ ρ
∂v
y
∂y
+ ρ
∂v
z
∂z
+
v
x
∂ρ
∂x
+ v
y
∂ρ
∂y
+ v
z
∂ρ
∂z
, (1.18c)
where grad(ρ)or∇ρ is the gradient of the density ρ. The gradient is a vector
function of a scalar field defined as follows
in one dimension (1D): grad(ρ) =
∂ρ
∂x
, (1.19a)
in two dimensions (2D): grad(ρ) =
∂ρ
∂x
,
∂ρ
∂y
, (1.19b)
in three dimensions (3D): grad(ρ) =
∂ρ
∂x
,
∂ρ
∂y
,
∂ρ
∂z
. (1.19c)
Therefore, both ∇ρ and v in Eq. (1.18) are vectors and scalar multiplying of these
two vectors (1.18c) gives the following result
in one dimension (1D): v grad(ρ) = v
x
∂ρ
∂x
, (1.20a)
in two dimensions (2D): v grad(ρ) = v
x
∂ρ
∂x
+ v
y
∂ρ
∂y
, (1.20b)
in three dimensions (3D): v grad(ρ) = v
x
∂ρ
∂x
+ v
y
∂ρ
∂y
+ v
z
∂ρ
∂z
. (1.20c)
Now by comparing Equations (1.1), (1.3) and (1.18) we can establish the relation-
ship between the Eulerian (
∂ρ
∂t
) and Lagrangian (
Dρ
Dt
) time derivatives of density
1.6 Comparing Eulerian and Lagrangian equations 21
Fig. 1.7 Schematic representation of the advective transport in the case of uniform
1D movement of a fluid with a linear density distribution. The dashed and solid
thick lines correspond to the density distribution for the moments of time t
0
and
t
1
, respectively. Circles A and B denote the density for two different Lagrangian
points passing through an Eulerian observation point (solid rectangle C) at the
different moments of time (t
0
and t
1
, respectively).
as
Dρ
Dt
=
∂ρ
∂t
+v grad(ρ). (1.21)
The extra term, v grad(ρ) in the Eulerian continuity equation is an advective trans-
port term that reflects changes of density in an immobile (Eulerian) point, due
to the movement of an inhomogeneous medium with existing density gradients
relative to this point (Fig. 1.7). Obviously, the advective transport terms in the
Eulerian continuity equation are only relevant (i.e. nonzero) in situations when
both the velocity of the medium and density gradient are nonzero. On the other
hand, substantive changes of density (Dρ/Dt) in the moving Lagrangian point do
not depend on density gradients and the Lagrangian continuity equation (1.3) thus
does not contain advective terms.
When the density in all moving material points does not change with time
(i.e.
Dρ
Dt
= 0) the Eulerian continuity equation reduces to the Eulerian advective
transport equation
∂ρ
∂t
=−v grad(ρ). (1.22)
The minus sign in the right-hand side of equation (1.22) reflects the relation between
the density gradient and the direction of motion (Fig. 1.7): if a medium is moving
in the direction of decreasing density (i.e. v grad(ρ) < 0), then the density in an
immobile observation point increases with time (i.e.,
∂ρ
∂t
> 0).
Let’s now derive the advective transport relation (Eq. (1.22)) for the simple 1D
case shown in Fig. 1.7. A fluid with a linear density distribution moves at a constant
velocity v
x
. At the moment of time t
0
the density (ρ
A
) in the Eulerian observation
22 The continuity equation
point C (solid rectangle) corresponds to the Lagrangian point A (solid circle). Due
to fluid motion, the density profile shifts to the right with time. Therefore, at a later
moment of time t
1
, the density (ρ
B
) at the same Eulerian point C will correspond to
another Lagrangian point B (open circle). Under the assumption that the difference
in time between the two moments (t =t
1
−t
0
) is small, we can approximate the
time derivative of density at the Eulerian point C as
∂ρ
∂t
≈
ρ
t
t
=
ρ
B
− ρ
A
t
1
− t
0
, (1.23)
where ρ
t
corresponds to changes in density with time at the Eulerian point C. If
on the other hand, the displacement of the density profile (i.e. the distance between
two Lagrangian points A and B, x =x
1
−x
0
) is also small, we can approximate
the local derivative of density with respect to x, at the moment of time t
1
as
∂ρ
∂x
≈
ρ
x
x
=
ρ
A
− ρ
B
x
1
− x
0
, (1.24)
where ρ
x
corresponds to a change in density with coordinate in the proximity of
the Eulerian point C at the instant when time equals t
1
. Note that ρ
t
and ρ
x
are different quantities, and that in the case considered they are similar in absolute
value but are different in sign (cf. Eq. (1.23) and (1.24)). Taking into account the
relationships
x = v
x
t
and
ρ
x
=−ρ
t
,
the following expression can be obtained (derive as an exercise)
ρ
t
t
=−v
x
ρ
x
x
. (1.25)
When t and respectively x tend to zero, the differences in Eq. (1.25) can be
replaced by derivatives and we obtain the 1D advection equation
∂ρ
∂t
=−v
x
∂ρ
∂x
. (1.26)
The advective transport equation is very important for geodynamic modelling since
it describes changes in the distribution of transport properties (such as density,
temperature, composition, etc.) in non-homogeneous, deforming media. We will
come back to this issue in Chapter 8.
Analytical exercise 23
1.7 Incompressible continuity equation
For many geological media (such as the Earth’s crust and mantle), one may assume
an incompressibility condition (i.e. density of material points does not change with
time), which is written as:
Dρ
Dt
=
∂ρ
∂t
+v grad (ρ) = 0. (1.27)
It is valid in cases when pressure and temperature changes are not very large and
no phase transformations leading to volume changes occur in the medium. In this
situation, one can use an incompressible continuity equation which is the same in
both Eulerian and Lagrangian forms
div(v) = 0. (1.28)
The incompressible continuity equation is broadly used in numerical geodynamic
modelling, although in many cases it is a rather big simplification (e.g. in the case of
the whole Earth mantle convection, e.g., Tackley, 2008). Typical examples of geo-
dynamic settings where deformations are strongly defined by the incompressibility
condition are, for example, corner flow in the mantle wedge above subducting
slabs (e.g. Turcotte and Schubert, 2002) and circulation of tectonic melanges in
subduction channels (e.g. Cloos, 1982).
Analytical exercise
Exercise 1.1
In a region of the Earth’s mantle, the velocity field is given by
v
x
= 10
−10
+ x · 10
−13
+ y · 10
−13
+ z · 10
−13
v
y
= 10
−10
− x · 10
−13
+ y · 2 × 10
−13
+ z · 3 ×10
−13
v
z
= 10
−10
− x · 10
−13
− y · 10
−13
− z · 2 ×10
−13
.
The mantle density field in the same region is given by
ρ = 3300 + x · 0.001 − y · 0.002 + z · 0.001.
Calculate ρ,div(v),
∂ρ
∂t
and
Dρ
Dt
for the point with coordinates x =1000, y =1000,
z =1000.
24 The continuity equation
Programming exercise and homework
Exercise 1.2
Write a MATLAB code for computing and visualising a 2D velocity field and its
divergence. Model design: an area of the mantle (1000 ×1500 km) is convecting
with one central upwelling in the middle of the model box and two downwellings
at the sides. The velocity field is given by the following equations
v
x
=−v
x0
sin
2π
x
W
cos
π
y
H
,
v
y
= v
y0
cos
2π
x
W
sin
π
y
H
,
where x and y are respectively horizontal and vertical coordinates inside the box in
m; W =1 000 000 m and H =1 500 000 m are the width and height of the model,
respectively (i.e., 1000 ×1500 km); v
x 0
=10
−9
m/s and v
y 0
=10
−9
m/s are scaling
values for respectively horizontal and vertical velocity components (10
−9
m/s ≈
3 cm/year). Compute (analytically) v
x
, v
y
,
∂v
x
∂x
,
∂v
y
∂y
and div( ¯v) =
∂v
x
∂x
+
∂v
y
∂y
on
a 2D grid of points (e.g. 31 ×31) which are regularly distributed inside the model
and visualise these parameters separately as colourmaps (pcolor) in order to see
how they are distributed relative to each other. Visualise the velocity as an arrow
field (quiver). Try to tune the scaling velocity values v
x 0
and v
y 0
, such that the
divergence of velocity goes to (almost) zero in the entire model (i.e. absolute values
of div(¯v) should be many orders of magnitude less then that of
∂v
x
∂x
and
∂v
y
∂y
). An
example is in Divergence.m.
2
Density and gravity
Theory: Density of rocks and minerals. Thermal expansion and com-
pressibility. Dependence of density on pressure and temperature. Equa-
tions of state. Poisson equation for gravitational potential and its
derivation.
Exercises: Computing and visualising density, thermal expansion and
compressibility.
2.1 Density of rocks and minerals. Equations of state
Many geodynamic processes are either directly or indirectly driven by the gravity
force due to the spatial variation of rock density inside the Earth. The density of
rocks (ρ) depends on pressure (P), temperature (T), chemical composition (C)and
mineralogical composition (M)
ρ = f (P,T,C,M). (2.1)
These factors are not fully independent. The mineralogical composition of a rock
with a constant composition may for example change due to changes in pressure and
temperature. Easy-to-remember densities of major rock types used in geodynamics
are:
felsic rocks (e.g., granites) ∼2700 kg/m
3
,
mafic rocks (e.g. basalts) ∼3000 kg/m
3
,
ultramafic rocks (e.g., peridotites) ∼3300 kg/m
3
.
Other relevant rock densities are given in Table 17.2 and in Appendix 2 of Turcotte
and Schubert (2002).
25
26 Density and gravity
Variations in the density of minerals and rocks with T and P are often charac-
terised by thermal expansion (α) and compressibility (β)
α =−
1
ρ
∂ρ
∂T
, (2.2)
β =
1
ρ
∂ρ
∂P
. (2.3)
The thermal expansion of rocks and minerals is typically positive (α>0), i.e.
density decreases with increasing temperature (cf. minus in the right part of
Eq. (2.2)). There are, indeed, exceptions, for example, beta-quartz possesses a
density which remains almost constant with increasing temperature at room pres-
sure (α ≈0). Typical values of α are on the order of n ×10
−5
K
−1
. Compressibility
is always positive (β>0) and density always increases with increasing pressure.
Typical values of β are on the order of n ×10
−2
GPa
−1
(or n × 10
−11
Pa
−1
). In
cases of constant α and β, integration of Equations (2.2)and(2.3) versus T and P
gives
ρ = ρ
r
e
β(P −P
r
)−α(T −T
r
)
, (2.4a)
where ρ
r
is the density of a given material at reference pressure P
r
(typically 10
5
Pa = 1 bar) and temperature T
r
(typically 298.15 K = 25
◦
C).
Since both α(T − T
r
)andβ(P − P
r
) are typically very small (much less then
unity), the equations can be simplified either to (using the rules that
e
a
≈ 1 + a
and
e
−a
≈ 1 − a when a 1)
ρ = ρ
r
[
1 + β(P −P
r
)
]
×
[
1 − α(T − T
r
)
]
, (2.4b)
or (using the rules that
e
a
≈ 1 + a and
e
−a
=
1
e
a
≈
1
1 + a
when a 1) to
ρ = ρ
r
1 + β(P −P
r
)
1 + α(T − T
r
)
. (2.4c)
Equations (2.4a)–(2.4c), however, do not account for changes in density due to
changes in mineralogical composition with changing pressure and temperature.
These equations of state (EOS) are much simplified since they are based on the
assumption that α and β remain constant with pressure and temperature. For real
rocks and minerals, these two parameters however strongly depend on both P and
T (Fig. 2.1) and therefore more realistic EOS must be used for describing the
density changes (e.g. Anderson, 1995). If compressibility of the mineral depends
on pressure, then the following well-known equations can be used: