See discussions, stats, and author profiles for this publication at: https://www.researchgate.

net/publication/328655704

Introduction to Tectonophysics

Book · November 2018

CITATIONS

1 author:

Patrice F Rey
The University of Sydney
155 PUBLICATIONS 2,347 CITATIONS

SEE PROFILE

Some of the authors of this publication are also working on these related projects:

BASIN GENESIS HUB – GEODYNAMICS AND EVOLUTION OF SEDIMENTARY SYSTEMS View project

AUGURY View project

All content following this page was uploaded by Patrice F Rey on 01 November 2018.

The user has requested enhancement of the downloaded file.

 Introduction to
Tectonophysics

Patrice F. Rey
Preface

i
PREFACE
This eBook, in a permanent state of unfinishedness, focusses mechanical coupling between the lithosphere and the astheno-
on tectonic processes, from the forces that drive them to the pa- sphere.
rameters that control the mechanical behavior of the Earth’s
How the Earth’s lithosphere deforms depends on the relation-
lithosphere.
ship between applied forces and the mechanical response
We start in Chapter 1 by looking at the Earth’s geotherm be- measured by the accumulated strain and the strain rate, hence
cause the mechanical behavior of rocks - and therefore that of Chapter 4 tackles the rheology of rocks. This chapter intro-
the lithosphere - is very sensitive to temperature and because duces the concepts of elastic, plastic and viscous strain and ex-
tectonic processes can significantly bring the geotherm out of plores the influence of parameters such as temperature, pres-
equilibrium forcing the cooling or heating of rocks. sure, and fluid, before presenting the important concept of
lithospheric strength envelopes.
The second chapter introduces the notion of isostatic equilib-
rium, critical to understand the topography of the Earth’s sur- Chapter 5 introduces the new discipline of computational tec-
face when its lithosphere is at rest. This chapter explains how tonics and geodynamics. Equipped with a reasonable under-
isostatic equilibrium differs from the notion of mechanical equi- standing of the thermal and mechanical properties of the litho-
librium. Isostasy introduces gravitational stresses because den- sphere, one can tap on the power of high-performance comput-
sity interfaces are no longer parallel to gravitational equipoten- ing to explore through numerically experiments tectonic and
tial surfaces. The notion of gravitational stress is important to geodynamic processes.
understand how the Earth’s lithosphere can deform in the ab-
Chapter 6 explores through numerical experiments a range of
sence of plate boundary forces.
tectonic processes at lithospheric scale, in a variety of tectonic
Plate boundary forces are the focus of Chapter 3. Slab-pull, settings from continental collisions, to extensional tectonics
ridge push, and drag from asthenospheric flow contribute to and transcurrent tectonics.
drive the motion of lithospheric plates at the Earth’s surface,
Finally, Chapter 7 focuses on the dynamics of mantle convec-
and their deformation. Ridge-push can be evaluated analyti-
tion.
cally using the concept of the gravitational forces. The evalua-
tion of the asthenospheric drag and that of the slab-pull is
more difficult because it requires some understanding of the copyright © Patrice F. Rey, 2018

ii
CHAPTER 1

The Earth’s Geotherm

The geotherm, i.e. the distribution of temperature with depth, is an important characteristic of the Earth's lithosphere
because temperature impacts on all physical properties of rocks (e.g. density, viscosity, conductivity, elasticity,
magnetism etc). In particular, temperature controls the rheology of rocks and therefore how they deform in response to
applied deviatoric stresses, and how the Earth's lithosphere reacts to tectonic forces. In this chapter we derive, from first
principles, a simple expression for the steady-state geotherm. We then consider the notion of transient geotherm.

3
SECTION 1

Heat Transfer in the Earth’s lithosphere

In this section

1. Heat transfer in the Earth’s

Lithosphere
2. Heat energy & temperature
3. Heat conduction
4. Heat advection
5. Radiogenic heat production

The continental geotherm is a function of the i/ rate at which given depth does not change through time). In contrast, when
heat is produced or consumed within the lithosphere, ii/ the the lithosphere has a net gain or a net loss of heat, the geotherm
rate at which the lithosphere looses heat to the atmosphere/ is said to be transient (i.e the temperature changes through
ocean system, and iii/ the rate at which the lithosphere gains time) until a new equilibrium is reached between heat lost and
heat from the hot convective mantle. When the heat lost bal- gained. On a billion year time scale, the geotherm is always tran-
ances the heat gained, an equilibrium is reached and the geo- sient because the primordial accretionary heat and the Earth’s
therm is said to be in steady state (i.e. the temperature at any supply in radiogenic isotopes progressively decrease, however,

4
on a scale of ~100 myr, and in the absence of tectonic or magmatic The main processes able to change the amount of heat energy in
activity, the geotherm can approach an equilibrium state. The litho- the lithosphere are:
spheric geotherm varies laterally. The base of the lithosphere is de-
•Heat conduction (transfer of kinetic energy between molecules
fined by the isotherm ~1300ºC. In cratonic regions, this isotherm
or atoms from a hot to a less hot region)
can stand at a depth between 200 and 300 km. In other continental
•Heat advection (replacement of a volume of rock at tempera-
regions it can be as shallow as 50 km, and at mid-oceanic ridge it is
ture T1 with an equivalent volume at temperature T2)
met at a few kilometers underneath the oceanic floor.
•Heat production (heat produced by radioactive isotopes, vis-
Here, we first review the processes involved in heat generation and cous heating, exothermic metamorphic reactions)
heat transfer, and we derive from the rate of these processes a gen- • Heat consumption (heat consumed by endothermic metamor-
eral equation which describes the change in temperature with phic reactions, in particular partial melting)
depth and through time. From this general equation we derive a
particular solution for the so called "steady state" continental geo-
The variation of temperature dT over an increment of time dt de-
therm (temperature changes with depth but not with time, i.e. zero
pends on the sum of heat variations dE due to each process. In
net heat gain or loss). In a second part, we discuss how the steady
what follows, we derive three expressions for i/ the rate of heat
state continental geotherm is affected by a number of geological
conduction, ii/ for the rate of heat advection, and iii/ for the rate
processes including, lithospheric thinning and thickening, burial
of radiogenic heating. From these, we derive the 1D conduction-
via sedimentary or volcanic processes, and basal heating via the
advection heat transfer equation from which an expression for the
spreading of mantle plumes at the base of the Earth's lithosphere.
steady state geotherm can be derived. Sounds more complicated
than it really is. So bear with me ...

Temperature and Heat

The temperature (degree of hotness or coldest) of a small volume Heat conduction

of rock somewhere in the lithosphere varies if heat energy (a form
Conduction transports heat from hot to cold regions. The flow of
of kinetic energy) is gained or lost. The relationship that gives the
heat (Q) is proportional to the negative temperature gradient (dT/
variation of temperature dT as a function of a variation of heat dE
dz) between the cold and the hot region, with the coefficient of pro-
is:
portionality being the conductivity (K). Mathematically this trans-
dT = dE / (Cp . m) lates into the Fourier's law where Q is in W.m-2 and K is in
W.m-1.K-1. In our geological reference frame, z increases downward
with Cp the heat capacity, and m the mass.

5
(T(z+dz)>T(z)). Conduction occurs in the direction of decreasing Radiogenic heat production
temperature (i.e. dEc is a gain for upward conduction) hence the z
Radiogenic disintegration of radioac-
sign "-" insures that Q is positive upward (dEc is positive when T
tive isotopes (238U, 235U, 232Th, and
increases downward).
40K) releases heat. The increment of

dT radiogenic heat (dEr) produced in a A

Q =−K
z
E1=-Q(z).a.dt
Q(z) dz small cylinder of rock of section a
dEr = A.a.dz.dt
dEc = E2 - E1 and length dz over an increment of
Let's consider a small cylinder of rock of section
dEc = (dQ/dz).dz.a.dt
time dt is:
a a (area in m2). If the incoming and outgoing z+dz a
z+dz
E2=-Q(z+dz).a.dt heat at both ends of the cylinder are the same, A⋅a⋅dz⋅dt = A⋅dV⋅dt
there is no net heat gain or loss, and the tem-
Q(z+dz)
perature remains unchanged. Temperature where A is the rate of radiogenic heat production. Radioactive heat
changes when the heat E1 leaving the volume over an increment of is the main internal heat source for the earth as a whole (it is meas-
time (E1=Q(z) ⋅ a ⋅ dt) is different to the heat E2 entering it ured in W.m-3).
(E2=Q(z + dz)⋅a⋅ dt).

The entering heat flow Q(z + dz) can be approximated with a Tay- Heat Advection
lor series in which only the two first terms are of significance. (n.b.:
mathematically f (xn + d x) can be approximated from f (xn) and the Advection of heat implies that
derivatives at location xn : f ′(xn) , f ′′(xn) , etc: an mass of material at tem-
E1 =Cp.ρ.a.uz .dt.T
a.uz . dt
perature T (in yellow on the z
sketch) is being pushed out of
dQ dz 2 d 2Q our cylinder and replaced by dEu = E2 - E1
Q(z + dz) = Q(z) + dz⋅ + ⋅ + ⋅⋅⋅ dEu = C .ρ.a.u .dt.dT
dz 2 dz 2 an equivalent mass of mate- p z
rial at temperature T+dT. The
Therefore the increment of heat (dEc) gained or lost in an increment increment of heat gained or
of time dt is dEc = E2 − E1 : a
lost (dEu) over an increment of a.uz . dt
E )
dQ d 2T time dt is proportional to the z+dz 2 =Cp.ρ.a.uz .dt.( T+dT
a⋅Q(z)⋅dt − a⋅Q(z + dz)⋅dt = − a⋅dz⋅ ⋅dt = K⋅ 2 ⋅dV⋅dt mass of material displaced (
dz dz
ρ ⋅a⋅uz ⋅dt ), the heat capacity

6
of the material (Cp in J . kg-1 . K-1: is the amount of energy required left to right, the conduction term, the term of radiogenic heat pro-
by 1 kg of the material to increase its temperature by 1 K), and the duction, and the advective term. In 3D, this equation becomes:
temperature contrast of the incoming and outgoing masses. (n.b.:
∂T K ∂2T ∂2T ∂2T A ∂T
because z increasing with depth, uz is negative for upward convec- = ⋅[ 2 + + ] + − uz ⋅
tion, however, dEu is a gain for upward convection hence the sign - ∂t ρ⋅Cp ∂x ∂y 2 ∂z 2 ρ⋅Cp ∂z
insures that when uz is negative (i.e. upward motion) the advective
heat gained is positive). We get that the rate of heat advection is:
−Cp ⋅ρ⋅a⋅uz ⋅dt ⋅dT The geotherm in the continental crust
Total heat gained or lost We derive here the equation for the steady state continental geo-
therm. For a steady state geotherm (no variation of temperature
Adding up the rate of conductive heat (dEc), the rate of convective
through time) we have by definition that dT/dt=0 and uz = 0 (no ad-
heat (dEu), and the rate of radiogenic heat (dEr) gives the total rate
vection or erosion or sedimentation); therefore the 1D heat transfer
of heat gained or lost. This variation of heat triggers a change in
d 2T A
temperature dT = dE/(Cp.m) therefore : equation simplifies to: " = −
dz 2 K
Cp.m. dT= dEr + dEu + dEc and: This is a second order differential equation. It says that the gradi-
ent of the temperature gradient is a constant equal to the negative
d 2T
Cp ⋅ρ⋅dV⋅dT = A⋅dV⋅dt − Cp ⋅ρ⋅a⋅uz ⋅dt ⋅dT + K⋅ 2 ⋅dV⋅dt ratio between the radiogenic heat production and conductivity.
dz This kind of equation is solved by integrating twice and by using
two boundary conditions. For example we may know the tempera-
This leads to the 1D conduction-advection heat transfer equation: ture at the surface let's say: T=To at z=0; and we may know the sur-
face heat flow for instance at z=0 the heat flow is -Qo (remember
dT K d 2T A dT
= ⋅ 2 + − uz ⋅ this is a positive number). Assuming that A is constant with depth,
dt ρ⋅Cp dz ρ⋅Cp dz the first integration led to the temperature gradient dT/dz:

This equation describes the variation of temperature with depth dT A

= − ⋅z + C1
and through time due to heat conduction, radiogenic heat and dz K
heat advection. This equation assumes no lateral heat flow (hence
This gives the gradient of the geotherm as a function of depth.
1D). On the right end side of the equation the three terms are, from
From this function we get that at the surface (i.e. z=0): dt/dz=C1

7
Because of the Fourier's law, the second boundary condition im- relationships, one for the crust (the very last equation), and one for
poses that at z=0: the lithospheric mantle.

dT Q0 If the production of radiogenic heat is zero in the mantle then the

=
dz K 1D conduction-advection heat transfer equation simplifies to:

by combining the last two equations for z=0 we get C1 = Qo/K, we d 2T

=0
therefore conclude that the thermal gradient at any depth z is ... dz 2
dT A Q
= − ⋅z + 0 Integrating twice we get:
dz K K
T(z) = C2 ⋅z + C3
Integrating a second time led to:
With the two following boundary conditions: T=Tc at z = zc = Moho
A 2 Q0
T(z) = − ⋅z + z + C2 depth, and Q = -Qm (the basal heat flow), we get that:
2K K Qm
C3 = Tc − ⋅z
The first boundary condition says that T=To at z=0, therefore we K c
get that: Qm
and that C2 =
A 2 Q0 K
T(z) = − ⋅z + z + T0
2K K The geotherm in the lithospheric mantle is therefore:

Qm
⋅ (z − zc) + Tc
This relationship is the steady state geotherm. It gives the distribu-
T(z) =
tion of temperature with depth in a layer with homogeneous radio- K
genic production A, conductivity K, with a surface temperature of
Hence, the geotherm in the lithosphere is defined by a two steps
To and a surface heat flow of Qo.
function:
The geotherm in the continental lithosphere and beyond
A 2 Q0
T(z) = − ⋅z + z + T0, for 0 < z ≤ zc
The continental lithosphere consists in two layers with contrasted 2K K
thermal properties. In particular the radiogenic heat production in
Qm
the mantle is negligible compared to that of the crust. The geo- T(z) = ⋅ (z − zc) + Tc, for z > zc
therm in continental lithosphere is therefore best described by two K

8
alternatively, if Qm instead of Qo is known the crustal geotherm (i.e. dT dT d 2T
A “differential equation” links together the derivatives ( , , 2 )
when 0 < z < zc): dt dz dz
dT dT
A 2 Q A of a function T(z, t). Here,
dt
is the rate of temperature change;
dz
T(z) = − ⋅z + [ m + ⋅zc]⋅z + T0, for 0 < z ≤ zc
2K K K is the gradient of temperature (i.e. the slope of the geotherm); and
d 2T
In the asthenospheric mantle, convective motion is such that the is the gradient of the gradient of temperature.
dz 2
temperature shows relatively little variation as convection acts as
an efficient mixing process. For z > zl (zl being the base of the litho- From this differential equation, one can derive a solution for the
sphere) the geotherm follows the adiabatic gradient ~0.3 K per kilo- particular case when thermal equilibrium is reached (i.e. there is no
dT
meter. For instance, from the top of the convective mantle (i.e. the change of temperature through time: = 0 and necessarily there
dt
base of the lithosphere) to the core/mantle boundary (i.e. a dis- is no advection of heat, hence uz = 0). This corresponds to the
tance of ca. 2700-2800km) temperature increases by only ~3000ºC. “steady geotherm” for which the differential equation simplifies to:
This contrast with the evolution of the temperature within the litho-
sphere where temperature increases from ca. 0ºC at the surface to d 2T A
T’’(z) = 2 = −
1300ºC at depth of 100-200km. dz K

This tells us that the change of the temperature gradient with

d dT
depth ( )( ) is equal to a constant -A/K (we assume here that
Appendix 1: Solving differential equations dz dz
the radiogenic production A and conductivity K are constant with
Very often the rate of natural processes (i.e. rate of heat conduction, depth). The negative sign implies that the rate of temperature in-
rate of heat advection, the rate of radiogenic heating etc) can easily crease decreases with depth.
be calculated or measured. In the context of the Earth’s geotherm,
we have seen that, from the rate of individual heat transfer mecha- Lets plot T’’(z):
nisms, we can express the overall rate of temperature change T''(z)
through time (t) at any depth (z) via the following differential equa-
tion:

dT K d 2T A dT z
= ⋅ 2 + − uz ⋅ -A/K
dt ρ ⋅Cp dz ρ ⋅Cp dz

9
 d 2T A The curve in blue is the geotherm, bearing in mind that this curve
The differential equation 2 = − can be solved to find out the
dz K is only valid for 0 < z < zmoho.
geotherm T(z). For this the differential equation is integrated twice:

The first integration leads to the temperature gradient :

T(z)
dT A
T’(z) = = − ⋅z + C1
dz K
T(0)
Graphically we get that:
T'(z)
zmoho z
C1 is therefore the temperature
c1
gradient at z=0 (Earth’s surface), -A/K
and -A/K is the slope of the gradient. We have recognized that the integration constants C1 and C2 have
the following meanings: C1 is the temperature gradient (the slope
of the geotherm) at the surface (i.e. at z=0); and C2 is the tempera-
z
-A/K ture at the surface. The temperature at the surface is typically in
T''(z) the range of 0 to 30ºC, and the temperature gradient at the surface
can be determined by measuring the temperature at the surface
The second integration leads to the geotherm: and at the bottom of a well a few 100s to a few 1000s of meter deep.
A 2
T(z) = − ⋅z + C1 ⋅z + C2 Appendix 2: Steady state geotherm with decreasing RHP
2K
The steady state crustal geotherm we derived earlier assumed that
Graphically we get that:
the radiogenic heat production is constant with depth. This is usu-
C2 is therefore the temperature ally not the case as cycles of partial melting and crustal differentia-
T'(z)
at z=0. tion through upward flow of melt tend to deplete the lower crust
c1 and enrich the upper crust in radiogenic elements. The decrease
-A/K
with depth of the radiogenic heat production is often modeled by
c2 an exponential law:
T(z)
−z
z A(z) = A0 . Exp( )
-A/K
hc
T''(z)

10
This means that the radiogenic heat production measured at the After substituting C1 and C2 into T(z), expanding and simplifying,
surface ( A0 ) in divided by e (e=2.71) every hc meters (typically one gets:
5000 m< hc <25000 m). Therefore the 1D conduction-advection
A0 2 −z Q A −z
dT K d 2T A(z) dT T(z) = ⋅hc ⋅(1 − Exp( )) + ( m − 0 ⋅hc ⋅Exp( c ))⋅z + T0
heat transfer equation is: = ⋅ 2 + − uz ⋅ K hc K K hc
dt ρ⋅Cp dz ρ⋅Cp dz
If one knows the surface heat flow instead of the mantle heat flow
In the case of a steady state geotherm, this equation becomes: then:
d 2T −A(z) A0 2 −z Q
= T(z) = ⋅hc ⋅(1 − Exp( )) + 0 ⋅z + T0
dz 2 K K hc K
d 2T A0 −z
replacing A(z) we get that = − ⋅Exp( )
dz 2 K hc
Following two successive integrations one gets:

dT A −z
= 0 ⋅hc ⋅Exp( ) + C1 and
dz K hc
A0 2 −z
T(z) = − ⋅hc ⋅Exp( ) + C1z + C2
K hc
To find out the two integration constants C1 and C2 we need to call
two boundary conditions. Lets say that we know the temperature
at the Earth' s surface T0 and the heat flow Qm entering the base of
the lithosphere.
A0 2
The first boundary condition leads to: C2 = T0 + ⋅hc
K
Qm A0 −z
... and second to: C1 = − ⋅hc ⋅Exp( c )
K K hc

11
SECTION 2

Earth’s Geotherm
In this section
TºC TºC TºC
0 400 800 1200 1600 0 400 800 1200 1600 0 400 800 1200 1600
1. Steady state geotherm
2. Transient geotherm 20 20 20

40 Moho 40 Moho 40 Moho to

t1
t5
60 60 Mechanical boundary layer to 60
t10
t100
Thermal boundary layer t1 t20
t5
80 80 t100 80
t10
t
t20

8
100 100 t 100

8
to

120 120 120

140 t 140 140

8
Depth (km)

Depth (km)

Depth (km)
Plume spreads at the top of the
Plume spreads under the lithosphere thermal boundary layer Plume spreads at the Moho

Steady state geotherm

tion). We then look at how sensitive is the geotherm to change
In the previous section we have derived the 1D conduction- in mantle heat flow, conductivity and thickness of the radio-
advection heat transfer equation, and we have derived an equa- genic crust.
tion describing the equilibrium geotherm, also called the
Distribution of radiogenic heat production
“steady state geotherm” because the temperature at every
depth doesn’t change through time. Here we explore how sensi- We have assumed so far that the volumetric radiogenic heat pro-
tive is the geotherm to the radiogenic heat production when it is duction A was depth-independent ( A remains constant with
not constant with depth (as we have assumed in the first sec- depth). However, because the upper crust is enriched in incom-

12
patible elements* the radiogenic heat production decreases with partial melting, they tend to concentrate into the melt phase. Due to the
depth. A common model assumes that A is divided by e (e=2.71) buoyancy of the melt, incompatible elements concentrate over time into
every h meters, h being the length scale of the exponential law. This the upper part of the crust. Through repeated cycles of melting, the con-
model of distribution is given by: centration of radiogenic elements contributes to the cooling of the
Earth’s geotherm.
A(z) = A0 ⋅Exp( − z /h)
Basal Heat Flow (Qm)
... in which A0 is the radiogenic heat production at z = 0 m.
The graph below illustrates the sensitivity of the geotherm to the
The graphs shows mantle heat flow (Qm), also
the geotherm for Radiogenic Heat (mW.m-3) Temperature (TºC)
2 4 6 8 500 1000
called in some textbook ba-
TºC
A(z) decreasing sal heat flow. This graph 0 500 1000 1500
0
exponentially shows three geotherms cal-
with various culated assuming same
25 25 20
length scale h (10, rate of radiogenic heat pro-
25, 50km). Ao is duction, same conductivity Moho
Moho 40
adjusted so that and same crustal thickness.
Depth (km)
50 50
Depth (km)

the total radio- TMoho: Only the mantle heat flow

724, 60
genic heat produc- varies.

Qm
935,
10 km
1284

=36
tion (R.H.P which 75 75
25 km

Qm
is given by the in- 50 km A basal heat flow of 12 x 80

=24
tegration of the R.H.P.: 0.078 W.m-2 10-3 W.m-2 (yellow geo-

Qm=12
100 100 therm) is characteristic of
radiogenic heat 100
Archaean (3.5-3 Ga)
profile with cold and thick cratonic

Depth (km)
depth) is the same lithospheres, whereas a ba-
in all models. The larger is h (i.e. the deeper are the radiogenic ele- sal heat flow of 36 x 10-3
ments) the hotter is the geotherm. This suggests that at an early W.m-2 is representative of
stage of the Earth evolution (before the extraction of the crust) the thinned lithosphere. In-
Earth had a warmer geotherm. creasing the mantle heat
flow from 12 to 36 x 10-3 W.m-2 increases the temperature at the
* n.b.: Incompatible elements such as U, Th and K have a large radius and Moho from ~320 ºC to 720 ºC.
therefore do not fit easily into crystals lattice. Hence, during episodes of

13
Thermal conductivity (K) thinning or fast thickening. During thinning the geotherm in-
creases as warmer temperatures are reached at lesser depths (the
The graph below illustrates the sensitivity of the geotherm to the
slope of the geotherm increases). During thickening the geotherm
thermal conductivity K (W.m-1.K-1). It shows three geotherms calcu-
decreases as cooler conditions are met deeper in the crust (the
lated assuming same rate of
slope of the geotherm decreases). Second, as the thickness of the
radiogenic heat production,
TºC crust changes, the steady state geotherm also changes because a
same same crustal thickness, 0 500 1000 1500
0 thicker crust will produce more radiogenic heat, therefore crustal
and same mantle heat flow.
thickening leads to warmer geotherm, whereas crustal thinning
The conductivity is known to
20 leads to cooler geotherm.
vary with temperature. Here
however, we assume that K is Moho Lithospheric thickening
40
constant through depth. K is
Thickening produces heat advection resulting in a rapid cooling of
proportional to the ability of
60 the geotherm. Thickening also increases the thickness of the radio-

k=
material to conduct heat

1.7
genic layer therefore increasing the production of radiogenic heat

k=2

5
away. Therefore, the larger 80

.25
in the lithosphere. Isostasy (cf. section Isostasy and Gravitational
the thermal conductivity of
Forces) produces uplift leading to erosion which in turn affects the
k=2
rocks, the lesser the capacity 100
amount and distribution of the radiogenic heat elements in the
of the crust to store heat, and .75
crust which cools geotherm. The interplay between the mode of the
Depth (km)

the cooler is the geotherm.

thickening (heterogeneous via thrusting vs homogeneous), the
Increasing the thermal con-
thickening rate, and the erosion rate leads to contrasted thermal his-
ductivity from 1.75 to 2.75
tories.
W.m-1.K-1 decreases the tem-
The graphs on the next page show transient geotherms (at 0, 0.5, 2,
perature at the Moho from
10...Myr) following heterogeneous thickening (thickening is
750 to 500ºC.
achieved by doubling the thickness of the crust via a single thrust,
Thickness of the radiogenic crust which explains the temperature discontinuity at t=0), and homoge-
neous thickening (the thicknesses of the crust and the lithospheric
Crustal thickening or thinning modify the geotherm in two differ-
mantle are doubled via pure shear deformation). Erosion is dis-
ent ways. First, during deformation, heat is advected mainly up-
carded here. The discontinuity in the case of heterogeneous thicken-
ward (in the case of thinning) and mainly downward (in the case of
ing is smoothed out in a few Ma. Transient geotherms in both cases
thickening) as rocks carry their temperature with them during fast
are similar for time > 10 Myr.

14
 stasy leads the subsidence of the surface of the lithosphere leading
0 200 400 600 800 1000 1200 TºC 0 200 400 600 800 1000 1200 TºC
0 0 to sedimentation which in turn affects the amount and distribution
of the radiogenic heat elements in the crust, which also impacts on
the geotherm. The interplay between the geometry of the thinning,
the thinning rate, and the sedimentation rate lead to contrasted

10
50 50

m
yr
0. thermal histories.

50
5m
2

250
10

yr

m
m

50

25
yr

yr
The graph on the bottom left shows transient geotherms (0, 10, 50,
my

0m

my
m

Moho Moho
r

yr

r
yr 100, 250... myr) following homogeneous thinning (the thicknesses
100 100
of the crust and the lithospheric mantle are halved by pure shear
deformation). Sedimentation is discarded here. Following the in-
crease of the geothermal gradient due to extensional deformation,
150 150
thermal relaxation leads to cooling and therefore the thickening of
Depth (km)
Depth (km)

Homogeneous thickening via pure shear the lithosphere. Slowly, the transient geotherm approaches the
Heterogeneous thickening via thrusting doubling the thickness of the lithosphere
and doubling the thickness of the crust steady state geotherm.

Sedimentation of burial of the crust

Sedimentation and burial of the continental crust under a few kilo-

Lithospheric thinning meter of volcano-sedimentary rocks have a significant long-term
0 200 400 600 800 1000 1200 TºC
0
Thinning drives heat advec- effect on the geotherm. The newly deposited layers effectively insu-
tion as rocks, and the heat at- late the buried heat producing layer. If the conductivity of the up-
Moho
tached to them, are displaced per layer is lower or equal to that of the buried layer, then heat ac-
vertically mainly upward re- cumulates increasing the geotherm.
50
sulting in a rapid warming of This effect could have played a major role in the differentiation of
the geotherm (the geothermal 10 the continental crust in the Archaean. At that time, radiogenic heat
50 my
gradient increases). However, m r production in the crust was 2 to 6 time larger that of present day
yr
25
0

rate of radiogenic heat production. Furthermore it was a time

10

thinning also decreases the

m

0
yr

100
when 5 to 20 km thick continental flood basalts (the so called green-
yr

thickness of the radiogenic

Depth (km)

crust, and therefore reduces stones) where deposited at the surface of the Earth, insulating the
60% homogeneous thinning of
the production of radiogenic the lithosphere via pure shear
heat producing crust.
heat in the lithosphere. Iso- The graphs show transient geotherms following the emplacement

15
 thick layer with a temperature of 1700ºC. The three graphs show
Archaean continental lithosphere Archaean continental lithosphere
TºC TºC
Greenstone 200 400 600 800 1000 1200 1400 Greenstone 200 400 600 800 1000 1200 1400
the results for various depth of emplacement.
Crust Crust
Crust
25 25 TºC TºC TºC
Geotherm to Geotherm to Geotherm t

8
0 400 800 1200 1600 0 400 800 1200 1600 0 400 800 1200 1600

to+50Ma to+10Ma
50 50
Geotherm t to+20Ma 20 20 20
to+100Ma

8
to+30Ma
75 to+200Ma 75 40 Moho 40 Moho 40 Moho to
to+40Ma
Mantle Mantle
t1
to+400Ma t5
to+50Ma 60 60 60
100 100 Mechanical boundary layer to
t100
t10
Thermal boundary layer t1 t20
to+100Ma t5

Depth (km)
Depth (km)

80 80 t100 80
t10
125 125 to+200Ma t
t20

8
to+1Ga
100 100 t 100

8
to
150 150
120 120 120

175 175
140 t 140 140

8
Depth (km)

Depth (km)

Depth (km)
Plume spreads at the top of the
Plume spreads under the lithosphere thermal boundary layer Plume spreads at the Moho

of a 6 km (graph on the left) or 12 km thick greenstone covers

(graph on the right), with no radiogenic heat production in them.
Transient geotherms
The temperature increase is large enough to lead to significant par-
tial melting of the crust. Calculation of the steady state geotherm is relatively straightfor-
ward. In contrast, the calculation of transient geotherms, such as
Emplacement of a mantle plume those displayed above, requires computational tricks as the tem-
perature changes in both space and time. Many analytical solutions
Mantle plumes initiate at the core-mantle boundary and rise fast
have been proposed for a range of problems. The well-known book
through the convective mantle. Upon approaching the more rigid
from Carslaw and Jagger, first published in 1946 "Conduction of
lithosphere, they spread laterally under the lithosphere. They may
Heat in Solids", provides hundreds of analytical solutions to many
also thermally erode to base of the lithosphere and spread higher
problems. Here, we give a couple of them.
up. This is equivalent to put a hot layer, a few tens of kilometer
thick, underneath or within the colder lithosphere. In the Ar- Progressive cooling of the oceanic lithosphere
chaean, the Earth was warmer and plumes were most likely more
numerous. The next graphs document the thermal impact of a man- The formation of oceanic lithosphere at mid-oceanic-ridge is a ther-
tle plume assuming that the plume's head spreads into a 50 km mal problem involving the progressive cooling of the astheno-
sphere. The temperature at the seafloor Ts is maintained constant,

16
so is the temperature in the asthenosphere Tm. The lithosphere be- (he assumed Tm =2000ºC) and that Earth had
comes cooler and therefore thicker as cooling proceeds. Neglecting cooled by conduction to its present surface
radiogenic heat in the basaltic crust and the mantle, and assuming heat flow (he assumed Q0= -30 10-3 W.m-2). He
no sedimentation, the 1D heat transfer equation becomes: assumed that all heat was lost at the surface
by conduction (with κ = 10-6.m-2.s-1). Thus the
dT d 2T
= κ⋅ 2 problem is reduced to that of finding the tem-
dt dz perature within a cooling half-space of infi-
nite extent as a function of time after the half
with κ the thermal diffusivity which can be expressed in terms of
space is set at a specific temperature, a prob-
conductivity (K), density ( ρ ) and heat capacity (Cp): κ = K⋅ρ⋅Cp lem for which the equations above apply. Kelvin found that the
Earth was 50 myr old. This 2 order of magnitude off the mark, but
The analytical solution of this differential equation is: still much better than the biblical 6000 yr.
z
T(z, t) = Ts + (Tm − Ts)⋅erf( ) Cooling of a dike
2 κ⋅t
The cooling of a dike with half width w is another problem which
By differentiating with respect to z one get the temperature gradi- has an analytical solution of the following form:
ent ...
Ts w−x w+x
dT Tm − Ts z2 T(x, t) = ⋅(erf( ) + erf( ))
= ⋅Exp( − ) 2 2 κ⋅t 2 κ⋅t
dz π⋅κ⋅t 4⋅κ⋅t
If the dike has a width of 2 m, i.e. w=1 m, and its initial tempera-
...from which the variation of the surface heat flow through time ture was Ts = 1000 °C, and κ = 10-6 m2s-1, then the temperature at
can get extracted as: the center of the dike would be about 640°C after one week, 340°C
after one month, and only 100 °C after one year.
dT T − Ts
= m
dz π⋅κ⋅t
Some swarms, such as the giant McKenzie swarm in Canada, ap-
pear to radiate from a point, commonly interpreted as a plume
source for the magmas. The mid-Proterozoic Coppermine River
In the 19th century, Lord Kelvin used this equation to find out the
flood basalts were erupted at the same time near the plume head.
age of the Earth. Kelvin made the assumption that the Earth had
Individual dykes range from 10 to 50 m in width, with some up to
formed as a molten body at the temperature at which basalt melts
200 m wide. Some dykes can be traced for up to 2000 km.

17
Numerical solutions and therefore:

Since the arrival of com- ∂T L ∂ 2T A ∂T

puters, which allow for (1) − (1 − ) = κ⋅ 2 + −U
∂t Cp ⋅(Tliq − Tsol) ∂z ρ⋅Cp ∂z
large numbers of calcula-
tion to be performed Boundary conditions: T(z, 0)=0; T(0, t)=T0; dT/dz(zl, t)=-Qm/K
very fast, new tech-
niques based on numeri- An approximation of this differential equation can be obtained by
cal algorithms have been replacing ∂t and ∂z with finite differences ∆t and h. The central
designed to solve differ- tenet of this computational technique is that for a given time t=n,
ential equations. These the temperature at a depth z = i (Tin) can be calculated from knowl-
numerical recipes are based on the "discretization" of differential edge of the temperature at depth z = i-h and z = i+h. In a similar
equations ... fashion, for a given depth z = i, temperature at time t = n can be cal-
culated from knowledge or the temperature at time t = n-∆t and
• Consider the 1D heat conduction-advection equation in a slab zl temperature at time t = n+ ∆t. The level a accuracy depend of the
thick. The upper surface of the slab has at temperature To. The tran- size of the time and space finite differences. There are many ways
sient temperature is then described by the conductive-advective to express T(z, t) as a function of T(z-h, t), T(z+h, t), T(z, t- ∆), T(z, t+
equation of heat balance. Here we consider a situation where inter- ∆). Here we present the Crank-Nicholson scheme where, equation
nal radiogenic heat is creation at a rate A, heat is lost or gain by ad- (1) is rewritten as:
vection at speed U (since z increases downward, if erosion U < 0, if
sedimentation U > 0), and heat is lost through phase transition (2)
solid-liquid (L < 0 for partial melting, L > 0 for crystallization) X
Tin+1 − Tin
( Cp ⋅(Tliq − Tsol) )
L
the melt fraction is a function of T and therefore z. − 1−
Δt
∂T ∂ 2T A ∂T ∂X L
= κ⋅ 2 + −U + ⋅ n
− 2Tin + Ti−1
n n+1
− 2Tin+1 + Ti−1
n+1
Tin+1 − Tin
2 ( )
∂t ∂z ρ⋅Cp ∂z ∂t Cp 1 Ti+1 Ti+1 A
= κ⋅ ⋅ + −U +
h2 h2 h ρ⋅Cp
re-arranging we get:
where: Tin = T(zi, tn), T n+1
i = T(zi, tn+1), etc
∂T 1 ∂T ∂ 2T A ∂T
−( ⋅ ) = κ⋅ 2 + −U
∂t Tliq − Tsol ∂t ∂z ρ⋅Cp ∂z

18
The finite difference equation (2) is known as the Crank-Nicholson This formulation results in N equations and N+4 unknowns, the
scheme. It is based on central differences for the spatial derivatives aver- four additional unknowns being u0n and uN+1n, u0n+1, uN+1n+1. How-
aged forward in time over time steps n and n+1. ever, these four additional unknowns lie outside the computational
grid (time, space). Also, from the boundary conditions the follow-
Expanding (2) and re-arranging to express temperature at time n+1
ing conditions must be satisfied at the nth and n+1 time steps:
as a function of the temperature at time n we get equation (3):
T1n = 1, TNn = 0, T1n+1 = 1, TNn+1 = 0

( h ) i ( 2h 2 )
−κ ⋅ Δt 1+n L κ ⋅ Δt UΔt 1+n −κ ⋅ Δt 1+n
⋅ T−1+i + 1 − + + ⋅T + ⋅ T1+i
2h 2 Cp ⋅(Tliq − Tsol) h2 This gives the following set of N equations (5) for the temperature

( h ) ( 2h 2 ) ρ ⋅ Cp
κ ⋅Δt n L κ ⋅Δt UΔt n κ ⋅ Δt n A ⋅Δt at the n+1 time step :
= ⋅T −1+i + 1 − − + ⋅T i + ⋅T 1+i +
2h 2 Cp ⋅ (Tliq − Tsol) h2
A ⋅ Δt
replacing and substituting: a11 ⋅ T01+n + b11 ⋅ Ti1+n + c11 ⋅ T21+n = a1 ⋅ T0n + b1 ⋅ T1n + c1 ⋅ T2n + , i=1
k
κ ⋅Δt U⋅Δt L 1+n A ⋅ Δt
Rd = , V = , Lt = a1i ⋅ Ti−1 + b1i ⋅ Ti1+n + c1i ⋅ Ti+1
1+n n
= a1 ⋅ Ti−1 + bi ⋅ Tin + c1 ⋅ Ti+1
n
, i = 2 to N − 1
+
2h 2 h Cp(Tliq − Tsol) k
1+n A ⋅ Δt
a1N ⋅ TN−1 + b1N ⋅ TN1+n + c1N ⋅ TN+1
1+n n
= aN ⋅ TN−1 + bN ⋅ TNn + cN ⋅ TN+1
n
+ , i=N
k
we get...
1+n
−Rd ⋅T−1+i + (1 + 2Rd + V − Lt) ⋅ Ti1+n − Rd ⋅T1+i
1+n
with (6):
A⋅Δt
n
= Rd ⋅T−1+i + (1 − 2Rd − Lt + V ) ⋅Tin + (Rd ⋅T1+i
n
)+ b11 = 1, c11 = 0, a1N = 0, b1N = 1
ρ ⋅Cp
b1 = 1, c1 = 0, aN = 0, bN = 0
further replacement...
These latter conditions ensure that the boundary conditions are al-
a1i = − Rd, c1i = −Rd, b1i = (1 + 2Rd + V − Lt)
ways satisfied.
and ai = Rd, ci = Rd, bi = (1 − 2Rd + V − Lt)
The N equations above can be rearranged into the following matrix
Rearranging terms, we obtain the following set of implicit equa- equations:
tions for the temperature at the n+1 time step in terms of the tem-
perature at the nth time step (4):
1+n A⋅ Δt
a1i ⋅T−1+i + b1i ⋅Ti1+n + c1i ⋅T1+i
1+n n
= ai ⋅ T−1+i + bi ⋅ Tin + ci ⋅ T1+i
n
+ = din
k

19
 b11 c11 0 0 .. .. 0 Numerical solutions an example
T1n+1 d1n
a12 b12 c12 0 .. .. 0 500 1000 TºC
T2n+1 d2n • We illustrate here the problem
0 a13 b13 c13 0 .. 0
T3n+1 d3n where a continental crust with initial
0 .. .. .. .. .. 0 ⋅ . = . (7)
.. .. .. .. .. . . thickness zc is thickened by a factor 2
0 0 n
n+1
TN−1 dN−1 via the emplacement of one single zt
0 .. .. .. a 1N−1 b1N−1 c1N−1
50

T n+1 dNn thick thrust. zc = 40 km, zt = 40 km.

0 .. .. .. .. a1N b1N N
Moho
The figure on the right shows the in-
100
with: stantaneous and steady state geo-
therm (potential geotherm reached
b1 c1 0 0 .. .. 0

Depth (km)
after infinite time). The aims is to dis-
d1n 1
a2 b2 c2 0 .. .. 0 play transient geotherms down to a
T2n
Heterogeneous thickening via thrusting
d2n and doubling the thickness of the crust
0 a3 b3 c3 0 .. 0 depth z = zl + zt at various time inter-
d3n T3n A⋅Δt
. = 0 .. .. .. .. .. 0 ⋅ .. ⋅ k vals, where zl is the thickness of the
.
n 0 .. .. .. .. .. 0 n
TN−1 lithosphere before thickening zl = 120 km.
dN−1
0 .. .. .. aN−1 bN−1 cN−1
TNn
dNn • We choose a spatial finite difference h = 4000 m (spatial grid).
0 .. .. .. .. aN bN
Therefore the number of column in the tridiagonal matrixes will be
zl/h. The Crank Nicholson scheme imposes the maximum time step
Thus the RHS of (7) is determined explicitly from the solution at
of h2/(2.κ) = 253,000 yr.
the nth time step. Note that the boundary conditions are satisfied by
ensuring that the coefficients given in (6) are satisfied. Hence in or- • We get that Rd = (κ.∆t)/(2h2) = 0.25. We assume no erosion: v = 0;
der to find the solution at the time step, one must solve N linear and we allow for partial melting:
equations. Since the coefficient matrixes in (7) and (8) are tridiago- Lt=Latent_heat/(Cp.(Tliquidus-Tsolidus)). With Rd, v and Lt we can
nal, one can make use of efficient algorithms (e.g. Thomas algo- determine ai1, bi1, ci1, the coefficient of the tridiagonal matrix at
rithm) to find the solution. time n+1, and ai, bi, and ci the coefficient of the tridiagonal matrix at
time n.

We choose to get the transient geotherm at 20 Myr interval up to

400 myr. The depth of the tridiagonal matrixes (number of row) is
therefore 400/20 = 20.

20
• With this, we construct the tridiagonal coefficient matrixes (7)
and (8) , the equation system (7) is solved via matrix inversion
(Numpy, R, MatLab or Mathematica are well equipped to do the 0 200 400 600 800 1000 1200 TºC
0
dirty work...). Graph on the right show transient geotherms from
0.5 to 250 myr and to infinite time.

50
0.
5m

10
yr

50

25
yr
m

0m
my
yr
Moho

yr
(7)

100

150

Depth (km)
Heterogeneous thickening via thrusting
(8) and doubling the thickness of the crust

21
Review 1.1 Heat transfer & the Geotherm Review 1.2 Heat transfer & the Geotherm Review 1.3 Heat transfer & the Geotherm

The relationship between heat and Using the Taylor series and bearing in mind the Fou- List five processes consuming or re-
temperature is: rier’s law: Q = -K . ( dT/dz), which of the following leasing heat into a volume of rock
is/are correct: somewhere in the crust.

A. T(z+dz) = T(z) + dz x dT A.
A. T = E x Cp x m
B. T(z+dz) = T(z) + dz x (dT/dz) B.
B. dE = dT / (Cp x m)
C. T(z+dz) = T(z) x (dT/dz) C.
C. dE = dT x (Cp x m)
D. T(z+dz) = T(z) - dz x (Q/K) D.
D. dT=dE / (Cp x m)
E. T(z+dz) = T(z) + (dT/dz) x (Q/K) E.

Check Answer Check Answer Check Answer

In a mathematical equation involving dimensional parameters, physical dimensions such as meter (m), second (s), kilogram (kg), Joule (J), Kel-
vin (K) or degree centigrade (ºC) can be balanced to verify whether or not a proposed relationship between physical quantities is correct. For in-
stance lets propose that the velocity (v) is a relationship between time (t) and distance (l) such that : v = t / l2.
Knowing that the unit of velocity is m/s, it is clear that the proposed relationship is wrong since the unit of the right hand side of the proposed
equation is s.m-2 . Considering that velocity is in m/s, one can easily figure out that the correction relationship must be: v = l/t.
Review 1.4 Heat transfer and Geotherm Expected learning outcomes

Which of the following is/are correct: At this stage you should be able to:

• Explain why and how temperature increases with

depth. This means: 1/ being able to comment on the
A. The long term consequence of crustal various terms of the equation describing the steady
state geotherm, and the meaning for the various pa-
thinning is a hotter geotherm. rameters; and 2/ being able to explain the shape the
geotherm in a temperature versus depth graph. This
B. The long term consequence of crustal should involve concepts like conductivity, advection,
radiogenic heat production etc.
thickening is a hotter geotherm.
• Explain the relationship between heat flow, the tem-
C. The long term consequence of deposi- perature gradient and thermal conductivity.
tion of a thick layer of sediments is a
• Explain how tectonic and other processes can
colder crustal geotherm. change the lithospheric geotherm i/ on the short
terms through advective heat transfer, and ii/ on the
D. Higher rock conductivity results in long term through head conduction.

hotter geotherm. • At a more advanced level you should be able to de-

rive an expression of the steady state geotherm from
E. For the same total radiogenic iso- the 1-D conduction-advection heat transfer equation
and knowledge of two boundary conditions.
topes content, a colder syeady steady
state geotherm is reach when radio-
genic isotopes concentrate in the up-
per crust.

Check Answer
CHAPTER 2

Isostasy & Gravitational Forces

The Tibetan plateau stands 5000 m
above sea level. It is the surface
expression of a 75 km thick
continental crust, which reduces the
weight of the continental
lithosphere. It imparts a significant
horizontal gravitational push force
on its surrounding.
Isostasy is the physical process that
explains the surface elevation of the
lithosphere at rest. It is based on the
Archimedes’ principle which relates
the elevation of a floating body
above the surrounding fluid i/ to
the density contrast between the
floating body and the supporting
fluid, and ii/ the vertical length of
the floating body.

24
SECTION 1

Isostasy, elevation, gravitational potential energy &

gravitational force

Photo: DscobiePhotography.co.uk/creativecommon

At rest, a lithosphere with a thick crust stands above equipotential surfaces, introduces gravitational forces
adjacent non-thicken regions, and it is said to be in iso- driving the flow of rocks from regions of high pressure
static equilibrium. The notion of isostatic equilibrium to regions of low pressure. In this section, we docu-
must not be confused with mechanical equilibrium. In- ment the concept of isostasy and its relationship to sur-
deed, a horizontal gravitational force (a volume force) face elevation. We then introduce the notion of gravita-
appears when density interfaces are no longer parallel tional potential energy and that of horizontal gravita-
to gravity equipotential surfaces. Hence, the action of tional force.
isostasy, which forces density interfaces away from

25
The sketch below represents a schematic cross-section through a lithospheric plate. Although the thicknesses of both the crust (in pink) and
the sub-continental lithospheric mantle (in green) vary, isostatic equilibrium is maintained which explains the lateral variation of the sur-
face topography.
A B C
• Isostatic equilibrium implies that at - and below - a
particular depth called the compensation depth, the
pressure becomes hydrostatic (i.e. the pressure
shows no lateral variation). In other terms, at or be-
low the compensation level the weight of columns
with the same cross-sectional areas standing on the
same gravitational equipotential surface are the
same. The minimum depth for the compensation
level is that of the base of the thickest lithospheric col-
umn, in other terms it is the minimum depth as
Compensation level
which there is no lateral change in density.
Isostatic equilibrium=>Lithostatic pressure@z = constant for z > Compensation level
• Isostasy controls the elevation of the Earth's sur-
face. Lithospheres with thin/thick continental crusts have lower/higher surface elevation compared to the average lithosphere. One can
easily determine the surface elevation of column A and C with respect to the elevation of the undeformed column B, the surface of which is
assumed here to be at sea level. When this heterogeneous lithosphere is at rest then isostatic equilibrium is reached via differential uplift
and subsidence, and the pressure at the base of columns A, B and C, all standing on the compensation level, is the same. The pressure at
the base of column B is: PB = ρcrust ⋅g⋅zc + ρsclm ⋅g⋅(zl − zc) + ρasth ⋅g⋅( fl ⋅zl + h − zl ), and the pressure at the base of column C is:
PC = ρcrust ⋅g⋅fc ⋅zc + ρsclm ⋅g⋅ ( fl ⋅zl − fc ⋅zc)

where zc is the thickness of the crust, zl is the thickness of the lithosphere (i.e. crust plus sub-continental lithospheric mantle), h the eleva-
tion of column C, and ρcrust , ρsclm and ρasth are the density of the crust, the sub-continental lithospheric mantle and the asthenosphere respec-
tively. To calculate elevation of the high plateau h one has to solve PB = PC . The calculate the bathymetry of column A one has to solve
PA = PB .

26
The sketch on the right illustrates the variation with depth of the lithostatic pressure
B C σzz = ρ.g.z(MPa)
(P(z) = ρ ⋅g⋅ z) along two lithospheric columns (B and C) in isostatic equilibrium. We
assumed here that the density of the crust and that of the mantle are are the same in
both columns, and that there is no density difference between the lithospheric man-

ρg
tle and the asthenosphere. Because it has a thicker crust, column C has a surface ele-

ρg

z
z

( C)
vation higher than that of column B. Consequently, the lithostatic pressure at any

(B)
GPE (C)
depth within the crust of column C is a higher than the lithostatic pressure at the
same depth in column B. However, at or below the compensation depth (here the
base of crust in column C) the lithostatic pressure measured at a particular depth is
GPE (B)
the same in both column. This is because both column are in isostatic equilibrium.

Isostatic equilibrium does not mean mechanical equilibrium. Indeed column C has a
larger Gravitational Potential Energy (GPE in Pa) than column B. On the graph the
GPE of each column is equal to the area under their respective lithostatic pressure
GPE (B) < GPE (C) => Horizontal force acting from C to B
profile and down to the compensation level (i.e. orange area for column B, yellow
area for column C which extends underneath the orange area).
What is the implication of this?

The graph on the lower right shows water jets driven by the hydrostatic pressure in a leaking tank. Water tank
The horizontal outreach of each jet is proportional to the lithostatic pressure at the outlet. If we con- free surface

nect two tanks of same height but different water depth, the difference in water level - and therefore Low pressure

in pressure - drives an outflow has shown in the sketch on the

Increasing pressure
dropping free surface left. This flow stops when the water level, and therefore the pres-
sure, in both tanks equalizes. The same physics applies to our two
rising free surface lithospheric columns. Hence, lateral variations of density produce
Outflow horizontal gravitational forces (body forces) acting inside the
lithosphere and pushing material laterally from high pressure to
High pressure
low pressure regions. The flow stops when all the density inter-
faces are parallel to equipotential gravitational surfaces, or when
the gravitational deviatoric stresses are no longer strong enough to overcome the resistance to flow.

27
The gravitational force (Fg in N.m-1) acting between C and B is equal to the difference of their B C σzz = ρ.g.z(MPa)
respective GPE (difference between the two colored surface areas, i.e. the thin yellow wedge
on the graph below). The excess in GPE in column C (Fg > 0 N.m-1) will drive extension,

ρg
ρg
whereas a deficit (Fg < 0 N.m-1) will drive contraction. Mathematically:

z(
z(

C)
B)
GPE (C)
Top Top

∫ ∫
GPEC = ρC(z)⋅g⋅z dz and GPEB = ρB(z)⋅g⋅z dz
GPE (B)
Bottom Bottom

Fg = GPEC − GPEB
GPE (B) < GPE (C) => Horizontal force acting from C to B

This movie below is the result of a numerical experiment documenting the gravitational collapse of an orogenic plateau. This model shows
a 65 km thick crust adjacent to a 40 km thick crust (brown), with the upper mantle in green. The experiment starts at a stage when thermal
relaxation has produced 20 to 25% of melt (dark pink) in the plateau’s lower crust. The only force acting is the gravitational force (i.e. the
vertical sides are fixed). The flow of the orogenic plateau is powered by the same physics responsible of the spreading of camembert.

Movie 2.1 Collapse of orogenic plateaux

28
THE CASE OF LITHOSPHERIC THICKENING AND DEBLOBBING Δσzz (MPa)
Compression 50 100 Tension

The figures on the right show the gravitational force (Fg) due to instantaneous homogeneous thick-
ening (at time to), and following convective thinning (to + 30 myr). Convective thinning is the gravi- 50
tational process upon which the cold and dense lithospheric keel is dragged into the convective

Depth (km)
mantle. Following homogeneous thickening, the excess in gravitational potential energy in the up-
100
per part of the thickened lithosphere is more or less balanced by a deficit in gravitational potential
energy in the deeper part of the lithosphere. Indeed the integration of σzz with depth gives a very
small gravitational force (Fg = 0.77 1012 N.m-1). 150

Following convective thinning (the cold and heavy lithospheric keel has been removed), the de- to Fg=0.77xTNm
-1

formed lithosphere has a larger excess in gravitational potential energy. This excess in GPE gives a
Δσzz (MPa)
gravitational force Fg = 6.7 1012 N.m-1. Such an Fg could easily balance or overcome tectonic forces, -50 50 100 150

and in some circumstances could be sufficient to drive extensional collapse of the mountain belt.
20

40

60

THE CASE OF LITHOSPHERIC THICKENING AND DEBLOBBING 80

-1
100 Fg=6.7xTNm
Gravitational force (Fg) due to instantaneous homogeneous stretching (at time to) and following to+30Ma Δσzz (MPa)
Compression -50 50 Tension
thermal relaxation (at time to + 200 myr). The yellow shaded areas represent the gravitational force
acting on the lithosphere. Fg is given by the integration with depth of the difference in lithostatic
pressure (σzz) between the thinned and the undeformed lithosphere. 50

Depth (km)
Following homogeneous thinning by a factor of 2 (50% reduction of the thickness of the crust and 100 -1
to Fg=0.17xTNm
that the whole lithosphere), Fg is rather small (< 1012 N.m-1).
Δσzz (MPa)
-150 -100 -50 50 100 150
As thermal relaxation and cooling proceeds, the thickness of the lithospheric mantle increases. This
20
leads to a decrease in the extensional gravitational force acting on the lower part of the lithosphere
40
and an increase in the compressional gravitational force acting in its upper part. The integration of

Depth (km)
60

Fg with depth gives a value > 1012 N.m-1. For contraction structures to develop Fg must be larger 80
-1
that the integrated strength of the thinned lithosphere. to+200Ma
100 Fg=1.2xTNm

29
GLOBAL VS LOCAL GRAVITATIONAL STRESS FIELDS

The plate average gravitational potential energy (GPE) defines a global reference level. A lithospheric column with an excess of GPE with
respect to this reference will be under extensional stresses, or under compressional stresses if it has a deficit. Continental and oceanic col-
umns with an elevation of ~70 m and -4.3 km respectively are in mechanical equilibrium with the global reference GPE. Mid-oceanic ridges
(MOR) have an excess of GPE which explains why they are in extension. Let's focus on the column B on the graph below. Depending on
its elevation, B has either an excess (B2) or a deficit (B1) of GPE with respect to the GPE of MOR. What would be evolution of the gravity-
related stress acting on B during thinning, which will shift the geometry of B toward that of C?

We consider first the case where the potential energy of column B is B1< GPE of MOR. As B becomes thinner, its GPE must increases to
reach that of MOR (path B1->C1). Therefore, column C is in extension with respect to its surrounding (column B). It is also in extension
with respect to the global Earth's litho-
sphere since C has a GPE C1 greater A B C
than that of the global mean potential. Thick Normal Thin
Hence, the regional gravity-related Potential Energy
Extension Pe>0 B2 Elevation B2>~70m
state of stress enhances the global M.O.R. C2
gravity-related state of stress. -4.3km C1 Plate's mean potential energy
In contrast, if the column B has a GPE B1 Elevation B1~70m
Compression Pe<0
B2 then thinning of B will lead to a de-
crease in GPE (path B2->C2) . With a M.O.R. A M.O.R.
B C
GPE C2, column C is in compression Crust
since it has a deficit of GPE with its im- Oceanic lithosphere
Lithospheric
mediate surroundings (B2). This re- Mantle
gional compression opposes the global Not to scale
extensional state of stress related to the
∆GPE between the thinned lithosphere (C2) and the global mean potential.

From this analysis we can define two gravitational stress field components: The Global Gravitational Stress Field (GGSF), and the Regional
Gravitational Stress Field (RGSF). The GGSF results from ∆GPE between a lithospheric column and the global mean potential energy,
whereas the RGSF results from the contrast in ∆GPE between a deformed lithospheric column and its immediate surrounding. The Effec-
tive Gravitational Stress Field (EGSF) is the superimposition of both gravitational stress fields.

30
SECTION 2

Consequences for deformation of the lithosphere

Horizontal gravitational forces arises when density inter- strength of rocks. Under the action of gravitational force,
faces are no longer parallel to equipotential surfaces. hot rocks flow from regions of high-pressure to regions of
These volume forces interfere in the development of oro- low pressure until all density interfaces are parallel to
gens and possibly in the development of rift basins. In equipotential surfaces. The flow can last for millions of
both cases gravitational forces tend to oppose the tectonic years, and occur in the absence of plate tectonic forces.
force driving thickening or thinning. When the tectonic
force wane, gravitational forces are only resisted by the

31
GRAVITATIONAL FORCE AND THE FORMATION OF CONTINENTAL PLATEAUX

Depending on its sign, GPE may oppose or enhance thickening. This exam-
ple illustrates the transition from a narrow mountain belt to an broad oro- 0 myr
genic plateau due to the build up of GPE. The upper and lower continen-
tal crust are in yellow and orange respectively. The lithospheric mantle is
in green and the asthenosphere in pale yellow. The dark shading is propor-
tional to stain rate. 5 myr

Following continental extension and rifting a narrow oceanic basin has

formed (0 myr). This basin is then inverted (i.e. the tectonic regime
changes from extensional to contractional) and the thinned continental
margins collide forming a narrow mountain belt via localized crustal scale 10 myr
reverse faults (5 to 10 myr).

As the crust becomes thicker, the excess in GPE produces a gravitational

force (Fg) that opposes the tectonic driving force (Fd). At one stage the
sum of Fg and the viscous strength of the thick crust (also a force to be
15 myr
overcomed) balances Fd. As a result deformation migrates into adjacent
weaker regions (15 to 25 myr). The mountain grows laterally and a pla-
teau develops.

Should the tectonic driving force Fd diminish or stop, there would be no 20 myr
forces to oppose Fg (beside the strength of rocks). This will led to gravita-
tional collapse of the orogenic plateau.

25 myr

32
GRAVITATIONAL FORCE AND CONVECTIVE THINNING Expected learning outcomes

Convective thinning describes the drag into the convective At this stage you should be able to:
mantle of the lower part of the lithospheric mantle. This
mechanism corresponds to the development of a Rayleigh- • Use your understanding of the Archimedes’ principle
Taylor gravitational instability (driven by a density inver- to explain the concept of isostatic equilibrium and de-
rive an expression for the pressure balance between
sion): the heavy and weak lower part of the lithospheric man-
two lithospheric columns in isostatic equilibrium.
tle is gravitationally unstable with respect to its surrounding.
• Use your understanding of the the gravitational po-
Depending on the rheology of the lithospheric mantle House-
tential energy of an object of mass m standing z meter
man and Molnar (1997, Geophys. J. Int. ,128,125–150.) esti- above ground to explain the concept of gravitational
mated that the part comprised in between the isotherm 900ºC potential energy of an entire lithospheric column.
and the isotherm 1300ºC is gravitationally unstable. The de-
tachment of this heavy keel from the rest of the lithosphere • Explain the concept of gravitational force acting be-
tween two lithospheric columns of contrasting poten-
dramatically changes the balance of forces. Convective thin-
tial energies.
ning results in a sudden increase in GPE stored in the thick-
ened crust, and therefore promotes extensional collapse. • Comment on the fact that isostatic equilibrium must
not be confused with mechanical equilibrium.

• Use the concept of gravitational force to explain the

formation of orogenic plateaux, and the concept of
gravitational collapse.

33
CHAPTER 3

Tectonic Forces
M.O.R.
Fg
Frp Oceanic
dFf Fo
dFr Continental lithosphere
Fox

100 KM
dFr
Fss dFaz dFa Fss
Mantle flow Mantle flow
Mantle flow

dFr

Ma
Ma

ntl
x dFp

ntl

e
flo
e
flo

w
w
Fsp
z

Ever since plate tectonic theory took hold in the late 1960s, geoscientists have argued over what
drives plates. Mantle upwelling at ridges that pushes plates apart, mantle circulation that drags plates
along, or slab pull have been, independently or in combination, proposed as driving forces for plate
motion. In this chapter, we explore the notion of tectonic forces.

34
FORCES DRIVING PLATES MOTION M.O.R.

Based on a simple model of lithosphere and man- Fg

Frp Oceanic
dFf Fo
tle interactions, geophysicists have found that dFr Continental lithosphere
Fox
three major tectonic forces — ridge push, slab pull

100 KM
dFr
and asthenospheric drive — can interact to explain
Fss dFaz dFa Fss
most observed plates motion at the Earth sur- Mantle flow Mantle flow

Mantle flow
face. The gravitational force acting between the
high-standing mid-oceanic-ridge and the distant, dFr

Ma
low-standing oceanic abyssal plain is called the

Ma

ntl
x dFp

ntl

e flo
ridge push force (Frp). It is of second order impor-

e
flo

w
w
tance when compared to the slab pull force (Fsp), Fsp
z
which is due to the pull exerted by the subduct-
ing slab onto the attached oceanic plate. Accord-
ing to Don Anderson, “... slabs drive tectonics.
There is no need for other driving mechanisms such as plumes or mantle convection that are independent of plate tectonics”. However, which
forces keep India moving northward since the detachment of the slab? The viscous drag imposed by the flowing asthenosphere
can drive plate tectonics. This force is called the asthenospheric-drive (Fad). The flow field in the asthenosphere is the product of con-
vective motion, as well as the flow induced by subducting slabs (slab suction force Fss a component of Fad), upper mantle avalanches
into the deeper mantle, and rising plumes (plume tectonics) pushing or pulling on the surrounding mantle. The
mantle flow can also drive the dynamic subsidence or uplift of the lithosphere above.

FORCES DRIVING PLATES DEFORMATION

FRP
phere
The dynamic of continental margins is under the influence of global plate tectonic Gravit
ational
Force Ocea
nic li
thos

forces, and local gravitational forces (Fg) acting between the high-standing con- FG

Asthenos
phere
here
lithosp
ental
tinental plate and the low-standing adjacent abyssal plain. The buoyant mantle Conti
n

Gravit
atio nal Fo
rce

wedge can also play an important dynamic role. The interplay between Frp , Fad, Fsp can Fg, FSP
FSS

can result in contrasting tectonic regimes. FSS

35
In the context of an active continental margin, the tectonic regime expe- M.O.R. Trench Plate

rienced by the upper plate (i.e. that standing above the subduction
Frp
Oceanic lith
zone) can vary from extensional tectonics, to contractional tectonics, osph
ere
Fg
Continental lithosphere

and to transcurrent tectonics. The variability of this tectonic regime ex-

presses the very changing balance of forces acting on the continental
margin.

MOR-Trench
Extensional tectonics is the tectonic regime which unfolds when the
Fsp
subduction trench moves away from the upper continental plate (a). a/ Trench-Plate

The trench retreats from the continent in response to slab rollback

M.O.R. Trench Plate
driven by the slab pull force, and/or under the push from the gravita-
Fg
tional force acting at the transition between of the oceanic plate and the Frp Oceanic lith
osph
ere Continental lithosphere
continent. This tectonic regime leads to the detachment of continental
fragments (or continental ribbons) from the main continental plate.

When the distance between the trench and a reference location inside

MOR-Trench
the continental plate decreases, the continental margin is shortened (b). Fsp

The tectonic regime is contractional , leading to the development of a

b/ Trench-Plate

cordillera orogen. Through times, the interplay between ridge-push,

M.O.R. Trench Plate
slab-pull, slab-suction and gravitational force results in a changing tec-
Fg
tonic regime with alternating periods of contraction, extension through Frp
Oceanic lith
osph
ere Continental lithosphere
gravitational collapse and transcurrent tectonic regime. These changes
can be triggered by small variations in the magnitude of plate velocity,
in the direction of plates motion, or asthenospheric flow. MOR-Trench

An active yet tectonically stable continental margin indicates that the Fsp
resulting effective force acting on the margin is very small and unable
Trench-Plate
to deform the margin (c). The margin is in mechanical equilibrium. c/

36
MANTLE FLOW: THE CASE OF STRONG COUPLING

Underneath the lithosphere, the flow of the asthenosphere exerts a

dynamic influence on the motion of overriding plates through the v
asthenospheric drive. The true extent of this influence is strongly de-
bated as it depends on the level of mechanical coupling between the
lithosphere and the convecting asthenosphere.

The sketches on the right illustrate two cases where the mantle ex-
erts a strong but opposite influence on the plates. The red arrows
show the velocity profile across the subducting lithosphere and into
the adjacent asthenospheric mantle.

A/ The subducting slab drags the surrounding asthenosphere into A

the subduction zone (this is the slab suction force). The astheno-
sphere resists and slow down the subducting slab. In this case, the
sub-lithospheric mantle flow is the result of plate tectonic proc-
esses, as the slab drags the nearby asthenospheric mantle into the
subduction zone. v

B/ The flow in the asthenosphere drives plate motion by dragging

the oceanic lithosphere into a downwelling zone (suction force
through shear traction). Here, it is the flow in the sub-lithospheric
mantle that powers plate-tectonic, as mantle convection drags the
slab into the subduction zone.

37
MANTLE FLOW: THE CASE OF WEAK COUPLING

If the strength of the asthenosphere/lithosphere coupling is not that

significant, the convective asthenosphere can still interfere with
plate motion. In particular, when the flow lines in the asthenosphere
are at an angle to the subducting slab, the “mantle wind” can force
the steepening or shallowing of the slab, hence promoting trench re-
treat (therefore extensional tectonics in the upper plate) or advance
(therefore contractional tectonics in the upper plate), as shown on
the sketches on the right.

C/ The asthenospheric wind under the subducting slab forces the

slab to bend downward, augmenting the slab-pull. This effect is
here enhanced by the downward asthenospheric flow above the C
slab forcing the steepening of the slab and promoting extensional
tectonics in the upper plate.

D/ In this case, the asthenospheric wind supports the subducting

slab, opposing the slab-pull. The slab is pushed upward forcing the
shallowing of the slab, promoting contractional tectonics in the up-
per plate.

Overall, the forces Frp , Fg, Fsp and Fad acting on the plates, the sub-
lithospheric flow, and the flow of the mantle in the vicinity of the
slab influence plate motion, and the tectonic regime at continental
plate margins.
D

38
RIDGE PUSH AND OTHER HORIZONTAL GRAVITA-
TIONAL FORCES

Gravitational force at mid-ocean ridge

Mid-oceanic-ridges (MOR) form a ~55,000 km long
mountain belt in the middle of the oceans. This moun-
tain belt stands ~2500 m above the average depth of oce-
anic abyssal plains (-5000 m), and the crest of the MOR
stands in average 2500 m below sea level. MOR, includ-
ing the region of hot mantle underneath, applies an hori-
zontal force directed from the MOR to the adjacent abys-
sal plains and adjacent continents. MORs are in iso-
static equilibrium, but not in mechanical equilibrium
since they are the locus of important and sustained ex-
tensional deformation, which accommodates the con-
tinuous production of new oceanic crust through decom-
pression partial melting of the asthenospheric mantle.

The map above shows the bathymetry of the north Atlantic ocean between North Africa and Newfoundland. The bathymetric pro-
file (along the thin red line) across the mid-oceanic ridge reveals the bathymetry across the MOR, with its crest at around 1000 m
depth, and the abyssal plains between -4500 and -5000 m.

Assuming temperature independent densities, an estimate of the total ridge-push per unit length parallel to the ridge axis is: "

(3 2)
L e
Frp = g⋅e⋅(ρm − ρw)⋅ +

where e is the elevation of the mid-oceanic ridge above the cooling lithosphere, ρm is the density of asthenosphere (3320 kg m-3), ρw
is the density of sea water (1030 kg m-3), L is the thickness of the oceanic plate (85 km), and g is the acceleration of the gravity field
(10 m s-2). With these values we get Frp = 2 × 1012 Nm-1. This estimate is an order of magnitude less than the slab-pull force. How-
ever, the value of the ridge-push force may increase up to 6.2 × 1012 Nm-1 when the ridge is underlain by a hot spot. A more thor-
ough formulation must integrate the temperature dependence of rocks density. Using the concept of gravitational force, the total

39
ridge push can be expressed as the integration with depth of
the difference in lithostatic pressure between a lithospheric col-
umn at the mid-oceanic ridge and a lithospheric column far
away from the mid-oceanic ridge. In other terms, the ridge
push corresponds to the difference between the gravitational
potential energy of the MOR and that of the adjacent abyssal
plain:

Fg = GPEMOR − GPEPLAIN
Continental margins are also regions with large contrasts in ele-
vation, in particular when these margins support mountain
belts. For instance, the South American Andes stands in aver-
age ca.8000 m above the adjacent abyssal plain of the east Pa-
cific ocean. The large contrast in gravitational potential energy
produces a gravitational force directed from the continent to-
wards the oceanic basin. The Andes are also in dynamic equilib-
rium supported by the convergence of the South American and
Nazca plates.

Gravitational force at continent-ocean transition

Even in the absence of cordillera, the transition between continents, standing in average a
few hundred meters above sea level, and adjacent abyssal plains are regions of significant
gravitational stresses. Not surprisingly, earthquakes tend to concentrate at plate bounda-
ries even when they involve passive margins. For instance, earthquakes tend to concen-
trate in a broad zone at the edge of the Australian continent. This zone shows a significant
gradient in gravitational potential energy and therefore a strong gravitational force.

40
SLAB PULL FORCE
M.O.R. Trench Plate

A simple formulation for the slab-pull per unit length parallel to

Frp zl
the trench is given by : zm. (zl. ∆ρ).g. This formulation assumes that Oceanic lith
osph
Fg
ere Continental lithosphere
the density of the plate and that of the asthenosphere are
temperature-independent. Assuming an average density contrast of
about 60 kg.m-3, a depth of the subducting slab zm = 660 km, a thick-
zm
ness of the slab zl = 100 km and g = 10 m s-2, we can calculate an abso-
lute maximum for the slab pull Fsp of ~4 × 1013 N m-1.
Fsp
Note: the value would be an order of magnitude greater (Fsp ~ 2 ×
1014 N m-1) for a slab going all the way down to the core-mantle Asthenosphere

boundary.

A more accurate formulation takes into account the temperature dependence of the density of rocks, the diffusion of heat, and the
velocity of the subducting slab. An estimate of the slab-pull force per unit length of subduction zone, Fsp(z), acting at depth Fsp and
caused by the density contrast between the cold oceanic plate and the mantle is given by:

8⋅g⋅α⋅ρm ⋅Ta ⋅L 2 ⋅Re

( 2⋅Re⋅L )
π 2 ⋅z π 2 ⋅d
Fsp = ⋅ Exp( − ) − Exp( − )
π 4 2⋅Re⋅L

ρm ⋅Cp ⋅ν⋅L
with: Re =
2⋅k
where Re is the thermal Reynolds number (the ratio of heat convection to heat conduction), z is the depth beneath the base of the
oceanic plate, α is the coefficient of thermal expansion (3 × 10-5 K-1), Ta is the temperature of the asthenosphere (1350ºC), k is the
thermal conductivity (2 W m-1 K-1), L is the thickness of the plate (85 km) and d is the depth of the upper mantle (d+L = 660 km), Cp
is the specific heat (1.17 × 103 J kg-1 K-1), and v is the rate at which the oceanic slab sinks into the mantle (~ 10 cm yr-1). Using the
above values we get Fsp = 2.5 × 1013 N m-1.

41
However, it is worth noting that Fsp linearly depends on the ill-defined coefficient of thermal expansion. With a coefficient of ther-
mal expansion varying in the range of 2×10-5 to 4×10-5 K-1 the slab-pull force varies between 1.7×1013 and 3.4×1013 N m-1. In addi-
tion, the olivine-spinel phase change, which occurs at around 350 to 420 km, increases the density in the subducting slab, provid-
ing an extra pull. Note that this extra pull depends on the thermodynamical characteristics of the phase transition, notably the
slope of its Clapeyron curve (dP/dT). The value of the slope of the Clapeyron curve ranges between 3 and 4 MPa K-1. The resulting
extra pull ranges between 1.2×1013 and 1.6×1013 N m-1.

The slab-pull force is resisted by the friction between the slab and the lithosphere, and the viscosity of the asthenospheric mantle.
This force is proportional to the velocity of the subducting plate and to the viscosity of the asthenosphere, that is also poorly con-
strained. Numerical calculations based on the differential equations for the flow of a viscous fluid suggest that the resistive force is
of the order of 1013 N m-1.

MAXIMUM EFFECTIVE TECTONIC FORCE

Overall, the maximum effective tectonic force available for lithospheric deformation is of the order of a few 1013 Nm-1. This force
must be able to double the thickness of the continental crust, to account for Tibet, the highest and largest plateau on Earth, and to
sustain the gravitational force associated with this thickening and with mountain belts at the edge of continents (e.g. the Andes).

Averaged over a lithospheric thickness of 100 to 200 km, the maximum differential stress available to deform the Earth litho-
spheres is in the range 100-500 MPa. This places some constraints on the rheology of Earth’s material.

re
sphe
ic litho
Ocean

b2/ lithosp
here Asthenosp
here
ental
Contin
a/

FRP
sphe
re b1/
e ic litho
l Forc Ocean FG
tiona FRP
Gravita re
sphe
FG sphere ic litho
FRP ental litho Ocean
here Contin
e sp rce
here Asthenospher ic litho tional Fo
lithosp Ocean Gravita here
ental Asthenosp
Contin l Forc
e FG
tiona
Gravita
e
here henospher
lithosp Ast
FSS ental
Contin
FSP e
l Forc
tiona
Gravita FSS
FSS

FSS FSS

FSS

42
# Load the marmap library
require(marmap)

#setwd stands for set working directory

setwd('/Users/Patrice/Documents/GIS/')

# Grab data from National Oceanographic and Atmospheric Administration (National Geophysical Data Center)
# More data available here: https://www.bodc.ac.uk/data/online_delivery/gebco/ APPENDIX: This R
map_Altiplano<-getNOAA.bathy(lon1=-60, lon2=-80, lat1=-30.00, lat2=-10, resolution=2)
script grabs bathymetry
# Hypsometry of the region and topography data
# zbreaks = quantile(map_Altiplano, seq(0, 1, length.out=256))
zbreaks = zbreaks <- seq(-8000, 8000, by=100) from the US National
plot(zbreaks) Oceanographic and At-
# Make colour palette mospheric Administra-
ocean.pal <- colorRampPalette(c("#000000", "#000413", "#000728", "#002650", "#005E8C", "#0096C8", "#45BCBB", "#8AE2AE",
"#BCF8B9", "#DBFBDC")) tion (NOAA), creates a
land.pal <- colorRampPalette(c("#336600", "#F3CA89", "#D9A627", "#A49019", "#9F7B0D", "#996600", "#B27676", "#C2B0B0", topographic map, and
"#E5E5E5", "#FFFFFF"))
#land.pal <- colorRampPalette(c("#467832", "#887438", "#B19D48", "#DBC758", "#FAE769", "#FAEB7E", "#FCED93", "#FCF1A7", plot a topographic profile
"#FCF6C1", "#FDFAE0"))
mypalette <-c(ocean.pal(sum(zbreaks<=0)-1), land.pal(sum(zbreaks>0))) (see p.39 and 40).
# Plot the map
plot(map_Altiplano, image=T, bpal=mypalette, land=T, deep=c(-6000, -3000, 0), shallow=c(-3000, 500, 0), step=c(1000, 500, 0),
col=c("grey10", "grey40", "black"), lty=c(1, 1, 1), drawlabel=c(T, T, T), asp=0, ylim=c(-30,-10), xlim=c(-80,-60))

# Topographic profile
points(c(-78, -61),c(-20, -20), type="o", col=2)
profile <- get.transect(map_Altiplano, x1 = -78, y1 = -20, x2 = -61, y2 = -20, locator=F, distance = TRUE)
plotProfile(profile, ylim=c(-9000,5000), xlim=c(0,1800))

# The following organize map, legend and profile into one single image
def.par <- par(no.readonly = TRUE)
nf<-layout(matrix(1:4,nc=2), height=c(3,1), width=c(5,1))
par(mar=c(4,4,1,1))
layout.show(nf)

plot(map_Altiplano, image=T, bpal=mypalette, land=T, deep=c(-6000, -3000, 0), shallow=c(-3000, 500, 0), step=c(1000, 500, 0),
col=c("grey10", "grey40", "black"), lty=c(1, 1, 1), drawlabel=c(T, T, T), asp=1, ylim=c(-30,-10), xlim=c(-80,-60))
points(c(-78, -61),c(-20, -20), type="o", col=2)
plotProfile(profile, ylim=c(-9000,5000), xlim=c(0,1800))
image(x=0, y=zbreaks, z=matrix(zbreaks, 1, length(zbreaks)), col=mypalette, breaks=zbreaks, useRaster=FALSE, xlab="", ylab="",
axes=FALSE)
axis(2, at=seq(-9000, 7000, 1000), las=2)
mtext("[meters]", side=2, line=3)
box()

43
 Expected learning outcomes

At this stage you should be able to:

• List and describe the three types of tectonic forces

acting on the Earth’s lithosphere, and tell which are
plate-boundary forces?

• Calculate the value of the ridge-push force from

knowledge of the principle of isostasy and the con-
cept of gravitational force.

• Explain the role of the strength of the mechanical

coupling between the convective mantle and the litho-
sphere on the dynamics of active continental margin.

• Comment on the changing dynamics (i.e. tectonic

regime) of active continental margins.

• Comment on the maximum value of the effective tec-

tonic force, and that of individual tectonic forces.

44
CHAPTER 4

Introduction to Rheology

Rheology is the study of flow, the mechanical response of material to applied deviatoric stresses. “Flow” is
used here in its broader meaning, which includes both viscous flow and frictional flow. The relationships
between, on one hand applied deviatoric stress and resulting strain, and on the other hand deviatoric stress
and strain rate characterize the macroscopic mechanical behavior of rocks. These relationships lead to the
constitutive equations linking deviatoric stress and strain rate; "constitutive" as they depend on the
constitution of the material.

45
SECTION 1

Elastic, Plastic and Viscous Flows

In this section

1. Brittle vs ductile Elastic Viscoelastic Viscoplastic

deformation
σ ε = σ (1 - e-t/τ)
2. Elastic, plastic, and ε=
viscous flow curves E E
3. Rheology of
polycrystalline rocks
4. Steady-state flow laws
5. Sensitivity of flow laws

Primary Creep
Secondary Creep

Rocks display a very large range of behavior which permanent strain accumulates. Rocks can
when submitted to deviatoric stresses. Over mil- deform elastically since they transmit mechanical
lions of years, hot rocks can flow under small de- waves (sound waves, seismic waves). The rheol-
viatoric stress, like a very low viscosity fluid do. ogy of rocks is therefore complex and covers a
Rocks can break under a sudden deviatoric stress broad range of mechanical behavior from elastic,
load (e.g. hammer). Rocks can resist moderate plastic and viscous.
deviatoric stress, until it reaches a threshold from

46
BRITTLE VS DUCTILE DEFORMATION: The role of temperature, composition and strain rate...

Objects deform as a continuum (i.e. continuous deformation) via elastic deformation or via viscous or plastic flow, when the flow unit are
atoms moving in the lattice of crystal. In contrast, when the flow unit are fragments of crystals or fragments of rocks moving relatively to
each other because of fractures or faults, deformation is said to be brittle (i.e discontinuous deformation via frictional flow).

Discontinuous deformation:
• Case 1: The media is pre-fractured, its strength depends on the frictional strength of pre-existing fractures and faults.
• Case 2: The media has no pre-existing fractures, its strength depends on its constitution in particular on the cohesion of its grains.
This cohesion must be overcome before fractures and faults can accommodate frictional sliding.

Continuous deformation:
• Elastic deformation
• Plastic deformation
• Viscous deformation

Ductile deformation. Photograph P. Rey Brittle deformation. Photograph Pui-Leng, flickr

47
CONTINUUM MODELS OF ROCK’S MECHANICAL BEHAVIOR

Each type of flow has a characteristic flow curve in the deviatoric stress vs strain space. The graph below shows the flow curves for three ele-
mentary deformation mechanisms, with for each their mechanical analogue:

• Linear elastic flow (in blue): For small stresses, most mate-
rial are elastic. The characteristics of elastic flow are: 1/ Elas- σ
tic flow occurs as soon as stress is applied. 2/ Strain increases
σ
Ideal Plastic
as long as stress keeps increasing. The elastic flow curve is lin-
ε
ear, its slope is 1/E with E the Young modulus, a physical σ*
ε= . f(t)
properties characteristics of elastic flow. 3/ Strain does not Elastic E
accumulate if the stress is maintained constant. 4/The mate- A Plastic B
rial recovers its original shape when stress is relaxed. σ*
The mechanical analogue for elastic deformation is a spring. σ
ε= Viscous

Deviatoric stress
E
• Plastic flow (in green): In most material elastic flow is lim- ε= σ . t
η
ited to a certain level of stress beyond which the flow
switches from elastic to plastic. This limit is called the yield
stress (σ ∗). The characteristics of plastic flow are: 1/ Plastic D
flow occurs as soon as the yield stress is reached (i.e. from
point A onwards). 2/ Plastic strain accumulates at a level of 1/E C
stress which is constant and equal to the yield stress. Hence, Plastic Elastic
the amount of strain depends on the duration t over which Cumulative strain (ε )
stress is applied. 3/ The plastic component of strain is perma-
nent. When the imposed deviatoric stress is removed (B>C),
the material recovers the component of elastic deformation only. 4/ For and ideal plastic material there is no elastic component.
The mechanical analogue of plastic deformation is a rigid block sliding on a rough planar surface.

• Linear viscous behavior (in purple): Linear viscous material (newtonian viscosity) have no yield stress. The flow curve is characterized
by a linear relationship between stress and strain. Strain accumulates at varying and also constant stress (dashed purple line from D) and
when the stress is removed the flow stops but the material does not return to its undeformed state.
The mechanical analogue of viscous deformation is a dashpot.

48
In a graph strain vs time, in experiments where the deviatoric stress is applied at t = 0 and maintained constant for t > 0, elastic, plastic and
viscous flows are characterized by three contrasting relationships.

Elastic strain: When the deviatoric stress is imposed, the material records a fi-
nite amount of elastic strain. This amount of elastic strain remains constant σ = Constant
σ Ideal Plastic
through time (and equal to ) as long as the deviatoric stress is maintained.
E σ .t
ε = σ*E-1 + α . t ε=
Plastic strain: The deviatoric stress is constant and equal to the yield stress. At η

Cumulative strain (ε )
σ∗ Viscous
t = 0, a component of elastic strain ( ) is instantaneously recorded before the
E
onset of accumulation of plastic strain. Plastic strain accumulates as long as
Elastic
the deviatoric stress is maintained. For ideal plastic flow, the rate of accumula-
σ
tion (α) is constant. ε=
E
Viscous strain: Viscous strain accumulates even for a very small deviatoric
σ Time
stress. This strain accumulate at a constant rate ( ), which depends of the
η
shear viscosity (η) of the material. Viscous strain accumulates as long as the de-
viatoric stress is maintained.

σ Elastic Plastic
RHEOLOGY OF POLYCRYSTALLINE ROCKS Hardening Creep
C
B
Crystalline rocks display a mechanical behavior that incorporates the three ele- σ*' C''
Failure
mentary flows; rocks are elasto-visco-plastic material. This graph shows a char-
acteristic flow curve for a polycrystalline material. At stresses below the yield σ*
A
stress (σ ∗), polycrystalline material behave elastically (blue curve). Above the

Deviatoric stress
yield stress, the material behave plastically. At low level of strain (green curve
from A to B) the material becomes stronger (hardening) as the applied stress
must increases in order to keep the material deforming (i.e. the yield stress
must be exceeded for strain to accumulate). At high level of strain (red line be- E
yond B) the material flows under a constant stress. Under the hardening plastic C' Cumulative strain (ε )

49
regime, the removal of the driving stress leads to the removal of the elastic component of the deformation (curve CC'). If the sample is re-
stressed, elastic deformation occurs under an extended domain as the yield stress has increased (curve C'C''). This means that the material
has become stronger.

To understand how hardening occurs one must understand how strain is achieved at the microscopic scale. Plastic strain affects the ar-
rangement of atoms in the lattice of minerals. During strain, this arrangement is perturbed by the introduction of defaults or gaps in the lat-
tice called dislocations. Each dislocation introduces a local elastic strain. Under stress, these defaults move around leading to permanent de-
formation. In the hardening plastic regime, the density of defaults increases but low stress levels impede their displacement, hence the elas-
tic energy in the crystal lattice increases making the material stronger. The animations show how dislocation moves about:

Dislocation wall
Motion of dislocations is driven by internal elastic
energy. Dislocations organize themselves into ar-
rays such as dislocation walls to minimize internal
elastic energy. This results in the multiplication of
sub-grains and therefore permanent plastic strain.

sub-grain

Movie 4.1 Dislocation glide Movie 4.2 Dislocation climb Movie 4.3 Dislocation glide and climb

50
The viscous flow of polycrystalline material can be illustrated on a
space strain vs time for experiment performed at constant stress. In Tertiary creep
this graph, the plastic flow curve can be divided into three regimes σ = Constant

Cumulative strain (ε )
Secondary creep A'
called: Primary, Secondary and Tertiary creep.
Primary ra i n Rate Failure
nt S t
creep ta
Primary creep corresponds to a reversible flow for which elastic defor- Cons
A
mation is instantaneously removed following unloading of the sam-

Elastic
ple (at time t1), whereas another component of strain called viscoelas-
tic deformation is also recovered but over of a time window t1 - t2.

Elastic

V isc
Secondary creep is characterized by a linear relationship between σ*
la

oe
strain and time implying that the material is deforming at constant stic

Elastic

Plastic
Vis
strain rate. Upon unloading elastic deformation is instantaneously re- oe

c
la s
tic
covered whereas the viscoelastic deformation is recovered over a pe-
riod of time (t'1 - t'2). The sample, however, records a permanent plas- t1 t2 t'1 t'2 Time
tic deformation. It is within the secondary creep regime that the creep
parameters that govern the rheology of rocks in the ductile regime are determined. The tertiary creep corresponds to the development of a
mechanical instability in which an increase in strain rate leads to the mechanical failure of the stressed sample. Real rocks display a com-
plex elasto-visco-plastic behavior. This mechanical behavior can be represented by a
combination of springs (elastic component), blocks on a rough surface (plastic com-
Elastic Viscoelastic Viscoplastic
ponent), and dashpots (viscous component) connected in series or in parallel to fit a
σ ε = σ (1 - e-t/τ)
real flow curve. ε=
E E
The rock analogue on the right has the same flow curve than than presented above.
At depth > 10 to 20 km, the deformation of rocks is characterized by slow steady-
state creep at constant strain rate, that can accommodates large amounts of ductile
deformation. It is generally assumed that steady-state constitutive equations of
rocks, or flow laws, can be used to characterize the large-strain high-temperature Primary Creep
ductile deformation that occurs in the Earth. Secondary Creep

51
SECTION 2

Strength envelopes
In this section

Yield
1. Low-Moderate stress (MPa)
stress Yield stress (MPa) Yield stress (MPa) Yield stress (MPa)
2. High-stress3000 2000 1000
regime 0 1000 2000 1000 0 1000 2000 1000 0 1000 2000 1000 0 1000
3. Brittle regime
4. Sensitivity of flow curves 20 20 20 20
5. Rheological profiles
Moho
6. Sensitivity of profiles 40 40 40 40

60 60 60 60

80 80 80 80
789MPa 255MPa 598MPa 193MPa 584MPa 213MPa 582MPa 212MPa
[12] 100 [11] 100 [16] 100 [1] 100

The notion of strength envelope is relative to sure, the accumulated strain and the tectonic
the evolution with depth of the deviatoric regime. At higher temperature where viscous
stress necessary either to exceed the yield creep dominates, the strength envelope is sensi-
stress, or to achieved a pre-defined strain rate tive to rock types, strain rate, and the magni-
through viscous flow. At low temperatures, tude of the deviatoric stress. Here we start by
where plastic strain dominates, the strength en- giving the relationship deviatoric stress vs tem-
velope is mainly a function of pore fluid pres- perature (i.e. depth).

52
FLOW LAWS FOR STEADY-STATE CREEP: LOW TO MODERATE-STRESS REGIME

The constitutive equation that accounts for most low- to moderate stress steady-state strain is the power-law creep equation, so
called because the steady-state strain rate is proportional to the differential stress raised to a power n. The following equations
•
give the strain-rate ( ϵ ) as a function of the differential stress, and the differential stress as a function of the strain rate:

• 1/n

(R.T) (A) ( n⋅R⋅T )

• −E ϵ E
ϵ = A⋅σ n ⋅Exp rearranging we get >> σ= ⋅Exp

where A is a constant (Pa-n s-1), n is the stress exponent that characterizes the sensitivity of strain rate on the differential stress (n is
dimensionless), E is the activation energy per mole for the creep process (J mol-1), it is the energy barrier that inhibits the creep
mechanism, R is the Boltzmann constant (8.3144 J mol-1 K-1), and T is the temperature (K). The constants A, E and n are characteris-
tic of material. The power law creep shows that both temperature and differential stress have a large effect on the strain rate. Thus
an increase in temperature increases the strain rate for a constant stress, or lowers the
stress required to produce a given strain rate.
n=1
This effect is accounted for by the rapid increase, with increasing temperature, of the expo- n=3

σ1−σ3
nential term in the equations above. The graph shows that as n increases from 1 to large n=

8
values, power-law material evolves from viscous material (for n=1 the power law becomes
ε
a linear relationship between stress and strain rate) to near ideal plastic material. For most
rocks at moderate stress level 2<n <5, whereas at low level stress, 1<n<2.

FLOW LAWS FOR STEADY-STATE CREEP: HIGH-STRESS REGIME

•

(ϵ)
R⋅T ϵ
Power-law creep implies that at 500ºC olivine would only deform at unrealistically high
σ = σd ⋅ 1 − ⋅ln •d
stresses. A better description of the behavior of olivine at high-stress regime (> 200MPa) is Ed
given by the Dorn's law, a relationship that is not as temperature dependent as the power
law, in which Qd is an activation energy, σd a critical stress that must be exceeded, and εd
is the critical strain rate.

53
FLOW CURVES IN THE BRITTLE REGIME

With pre-existing fractures: At low temperature or at high strain rate or under high pore-pressure, but mainly in the upper crust and
the upper mantle, the failure mechanism is modeled as frictional sliding:

σ = β⋅(1 − λ)⋅(ρ⋅g⋅z)
where g is the gravitational acceleration, λ is the ratio of fluid pore pressure to the normal stress, ρ the density. β is a parameter de-
pendent on the type of faulting given by:

( )
−2
R−1
β= with: R= 1 + μ2 − μ
1+(σ ) ⋅(R − 1)
σzz − σ3
1 − σ3

and μ the coefficient of internal friction that characterizes the roughness of the fracture plane.

Because the coefficient of internal friction varies little around 0.75 (pretty much independent on the rock composition), we get that
R=4, and therefore β varies continuously from 0.75 for normal dip slip faults, to 1.2 for strike-slip fault, to 3 for reverse dip slip
faults depending of the exact value of the principal stresses ratio. The frictional sliding equation shows that it takes less differen-
tial stress to achieve brittle failure under extensional stress regime (maximum principal stress vertical) than under a compressional
stress regime (maximum principal stress horizontal). This makes sense since in a compressional stress regime gravity acts against
reverse faulting but enhances the principal stress in extension.

Without pre-existing fractures: In this case, the strength of rocks includes the cohesion (C) holding the grains together at atmospheric
pressure. The failure mechanism is modeled by the Coulomb's criteria linking the shear stress ( τ ) acting on a fault plane and the
normal stress ( σn ) acting perpendicular to it: τ = C + μ⋅σn

The normal stress and the shear stress are related to the orientation of the fault plane with respect to the principal stress axes (
σ1, σ2 and σ3):
σ1 + σ3 σ1 − σ3 σ1 − σ3
σn = − ⋅cos(2α) and τ= ⋅sin(2α) with α the angle between the fault plane and σ1.
2 2 2

54
SENSITIVITY OF FLOW CURVES TO THERMODYNAMIC PARAMETERS AND CHEMICAL ENVIRONMENT

Confining Pressure and Temperature

Confining pressure prevents rocks from falling apart. It is therefore not surpris-
250 Wombeyan Marble
ing that increasing confining pressure increases the amount of deformation a σ3=100 MPa
sample can accumulate before failure. For instance, in the example of the Wom- 200
beyan marble brittle failure is either reached for a larger amount of accumu-
150 σ3=35 MPa
lated strain or even impeded when a higher confining pressure is applied. This
higher confining pressure allows the sample to sustain a larger differential 100
stress. Temperature enhances the ductility of material. For instance, in the exam-

σ1−σ3 (MPa)
σ3=10 MPa
ple of the Solenhofen limestone, raising temperature decreases the magnitude 50
of the deviatoric stress that can be supported, and increases the amount of accu- Failure
Failure
mulated strain before brittle failure . Under high temperature material can there- ε
1 2 3
fore accumulate more strain but they can do so for a smaller amount of differen- Cumulative strain (ε )
tial stress.

Solenhofen Limestone
400 300ºC

300 Failure
500ºC
Failure
200 600ºC

Failure

σ1−σ3 (MPa)
100

Pc~40MPa

2 4 6 8 10 ε
Cumulative strain (ε )

Paterson and Wong, 2006. Experimental Rock Deformation – The Brittle Field. Surveys
in Geophysics, 27, 4, pp 487-488.

55
H2O and Pore Fluid Pressure
Water in rocks acts as a softening agent. Indeed, compared to a dry sample, a wet sample deforms at lower deviatoric stress (bot-
tom left figure). The presence of fluids enhance mechanism of deformation controlled by coupled dissolution (at high stress re-
gion) and precipitation (in low stress regions). Pore fluid pressure acts against the confining pressure and therefore promotes fail-
ure at lower differential stress and lower strain (figure on the right). This process is called “hydraulic fracturing”, this is the proc-
ess explaining micro-seismicity following dam water loading.

Yule Marble Dry

300 600 Solenhofen Limestone
0.5% H20

60 MPa
200 400

2.9% H20
Failure
σ1−σ3 (MPa)

100 200 Failure 65 MPa

σ1−σ3 (MPa)
Failure
93 MPa 78 MPa
T=300ºC
Pc~500MPa

4 8 12 18 ε 3 6 9 12 ε
Cumulative strain (ε ) Cumulative strain (ε )

σ 1
@ Mervin Paterson, ANU

σ 3 σ 3

σ 1

Fluid assisted diffusive mass-transfer through the lattice (Nabarro-

Herring creep) or along the grain boundaries (Cobble creep)

56
RHEOLOGICAL PROFILE OF THE CONTINENTAL LITHOSPHERE

It is time now to synthesize what we have learn about the rheology of crystalline rocks to determine the strength profile (i.e. en-
velop) of the continental lithosphere.

The strength of the lithosphere can be defined as the vertical integration, from the top to the bottom of the lithosphere, of the differ-
ential stress required to trigger either brittle failure or the flow failure of rocks. Failure can occur by power-law creep or Dorn law
creep at high temperature and low strain-rate, or by frictional sliding at low temperature and high strain rate, in which case the dif-
ferential stress depends of the tectonic regime. However at any given depth, the failure mechanism is the one that requires the
minimum differential stress to operate.

The strength of rocks at any depth in the crust is the lowest differential stress to reach either brittle or flow failure:

• 1/nc

( Ac ) ( nc ⋅R⋅T(z) )
ϵ Ec
σ1 − σ3 = β⋅(1 − λ)⋅(ρc ⋅g⋅z) or σ1 − σ3 = ⋅Exp

The strength of rocks at any depth in the lithospheric mantle is the lowest differential stress to reach either brittle, or power law, or
Dorn law failure:

σ1 − σ3 = β⋅(1 − λ)⋅(ρm ⋅g⋅z)

• 1/n

( Am ) ( )
ϵ Em
or σ1 − σ3 = ⋅Exp for (σ1 − σ3) < 200MPa
nm ⋅R⋅T(z)

(ϵ)
R⋅T(z) ϵd
or σ1 − σ3 = σd ⋅ 1 − ⋅ln • for (σ1 − σ3) > 200MPa
Ed

57
The figure below shows the rheological profiles for a compressional and an extensional tectonic regime. The rheological parame-
ters for the crust and the mantle are that of a granite (quartz > 40%) and a dunite (peridotite with > 80% olivine) respectively. The
•
power law and Dorn creep law are dependent on both the temperature T and the strain rate ϵ . Therefore, these parameters must
be defined first. The geotherm here is such that the temperature at the Moho is 460ºC and the chosen strain rate is 10-15 s-1.

The straight parts of the profiles (in the upper crust and the upper mantle) represent brittle failures. The red lines are the Dorn law
creep curves, whereas the black curved lines are the power law flow curves. One can see that the differential stress for both the
power law and Dorn law creep strongly decrease with T.

Strength in compression σ1-σ3 (MPa) 1000 0 1000

Crust
β=3 Ac = 5x10-6 MPa-ns-1 Compression Extension
nc = 3.0 20
Strength in extension
Ec = 190 kJ.mol-1 Moho
β=0.75 40

60
Mantle

Depth (km)
Am = 7x104 MPa-ns-1 Ed=540 kJ.mol-1 Dorn Law + 80
Power law
nm = 3.0 Stress Treshold=8500 MPa
Em = 520 kJ.mol-1 ed=3.05x1011s-1 100
TMoho=460ºC

A few things to note:

• With the Dorn law, brittle failure does not occur in compression, which significantly reduces the strength of the upper mantle.
• The brittle part of the crust is thicker in extension that it is in compression.
• For a normal geotherm (TMoho <650ºC), the upper crust and the upper mantle are the strongest layers of the lithosphere, the
lower crust and the lower lithospheric mantle are comparatively much weaker.
• Because quartz deforms by ductile flow at lower temperature (~300ºC) than olivine (ductile at ~600ºC) the upper mantle is much
stronger than the lower crust for a normal geotherm (TMoho <650ºC).

58
SENSITIVITY OF STRENGTH ENVELOPS TO TEMPERATURE
The dependence of the strength envelops on temperature is documented below. From left to right the temperature at the Moho in-
creases from 400ºC to 700ºC. The strength of the lithosphere, in both compression and extension (the surface area of the dark and
pale blue regions respectively) has been averaged over the lithospheric thickness assuming for the mantle either a power law
creep or a combination of power law and Dorn law. As temperature increases, the averaged strength of the lithosphere decreases.

σ1-σ3 (MPa) 2000 1000 0 1000 2000 1000 0 1000 2000 1000 0 1000 2000 1000 0 1000 1000 0 1000

20 20 20 20 20

Moho
40 40 40 40 40

60 60 60 60 60
Power law 737 MPa 237 MPa 570 MPa 195 MPa 410 MPa 151 MPa 236 MPa 108 MPa 106 MPa 61 MPa
80 80 80 80 80
Dorn law+ 438 MPa 221 MPa 333 MPa 183 MPa 239 MPa 142 MPa 158 MPa 99 MPa 101 MPa 57 MPa
Power law -40% 100 -7% -41% 100 -6% -41% 100 -6% -33% 100 -8% -5% 100 -7%

TMoho=400ºC TMoho=460ºC TMoho=540ºC TMoho=620ºC TMoho=700ºC

The integrated strength of the continental lithosphere (in Nm-1) is defined the integration
500ºC 700ºC
over depth of the rheological profile. Since the frictional sliding depends on the tectonics re- 40

Integrated Strength x 1012 N.m-1

gime, three "strengths" can be defined for compressional (Fedt), transcurrent (Feds), and ex- 30
20
tensional (Fedc) tectonic regime:
10
The figure on the right shows the integrated strength as a function of temperature at the 0
TMoho (ºC)

400 600 800 1000

Moho. For TMoho ~500ºC the upper mantle is the stronger layer of the lithosphere (in blue -10
Strength
-20
in the inset). Because the power law creep is exponentially dependent on temperature a
-30
small increase in temperature strongly reduces the integrated strength. Indeed, when -40

Depth
700ºC
TMoho is close to 700ºC the strength of the upper mantle drastically decreases. Over 700ºC, -50
500ºC
the stronger layer of the lithosphere is the brittle upper crust the rheology of which has no
dependence on temperature. The red and the blue curves are the integrated strength in ex-
tension and compression respectively.

59
SENSITIVITY OF STRENGTH ENVELOPS TO LITHOLOGIES

The continental crust displays a wide Yield stress (MPa) Yield stress (MPa) Yield stress (MPa) Yield stress (MPa)
3000 2000 1000 0 1000 2000 1000 0 1000 2000 1000 0 1000 2000 1000 0 1000 Upper crust
range of compositional variation from
mafic rocks (gabbros, amphibolites) to 20 20 20 20
Lower crust
quartz dominated composition (gran- Moho
40 40 40 40
ites). This contrast with the relative ho-
60 60 60 60 Moho
mogeneity the mantle (peridotite with
Mantle
60%<olivine<95%). The figure on the 80
789MPa 255MPa
80
598MPa 193MPa
80
584MPa 213MPa 582MPa 212MPa
80

right illustrates the dependence of the [12] 100 [11] 100 [16] 100 [1] 100
Strength in
rheology of the lithosphere on its litho- 3000 2000 1000 0 1000 2000 1000 0 1000 2000 1000 0 1000 2000 1000 0 1000
compression
logical composition. All rheological pro- 20 20 20 20

files involve the same geotherm with a Moho Strength in

40 40 40 40 extension
temperature at the Moho of 540ºC (thin
60 60 60 60
dashed line TMoho = 600ºC).
Depending of the composition and the 80 80 80 80 TMoho=540ºC
460MPa 178MPa 433MPa 166 MPa 425MPa 155 MPa 406MPa 163MPa
rheological parameters the integrated [6] 100 [4] 100 [14-15] 100 [13] 100 TMoho=600ºC
strength in contraction of the continental 3000 2000 1000 0 1000 2000 1000 0 1000 2000 1000 0 1000 2000 1000 0 1000
lithosphere varies over one order of mag-
nitude from 790 MPa down to 78 MPa. 20 20 20 20 20

Moho
In extension the integrated strength 40 40 40 40 40

ranges from 255 MPa down to 50 MPa. 60 60 60 60 60

The large variability of the rheological 80 80 80 80 80

properties of common crustal and man- 404MPa 160 MPa 379MPa 149MPa 365MPa 153MPa 355MPa 147MPa 78MPa 49MPa
[5] 100 [7] 100 [9] 100 [10] 100 [2wet] 100
tellic rocks, prevents the definition of a
“standard” rheology for the continental
lithosphere.

60
#===============================================================
# This R script calculates and plot the continental geotherm and rheological profiles
#===============================================================
# ====== Geotherm =========
library(mosaic)
library(mosaicCalc)
# SI unit: m, sec, kg

# Radiogenic heat production (W.m-3)

H<-0.7*10^-6
# Thermal conductivity (diffusivity x density x heat capacity), (W.m-1.K-1)
K<-(10^-6)*2650*1000
#Mantle heat flow (W.m-2)
Qm<-0.020
#Surface Temperature (K)
To<-293
#Thickness of the continental crust (m)
zc<-42000
#Temperature defining the base of the lithosphere (K)
Tl<-1603

#Crustal geotherm
Temp_crust<-function(z){(-H/(2*K))*z^2+(((Qm+H*zc)/K))*z+To}
#Geotherm in the lithospheric mantle
Temp_lithos_mantle<-function(z){(Qm/K)*(z-zc)+Temp_crust(zc)}
#Depth of the base of the lithosphere
zl<-findZeros(Temp_lithos_mantle(z) - Tl ~ z, z.lim=range(100000, 200000))
#Geotherm in the asthenosphere (called the adiabat because it follows the adiabatic gradient)
Temp_asthenos<-function(z){Tl+0.0003*(z-zl)} The rheological strength profiles above were
#Geotherm in the first top 200 km … derived using three different reference
Geotherm<-function(z){ifelse(z<0, To, ifelse(z<=zc, Temp_crust(z), ifelse(z<=zl, Temp_l- strain rates, and the R script on the left.
ithos_mantle(z), Tl)))}
#Plot of the geotherm above
plotFun(Geotherm(z) ~ z, z.lim=range(0, 200000))

# ====== Rheology =========

61
 Expected learning outcomes

At this stage you should be able to:

• Explain the concepts of continuous and discontinuous deformation, and how

elastic, plastic and viscous deformation fit in that framework.

• Explain the concepts of elastic, plastic and viscous deformation, and the notion
of yield stress.

• Using appropriate mechanical analogues, explain the concepts of visco-elastic,

visco-plastic and visco-elasto-plastic deformation.

• Document in a space stress vs strain the meaning of elastic, plastic and viscous
strain.

• Document in a space strain vs time the meaning of elastic, plastic and viscous
strain.

• Document in a space stress vs strain the rheology of polycrystalline rocks, and

use the concept of dislocation.

• Comment on the flow law in the moderate and high stress regimes, and in the
brittle regime.
n.b.: We do not expect you to know these laws by heart, just to comment on their vari-
ous terms and parameters.

• Comment on the effect of pressure, temperature and fluids on the rheology of

rocks.

• Combine the various flow laws to graph a strength profile for the continental litho-
sphere, and comment and the rheological layering of the lithosphere.

62
CHAPTER 5
Computational Tectonics & Geodynamics

Computational tectonics/geodynamics provides a robust platform

for hypotheses testing and exploring coupled tectonic and geody-
namic processes. It delivers the unexpected by revealing previously
unknown behaviors.

63
SECTION 1

The rise of computer science

During the mid-19th century, the laws of thermodynam- tonic processes which also involves the flow - brittle or
ics - which describe the relation between heat and forces ductile - of rocks. 21st century computers have now the
between contiguous bodies - and fluid dynamics - which capability to compute the flow of fluid with complex,
describes the flow of fluids in relation to pressure gradi- temperature-dependent rheologies, and contrasting physi-
ents - reached maturity. Both theories underpin our un- cal properties we typically encounter on Earth. For geosci-
derstanding of many natural processes from atmospheric ences, this is a significant shift.
and oceanic circulation, mantle convection, and plate tec-

64
Computational tectonics and geodynamics: There is a revolution Navier-Stokes equations: Written in Georges G. Stokes (1819-1903)
unfolding at the moment in all corners of science, engineering, and the 1850’s, the Navier-Stokes equations
other disciplines like economy and social sciences. This revolution are the set of differential equations that
is driven by the growing availability of increasingly powerful describe the motion of fluid and associ-
high-performance computers (HPC), and the growing availability ated heat exchanges. They relate veloc-
of free global datasets and open source data processing tools (Py- ity gradients to pressure gradients and
thon, R, Paraview, QGIS-GRASS, etc). Here in Australia, supercom- can be analytically solved only when
puters such as Magnus (36,000 cores) at the Pawsey Supercomput- considering steady laminar flows.
ing Center in Perth, and Raijin (58,000 cores) at NCI (National Com-
putational Infrastructure) in Canberra, have enabled and democra-
tized the art of numerical modeling.

Broadly speaking, numerical modeling is a discipline which en-

ables the exploration of the behavior of complex systems. It also
gives scientists, engineers and economists the capacity to extract
new knowledge and understanding from big data. In the context of
geology and geophysics for instance, numerical modeling allows
us to build a model of lithosphere and explore how this lithosphere
behaves when submitted to tectonic forces. It enables geoscientists
to build spherical models of the Earth to explore mantle convec-
tion. These models can account for a very broad range of petro-
physical properties including radiogenic heat, heat diffusivity, den-
sity, heat capacity, rheology (brittle and ductile), solidus and liqui-
dus, etc. The behavior of the lithosphere and that of the convective
mantle is governed by the laws of thermodynamics which de-
scribes exchange of energies within the system, and fluid dynamics
which relates deformation (i.e. flow) to pressure gradients. Both
thermodynamics and fluid dynamics are fully developed theories
from 19th century physics. Yet, it is only over the past decade that
HPC became powerful enough to efficiently solve in 3D the rele-
vant equations at a reasonable spatial and temporal resolution.

65
Let’s try to make sense of all this. Newton’s second law of motion one component must be gained by another component. Hence,
states that the sum of all forces applied to an object must be equal heat exchanges are described via the equation of conservation of
to its mass (m) times its acceleration (a): m⋅a = ∑F energy.

This law applies to any object including a parcel of fluid. The right Remembering the 1D conduction-advection heat transfer equation
hand side of the equation includes the force of gravity, the viscous (see Chapter 1)...
forces due to friction with the surrounding fluid, pressure gradient
dT K d 2T A dT
from the surrounding fluid, surface tension, etc. The left hand side = ⋅ 2 + − uz ⋅
dt ρ ⋅ Cp dz ρ ⋅ Cp dz
has by necessity the dimension of a force, it is the inertia force. The
acceleration a can be splitted into a local acceleration (also called ... re-arrange so it looks similar to that given on page 65 ...
the unsteady acceleration because the velocity changes with re-
ρ ⋅Cp ⋅dT dT d 2T
spect to time), and the convective acceleration when a fluid parti- + ρ ⋅Cp ⋅uz ⋅ = K⋅ 2 + A
cle moves through a non-uniform flow field (i.e. field with spatial dt dz dz
gradients of velocity that may not vary with time), hence: This equation assumes that heat capacity Cp and conductivity K do

( ∂t )
∂V not vary with temperature, which is not the case in nature. Hence
m ⋅a = m⋅ + v ⋅ ∇v we re-write the equation above for a more general case in which Cp
and K depend on T:
The conservation of mass m is described by stating that the diver-
ρ ⋅ dCp ⋅dT dCp
dz ( dz )
gence of the velocity field is zero. The divergence is a scalar field dT d2 d 2T
+ ρ⋅ ⋅u ⋅ = 2 ⋅ K⋅ 2 + A
which measures how much the flow is expanding (positive diver- dt dz z dz
gence) or contracting (negative divergence) at any given point.
Zero divergence implies that the fluid is incompressible, a reason- The 3-D generalization of this equation is the equivalent per cubic
nable approximation for rocks (known as the Boussinesq approxi- meter of the equation on page 65:
mation). ρ ⋅∂Cp ⋅∂T
( )
+ ρ ⋅uz ⋅∇Cp ⋅∇T = ∇⋅ K⋅∇T + A
As a fluid parcel flows through a surrounding fluid, it gains or ∂t
looses energy via conduction, advection, shear heating, radiogenic
The Reynolds’ number is the ratio inertia to viscous force. Laminar
heat production, partial melting, phase changes etc. As its tempera-
flow occurs at low Reynolds numbers, where viscous forces domi-
ture varies its density changes, and therefore its inertia as well as
nantes. Flow at low Reynolds numbers is smooth and uniform (no
the gravitational force acting on it. At the scale of an entire system,
turbulence). This mode of flow best describes geological flows.
the variation of energy is a zero sum game as the energy lost by

66
Some examples of simple laminar flow for which analytical solu- In short, in the second half of the 20th century, computers were not
tions exist include the free fall of a spherical object into a newto- powerful enough to take advantage of these new numerical meth-
nian fluid (e.g. the settling of crystals into a magma), the flow in- ods. As computers grew in
duced by the motion of a rigid plate power, so did the complexity
above a newtonian fluid (Couette flow of fluid flow problems that
in relation to the motion of tectonic plate one could tackle. However,
above the asthenosphere), and the flow for a long time, computational
of newtonian fluid between two static tectonics involving a layered
plates (Poiseuille flow, for instance the lithosphere made of stronger
1950’s mainframe computers
flow of the lower crust in orogenic pla- layers (upper crust and upper
teaux). mantle) and weaker layers (lower crust and lower lithospheric
mantle) was limited to 2 dimensional models in which the node of
the computational grid had to conform to the deformation of the
model (Eulerian grid, see image below). In contrast, computational
mantle geodynamics could afford 3 dimensional models because
models of mantle convection could use a more efficient fixed grid
(Lagrangian grid) as long as the convective mantle was made of
one single newtonian fluid.

In the mid 1950’s, with the advent of mainframe computers, the

Navier-Stokes equations could be discretized and solved at the
nodes of a numerical grid to explore the type of complex, time-
dependent fluid flow we encounter in nature. At the same time, nu-
merical methods and computer algorithms were progressing so @ Yuen, Minneapolis
@ Beaumont, Halifax
fast that they quickly overtook the computer capability of the time.

67
Particle-in-cell numerical methods. In the late 1990’s, computers Over the past decade, the growing availability of powerful high-
become powerful enough to allow the implementation of a 1950’s performance computers has unleashed the power of computer-
numerical method called particle-in-cell (PIC, developed at Los Ala- based modeling in all branches of science. Geoscientists are now
mos National Laboratory) in which individual fluid elements, car- able to simulate in four dimensions, using realistic coupled ther-
rying material properties and flow history, are advected through a mal and mechanical properties, lithospheric-scale deformation and
fixed computational grid. For geologists, this progress meant that mantle geodynamics. 21st century HPC have finally caught-up
with 19th century physics, and numerical recipes from the mid-50s.

Numerical experiments vs numerical modeling vs numerical

simulation: Before going further we need to clarify the difference
between experiment, modeling and simulation. For most people these
concepts are interchangeable, but not for the experts.

Numerical experiments
PIC = Eulerian mesh + Lagrangian particles
like physical experi-
one could simulate the deformation of mechanically layered sys- ments do not pretend to
tems such as the Earth’s lithosphere; albeit in two dimensions. Ellip- reproduce a natural
sis developed by Moresi, first at Caltech then at CSIRO, is one of process in a realistic
the earliest and most robust codes (along with Citcom from which manner. The aim is to
Ellipsis has emerged) implementing an efficient PIC method on the illustrate a concept, or
back of a robust multigrid solver. to try to understand the
few most important pa-
rameters involved in a
particular process. The famous analogue experiment from Paul
Tapponnier, performed for the first time in the very late 70‘s, per-
Rey, Coltice, Flament, fectly illustrates the concept of escape tectonics, a process which
Nature 2014 accommodates convergence via the lateral escape of continental
blocks in front of a rigid indenter. Clearly, this experiment bears
Ellipsis has been used to show that early proto-continents and thick oceanic plateaux very little resemblance with tectonics in Asia and South East Asia,
had enough gravitational power to slowly force adjacent oceanic lithospheres to sub- (no orogenic plateau ...) but it does a nice job in illustrating a con-
duct. The model above includes the continental crust (red), lithospheric mantle (pink),
partially molten mantle (bright blue), the rest is the mantle asthenosphere and litho- cept which has changed our understanding of collisional tectonics.
sphere).

68
 Numerical modeling aims at under- in the context of geosciences, assumes that we are able to imple-
standing a particular process within ment a very broad range of coupled natural mechanical, petrologi-
a particular tectonic context. The ini- cal, geochemical processes including the migration of partial melt
tial and boundary conditions are and associated advection of heat, the release and consumption of
carefully thought through to describe water and other volatiles, evolving anisotropic rheologies due to
a geologically realistic setting. Re- grain reduction and crystallographic fabrics, dilatancy due to
sults of the modeling are detailed micro-cracking and other softening mechanisms. It will take dec-
enough to be able to be compared - to ades before we can properly engage with numerical simulation in a
the first order - to natural meaningful manner.
geological exam-
ples, without FEM and Multigrid solver: The Finite Element Method
trying to (FEM) consists in discretizing of a set of partial differential
match accu- equations across a 2D or 3D domain (i.e. ∂# are transformed
rately any- into ∆ # ). This allows solving efficiently the Navier-Stokes
one of them. equations at the nodes of a grid covering the model, but at
the price of a small error. To minimize this error, users can
use a very fine grid. The finer the grid the higher the resolu-
tion but the longer the compute time. Hence, the need to bal-
ance accuracy and compute time.
To solve the set of equations efficiently, computer scientists
Numerical simulations aim at reproducing with the greatest level of design methodologies taking advantage of parallel comput-
details and the greatest level of realism a particular process in a par- ing, as well as other techniques. One technique uses a multri-
ticular region. For most people this is what modeling is all about, grid solver. It involves the use of a stack of computational
and therefore they are quick to point to the many shortcomings of grids of increasing resolution. Instead of solving the set of
Paul Tapponnier’s experiments and dismiss its relevance since it equations directly at the node of one single high-resolution
neither account for the presence of the highest mountain chain nor grid, a long wavelength solution is first computed on a
the highest plateau on Earth. Numerical simulation assumes that coarser grid, and this solution is iteratively refined at the
our models are able to include parameters with realistic value dis- nodes of grids of increasing resolution taking into account
tributions and time-dependencies, as well as nested heterogenei- short wavelength pressure gradients.
ties on a broad range of scales. In addition, numerical simulations,

69
SECTION 2

Ellipsis: a FEM-PIC code

Ellipsis and Underworld apps solve the equations of Stokes rate, temperature etc). Lithologies are represented by particles
flow (Navier-Stokes equations in the limit of small Reynolds distributed inside the grid’s cells (hence the acronym PIC for
number) at the nodes of a computational grid using a Finite Particle In Cell). In this section, you will learn to compile Ellip-
Element Methods (FEM). Solving the Stokes equations pro- sis on your own Mac or Linux computer, and you will learn
vides the velocity and pressure fields from which a range of the command to launch and stop Ellipsis from a terminal win-
other field properties can be computed (viscosity, stress, strain dow.

70
Introduction: Ellipsis has a long history of development starting in the late 90‘s at Caltech by Prof. Moresi (now Melbourne Univer-
sity). It is derived from Citcom, a code dedicated to mantle convection modeling still under development today at Caltech. Ellipsis
was designed to explore 2-D coupled thermal and mechanical tectonic processes on a single CPU. It can simulate both brittle and
ductile deformation where brittle rheologies evolve with accumulated strain, and ductile rheologies evolve with temperature,
stress, strain rate and melt fraction. Underworld is a parallel version of Ellipsis with much expanded capacity and flexibility. It can
handle 2-D models on a laptop and 3-D models on high-performance computers. It is developed by a team led by Prof. Moresi,
with support from the Australian government through its Auscope NCRIS program. Its installation from its source code requires
an appropriate level of expertise. However, recent developments on virtual machine allow to deploy both Ellipsis and Underworld
in a Docker container on a local computer for the purpose of simple 2-D modeling.

Installation: Both codes run on unix-based OS (including Linux and Mac OSX), and it can also be installed on a PC running Win-
dows OS. Their source codes can be downloaded from geodynamics.org and underworldcode.org where detailed instructions can
be found to compile Underworld. To compile Ellipsis source one needs access to an appropriate compiler and a terminal window.
On a Mac computer, both come with the Xcode.app (free from the Mac App store). If you need to install a compiler go to
http://hpc.sourceforge.net/ and download the gcc compiler for your OS, then open a terminal window and execute: sudo tar -xvf
/path_to/gcc-5.1-bin.tar.gz -C /. On Linux, simply download a gcc compiler for your version of Linux and the 32 bit library (e.g., apt-
get install -y gcc-multilib lib32gcc-5.1). Then in a terminal navigate to the Ellipsis source and exe-
cute the three lines on the right (you may have to add the path to your compiler by adding ./configure CFLAGS='-m32'
CC=/path_to/cc_compiler after CFLAGS=’-m32’ ): make clean
An executable (only a few 100s bytes) called ellipsis3d will be created in the source directory. make
Copy this executable into the folder hosting your input files. A good practice is to keep a copy
of the Ellipsis executable with each model, so you can keep track of which version of Ellipsis you have been running. You are read-
ing to go ...
To run (i.e. execute) an input file, simply open a terminal window, navigate to the folder holding ellipsis3d and the input file (here
mycoolmodel.input) and enter the following command: ./ellipsis3d mycoolmodel.input
To stop the model, click the sequence control-C. More often than not, the reference model described in the input file requires access
to a file laying out the temperature at each node of the computational grid. This file is defined in the input file
(previous_temperature_file=”mytempfile.node_data”). We will see later how this file can be generated. This temperature file must be

71
present in the folder in which you are running the experiment, as
shown on the figure on the right.

As Ellipsis runs, it dumps at each time step a bunch of files into the
folder. These files, specified by the users in the input script, store
grid and particles information (temperature, pressure, velocity, stress, strain rate, viscosity, ...).
#.ppm0: These files are graphics ouput showing various lithologies, melt fraction, stress, viscosity etc.
#.node_data: These files store parameters such as temperature and pressure at the node of the computational grid. The informa-
tion in node_data can be used for data mapping (mapping of isotherms, super-solidus temperature, melt fraction, ...). For this you
will need to process the data stored in .node_data using Matlab, R, Python or any other appropriate tools.
#.profiles: These files give access to profiles (horizontal or vertical) passing through "sampling tracers". In the input file one can
specify the initial position of up to 50 sampling tracers (this position can remain fixed to the computational grid, or can move with
the rocks through the grid), the direction of the profile (horizontal or vertical), and the type of data to be
profiled (temperature, melt, depletion, pressure, velocity, strain rate, stress ... )
#.sample : These files are updated at each time step to add the new value of parameters attached to
one of the sampling tracers. Hence, there is one .sample file per sampling tracer. These files record
the evolution through time of data attached to a tracer (x, z, Data), where data can be temperature,
pressure, etc. These can be used to collect data to plot PTt paths.

Resolution of the computational grid and particle density: In Ellipsis, solving the Stokes equa-
tions is handled via a multigrid solver which uses a stack of increasingly finer computa-
tional grids. The pressure at a particular location is influenced by local/short wavelength
and distant/long wavelength sources. A coarse grid can effectively take care of distant
sources, the solution is then interpolated onto a finer grid to account for local sources, the
new solution is then averaged back onto a coarser grid. The larger the number of grids the
better the resolution, but the longer the running time. The pressure field is calculated using
a subset of the multigrid stack (pmg_levels), the maximum depth of which is levels-1.

72
The finer grid also contains Lagrangian tracers carrying information about material properties. The more tracers there are, the bet-
ter defined are the density interfaces. There is no simple rule for how many grid levels and tracers one should use. Consequently,
a model may take any time between a few minutes to several months. In most cases, levels=4 is sufficient when testing
lithospheric-scale models. Once the model has been properly tested at levels=4 to ensure that geometry and boundary conditions
are properly set, it is a good practice to let it run at levels≥4.

The figure on the right shows a grid stack (Levels=3: grid 0, grid 1 and grid 2) for
the multigrid solver. The number of cells in the finest grid is determined by the
Levels (here 3) and the parameters mgunitx (number of cells along the x direction
in grid 0, here 5) and mgunitz (number of cells along the z direction in grid 0,
here 3). The tracer_density is defined on the finest grid level. In this example trac-
er_density=4, meaning there are 4x4=16 tracers per element.

Initial temperature field: The temperature field is computed by running the in-
put script without any kinematic or dynamic boundary conditions (i.e. the model
is not extended or shortened). For maximum computational efficiency one can
increase all viscosities and turn plastic rheologies off so the model doesn’t deform while the geotherm evolves. This is achieved us-
ing the following overriding toggles (at the end of the Ellipsis input script):
TDEPV=off
" visc_max=5e31' '
' visc_min=1e31
' YIELD=off

… and increasing the fixed_timestep to 1 million years. Since we use second for unit of time ...:
" fixed_timestep=3.1536e13
Then, any #.node_data output can be used as a previous_temperature file by setting up in the Steady-States and Boundary Conditions
section of the input file:
previous_temperature_file="A_My_Previous_Temperature_File.node_data"

73
Time step: Unless when running thermal equilibration, time step should be kept below 6.3072e12 sec (that’s 200,000 years). If
weird things happen (i.e. unexpected short-term fluctuation of some parameters like viscosity) try to reduce the time step.

Ellipsis can also be delivered to a local machine via a Docker container (a kind of virtual machine) in which an Ellipsis executable is
ready to run. The next section explains how it works.

Using Ellipsis on Jupyter environment within a Docker container:

“Docker is an open platform for developing, shipping, and running applications.” A Docker container con-
tains a lean Linux operating system with all necessary libraries to compile codes and run apps. A
Docker container can run on any machine and any operating system as long as it as the Docker app
installed (https://www.docker.com).
We have created a Docker container for Ellipsis which is available for download through Docker
Store and/or Kitematic both web-based community libraries for Docker containers, available from
the Docker drop-down menu. nb1: A Docker image is the container in which your code run. The
command docker push publishes a Docker image. nb2: A Dockerfile is the file that is used to build (i.e. compile) the Docker image via
the command docker build.

The Ellipsis Docker contains the Ellipsis app, and a series of input files. It also includes Jupyter a browser-based integrated develop-
ment environment from which Ellipsis input files can be created and edited, Ellipsis executable can be launched, and Ellipsis out-
puts can be processed.

To run Docker containers users must first download the Docker app onto their machine.
This can be done via: https://www.docker.com/community-edition

Ellipsis can be run from the jupyter environment. Jupyter offers a collaborative environ-
ment in which one can run codes, scripts and calculations, and process data using the
python ecosystem, as well as other computer languages (>50 of them), to produce new data and visual outputs.

74
SECTION 3

Underworld Growing access to high-performance computers and

their relentless increase in power are changing the way
we do Science. We have now the numerical capability
to explore complex non-linear geological processes in a
well constrained physical framework. Underworld is to
tectonicists what a laboratory is to physicists, a space
where hypotheses can be tested not through physical
experiments, but through numerical experiments.

Underworld is an open source, modular, hierarchical,

computational framework, designed to perform parallel
2D or 3D time-dependent coupled thermal and mechani-
cal tectonic and geodynamic experiments. It uses an
expandable range of rheologies covering elastic, plastic
and viscous behaviors. It accommodates radiogenic
(i.e. internal) heating, partial melting - and other meta-
morphic reactions - with feedback on rheologies, densi-
ties and temperature.

Underworld development is led by Prof. L. Moresi (ANU)

and a team of developers led by John Mansour (at
Monash University) and Julian Giordani (Melbourne Uni-
versity). Underworld started in 2005 with funding from
the Victorian Partnership of Advanced Computing, and
carried on with support from the Australian NCRIS/
Auscope program.

75
What does Underworld do?

Underworld solves the equations governing Stokes flow (after George Gabriel Stokes 1819-1903)
a class of flow characterized by a small Reynolds number (after Osborne Reynolds, 1842-1912)
a dimensionless number which expresses the ratio between inertial forces to viscous forces (cf.
page 80). Stokes flow are characterized by smooth laminar fluid motion, in contrast to turbu-
lent and chaotic flow occurring at high Reynolds number, for which the full Navier-Stokes
equations are necessary. The sketch on the right shows an example of Stokes flow (source:
wikipedia).

Underworld solves the Stokes equations on a cartesian 2-D or 3-D grid assuming the flowing ma-
terial is incompressible, which imposes a zero divergence condition on the velocity field, and
therefore the mass continuity equation ui,i = 0 which ensures conservation of mass. The
Stokes’ equation in a tensorial form (see page 84 to learn how to unpack this form) is:

" " " " σij = τij − pδij = fi

This equation expresses the conservation of momentum. It states that the total stress tensor σij (the sum of the deviatoric stress ten-
sor τij and the pressure tensor pδij ) balances the buoyancy force fi (fi = ραTg) due to density variations (also known as the gravita-
tional body force) .

The constitutive equation expresses the relationship between the deviatoric stress tensor, the viscosity tensor (or constitutive ten-
sor) and the strain rate tensor: τij = 2ηijkl . ϵ̇kl

Hence, the resulting Stokes’ equation can be writen as: 2 ηijkl . ϵ̇kl − pδij = fi

or in a vectorized form the momentum equation is: η ∇2 u − ∇p = ρg

In which ∇2 u and ∇p are the velocity and pressure gradient respectively, η is the viscosity, ρg the driving force and η ∇2 u is the
stress gradient. The left term is the divergence of stress, the right term is the body force.

76
By solving the Stokes’ equation Underworld retrieves the pressure and the velocity at the nodes of a 2-D or 3-D mesh. Material prop-
erties including density, viscosity, strain-rate, etc are updated and stored into particles which are advected through the grid.

Keeping track of the temperature is critical as it impacts on both densities and viscosities. The conduction-advection heat transfer
equation takes into account the latent heat of fusion consumed during partial melting, a metamorphic reaction of particular impor-
tance because it consumes vast amount of latent heat and because a small fraction of melt can drastically reduces the viscosity of
partially melted rocks. The following equation describes the rate of change of temperature as a function of heat diffusion ( κ ), ra-
diogenic heat production (H), heat advection ( uz ), and partial melting (DF is the variation of melt fraction):

Cp ( )
∂T H ∆ S DF
" " " " " " " " " = κ ∇2 T + − uz ∙ ∇ T − T
∂t Cp Dt

The equation above ensures conservation of energy and allows the coupling of densities and viscosities to temperature through
the following equations:
" " " " " " " " " " ρ(T ) = gρ0(1 − α( ∆ T ) − βF)

( ) 2 (n ⋅ R ⋅ T)
1 E
" " " " " " " " " η T = A −1/n ⋅ Exp ⋅ ε̇(1−n)/n

Where A is the pre-stress factor, n is the stress exponent, E is the activation energy, and ϵ̇ is the strain rate.

With this the momentum equation becomes: η(T ) ∇2 u − ∇p = gρ(T )

Plastic deformation is introduced by considering that above a yield stress ( τ ) plastic materials become viscous. The yield stress is:

( )
" " " " " " " " " " τ= C0 + μeff ⋅ (ρgz) ⋅ f (ϵ)

Where f (ϵ) is strain weakening function (fault zones become weaker as strain accumulate), Co is the cohesion, μeff is the effective co-
efficient of friction (i.e. it takes into account the effect of pore pressure), ρgz is the confining pressure.

77
Once the material passes its yield stress, the material behave as a viscous fluid with a post-yielding viscosity described by:""
τ
" " " " " " " " " " ηyield =
(2 ∑ij ϵ̇ij ϵ̇ij)
1/2
1
2⋅

Where the denominator is twice the second invariant of the strain tensor (also called the effective strain rate).

∂uj
2 [ ∂xj ∂xi ]
1 ∂ui
" " " " " " " " " " ϵ̇ij = +

78
How does Underworld do it?

Underworld is one component in a modular numerical framework made of a stack of specialized modules (numerical libraries)
working together to solve efficiently the Stokes’ equations in the context of a tectonics and/geodynamics model. This framework
also provides gLucifer a visualization platform developed by the Underworld team.

•Underworld - Geophysical formulations, e.g. rheologies.

Underworld
•gLucifer - Visualisation package. Constitutive equations: Rheology, Radiogenic heating

•PICellerator - Particle in Cell methods (Lagrangian / Eule-

History
Particle-in-cellerator gLucifer
Handle Particles Parallel visualization
rian framework). sensitive
materials

•StgFEM - Finite Element Method (FEM) framework (Eule- StGFEM

Solve system of linear equation

rian). Includes input files and plugins stored in

StgFEM_Components or StgFEM_plugins directory. StGDomain
Math, geometry, mesh, particles

•StgDomain - Basic mathematical formulations, model ge-

ometry, mesh, shape and basic swarm framework.
StGermain
Parallel scientific applications

•StGermain - Core parallel programming architecture, mem-

ory management and Object Oriented frameworks.

Many of these components are built themselves upon other libraries (PETsc, OpenMPI or MPICH, HDF5 ...) which explains that
the deployment of the Underworld framework can be challenging.

Users access Underworld through “input files” the purpose of which is to describe the architecture of the model, the physical prop-
erties of materials that compose it, the time-dependent boundary conditions under which the model exists, and the optional and
user specified passive markers necessary to capture the thermal and mechanical evolution of the model.

79
Underworld 1.X - Up to the version 1.8 Underworld input file is distributed across a few xml files (xml: extensible markup lan-
guage). xml allows the description, and transport of data in a simple text format with a trees-like hierarchical structures of ele-
ments (elements, sub-elements, sub-sub-elements ...) that can be stich together using the tag “include”: <include>
mycrustRheology.xml </include>

Tags (<tag></tag>) in xml are not pre-defined like in HTML. xml allows for customizable case sensitive tags with quoted attrib-
utes: <density=”2700”>

Because tags are defined by users, xml is self descriptive and more information can be added using <!-- this is a comment -->:

Example1: Example2:
Here is a shape: Here is a shape involving the Intersection (!) function:
<struct name="my_diapir"> <struct name="my_diapir_host">
<param name="Type">Sphere</param> <param name="Type">Intersection</param>
<param name="CentreX" units=”km” > 0.5 </param> <list name="shapes">
<param name="CentreY" units=”km” > 0.5 </param> <param> myrectangle </param>
<param name="CentreZ" units=”km” > minZ </param> <param> ! my_diapir </param>
<param name="radius" units=”km” > 0.1 </param> </list>
</struct> </struct>

Example3: Example4:
Here we define a simple viscous rheology: Here we define a viscosity limiter function:
<struct name="my_diapirViscosity"> <struct name="viscosityLimiter">
<param name="Type">FrankKamenetskii</param> <param name="Type">ViscosityLimiter</param>
<param name="TemperatureField">TemperatureField</param> <param name="maxViscosity" units="Pa">1e23</param>
<param name="eta0" units="Pa*s">1.0e26</param> <param name="minViscosity" units="Pa">1e18</param>
<param name="theta">13.815510558</param> </struct>
</struct>

80
Here we define a plastic rheology: We put everything together adding a few bits to define the thermal and
Example 5 mechanical properties of materialOne
<struct name="myVonMises">
<param name="Type">VonMises</param>
Example 7
<param name="StrainRateField">StrainRateField</param>
<struct name="materialOne">
<param name="cohesion">50.0</param>
<param name="Shape"> my_diapir_host </param>
<param name="StrainWeakening">strainWeakening</param>
<param name="density" units=”kg.m-3”>2800</param>
</struct>
<param name="alpha" units=”K-1”>3e-5</param>
<param name="Diffusivity">1e-6</param>
Here we define a weakening function in which:
<list name="Rheologies">
initialDamageFraction is the fraction of particle with initial strain. <param> myVonMises </param>
initialDamageFactor is the maximum initial strain of a particle <param> strainWeakening </param>
initialSofteningStrain is the strain at which material starts to yield <param> my_diapirViscosity </param>
<param>viscosityLimiter</param>
finalSofteningStrain is the strain at which the material has fully yielded
</list>
initialDamageWaveNumberCosI (or ...SinI) is the inverse of the horizon- <list name="heatingElements">
tal length of the damaged box. <param name="Q" units=”W.kg-1”>1.8e-12</param>
<param name="Cp" units=”J.kg-1.K-1”>1000</param>
</list>
Example 6
</struct>
<struct name="strainWeakening">
<param name="Type">StrainWeakening</param>
<param name="TimeIntegrator">timeIntegrator_aux</param>
<param name="MaterialPointsSwarm">materialSwarm</param>
<param name="softeningStrain">0.1</param>
<param name="initialDamageFraction">0.0</param>
<param name="initialDamageWaveNumberCosI">0.5</param>
<param name="initialDamageFactor">0.5</param>
<param name="healingRate">0.0</param>
</struct>

81
Underworld 2.0 - Since the version 2.0, Underworld input files are written in python, and therefore can be assembled within a Ju-
pyter notebook. The Underworld2 code itself (UW2) can be either deployed and compiled as an executable application stored in the
user’s local machine, or deployed and compiled into a remote high-performance computer. This process demands a robust exper-
tise, which involves the installation of numerical libraries (MPI, PETSC, HDF5, etc). Once the input Jupyter notebook has been cre-
ated, it is saved as a python script via the terminal command jupyter nbconvert --to python myinputfile.ipynb or via File > Download
as .py. To execute this UW2 input script, open a terminal window and execute: mpirun -np 4 python myinputfile.py, where -np 4
specifies to run UW2 in parallel on 4 cpu (alternatively mpiexec 4 python myinputfile.py).

To simplify things a little, UW2 executable is made available inside a Docker container. After installing Docker, go to Kitematic (in
the Docker’s drop-down menu) and search for Underworld2. There is a number of versions including the robust Underworld2_latest.
nb: The installation of Docker and Jupyter are detailed at: https://geoslearn.github.io/LabDeploy.html.

UW2 jupyter notebooks can be run inside Docker using the host machines cpus.

Readers who want to experiment with Underworld2 should consult the following website:
http://www.underworldcode.org
https://github.com/underworldcode
https://github.com/rbeucher/UWGeodynamics.

82
Interactive 5.1 An example of UW-Geodynamics Jupyter input script (double click on the image below to scroll down)

83
INERTIA, VISCOUS FORCE AND THE REYNOLDS’ NUMBER Tensors describe linear relationships between scalars, vectors and
other tensors. The dimensionality of tensor, its rank, informs about the
Inertia is the resistance to a change of motion it is expressed by New- nature of the tensor, e.g.:
ton’s second law of motion F = m.a. Earth geodynamics (mantle flow a tensor of rank 0 is a scalar;
processes) and tectonics (plate tectonic processes) involve very small a tensor of rank 1 is a vector (1- dimen-
acceleration, hence very small inertia. sional array);
a tensor of rank 2 is a matrix (2-
dimensional array) ...
The viscosity of rocks is a measure of their capacity to resist flow.
Tensors provide a concise way for for-
Hence, the viscous force is the resistive force to flow. The viscosity of
mulating relationship between quanti-
solid rocks is very large typically >1e18 Pa.s. ties, the drawback is that they can be
hard to read.
The Reynolds’ number is the ratio between inertia force and the vis-
cous force. Since on Earth the inertia force is very small and the vis- The stress tensor is a second order tensor which collects into col-
cous force is very large, then the Reynolds’ number for tectonic and umns the coordinates of 3 traction vectors (Tei) acting on 3 perpendicu-
geodynamic flow is very small. lar faces of a cube. The value of these coordinates depends on the
choice of the reference framework. Nevertheless, these coordinates
Quiz: Let’s consider a tectonic plate 4000 km x 2500 km x 150 km ac- fully describe the state of stress on the central point in the cube. By
celerating from 0 to 5 cm.yr-1 in 100,000 years. The density of the convention the component σij refers to coordinates along axis i of the
plate is 3000 kg.m-3, and the viscosity of the asthenosphere is 1e20 traction acting on face j. The number of indices relates to the order of
Pa.s. What is the Reynolds number of the mantle flow induced by this the tensor (order, degree, rank are synonymous).
acceleration, assuming that the plate motion causes shearing distrib-
Einstein’s summation convention implies summation when an index
uted over 10 km of asthenosphere?
variable appears twice in a single term:
n.b.: shear stress = ~ strain rate * viscosity 3

∑
ci xi ⟺ ci xi ⟺ c1x1 + c2 x2 + c3 x3
i=1
3 3

∑∑
hence: τij = 2ηijkl ⋅ ϵ˙kl = 2ηijkl ⋅ ϵ˙kl
k=1 l=1

84
CHAPTER 6

Lithospheric Deformation

The mechanical behavior of the Earth’s lithosphere under dynamic or kinematic boundary conditions is strongly
non-linear, and therefore difficult to predict. This is due to the exponential dependency of the rheology of rocks
to temperature, as well as weakening mechanisms including strain (which reduces the grain size of rocks and
helps fluids circulation) and partial melting. Hence, small variation in temperature, strain or melt fraction can
result in strong variations in mechanical properties and mechanical behaviors. In this chapter we turn to
computational tectonics to explore lithospheric deformation in a variety of tectonic settings.

85
SECTION 1

Collisional tectonics
Content

1. A generic 2-D thermo-

mechanical experiment of
collisional tectonics
2. The European Alps

This section shows that numeri-

cal experiments, involving realis-
tic mechanical and thermal prop-
erties, can produce collisional
structures observed in the Euro-
pean Alps.

Tectonic convergence between plates is accommodated continental collision. In what follows we explore colli-
by a range of processes including subduction (oceanic or sional processes through numerical experiments. In our
continental), coupled shortening and thickening, and coupled thermal-mechanical numerical experiments the
transcurrent tectonics involving one or more lithospheric- computational grid is 768 km x 1536 km, with a cell reso-
scale strike-slip faults. Following the closure of oceanic lution of 2 km x 2 km. The rheology of all rocks in the
domains, continents may subduct. However, continental model is visco-plastic, excepted for the asthenospheric
subduction is limited by the buoyancy of the continental mantle which is assumed iso-viscous (1e20 Pa.s).
crust, and continental subduction is rapidly followed by

86
Subduction and continental collision -

(a) - The converging continents are 160 km thick, and are ad-
jacent to oceanic lithospheres of thicknesses 110 km and 40 a
km. The continental crust consists in 20 km of upper crust
(lavender) above a basement 20 km thick (orange). The oce-
anic crust is 6 km thick (pink). A fault dipping left decouples the
thick and thin oceanic lithospheres. Convergence is achieved by
pushing the right continental lithosphere to the left, whereas the
left wall of the continental lithosphere to the left is fixed. Once a
~300 km slab has formed, slab pull becomes the only driving b
force. Erosion removes rocks above 5 km elevation, and sedimen-
tation fills all spaces below -2 km depth. Sediments are color
coded yellow to blue from older to younger sediments.

(b-c) - Initiation of subduction leads to extension and thinning of c

the margin of the continent on the left, and deposition of sedi-
ments above the subducting oceanic plate. The flexure of the sub-
ducting plate provides the space for extension.

(d-e) - The thin oceanic plate is forced to subduct, whereas the d

sedimentary basin is transferred from the original downgoing
plate onto the upper plate. As the continental crust from the up-
per plate enters the subduction zone its motion is blocked by the
colliding lithospheric mantles. A crustal-scale overturned fold
develops, and the structural vergence switch to the left.
e

87
(f) - The system has reached the stage of continen-
tal collision. The thinned margin of the left conti-
nent is being subducted underneath a thick
wedge made mainly of crustal material from the
right continent. Crustal rocks reach a depth > 100
km. Its only at this stage that high topography is
being generated.

(g-h) - The lithospheric mantle of the continent on

f
the right is peeling-off as the crust is being trans-
ferred onto the left continental plate, the margin
of which is been completely smeared. The upper
crust of the left continent and the sedimentary ba-
sin are being pushed to the left. A set of left verg-
ing folds develop. The upper part of the hinge of
the overturned fold is sheared by a normal fault,
and translated to the left above the sedimentary g
basin. This leads to the exhumation of deep crust.

88
The European Alps is an example of collisional orogen which follows the subduction of oceanic lithosphere. It shows, at various
scale, a range of structures resembling those observed in our numerical experiment. When considering similarities between nu-
merical experiments and natural examples, one needs to exert caution and keep in mind that nature is always far more compli-
cated that our numerical representations of it.

The Alps results from the convergence and collision of the Adriatic micro plate, an African promontory sometimes called Apulia
(in yellow on the map) and Europe (in blue), in between which the Briançonnais microcontinent (in green) was squeezed. Conver-
gence started at around 70 myr ago with the south dipping subduction of the Ligurian-Piemont ocean that separated the Briançon-
nais terrane and the African margin. By 35 myr ago this ocean was closed and the Briançonnais, with the narrow Valais ocean that
separated the Briançonnais
and Europe, were squeezed
in between the colliding
European and African mar-
gins. The deeper, internal
part of the Briançonnais ter-
rane was buried to great
depth (>100 km), under the
African margin (in orange in
the cross-section on the
next page), whereas the ex-
ternal part (zone houillère,
pale green) escaped high-
pressure metamorphism.
High-temperature metamor-
phism (up to 800ºC) was
limited in time (~30-15
myr) and space (subpen-
ninic nappes of the central
Alps).

89
The flysch sediments (dotted regions on the cross-sections)
that derived from erosion of the African plate and deposited
on the subducting Valais margin where translated onto the
European margin, as the Briançonnais was exhumed, trans-
lated over the European margin, and folded into recumbent
folds. A combination of folding and gravitational spreading
exhumed the Helvectic nappes (pale blue) in a tectonic win-
dow. These thin-skinned nappes are made of sedimentary
rocks deposited over the European rifted platform. They
were sheared-off from their pre-Mesozoic basement to form
a fold-and-thrust belt. The European pre-Mesozoic (i.e. Varis-
can) basement was also deformed into the so-called Subpen-
ninic nappes (in pink). Ongoing collision resulted on back-
thrusting of the Briançonnais over the Adriatic plate.

A the plate-scale, the Alps involves a complex and seg-

mented subduction zone system wrapping around Adria
(see figure next page). The Alps now stands above a suture
Transect NFP-20 West, Schmidt, et al., The TRANSMED Atlas, 2004
zone involving: to the west an east to south dipping subduc-
tion, and to the east a north to northeast dipping subduction. At the surface, the Alps presents itself as a curvilinear belt with
sharp bends: the Cantabrian arc, the arc of the western Alps, the Carpathian arc, and the Helenic arc. The orogen itself is rela-
tively symmetric and bivergent with opposite structural vergences.

Our simple 2-D numerical models cannot account for this sort of complexity. Nevertheless, it is interesting to note that to the first
order our numerical experiment reproduces:
• The burial of basement rocks down to 100 km.
• Their deformation into crustal-scale recumbent folds.
• Formation and denudation into a tectonic window of thin-skin Helvetic-style nappes.
• The emplacement of nappes of external (i.e. upper crustal) basement rocks above younger allochthonous flysch.

90
 # R script to produce the map on the next page from gebco data.

# Load the libraries

library(marmap)
library(raster)
library(RColorBrewer)
library(ggplot2)
library(devtools)

#setwd stands for set working directory

setwd('/path_to_your_directory/')

# Etopo color scheme, modified from

https://github.com/cran/marmap/blob/master/R/palette.etopo.R

etopoAlps <- read.csv(textConnection(

"altitudes,colours
Handy et al., Int. J. Earth Sci., 2014. 4500,#FBFBFB
2000,#864747
1950,#7E4B11
1000,#9B8411
950,#BD8D15
300,#F0CF9F
0,#307424
-1,#AFDCF4
-3000,#090B6A"
), stringsAsFactors=FALSE)

etopoAlps$altitudes01 <- scales::rescale(etopoAlps$altitudes)

etopoAlps.colors <- function(n) {

colorRampPalette(etopoAlps$colours)(n)
}

scale_fill_etopo <- function(...) {

ggplot2::scale_fill_gradientn(colours=etopoAlps$colours, values=etopoAlps$alt-
itudes01, limits=range(etopoAlps$altitudes), ...)
}

scale_colour_etopo <- function(...) {

Moho depth and suture zone around the Adria, Luth ggplot2::scale_colour_gradientn(colours=etopoAlps$colours, values=etopoAlps$alt-
et al., Tectonophysics, 2013. itudes01, limits=range(etopoAlps$altitudes), ...)
}

91
Topographic map of the Alps (see R script), with geology of the western Alps.

Zhao et al., Geology, 2015

92
SECTION 2

Transcurrent tectonics
Lorem Ipsum

1. A generic 3-D thermo-

mechanical experiment of
pull-apart basin
2. The Dead Sea basin

Here we document crustal flow in an orogenic pull-apart velop in the deep crust. A block of lithosphere
where localized extension and crustal thinning focus the 384x256x128 km is mapped over a computational mesh
exhumation of deep crust. Using numerical experiments consisting of 240x160x80 elements. It includes an 8 km
we show that an extensional strain regime, collinear with thick layer of air, a layer of continental crust 60 km thick
the direction of plate motion, is partitioned into the shal- (temperature at the Moho is 830ºC), and a layer of mantle.
low crust, whereas contractional structures and fabrics - Two non-overlapping, parallel and vertical faults extend-
at a high angle to the direction of imposed transport - de- ing from the surface to 20 km depth, and laterally over 64

93
km from the two opposite vertical walls. Kinematic bound-
View from the top showing both
ary conditions drive divergence at a total velocity of 2.56 cm
extensional and shear fractures.
yr-1 promoting the formation of a dextral pull-apart.

After 1 myr of displacement (i.e. total displacement of 25.6

km), a fractures network invades the brittle upper crust in
the region in between the two strike slip faults. This net-
work consists of i/short extensional fractures oriented per-
pendicularly to the direction of extension imposed by the ve-
locity boundary conditions, and ii/ longer shear fractures
that merge into the two faults zones. Damages from this frac-
ture network explain the localization of further deformation
into this region.
256 Brittle strain damage after 1 myr of
km
256 km
1.2

deformation
8c
m/
yr

Fault
1.2

Fault
8 cm

64 k

20
1.2
m
/yr

8c

km
m/
yr

64 km 16 k
m
1.2

st
8 cm

120
0 r cru km
p e 40
/yr

800 Up ust km
400 r cr 60
0 e
Low
le
m ant
ri c ºC
o s phe 1330
Lith
120 km

ere
o sph
hen
Ast

94
Meanwhile in the deeper crust ... As displacement increases the upper crust gets thinner, and pressure gradients in the
weak deeper crust drives the flow of the deep crust towards the pull-apart region, resulting in the formation of a dome.

The long axis of the dome joins the tips of the strike-slip faults. Strain markers show that the direction of stretching is co-
linear with the imposed velocity boundary condition. As a result, the stretching lineation plunges to the west at the west-
ern end of the dome, and to the east at the eastern end of the dome. Deeper in the dome the strain defines two sub-domes
separated by a median high-strain zone.

Gallery 6.1 Architecture of a dome in a pull-apart. Strain marker

0 km

N
16 km

60 km

Depth (km)
-40 -30 -20 -10

-43.8 -4.43

N
Shallow crust

Deep crust

Strain markers Solidus

initially at
Moho
42 km depth

View to the NW showing the dome underneath the upper brittle crust (dark
and pale blue layers). The red surface is the solidus.

95
Movie 6.1 Faults network in a pull-apart basin between two non-overlapping vertical faults.

96
Movie 6.2 Faults network in a pull-apart basin between two overlapping vertical faults.

97
The Dead Sea: an example of pull-apart basin

The Dead Sea transform fault system runs along the mar-
gin of the Arabian shield, and mechanically links the
opening of the Red Sea to the south and to the Taurus-
Zagros collision to the north. Over the past 15 to 20 myr, it
has accommodated ca. 100 km of sinistral motion be-
tween the African plate and the Arabian plate, motion
slowing down to 6 mm per year over the past 5 myr. The
Dead Sea basin forms in a narrow corridor no large than
15 km in between two left lateral strike-slip fault zones:
the Jericho (Jordan) fault to the west and the Arava fault
to the east. The Dead Sea itself is a 350 m deep hyper-
saline lake (dark blue in far right map), its water level stands 400 m below the Mediterranean sea level.

The sedimentary basin is 8.5 km deep with very steep subsidence Shalev et al., Tectonophysics, 2012
gradient across closely spaced faults, and an unusual flat basin floor.
It is segmented into sub-basins by transverse normal faults highly
oblique to the strike slip faults (Smit et al., Tectonophysics, 2008).
The flanks of the basin show different characteristics. The east flank
is higher and dissected by river gorges highly oblique to the basin.
The west flank shows NNE trending topographic ridges that con-
form with the left-lateral sense of shear, and suggests that the litho-
sphere on the west flank is relatively weaker compared to the Ara-
bian shield. Heat flow data reveals a negative thermal anomaly that
may corresponds to the downward advection of cold sediments.

98
Interestingly, seismic data shows that the lower crust from ca. 18 km down to the Moho is not affected by the basin devel-
opment.

In summary: The Dead Sea basin is characterized by a set on

particular features including:
• Two closely spaced bounding faults
• A flat basin floor
• No heat flow anomaly
• A deep crust that shows no thinning
• A Moho that shows no offset

To explains these features one can hypothesize that when at

the Earth’s surface the spacing between the two bounding
faults (here the Jericho fault and the Arava fault) is smaller
than that thickness of the brittle crust, the faults merge into
a ductile strike slip shear zone upon entering the deep
crust. Hence: extension and thinning is strongly partitioned
into the upper brittle crust; sedimentation keeps up with Bottom panel: E-W seismic profile across the Dead Sea basin
showing P-wave velocity contours. Interpreted geology: Basin
subsidence of the basin floor; no deep isostatic flow is re-
fill (dark yellow); pre-basin fill sediments (light yellow); upper
quired; and as a consequence the Moho remains flat and crust (green); lower crust (orange). Top panel: comparison be-
there is no enhanced heat flow. To test this hypothesis a 3- tween observed and calculated free-air gravity anomaly using a
D numerical experimentation can be set-up. density model derived from the P-wave velocity model.

99
SECTION 3

Extensional tectonics
Lorem Ipsum

1. Extensional collapse
2. Continental extension
3. Continental rifting
4. Mid-oceanic ridges

Luke Mondy, 2016

Continental extension occurs in a limited number of modes. Because continents have insulating properties, the tempera-
ture excess that buildup in the mantle underneath supercontinents lead to dynamic topography which generate exten-
sional stresses that drive continental breakup, deformation focussing on pre-existing zones of weakest such as ancient su-
ture zones. In a process similar to the calving of iceberg from the margin of polar ice caps, the recurrent detachment of
continental blocks at the margins of continents is another mode of continental extension. In this case, slab rock back
and/or the upward flow of the mantle wedge drive extension. Extension can also be driven by gravitational stresses in-
side the lithosphere. These stresses are the result of lateral variation of gravitational potential energy due to lateral

100
change in density. On a regional scale, this process explains extensional tectonics in orogenic systems. The north-south
trending normal faults in the Himalayas and in the Tibetan plateau, and the north-south trending normal fault array in the
Basin and Range (western USA) are the product of gravitational collapse. At a local scale, gravitational collapse can also ex-
plain the normal faults in a pull-apart basin (see the section on the Dead Sea Basin).

The collapse of orogenic plateaux- During their build-up phase, orogenic plateaux store significant excess of gravitational
potential energy. As the plateau region grow weaker through thermal relaxation and heating, this excess of energy drives
their gravitational collapse, and - in the case of fixed boundary collapse - contractional shortening of foreland regions. One
can use 2-D numerical experiments to explore the process of plateau collapse. The Underworld model below is 150 km deep
and 900 km long; and it is mapped over a grid of 2 km resolution. It uses appropriate thermo-mechanical properties (radio-
genic heat production, diffusivity, rheology, density etc), fixed lateral boundary conditions, constant temperature on top,
and an isostatic condition maintain the pressure constant at the base of the model. In this experiment, the plateau lower
crust flows horizontally at a velocity of a few cm per year into the foreland which becomes thicker while the orogenic crust
becomes thinner. In this model, the strength of the upper crust is such that there is no measurable deformation at the sur-
face.
Passive tracers

40 km Continental crust
Orogenic plateau 65 km

Mantle

101
Extension of orogenic lithosphere and the formation of metamorphic core complexes- Extension of thick and hot crust
leads to the development of metamorphic core complex (MCC). The Underworld experiment below represents a 120 km x
360 km lithosphere mapped over a computational grid with a resolution of 2 km. The crust is 60 km thick and the tempera-
ture at the Moho has reached ~800ºC. A ~2 km thick
Upper crust Fault
fault-shaped prism made of weak material is embedded
60 km
into the upper crust. Its role is to localize deformation Passive tracers
Lower crust
away from the vertical walls in order to mitigate bound-
ary effects. This experiment uses appropriate thermo- Mantle

mechanical properties. The boundary conditions involve a

kinematic condition pulling of both vertical
walls at a velocity of 1.8 cm per year, a con-
stant temperature on top, a constant heat
flow and isostatic condition at the base of Partially melted lower crust
the model. To assess the deformation of the
lower crust, and array of circular passive
tracers is embedded into the deep crust.

Extension leads to the thinning of the conti-

nental crust, and isothermal decompression
of the hot deep crust which results in decompression melting. The presence of melt in the lower crust reduces the viscos-
ity, and increases the mobility of the deep crust which flows easily in response to lateral pressure gradient. Extension and
thinning in the upper crust is localized around zones of brittle deformation. The necking of the brittle upper crust is ac-
companied by the formation of gneiss dome. The structure in the deeper part of the dome involves a sub-vertical high-
strain zone separating two colliding recumbent folds, which evolves into two sub-domes.

102
The “calving” of continental ribbons at passive margins - To the east of Aus-
tralia’s east coast, the Tasman sea separates Australia and a girdle of continen-
tal rafts with include the north and south islands of New Zealand. Most of
these rafts consist of submerged ribbons of thinned continental crust. They
were detached at around 85 myr ago from Gondwana’s Pacific margin at a
time when Antarctica and Australia were still attached.

A number of hypothesis have been proposed to explain the detachment of

continental blocks. These include:
• The slab pull force acting on a passive continental margin attached to a
oceanic lithosphere being subducted.
• The pull on a continental margin from a subducting slab retreating to-
ward the ocean (i.e. slab rollback).
• Extensional gravitational stresses due to the buoyant rise of the mantle
wedge above a subducting slab.
• The collapse of a cordillera due to a reduction of plate velocity or a
change in the direction of plate motion.

b1/

FRP
re
sphe
ic litho
Ocean sphe
re FSP
FG ic litho
Ocean

here
lithosphere Asthenosp ph ere
ental enos
Contin Asth
l Force
tiona
Gravita sphe
re
tal litho l Force
inen tiona
Cont ta
Gravi
FSS

FSS

FG
FRP
re
sphe
here ic litho
lithosp Ocean
ental
Contin
l Force
tiona
Gravita here
Asthenosp

FSS

FSS

103
Continental blocks and continental ribbons are a common occurrence in particular along the west margin of the Pacific
ocean. Numerical experiments can help understand the range of conditions amenable to their formation. The model be-
low is a 2-D representation of a passive continental margin. Its dimensions are 768 km x 192 km, and it involves an oceanic
lithosphere, a continental lithosphere, and a transitional domain in between which represents the continental margin.
Each domain is 256 km long. The model is extended through the imposing of a velocity condition. We explore the role of
the extensional velocity and the age of the oceanic lithosphere.
Ocean Margin Continent
When pulling at 3 cm/yr at both sides, our ex-
periment suggests that a 10 myr old oceanic
Mantle
lithosphere is weaker than the continental mar-
gins, in which case rifting occurs at the junction
between the ocean and the continental margin.
A 20 myr old oceanic lithosphere is stronger than
the continental margin. In this case extension and 10 myr old oceanic lithosphere
thinning occurs at the junction between the conti-
nent and the margin. Interestingly, extension gets
Velocity: 3 cm/yr at each vertical wall
localized in two rifts.

Reducing the total extensional velocity from 6 cm/

20 myr old oceanic lithosphere
yr down to 2 cm/yr leads to one single rift. One can
concludes that a single rift can accommodate exten-
Velocity: 1 cm/yr at each vertical wall
sion up to a given velocity. For faster extension,
the system evolves multiple rift zones to keep up
with the imposed velocity, in which case multiple
20 myr old oceanic lithosphere
ribbons can develop.
Nikita Golkin, Sydney Uni. Honours thesis.2015

104
The “calving” of continental ribbons at active margins - In this experi- cont. crust
cordillera lithospheric mantle
ments we simulate an active continental margin supporting a cordillera fault
orogen (crustal thickness 60 km). The model is 1536 km long and 768
km deep. The initial reference model at time t = to is produced by push-
convective mantle
ing an oceanic lithosphere toward the continent at 5 cm/yr. A fault-
like zone of weakness insures that the oceanic plate slides underneath
the continent. To maintain a constant mass into the model and insure
isostasy, an equivalent outflow is imposed to the right wall (0.791 cm/
yr). A zero horizontal velocity condition (ux = 0) is imposed on the left
wall. We turn off the push and the outflow once the slab reaches ~300
self-consistent experiment
km (time = 5.30 myr). From this time, that we take as our t = to, the de- starts here ....
formation is only driven by the slab-pull force (i.e. ux = 0 on both verti-
cal walls).

The oceanic is pulled towards the continent underneath which it sub-

ducts. The orogenic crust spreads towards the ocean, and a block of oro-
genic crust gets almost separated from the continent. However exten-
sion in the cordillera is short lived from ~5.3 to ~7.8 myr, and the oro- Arrows show
genic block fails to fully detach. The trench retreats towards the ocean, velocity field
however the slab itself is not retreating. From ~7.80 myr onwards the
thinned cordillera is in transient mechanical equilibrium.

An experiment including a buoyant mantle wedge shows that the de-

tachment of the orogenic block is enhanced, but still not completed.
This may be explained by the absence of diking in our numerical ex-
periment, which in nature could be a strong weakening mechanism.

105
Movie 6.3 Stretching of a cordillera orogen

106
 CHAPTER 7

Mantle Geodynamics

Hassan et al. 2016.

Nature, 533, 239-242.

The volume of rocks expands when they become warmer, and contracts when they cool down. This leads to a decrease and
a increase in density, respectively. The coefficient of thermal expansion of most rocks, as a volume fraction per degree, is
about 3e-5 ºC-1. This means that close to the core-mantle boundary, where temperature reaches ~4000ºC, rocks are less
dense by about 12% (3e-5x4000=0.12) compared to rocks closer to the surface. This drives mantle convection, i.e. the
buoyant rise of deep hot mantle rocks, and the sinking of an equivalent volume of colder, denser rocks.

107
SECTION 1

The Rayleigh number

Summary Figure 7.1 Example of mantle flow underneath an extending lithosphere.

1. Why mantle convection

2. Rayleigh number

3. Mantle heating and the

Rayleigh number

The purpose of this section is to explore the nal conditions. We will also illustrate what
dynamic of the convective mantle, to under- role lithospheric plates play on the pattern of
stand the contribution of this process to cool- mantle convection. We start here by introduc-
ing the interior of planets, and - through nu- ing the Rayleigh number, a dimensionless sca-
merical experiments - gain insights into how lar which represents the ratio of the forces that
mantle convection and plate tectonics relate to drive convection over the forces that resist con-
each other. We will test different modes of con- vective motion. In short, the Rayleigh num-
vection driven by various boundary and inter- ber is a measure of the vigor of convection.

108
 from the necessity to minimize internal energy. Convection
Temperature increase with depth, but why? forces deep, hot and less dense material up, and pushes
shallow, cold and denser material down. Because of the re-
In the Earth’s interior, temperature increases with depth for
lationship between density and temperature the Earth's
3 main reasons:
mantle is maintained in a state of permanent gravitational
1/The radioactive decay of radiogenic elements (mainly K40, unstability, as hot material rising toward the surface pro-
U235, U238, U232 and Th) involves exothermic reactions that lib- gressively cool down and become more dense. Heat (a
erate thermal energy. form of energy) is consumed by mantle convection which
2/ During its formation, the Earth has trapped a finite
acts as a giant heat exchanger, to cool down very efficiently
amount of accretional energy (gravitational energy trans-
formed into heat). the deep Earth’s interior.
3/ Last but not least, the Earth is cooled at the surface since
it is suspended into the chilling ~0 K Universe.
Rayleigh Number And The Strength Of Convection

The Rayleigh number measures the vigor of the convection.

Mantle Convection
It is a dimensionless number that expresses the ratio of the
Upon heating, rocks expand and consequently their den- forces driving mantle convection, and the forces opposing
sity decreases. Since temperature increases with depth one it.
concludes that the density of rock decreases with depth as Activity: Try to figure out which physical parameters are rele-
well. However, this is not ideal, as a system in which den- vant to the process of convection? Which of these parameters pro-
sity decreases with depth has a higher internal gravita- mote convection, and which resist convection?
tional potential energy compared to that of a system in
which density increases with depth. Indeed, a fundamental The buoyancy force - that arises from volumetric expan-
principle of physics states that a system naturally evolves sion (thermal expansion) - drives mantle convection. In
to minimize its internal energy. Mantle convection emerges the other hand, the effective viscous force (due to the vis-
cosity of the mantle) and the thermal diffusivity (which

109
smooth out thermal gradient) oppose mantle convection. where K is the thermal conductivity.
The ratio between the buoyancy force to the viscous force
RabT and Rabhf are identical as long as Qm.zm/K = T1 hence
is the Rayleigh number. For convection to occur, the Ray-
Qm = (T1.K)/zm.
leigh number must be larger than a critical value, which de-
pends on whether the surface is free of fixed (i.e. zero veloc-
ity at the surface). The Rayleigh Number driven by internal heating:
ρ0 .g.α.zm 3 H.zm 2
" " " Raih = .( −T0 )
The Rayleigh number takes slightly different forms depend- η.κ K
ing on how heat is delivered to the mantle (e.g. internal where H is the radiogenic heat production per cubic meter.
heating, basal heating or a mix of both). In details, things Rabhf and Raih are identical as long as H = Qm/zm and there-
are complicated but here we simplify and pose that: fore as long as: H = (T1.K)/zm2.

The Rayleigh number driven by a hot base (at temperature The Rayleigh number driven by both basal heat flow and in-
T1) and a cold surface (at temperature T0 at the top of the ternal heating:
convective layer): ρ0 .g.α.zm 3 Qm.zm+H.zm 2
ρ .g.α.zm 3 Ramh = .( −T0 )
" " " RabT = 0 .(T1−T0 ) η.κ K
η.κ
The second term in the Rayleigh numbers describes the heat-
Where ρ0 is the density are room temperature, g is the gravi- ing mode. As long as the first term of the heating mode (T1,
tational acceleration, α is the thermal expansivity, zm is the Qm.zm/K, H.zm2/K, ...) are equals, the Rayleigh numbers are
thickness of the convective layer, η is the average viscosity the same.
of the mantle, and κ is the thermal diffusivity.
For convection to occur, the Rayleigh number must be
The Rayleigh number driven by a hot base (at constant basal larger than a critical value, which varies from 600 to 3000
heat flow Qm) and a cold surface (at temperature T0 at the depending on the boundary conditions (e.g. stagnant lid,
top of the convective layer): plate tectonics, etc).
ρ0 .g.α.zm 3 Qm.zm
" " " Rabhf = .( −T0 )
η.κ K

110
SECTION 2

Exercise: Keeping it cool ...

Summary

1. Learning outcomes: To http://topography.geotheory.co.uk/

develop a deeper
understanding of mantle
dynamic and its dependence
to the mode of heating of
the convective mantle
2. Generic skills: Problem
solving ability,
computational skills,
analytical skills
3. Assumed knowledge: high-
school math, some notions
of tectonophysics (heat and
mass transfer).
4. Tools: R, Python, Matlab
5. Reading: Turcotte &
Schubert: Geodynamics
Exercise 1: Basal heating vs internal heat- To ensure we compare apple with apple we
ing vs mixed heating need to compare models that are at once
thermally and convectively equivalent. We will
Here we compare mantle convection pow-
consider that two models are: 1/ thermally
ered by various modes of mantle heating,
equivalent if in their non-convective state
i.e. basal heating (either constant basal tem-
(infinite viscosity) they have the same
perature or constant heat flow), internal ra-
steady-state average temperature, and 2/
diogenic heating, or a combination of both.

111
convectively equivalent if in their normal convective state spectively; H is the volumetric rate of radiogenic heat pro-
they have the same Rayleigh number. duction.

Find the relationships that sets the conditions for thermal You may need to know that:
equivalence and convective equivalence and comment on α β

( n+1 )
n+1
a.x
∫
your results. We give the following clues: a x n +b dz = +b.x+c
β α
• The steady state conductive geotherm for constant basal
temperature is:
(T −T )
TbT (z) = 1 0 .z+T0 Exercise 2: What’s the Earth’s Rayleigh number?
zm
• The steady state conductive geotherm for constant basal The Earth is currently releasing heat into space at the rate
heat flow is: of ~44 TW. About 20% of this amount (~8 TW) originates
Q from radiogenic isotopes trapped in the plates. This plate
Tbhf (z) = m .z+T0
K component of heat doesn’t contribute to power mantle con-
• The steady state conductive geotherm for internal heat- vection. Therefore, Earth is left with a convective heat flow
ing is: of 36 TW to drive mantle convection.
H.z 2 H.zm
Tih(z) = − + .z+T0 If mantle convection is driven by basal heating alone, then
2K K
the 36 TW must come from the Earth’s core in the form of a
• The steady state conductive geotherm for mixed basal
basal heat flow. In the other hand if mantle convection is
heat flow and internal heating is:
powered by internal heating alone then the 36 TW must be
H.z 2 Qm H.zm
Tmh(z) = − +( + ).z+T0 partitioned into the convective mantle as radiogenic heat
2K K K
production. In nature, we expect the convective heat flow
where Qm, K and T0 are the basal heat flow entering the
to be split into both internal and basal heating.
base of the convective mantle at the core-mantle boundary,
the thermal conductivity and the surface temperature re-

112
There is however a debate about the partitioning of the tion 2” W.m-3; Qm = “your value from question 1” W.m-2; η
heat flow between the core and the mantle. This partition- = 5e21 Pa.s; T0 = 273 K
ing is described by the Urey ratio, the ratio between internal
nb1: For the mixed heating mode, use 24 TW for the internal heat-
heating in the convective mantle to total convective heat
ing and 12 TW for basal heating.
flow (36 TW). Estimates varies from 0.10 to 0.90.
nb2: For the case involving a constant basal temperature one can
1/ Calculate the basal heat flow (in W.m-2) in the case calculate the temperature at the core-mantle boundary by extrapo-
where the convective heat flow is delivered to the convec- lating the linear conductive geotherm in the lithospheric plate
tive mantle in the form of basal heating at the core-mantle (273 K at the surface and 1603K at 150km depth).
boundary.
4/ Assuming no mantle convection, estimate the tempera-
2/ Calculate the rate in radiogenic heat production (in ture at the core-mantle boundary in both basal heating (con-
W.m-3) in the case where the convective heat flow is in the stant heat flow) and internal heating models using the con-
form of internal heating. ductive geotherm equations given earlier. Compare these
Please note that: The Earth’s radius is: 6371 km, the litho- values with the estimated range of temperature of the core-
sphere average thickness is 150 km, the depth of the core- mantle boundary in the convective Earth (3000 to 4000 K),
mantle boundary is at 2890 km. Volume of a sphere: and comment.
4
.π.R 3 ; surface of a sphere: 4.π.R 2 (R is the sphere radius).
3
3/ Calculate the Rayleigh Number for each heating modes
(internal heating, basal heating, and mixed heating) using
the following values:

ρ0 = 4000 kg.m-3; g = 10 m.s-2; α = 3e-5; zm = (2890-150) km;

K = 4 Wm-1 K-1; κ = 1e-6 m2 s-1; H = “your value from ques-

113
SECTION 3

Exercise: Convective patterns

Summary

1. Learning outcomes: To
develop a deeper
understanding of mantle
convection and its
relationship to plate
tectonics
2. Learn to use non-
dimensional approach
3. Generic skills: Problem
solving ability,
computational skills,
analytical skills
4. Assumed knowledge: high-
school math.
5. Tools: Ellipsis, R, Python
6. Reading: Turcotter &
Schubert: Geodynamics Ellipsis is a finite-elements code that solves multigrid solver (solutions are iteratively
complex 2D fluid dynamic problems. It refined by using an increasingly finer grid).
solves the relevant equations for incom- Particles, representing different materials,
pressible fluids (Navier Stokes equations, move through the grid carrying relevant in-
equations of conservation of mass and ener- formation to solve the equations. The com-
gies and constitutive equations) at the bination of a fixed grid (Eulerian grid) and
nodes of a fixed computational grid using a free particles (Lagrangian particles) allows

114
numerical experiments involving large deformation of # Initial conditions
gravacc=1
rheologically complex models. toptbcval=0 bottbcval=1

The Ellipsis input script (a simple text file) describes the # Embedded Temperature
Temp_rect=3 # Nb of rectangle
model including the boundary and internal conditions. El-
Temp_rect_x1=0.00,0.0,0.00 # Coordinates
lipsis then executes the script. Since mantle convection de- Temp_rect_x2=6.00,6.0,6.00
Temp_rect_z1=0.00,0.1,0.90
pends only on the Rayleigh number, the values of all pa-
Temp_rect_z2=0.10,0.9,1.00
rameters entering its definition are of little importance. Temp_rect_mag=0.25,0.5,0.75 # magnitude
Temp_rect_ovl=R,R,R # A/M/R = add/multiply/replace
Hence, in Ellipsis we will set gravity, diffusivity, viscosity,
density, mantle thickness, radiogenic heat production to 1, We fill our box with one material representing the mantle.
excepted for the coefficient of thermal expansion whose This box is defined by the coordinates of its top-left corner
value becomes the Raleigh number. (x1, z1) and bottom-right corner (x2, z2).
With this in mind, let’s build our first model in the case of nb: (0,0) is the top-left corner.
constant basal temperature. We first define the dimension # Material distributions
Material_rect=1 # number of rectangle regions
of the model. The model horizontal length must be a few
Material_rect_property=0 # Material id
times longer than it is depth. Since the depth of the con- Material_rect_x1=0.0 # Coordinates
Material_rect_x2=6.0 # Use columns to add rect.
vecting mantle is 1, we choose 6 for the length for the x di-
Material_rect_z1=0.0
mension. Material_rect_z2=1.0

# Dimensions of the model:

Then we define the boundary conditions. We use a con-
dimenx=6.0
dimenz=1.0 stant temperature of 0 at the top of the model, and a con-
stant temperature of 1 at its base.
For the initial conditions, we have gravitational accelera-
# Boundary conditions
tion set to 1 and the initial geotherm involves three layers
Temp_bc_rect=2 # number of bc
at temperature 0.25 (from 0 to 0.1), 0.5 (from 0.1 to 0.9) and Temp_bc_rect_norm=Z,Z # Z,X for horiz./vert. BC.
0.75 (from 0.9 to 1.0). Temp_bc_rect_icpt=0,1 # Coordinate of boundary plane

115
 Temp_bc_rect_aa1=0,0 # lateral coordinate extent # OUTPUT FILES
Temp_bc_rect_aa2=6,6 # lateral coordinate extent datafile="ConvMod1" # root name of output files
Temp_bc_rect_hw=0,0 # smoothing half-width datatypes="Temp,Pres" # output variables at nodes
Temp_bc_rect_mag=0,1 # magnitude of bc averages="Temp" # horizontally averaged values,
storage_timesteps=1 # writing interval
Heat_flux_z_bc_rect=0 # no basal heat flow checkpt_timesteps=1 # PPM writing interval
storage_timing=1.e-8# timestep writing interval
In what follow, we define the thermal and mechanical prop- checkpt_timing=1.0 # PPM writing interval (time)

erties of the mantle. Again, because all relevant parame- # Specifications for graphical output files
ters are included in the Rayleigh number, density, viscosity, PPM_files=1 # number of PPM files at each output step
PPM_height=256 # vertical size of output PPM file
diffusivity and heat capacity are equal to 1. Because we con- PPM_coloring=1 # variable upon which to base colouring
sider a basally heated model we set the internal heating to # 1=temperature, 2=viscosity, 3=stress, ...

0, and the coefficient of thermal expansion is the Rayleigh Finally, we may want to extract some temperature (or other
Number. variables) profiles across our model:
# MATERIAL PROPERTIES
# for extracting profiles/histories
Different_materials=1 # Nb of materials
Sampling_tracers=10 # number of sampling tracers

## Mantle
Material_0_density=1 Tracers (fixed to the grid or to the material) can record in-
# Rheological model formation such as Temperature, Pressure, etc. They can also
Material_0_rheol_T_type=1# visc(T)=N0*exp(-T1*T)
Material_0_viscN0=1 # N0 in viscosity model record parameter profiles (vertical or horizontal).
Material_0_viscT1=0 # T1 in viscosity model
# Thermal parameters One last thing, to run this model over 1600 time steps. #
Material_0_therm_exp=7.34e7 # coeff. thermal exp
# Rayleigh Nb
Material_0_therm_diff=1 # thermal diffusivity ADVECTION-DIFFUSION PARAMETERS
Material_0_Cp=1.0 # isobaric heat capacity
Material_0_Qt=0.0 # Rad. heating rate / kg maxstep=1200 # maximum number time steps

# Colouring (Try to experiment with this.) In a terminal navigate to the folder holding the input scritp
We now consider the output files we would like to get an executable and run the following command:
from Ellipsis: ./ellipis3d my_ellipsis_model.input

116
Exercise 3: 2D Experiments of Mantle Convection using. with:
Material_0_viscN0=5 # N0 in viscosity models
1/ Run over 1200 steps the previous script (convection
Material_0_viscT1=5 # T1 in viscosity models
with basal heating, constant temperature).
Verify that the average viscosity of the convective mantle is
2/ Build and run a model with internal heating
still 1.
(Ra=2.62e8). For this model, you will need to 1/ remove
the constant temperature at the base of the model, 2/ set to 4/ Analyze and comment your results. Give particular at-
zero the magnitude of the basal heat flow (insulation condi- tention to the shape of the averaged geotherm for both in-
tion, no heat is lost at the base of the model), and 3/ set the ternally heated and basal heated models. How do the geo-
internal heating rate to 1. therms compare?

5/ Design a model with both internal heating and basal

# Boundary conditions
Temp_bc_rect=1 # nb of bc. heating using constant basal heat flow instead of constant
Temp_bc_rect_norm=Z,Z # Z (horiz), X (vert.) BC
basal temperature (use the Rayleigh number from exercise
Temp_bc_rect_icpt=0,1 # Coordinate of boundary plane
Temp_bc_rect_aa1=0,0 # lateral coordinate extent 2). Run this model and compare it with a model in which a
Temp_bc_rect_aa2=6,6 # lateral coordinate extent
lithospheric plate is added. The plate has a thickness of 0.1
Temp_bc_rect_hw=0,0 # half-width of bc smoothing edge
Temp_bc_rect_mag=0,1 # magnitude of bc and extends over 50% of the model. Its density must be
20% lower than the mantle and its viscosity (constant with
Heat_flux_z_bc_rect=1 # nb bc ranges (surfaces)
Heat_flux_z_bc_rect_norm=Z # normal to boundary temperature) must be 1000 times larger than that of the
Heat_flux_z_bc_rect_icpt=1 # normal-axis intercept
Heat_flux_z_bc_rect_aa1=0 # lateral coord. extent
mantle. The plate can be maintained fixed or let to move
Heat_flux_z_bc_rect_aa2=6 around. You may want to experiment with both options.
Heat_flux_z_bc_rect_hw=0 # smoothing half-width
Heat_flux_z_bc_rect_mag=0 # magnitude of bc
For a fixed plate, use the following (change X1 and X2 by
3/ Run the models 1 and 2 using a temperature dependent the horizontal coordinates for your plate):
viscosity:
Velocity_z_bc_rect=1 # nb of bc (surfaces)
Velocity_z_bc_rect_norm=X # normal to bc pla

117
 Velocity_z_bc_rect_icpt=0 # coord. of the bc plane
Velocity_z_bc_rect_aa1=X1 # lateral coordinate
Velocity_z_bc_rect_aa2=X2 # lateral coordinate
Velocity_z_bc_rect_hw=0 # smoothing half-width
Velocity_z_bc_rect_mag=0 # magnitude of bc

6/ How plates influence the convection, the average tem-

perature and the temperature distribution in the mantle?

118
SECTION 4

The geodyamics of mantle melting

Basaltic volcanism is the result of partial melting of the Earth’s compression melting in small-scale convection cells along thick
mantle (Fig. 1A). Volcanic activity in arcs and along mid-ocean lithospheric keels has been proposed to explain minor intraplate
ridges and continental rifts is related to plate tectonics, through volcanism (Davies and Rawlinson, 2014).
decompression melting of the shallow sublithospheric mantle in
One of the main differences between plate-tectonics–related
regions of lithospheric extension and thinning, and/or hydra-
melting of the sublithospheric mantle and partial melting in hot
tion of the mantle wedge above subducting slabs. In contrast,
upwelling mantle material is the pressure (P) and temperature
intraplate volcanism, responsible for thick continental flood ba-
(T) at which the adiabat (i.e., geotherm in the convective man-
salts, has usually been linked to upwelling and decompression
tle) intersects the mantle solidus (Fig. 1B). These cannot be ob-
melting of large volumes of deep, hot mantle (e.g., mantle
tained directly from the composition of basalts because they are
plume; Morgan, 1983; Richards et al., 1989). More recently, de-
affected by fractional crystallization and interaction with sur-

119
rounding rocks, as well as the progressive evolution of the composi- A B
1200
Temperature (ºC)
1400 1600 1800 2000
1 2
0
tion of the source. The geochemistry of the first-formed batch of 0
Continental crust

Liquid
50

Depth (km)
melt as fertile mantle crosses its solidus (i.e., the primary melt) is a Lithospheric mantle a’

us
15
100 50

MgO
So

Melt
lid
function of the mantle potential temperature, Tp , as shown by for-

us
a

%
150 Con

adiab
ducti
Tectonic thinning and ve b’
ward thermodynamic modeling (Asimow, 1999; Asimow et al., 200 exhumation of 100 geo
the

at
rms
sublithospheric mantle Deep mantle 2
250
2001; Herzberg and O’Hara, 2002; Herzberg, 2004). Tp is the tem- upwelling

Adiabat
Depth (km)
1
Asthenosphere 150 3
perature measured at the Earth’s surface from extrapolation of the
Super continent

25
b
mantle adiabat (McKenzie and Bickle, 1988), a convenient way to

30
3 200 Sublithospheric Mantle
mantle warming

35
express the intersection temperature of adiabat and solidus. The Mantle
warming mantle
exhumation
ive
ect 250 Deep
composition of the elusive primary melt of a volcanic province de- nv

Adiabat
mantle

Co
upwelling
rives from processing the most magnesian liquids (i.e., parental 300
melts) inferred from a suite of basalts (Herzberg et al., 2007). Com-
parison of the calculated and modelled primary melt constrains the Fig. 1. A: Geodynamic settings of mantle melting. B: Phase diagram, contoured for
MgO, of mantle melting. 1—lithospheric thinning; 2— mantle upwelling; 3—mantle
Tp from which the geodynamic setting of the partial melting can be
warming. During (1) and (2), the convective mantle is exhumed adiabatically (blue
inferred. By targeting primary melts, petrologists mitigate geo- and red dashed lines in B) intersecting the solidus at ~70 km and 180 km, respec-
chemical alteration of parental magmas through fractional crystalli- tively, with Tp (1330 ºC and 1620 ºC) shown by the blue and red dots at Earth’s sur-
zation and interaction with surrounding rocks during melt ascent face. Following mantle warming (3), the adiabat (black dashed line) intersects the
(e.g., Herzberg et al., 2007). Keeping in mind the many pitfalls solidus at ~120 km depth, with Tp (~1470 °C) shown by the gray dot. Primary melts
are produced at solidus temperatures (blue, gray, and red stars). Should these melts
(Herzberg et al., 2007; Herzberg, 2011), data from present continen-
be extracted as soon as they are produced, they would follow their respective melt
tal rifting and spreading ridges derive a Tp of less than 1400 °C (Fal-
adiabat (solid blue, gray, and red arrows). In the case of continental rifting leading to
loon et al., 2007; Herzberg et al., 2007). Hole (2015) shows that data a mid-ocean ridge, melt is usually extracted at melt fraction of <1% (i.e., fractional
from plume-related flood basalts point to a Tp of more than 1500 melting). As latent heat escapes with the melt, the source rock cools and follows the
°C, in agreement with, for example, Thompson and Gibson (2000), solidus during exhumation (path a-a’). As the source becomes more residual, it pro-
Herzberg and O’Hara (2002), and Herzberg et al. (2007). The range duces parental melts with a range of compositions (empty blue arrows). In plumes,
melt continuously reacts with the residue, with which it remains in equilibrium. As
of adiabats during tectonic extension and exhumation of the shal-
latent heat remains with the melt in the source, the upwelling mantle keeps follow-
low mantle (blue region in Fig. 1B) and during deep mantle upwel- ing the adiabat (path b-b’). As the plume head slows down, melt segregates toward
ling (red region in Fig. 1B), do not overlap. The concept of Tp can the top of the partially melted column, before escaping to the surface (i.e., batch or
thus be used to constrain the geodynamic setting of mantle melt- equilibrium melting) following the parental melt adiabat (empty red arrow). The
ing. Let’s now consider the origin of the ca. 200 Ma Central Atlan- pressure at which melt is extracted is the final melting pressure. Both fractional melt-
ing and equilibrium melting can occur in succession depending on strain rate, perme-
tic magmatic province (CAMP).
ability, and the relative upwelling velocities of the residue and melt.

120
The CAMP extends over Africa, South America, Europe, and North pends on plate tectonics: continents amalgamate to create supercon-
America, covering ~107 km2 (Marzoli et al., 1999). It post-dates by tinents through the Wilson cycle, then break into smaller continen-
50 m.y. the final stages of Pangea’s assembly. The peak igneous ac- tal plates. Two- and three-dimensional numerical studies show that
tivity at 199–200 Ma pre-dates by ~10 m.y. the initial opening of the the Tp of the convective mantle can increase 10% to 15% (~130–200
central Atlantic Ocean (Marzoli et al., 1999; Sahabi et al., 2004). The °C) as plates aggregate into a supercontinent covering 20% to 35%
attribution of the CAMP to a mantle plume-head impinging the of Earth’s surface (Coltice et al., 2009), thus the Tp of partial melt-
base of the lithosphere (Hill, 1991; Wilson, 1997; Leitch et al., 1998; ing following the formation of supercontinents is between that of
Courtillot et al., 1999; Janney and Castillo, 2001; Ernst and Buchan, plate tectonic processes (<1400 ºC) and deep mantle upwelling
2002) has been challenged on the basis that chemical and isotopic (>1500 ºC). The potential temperature determined by Hole from
compositions of basalts point to shallow-mantle sources, including the CAMP’s primary magmas using the most advanced model
lithospheric ones, enriched by earlier subduction (Bertrand et al., (PRIMELT2), combined with the observation that magmatism fol-
1982; Bertrand, 1991; Puffer, 2001; Deckart et al., 2005; Verati et al., lowed supercontinent aggregation, confirms that the CAMP was
2005). However, the origin of the CAMP is still strongly debated be- the result of large-scale mantle melting following the aggregation
cause basaltic magmas change on their way to the surface. of the last supercontinent (Coltice et al., 2009).

Hole (2015) processes the composition of parental magmas of the REFERENCES CITED
CAMP to recover that of their primary melts, and concludes most Asimow, P.D., 1999, A model that reconciles major- and trace-element data
CAMP primary magmas point to a Tp between 1400 °C (olivine nor- from abyssal peridotites: Earth and Planetary Science Letters, v. 169, p. 303–
mative basalts) and 1500 °C (quartz normative basalts). This tem- 319, doi:10.1016/S0012-821X(99)00084-9.
perature range is too hot for sublithospheric decompression melt- Asimow, P.D., Hirschmann, M.M., and Stolper, E.M., 2001, Calculation of peri-
ing related to lithospheric thinning, and too cold for a source in a dotite partial melting from thermodynamic models of minerals and melts, IV.
hot upwelling mantle. Adiabatic decompression and the composition and mean properties of mid-
ocean ridge basalts: Journal of Petrology, v. 42, p. 963–998,
Continental lithosphere is thicker and relatively more stable than doi:10.1093/petrology/42.5.963.

oceanic lithosphere, and thus impedes the removal of heat from the Bertrand, H., 1991, The Mesozoic tholeiitic province of northwest Africa: A
hotter convective mantle. As continents aggregate into superconti- volcanotectonic record of the early opening of the central Atlantic, in Kam-
punzu, A.B., and Lubala, R.T., eds., Magmatism in Extensional Structural Set-
nents, they enlarge the wavelength of the convective pattern, fur-
tings: The Phanerozoic African Plate: Berlin, Springer-Verlag, p. 147–188.
ther impeding heat loss, thus forcing large-scale warming of the
sublithospheric mantle (Grigné et al., 2005). Modulation of mantle Bertrand, H., Dostal, J., and Dupuy, C., 1982, Geochemistry of early Mesozoic
tholeiites from Morocco: Earth and Planetary Science Letters, v. 58, p. 225–
temperature below continents depends on the size, number, and 239, doi:10.1016/0012-821X(82)90196-0.
distribution of continental plates (Coltice et al., 2007), and thus de-

121
Coltice, N., Phillips, B.R., Bertrand, H., Ricard, Y., and Rey, P.F., 2007, Global Herzberg, C., and O’Hara, M.J., 2002, Plume-associated ultramafic magmas of
warming of the mantle at the origin of flood basalts over supercontinents: Ge- Phanerozoic age: Journal of Petrology, v. 43, p. 1857–1883,
ology, v. 35, p. 391–394, doi:10.1130/G23240A.1. doi:10.1093/petrology/43.10.1857.

Coltice, N., Bertrand, H., Rey, P.F., Jourdan, F., Phillips, B.R., and Ricard, Y., Herzberg, C., Asimow, P.D., Arndt, N., Niu, Y., Fitton, J.G., Cheadle, M.J., and
2009, Global warming of the mantle beneath continents back to the Archaean: Saunders, A.D., 2007, Temperatures in ambient mantle and plumes: Con-
Gondwana Research, v. 15, p. 254–266, doi:10.1016/j.gr.2008.10.001. straints from basalts, picrites, and komatiites: Geochemistry Geophysics Geo-
systems, v. 8, doi:10.1029/2006GC001390.
Courtillot, V., Jaupart, C., Manighetti, I., Tapponnier, P., and Besse, J., 1999,
On causal links between flood basalts and continental breakup: Earth and Hill, R.I., 1991, Starting plumes and continental break-up: Earth and Plane-
Planetary Science Letters, v. 166, p. 177–195, tary Science Letters, v. 104, p. 398–416, doi:10.1016/0012-821X(91)90218-7.
doi:10.1016/S0012-821X(98)00282-9.
Hole, M.J., 2015, The generation of continental flood basalts by decompres-
Davies, R., and Rawlinson, N., 2014, On the origin of recent intraplate volcan- sion melting of internally heated mantle: Geology, v. 43, p. 311–314,
ism in Australia: Geology, v. 42, p. 1031–1034, doi:10.1130/G36093.1. doi:10.1130/G36442.1.

Deckart, K., Bertrand, H., and Liegeois, J.P., 2005, Geochemistry and Sr, Nd, Janney, P.E., and Castillo, P.R., 2001, Geochemistry of the oldest Atlantic oce-
Pb isotopic composition of the Central Atlantic Magmatic Province (CAMP) anic crust suggests mantle plume involvement in the early history of the cen-
in Guyana and Guinea: Lithos, v. 82, p. 289–314, tral Atlantic Ocean: Earth and Planetary Science Letters, v. 192, p. 291–302,
doi:10.1016/j.lithos.2004.09.023. doi:10.1016/S0012-821X(01)00452-6.

Ernst, R.E., and Buchan, K.L., 2002, Maximum size and distribution in time Leitch, A.M., Davies, G.F., and Wells, M., 1998, A plume head melting under a
and space of mantle plumes: Evidence from large igneous provinces: Journal rifting margin: Earth and Planetary Science Letters, v. 161, p. 161–177,
of Geodynamics, v. 34, p. 309–342, doi:10.1016/S0264-3707(02)00025-X. doi:10.1016/S0012-821X(98)00147-2.

Falloon, T.J., Danyushenvsky, L.V., Ariskin, A., Green, D.H., and Ford, C.E., Marzoli, A., Renne, P., Piccirillo, E., Ernesto, M., Bellieni, G., and De Min, A.,
2007, The application of olivine geothermometry to infer crystallization tem- 1999, Extensive 200-million-year-old continental fl ood basalts of the Central
peratures of parental liquids: Implications for the temperature of MORB mag- Atlantic Magmatic Province: Science, v. 284, p. 616–618,
mas: Chemical Geology, v. 241, p. 207–233, doi:10.1126/science.284.5414.616.
doi:10.1016/j.chemgeo.2007.01.015.
McKenzie, D., and Bickle, M.J., 1988, The volume and composition of melt
Grigné, C., Labrosse, S., and Tackley, P.J., 2005, Convective heat transfer as a generated by extension of the lithosphere: Journal of Petrology, v. 29, p. 625–
function of wavelength: Implications for the cooling of the Earth: Journal of 679, doi:10.1093/petrology/29.3.625.
Geophysical Research, v. 110, B03409, doi:10.1029/2004JB003376.
Morgan, W.J., 1983, Hotspot tracks and the early rifting of the Atlantic: Tec-
Herzberg, C., 2004, Geodynamic information in peridotite petrology: Journal tonophysics, v. 94, p. 123–139, doi:10.1016/0040-1951(83)90013-6.
of Petrology, v. 45, p. 2507–2530, doi:10.1093/petrology/egh039.
Puffer, J.H., 2001, Contrasting high field strength element contents of conti-
Herzberg, C., 2011, Basalts as temperature probes of Earth’s mantle: Geology, nental flood basalts from plume versus reactivated-arc sources: Geology,
v. 39, p. 1179–1180, doi:10.1130/focus122011.1. v. 29, p. 675–678, doi:10.1130/0091-7613(2001)029<0675:CHFSEC>2.0.CO;2.

122
Richards, M.A., Duncan, R.A., and Courtillot, V., 1989, Flood basalts and hot-
spot tracks, plume heads and tails: Science, v. 246, p. 103–107,
doi:10.1126/science.246.4926.103.

Sahabi, M., Aslanian, D., and Olivet, J.L., 2004, A new starting point for the
history of the central Atlantic: Comptes Rendus Geoscience, v. 336, p. 1041–
1052, doi:10.1016/j.crte.2004.03.017.

Thompson, R.N., and Gibson, S.A., 2000, Transient high temperatures in man-
tle plume heads inferred from magnesian olivines in Phanerozoic picrites: Na-
ture, v. 407, p. 502–506, doi:10.1038/35035058.

Verati, C., Bertrand, H., and Féraud, G., 2005, The farthest record of the Cen-
tral Atlantic Magmatic Province into West Africa craton: Precise 40Ar/39Ar
dating and geochemistry of Taoudenni basin intrusives (northern Mali):
Earth and Planetary Science Letters, v. 235, p. 391–407,
doi:10.1016/j.epsl.2005.04.012.

Wilson, M., 1997, Thermal evolution of the Central Atlantic passive margins:
continental break-up above a Mesozoic super-plume: Journal of the Geologi-
cal Society, v. 154, p. 491–495, doi:10.1144/gsjgs.154.3.0491.

123
SECTION 5

Modelling mantle melting

Summary

1. Solidus/liquidus
2. Basalts MgO content
3. Melting during fast
decompression
4. Melting during slow
decompression
5. Exercise: Melting in a mantle
plume

The mantle solidus and liquidus are relatively well constrained. For instance, according to Katz et al., (2003):
Tsol(z): P<15GPa = 1358.7 + 132.899 (ρ.g.z)/(1e9)- 5.104 ((ρ.g.z)/(1e9))2
Tsol(z): P≥15GPa = 1510.76 + 46.27 (ρ.g.z)/(1e9)- 0.8036 ((ρ.g.z)/(1e9))2
Tliq(z): P<15GPa = 2053 + 45 (ρ.g.z)/(1e9)- 2.00 ((ρ.g.z)/(1e9))2
Tliq(z): P≥15GPa = 1470.3025 + 55.53 (ρ.g.z)/(1e9)- 0.9084 ((ρ.g.z)/(1e9))2
Melt fraction F is calculated using the concept of super-solidus (Tss) and eq. 21 of McKenzie and Bickle (1988):
F(Tss) = 0.3936 + 0.8936 Tss + 0.4256 Tss 2 + 2.988 Tss 3
Basalts can be classified based on their MgO content: Komatiite (MgO>18%), komatiitic basalts (12%<MgO<18%), tholei-
itic basalts (6%<MgO<12%). MgO content depends on the pressure at which melt is extracted (Herzberg and Gazel, 2009):
Tc(MgO, z) = 935 + 33 MgO - 0.37 MgO2 + 54 (ρ.g.z) - 2 (ρ.g.z)2

124
Mantle melting during fast exhumation: The graph below shows a phase diagram with melt fraction mapped between the
solidus and the liquidus. This graph describes what happened during a phase of fast homogeneous continental rifting
back in the Archean when the mantle was ~200ºC hotter than present. The continental lithosphere (pale green) is thinned
rapidly (225 km to 110 km), which triggers the decompression of the mantle below the continent. When extension and thin-
ning is rapid each point on the geotherm follows the adiabatic gradient (i.e. no conductive cooling). As the geotherm be-
comes warmer, it intersects the solidus at two locations
Temperature (K)
(A0’ and A3’) which define the top and base of a par- 1400 1600 1800 2000 2200
0
tially molten mantle column. The base of the column re- Tpot = 1820 K

Liqui
So

0
lid

Melt
mains close to 150 km throughout the volcanic history.

us
20 40 60 80

dus
Geoth
erm (

adiab
Because extension is rapid, the partially molten column 50
t1 )

at
(i.e. solid matrix together with the melt) follows also A3’ A2’
Ge
the adiabat since there is no melt extraction, and there- oth

2
erm
(t0
fore no loss of heat. After ~100% extension, points on 100 )
LAB at t1

Pressure (GPa)
A1’
the initial geotherm located at A0 to A3 have moved up-
A3 melt %
ward to location A0’ to A3’. The partially molten man- A0’

4
150

Depth (km)
tle column consists of partially molten fertile sub- A2 Exhumation
SCLM of the LAB
lithospheric mantle (from ~150 to ~110 km), and par-

After extension
tially molten SCLM (from ~110 to ~55 km). Because the 200

6
melt accumulates in the column, the buoyancy of the
molten column increases and its viscosity decreases, un- A1 LAB at t0
Before extension

Adiabatic path
til melt escapes to the surface to form a basaltic crust. 250
Fertile
The maximum thickness of this crust is given by the Sub-Lithospheric mantle

8
(SLM)
area between the solidus and the geotherm (purple tri-
300
angle), here it is ~ 16 km (i.e. ~ 0.5 x (depth(A0’) - A0
depth(A3’)) * maximum melt fraction).

125
 Temperature (K)
The phase diagram on the right shows the MgO content between the solidus and 1400 1600 1800 2000 2200

the liquidus, according to Herzberg and Gazel (2009). The MgO content of melts ex-

Liqui
So
lid

Melt
tracted at various locations in the partially molten column ranges from 9 to 21 %.

u
Geoth 10 20 30

dus
s
erm (

adiab
t1 )

at
A ’
Mantle melting during slow exhumation: When extension is slow (graphs below), A3’ 2

the solid matrix flows upwards very slowly, and the melt - whose buoyancy in- 100
MgO %
A1’
creases as the pressure at which it is produced decreases - flows faster through the
porous source and ponds at the top of the partially molten column, before escap- A0’
150

Depth (km)
ing to the surface via diking. In the molten column, melt forms an interconnected
network for a critical melt fraction as low as 0.1 %, limiting the amount of melt in
200
the column to less than 1%. While the flowing melt follows the melt adiabat, the de-
pleted residue follows the solidus as latent heat escapes with the melt. The total
250
Temperature (K) Temperature (K) amount of
1400 1600 1800 2000 2200 1400 1600 1800 2000 2200
latent heat
Melt

Liqui
So
Liqui
So

lid
lid

300
consumed
us

10 20 30
us

adiab

40 60 80 Geoth

dus
dus

Geoth erm (
erm ( 20 t1 )
t1 )
at

50
depends only on the volume of rock advected
Melt
A3’ A3’

adiab
A2’

at MgO % above the solidus, not on the history of melt extrac-

100 LAB at t1
100

A1’ tion. Hence, the melt productivity, and therefore

melt %

150
A0’
150 A0’
the depletion of the source, is not significantly af-
Depth (km)
Depth (km)

A2 Exhumation
of the LAB fected by the residence time of the melt in the
SCLM
200 200 source before its extraction. Assuming 100% exten-
sion, the thickness of the basaltic crust also
A1 LAB at t0
reaches ~16 km. The MgO content of the melt var-
Adiabatic path

250 250
Fertile
Sub-Lithospheric mantle
(SLM)
ies with the pressure at which it is produced and
300
A
300
varies between 9 and 17%.
0

126
The small fraction of melt remaining in the source has only a modest impact on the buoyancy of the partially molten re-
gion. Its impact on the rheology is, however, more significant. On one hand, the viscosity of the partially molten mantle
strongly decreases because melt lubricates grain boundaries. On the other hand, partial melting drains water out of the
solid matrix and increases the olivine content. Although the impact of dehydration on the viscosity may not be as signifi-
cant as previously thought (e.g., Fei et al., 2013, Nature 498), both processes should contribute to increase the viscosity of
the depleted residue once its temperature drops below the solidus.

Exercise - Melting in a mantle plume: We model the rise of a mantle plume in a box 700 km deep and 4200 km wide (i.e.
simulating the upper mantle). The plume is modelled using a positive thermal anomaly 300ºC higher than the surround-
ing mantle (1347ºC). We take advantage of the symmetry of the problem and embed the anomaly on the lower left corner
of the box. The mantle is modelled as a visco-plastic material with temperature, stress, strain-rate and melt dependent vis-
cosity. The viscous rheology of the mantle is based on dry olivine (Brace & Kohlstedt, 1980). Its plastic rheology is approxi-
mated with a Coulomb failure criterion with cohesion 40 MPa, a coefficient of friction ((ref ) 0.268, and a weakening factor
dependent on the accumulated plastic strain (εp). The coefficient of friction (() evolves as a function of the accumulated
plastic strain as:

for εp < 0.15: ( = (ref ×(1-(1-0.0373)×( εp /0.15)0.25)

for εp > 0.15: ( = (ref × 0.0373

This weakening leads to a strong strain localisation simulating faults with nominal viscosity ~25 times weaker than that of
the surrounding rocks when plastic strain reaches 15%. This weakening, which is a key factor for the operation of plate tec-
tonics on present Earth, simulates the formation of phyllosilicates (serpentine, talc, micas) during the strain-induced hydra-
tion of mantle rocks and basalts. For semi-brittle deformation independent of pressure, we impose an upper limiting yield
stress of 300 MPa (Ord & Hobbs, 1980; Escartin & Hirth, 1997).

The top 15 km of the model consists of a thick layer of basalts, which simulates the oceanic crust. Because these basalts are
emplaced below sea level, they are strongly hydrothermally altered. For simplicity, we assume that this layer of basalt has

127
a nominal viscosity 1000 times weaker that the underlying mantle, a weak cohesion of 1 MPa, a reference coefficient of fric-
tion of 0.134 and the same weakening properties as the rest of the mantle.

In our modelling of a rising plume, we assume that melt in the partially molten column is extracted when the melt fraction
reaches 1%. In this case, latent heat escapes with the melt, and the ascending depleted residue follows the solidus closely.
We assume the fusion entropy to be 400 J kg-1 K-1.

The small fraction of melt (<1%) has only a modest impact on the buoyancy of the partially molten region. However, the
density of the residue decreases as it becomes more depleted. We assume a maximum density decrease of 1.5% upon full
depletion. Because even a small fraction of melt lubricates grain boundaries, it affects the viscosity of the partially molten
column. Hence, we impose a linear viscosity drop to a maximum of one order of magnitude when the melt fraction
reaches 1%. On the other hand, partial melting drains water out of the solid matrix and reduces the number of phases. Al-
though the impact of dehydration on the viscosity may not be as significant as previously thought (Fei et al., 2013), both
processes should contribute to increase the viscosity of the depleted residue once its temperature drops below the solidus.
In our experiment we impose an increase in viscosity proportional to depletion, assuming a viscosity increase of two or-
ders of magnitude for 100% depletion.

The geotherm is linear in the lithosphere (T(0km)=293 K, T(150km)=1620K). In the convective mantle, the geotherm fol-
lows the adiabatic gradient (~0.3 ºC per km).

Run the model (A_Plume100.input), analyze the result, and summarize your findings into a 4-page report (less than 3-4
Mbyte). For this you may want to consider the information stored in the following files (all readable with a simple text edi-
tor): .node_data, .profiles, .samples and .timelogs.

128
References
Brace, W. F., Kohlstedt, D. L. Limits on lithospheric stress imposed
by laboratory experiments. Journal of Geophysical Research 85,
6248-6252 (1980).

Escartin, J. & Hirth, G. Nondilatant brittle deformation of serpen-

tinites: Implications for Mohr-Coulomb theory and the strength of
faults. Journal of Geophysical Research 102, 2897-2913 (1997).

Fei, H., Wiedenbeck, M., Yamazaki, D. & Katsura, T. Small effect of

water on upper-mantle rheology based on silicon self-diffusion co-
efficients. Nature 498, 213-216 (2013).

Herzberg, C. & Gazel, E. Petrological evidence for secular cooling

in mantle plumes. Nature 458, 619-623 (2009).

Katz, R. F., Spiegelman, M. & Langmuir, C. H. A new parameteriza-

tion of hydrous mantle melting. Geochemistry, Geophysics Geosys-
tems 4 , doi:10.1029/2002GC000433 (2003).

McKenzie, D. & Bickle, M. J. The volume and composition of melt

generated by extension on the lithosphere. Journal of Petrology 29,
625-629 (1988).

Ord, A., Hobbs, B. E. The strength of the continental crust, detach-

ment zones and the development of plastic instabilities. Tectono-
physics 158, 269- 289 (1989).

129
 Temperature (K) Temperature (K)
1400 1600 1800 2000 2200 1400 1600 1800 2000 2200

0
0

Liqui
Li qui

So
So
lid

lid
us

u
20 40 60 80 10 20 30

dus

dus
s
50

2
Pressure (GPa)
melt %

Pressure (GPa)
100 MgO %
100

4
150
Depth (km)

150

Depth (km)
6
200

6
200

250
250
8

8
Solidus and liquidus according to Katz et al., Solidus and liquidus according to Katz et al.,
2003. Melt fraction from McKenzie & Bickle 2003. MgO content from Herzberg & Gazel
300 (1988) (2009)
300

130
View publication stats