The mantle, the core and magma oceans in the thermal
evolution of the EarthStéphane Labrosse1 , Adrien Morison1 , Daniela Bolrão2 , Roberto Agrusta1 , R. Deguen1 , T.
Alboussière1 , A. Rozel2 , M. Ballmer2 , P. J. Tackley2
1 LGL-TPE, Université Claude Bernard Lyon 1, Ens de Lyon, CNRS
2 ETH Zürich

FDEPS, Nov, 2019

Principle for thermal evolution models
QI Global energy balance
MCPdT= −Q(T ) + H(t)dtI Evolution time scale
τ =T
I The Urey number (time dependent)
Ur = 100
I Question:
I what is the function Q(T ) for mantle convection?

HQ

MCP TQ(T )Boundary layer scaling for thermal evolution model

I Boundary layer scaling⇒q=Ck∆T 1/3Rad

Ti
∆T

4/3

I Assuming the scaling applies to mantle convection, hence to the present time (t = 0): Q0 , T0 .
I Scaling of surface heat flow:

Q(T ) = Q0TT0

4/3 

η(T )
η0

−1/3

I Given the present-day Urey number: Ur = 100 × H/Q ∈ [20% − 50%].
I Solve the energy equation backward from the present time:MCdT= H(t) − Q(T )dtThe low Urey number ”paradox”
PowerHeat flow, Q

Q(T)

H(0)
Time, t01300˚C

dT= H(t) − Q(T );dt

Q(T ) = Q0-4.5Ga

MCTemperature, T

Urey number: Ur = 100H/Q

TT0

4/3 

η(T )
η0

−1/3

The low Urey number ”paradox”
PowerHeat flow, Q

Q(T)H(-dt)
Time, t-4.5Ga01300˚C

dT= H(t) − Q(T );dt

Q(T ) = Q0-dt

MCTemperature, T

Urey number: Ur = 100H/Q

TT0

4/3 

η(T )η0−1/3The low Urey number ”paradox”
PowerHeat flow, Q

Q(T)

H(-2dt)Time, t-4.5Ga-2dt

MCTemperature, T01300˚CdT= H(t) − Q(T );dt

Q(T ) = Q0-dt

Urey number: Ur = 100H/Q

TT0

4/3 

η(T )
η0

−1/3

Feedback with temperature-dependent viscosity

Temperature (K)

Hrad

300
1600
200
1400

1200

100

−4−3−2Age, Gyr−10−4−3−2Age, Gyr−10

Radiogenic Power, Qsurf (TW)

Ur=43

1800

400Qsurf

Feedback with temperature-dependent viscosity

Temperature (K)

HradUr=87.20606
Ur=87.2065

300

1600
200
1400

1200

100−4−3−2Age, Gyr−10−4−3−2Age, Gyr−10

Radiogenic Power, Qsurf (TW)

Ur=431800

400Qsurf

Proposed solutions

I Layered mantle convection =⇒ lower cooling efficiency (McKenzie & Richter, 1981). But theinferred lower mantle temperature is too large to keep it solid unless it is depleted in heatproducing elements (Schubert & Spohn, 1981; Spohn & Schubert, 1982).
I A smaller exponent β in the q = ARaβ scaling law (Christensen, 1985; Conrad & Hager, 1999;
Sleep, 2000; Korenaga, 2003).
I Differential core-mantle cooling (another type of layered model).

Proposed solutions

I Layered mantle convection =⇒ lower cooling efficiency (McKenzie & Richter, 1981). But theinferred lower mantle temperature is too large to keep it solid unless it is depleted in heatproducing elements (Schubert & Spohn, 1981; Spohn & Schubert, 1982).
I A smaller exponent β in the q = ARaβ scaling law (Christensen, 1985; Conrad & Hager, 1999;
Sleep, 2000; Korenaga, 2003).
I Differential core-mantle cooling (another type of layered model).

Proposed solutions

I Layered mantle convection =⇒ lower cooling efficiency (McKenzie & Richter, 1981). But theinferred lower mantle temperature is too large to keep it solid unless it is depleted in heatproducing elements (Schubert & Spohn, 1981; Spohn & Schubert, 1982).
I A smaller exponent β in the q = ARaβ scaling law (Christensen, 1985; Conrad & Hager, 1999;Sleep, 2000; Korenaga, 2003).
I Differential core-mantle cooling (another type of layered model).

A lower β exponent to decrease the feedback?
Christensen (1985):⇒ A decreased feedback between the mantle temperature and the surface heat flow.

Effect of the β exponent on thermal evolutionQ(T ) = Q0TT0

1+β η(T )η0−βNu=RaβQ
4343

β=0.0

2
1

β=0.0521T/°C
1700
1600

28=Urey#

1500
1400

432Ga BP01700
1600

43

1500
1400
1Q
43

β=0.1021

1700
1600

49

43

1500
1400

56
43

2410(Christensen, 1985)

β=0.3021T/°C
1700
1600

58

1500
1400

64
3210484.5
84.5
32Ga BPSame solution was advocated by Conrad & Hager (1999, 2001), Sleep (2000), Korenaga (2003).
Problem: No self-consistent dynamical model gives such low values of β10

Effect of the β exponent on thermal evolutionQ(T ) = Q0TT01+β 

η(T )
η0

−β

Nu=RaβQ
434
3

β=0.0

2
1

β=0.0521T/°C
1700
1600

28=Urey#

1500
1400

432Ga BP01700
1600

43

1500
1400
1Q
43

β=0.10211700
1600

49

43

1500
1400

56
432410(Christensen, 1985)

β=0.3021T/°C
1700
1600

5815001400

643210484.584.532Ga BPSame solution was advocated by Conrad & Hager (1999, 2001), Sleep (2000), Korenaga (2003).
Problem: No self-consistent dynamical model gives such low values of β

10Heat flow and plate size

1/q2

Plate "age" Subduction

Ridge

Subduction

Ridge

Subduction

Ridge

x=0

Subduction

0.1

x=L
T=T0

z=0

u

qtop

0.0
Surface heat flow
50w

utop

0
Surface velocity

T=Tiw

0

Temperature
1u

z=d2

468(Grigné et al, 2005 ; Labrosse & Jaupart, 2007)

1012

14

16

1820

2224

26

2830

32

T=T0+ ∆Tz(Grigné et al, 2005)

Heat flux coefficient, A0.80.6I Loop model: balance between buoyancy and viscous resistance

Nu = A(w) Ra1/3Tm4/3

4/3⇒ qtop = C (L)Ra1/3m Tm

0.4

0.2

0.001

⇒ classical scaling supported by convection models with
self-consistent plate tectonics (pseudo-plastic rheology).

λ=1.0
λ=0.5
λ=0.3

(Grigné et al, 2005)

23

4

Cell width, w

56x

Alternative scenario

I Standard approach:dT= H(t) − Q(T )dtparameterised by the mantle potential temperature only.
MC⇒ Core and mantle assumed to cool at the same pace.
I Assume instead that the core is cooling and not themantle:⇒ No feedback from temperature dependence of themantle viscosity!
dTM= H(t) − Q(TM ) + QCMB
dt
dTC
MC CC= −QCMB
dt

MM CM

Alternative scenario

I Standard approach:dT= H(t) − Q(T )dtparameterised by the mantle potential temperature only.
MC⇒ Core and mantle assumed to cool at the same pace.

Past
Slightly
warmerMuch
HotterPresent

I Assume instead that the core is cooling and not themantle:⇒ No feedback from temperature dependence of themantle viscosity!
dTM= H(t) − Q(TM ) + QCMB
dtdTC
MC CC
= −QCMBdt

MM CM

A modified Urey number
Labrosse, in AGU Mon. 207, 2016

I Core heat flow is a heat source to the mantle. Assume QCMB = 13.1 TW (Buffett, yesterday, veryprecise).
⇒ Modified Urey number:
H + QCMB
13 + 13.1== 0.67
Qsurface
39
I Arguments for differential cooling of the core and mantle?
Ur ? =

Total cooling of the mantle
1000

Temperature (K)
2000
3000

4000

Pressure (GPa)1020

I Total mantle cooling in 4.5 Gyr constrained bythe phase diagram of the upper mantle:∆Tm < 200K

3040(Abe, 1997)

Core cooling constrained by the need to sustain a geodynamo

1600
61000

4800

3600
2400

60504030
CT04C10120
1000.1

0.3

1Total dissipation (TW)30100

QCMB (TW) with SiO2 cryst.

-dTCMB / dt (K/Gyr)

705

1200

20080(Hirose et al, Nature, 2017)

30
20
100

QCMB (TW) without SiO2 cryst.

with SiO2 cryst., sub-isentropic
with SiO2 cryst., super-isentropic
without SiO2 cryst., super-isentropic

Rate of SiO2 cryst. (wt% /Gyr)

1400Thermodynamics of a convective
dynamoI For a given dissipation (> 0), one
can compute the CMB heat flow.
I Compute the cooling rate of thecore
I CMB heat flow much smaller when
compositional convection is active(IC growth and/or light element
extraction)
I Results strongly depend on thevalue of the thermal conductivity(debated)Example cooling histories

High k (Gomi et al, 2013), [K] = 0
Low k (Seagle et al), [K] = 0

5750

5250
Range of silicate melting

CMB Temperature, K

5500

5000
4750
4500
4250
4000
4.5

4.0

3.5

3.02.5 2.0
Time, Gyr

1.5

1.0

0.50.0I Model with QCMB = 1.15Qs at all time
I Current CMB temperature close to the solidus=⇒ lowermost mantle likely molten in thepast.

Dense partial melt pocket at the base of the mantle

I Large VS anomalies in the lower mantle → thermal and chemical heterogeneity.
I ULVZs at the edges of dense thermo-chemical piles. Interpreted as pockets of dense partial melt.
90˚0˚−90˚−3%0˚90˚180˚0˚270˚Lateral variations of shear wave velocity (Ritsema et al 1999)−3.0

−2.5

−2.0

−1.5

−1.0

−0.50.0(Adapted from McNamara et al, 2010)

0.5

1.0

1.5

ULVZ :
ObservedAbsent

2.0

2.5

3.0

Depth= 2850 km+3%

Mantle melting(Andrault et al, Science, 2014)

Ideas on melt processes

Temperature

L

S+LS(Mg,Fe)SiO3

CaSiO3

Earth?

Eutectic

Temperature

TimeL
S+LS

MgSiO3
Earth?

FeSiO3I For chemically complex systems, the composition ofliquid and solid are different.
I In particular: Fe partitions preferentially in the liquidsilicate rather than solid.
I Density change on melting: ∆ρ = ∆ρφ + ∆ρχ .
I ∆ρφ decreases with pressureI ∆ρχ < 0Composition of lower mantle melt

(Andrault et al 2012)(KD values, not DFe)I Precise value of partition coefficient still debated but,
I Fe content of the liquid is larger than that of the solid⇒ gets richer over time.

Crystallisation of a basal magma ocean (BMO)
4.5 Ga4.49 Ga~1.5 Ga

Present

Solid mantle

BMO: liquid
Core: liquid

ULVZs(Labrosse, Hernlund, Coltice, 2007)

I ULVZ: Dense partial melt at present
I Cooling of the core evidenced by the maintenance of the geodynamo for at least 3.5 Gyrs.
I ⇒ More melt in the past!
I Fractional crystallisation ⇒ compositional variations.

Phase change boundary conditions - 1
First developed for the inner core (Deguen, Alboussière, Cardin, et al)vertical coordinate, z

0.7

0.6

Turbulent oceanLatent
Meltingheata

Reference temperature
LiquidusndFheO

q < qrefFreezing
0.5

0.4

q > qref Temperature

Convecting solid0.3

I Viscous stress in the solid mantle ⇒ topography builds with timescale
I Heat transfer in the liquid erases topography with timescale
I Competition of the two processes controlled by Φ = ττφηητη = ∆ρgdτφ = −ρs L∂ Tm∂zρl cpl ulPhase change boundary conditions - 1
First developed for the inner core (Deguen, Alboussière, Cardin, et al)vertical coordinate, z

0.7

0.6

Turbulent ocean
LatentMeltingheata

Reference temperature
LiquidusndFheO

q < qrefFreezing
0.5

0.4

q > qref Temperature

Convecting solid0.3

I Viscous stress in the solid mantle ⇒ topography builds with timescale
I Heat transfer in the liquid erases topography with timescale
I Competition of the two processes controlled by Φ = ττφηητη = ∆ρgdτφ = −ρs L∂ Tm∂zρl cpl ulPhase change boundary conditions - 1
First developed for the inner core (Deguen, Alboussière, Cardin, et al)vertical coordinate, z

0.7

0.6

Turbulent oceanLatent
Meltingheata

Reference temperature
LiquidusndFheO

q < qrefFreezing
0.5

0.4

q > qref Temperature

Convecting solid

0.3

I Viscous stress in the solid mantle ⇒ topography builds with timescale
I Heat transfer in the liquid erases topography with timescale
I Competition of the two processes controlled by Φ = ττφηητη = ∆ρgdτφ = −ρs L∂ Tm∂zρl cpl ulPhase change boundary conditions - 2I Continuity of the normal stress:

Light upper ocean

vertical coordinate, z0.5h+±−p + (ρs − ρ±l )gh + 2η

0.3

Convecting solid

0.1

I Quasistatic relationship between thetopography and the solid velocity, with τφ aphase change time scale:

−0.1
−0.3
−0.5h−w=Dense lower oceanh±τφ±

I Dimensionless boundary condition for vertical velocity:±Φ± w + 2

∂w
= 0.
∂z

∂w
− p = 0,∂zwith

Φ± =τφ± |∆ρ± |gH
τφ±=
ητη

I Φ → ∞ ⇒ classical non-penetrative boundary condition (w = 0).
I Φ → 0 ⇒ permeable boundary condition (w 6= 0).
I Also, continuity of horizontal stress and fixed (phase change) temperatures.

Thermal structure with one boundary with Φ = 0.10Temperature

10z
10z
1-2.490x2.49

Nusselt number, Nu(a)Φ- = 0.1
Φ- = 1
Φ- = 10Φ- = 100
Φ- = 1000

2.6

2.2

1.8

1.4

1.0(b)

0.6

Temperature, 〈T〉Heat transfer and mean temperature - close to onset

0.55

Φ- = 0.1
Φ- = 1
Φ- = 10Φ- = 100
Φ- = 1000

0.5
200

400

600

800

Rayleigh number, Ra

1000

1200

200

400

600

800

Rayleigh number, Ra

1000

I Good match of the fully non–linear results (DNS) and the weakly non–linear ones for small
Ra/Rac .
I Deviation at high Ra, more rapidly for heat flow (Nusselt number) than average temperature.

1200

Φ- = 0.1
Φ- = 1
Φ- = 10Φ- = 100
Φ- = 1000

Φ- = 0.1
Φ- = 1
Φ- = 10Φ- = 100
Φ- = 1000

1.0Temperature, 〈T〉Nusselt number, Nu

Heat transfer and mean temperature - high Rayleigh number101

T = 1 - 2.64 Ra-0.33

0.75

Nu = 0.37 Ra0.33
Nu = 0.21 Ra0.33
100102103104

105

Rayleigh number, Ra

1060.5102

103

104

105

Rayleigh number, Ra

I At high Ra, Nu ∼ C Ra1/3 .
I Coefficient C larger for small Φ ⇒ heat flow about twice larger for a given Ra.
I Consistent with a dynamics controlled by the only active boundary layer.

106

Toward a fully dynamical evolutionary model

I A joint effort:
I Lyon: A. Morison, R. Agrusta, R. Deguen, T. Alboussière, S. Labrosse.
I Zürich: D. Bolrão, A. Rozel, M. Ballmer, P. J. Tackley.
I Ingredients (within StagYY):
I Convection in the solid mantle,
I Tracers to treat composition (FeO content),
I Moving boundary for the net motion of the interface,
I Boundary condition to model phase change at the interface,
I Compositional and thermal evolution of the BMO and core.

Temperature

Composition: conditions for the phase changeFreezingI Liquid of uniform composition ξl
LiquidInterface

Reaction
meltingSolidξsξl

FeO content

I Freezes at Tl (ξl ) ⇒ solid of massfraction ξs = Dξl .
I Liquid gets locally enriched.
I Solid with ξ < ξs melts at T < Ts (ξs ) byreacting with the liquid.
I Needs FeO transfer in the magma ocean
from regions of freezing to regions ofmelting.

Example of evolution

I Change of dynamical regime with time.
I Gradual stabilisation of a dense layer at the bottom.
I Timing strongly depends on parameters (Rayleigh number, Buoyancy number, rheology, internal
heating, etc). To be explored further!

Preliminary simulation with crystallizing BMO

0.0

0.2

0.4
0.6
Temperature

0.8

1.0

0.10

0.12

0.14

0.16
FeO content

0.18

0.20

0.22

Nusselt at bot
Nusselt at top

Solid thickness FeO content

Nu

100
50
0
0.14
0.12
0.10
1.0
0.8
0.00

0.01

0.02

0.03

0.04
Time

0.05

0.06

0.07

Explanation for broad plumes observed in the mantle?
abSamoa

Samoa

TahitiMarquesas

–2.0

0.0

0.2

0.4
0.6
Temperature

0.8

Galapagos+2.0
δV s /V s (%)

1.0 0.10(French & Romanowicz, 2015)

0.12

0.14

0.16
FeO content

0.18

0.20

0.22

Conclusions

I The heat budget of the Earth is rather well established (with known uncertainties).
I Thermal evolution models can match the observation if exotic behavior is assumed for the mantleor (preferred) a core cooling faster than the mantle.
I Implies long-term persistence of a basal magma ocean.
I Convection in the solid greatly influenced by the possibility of melting and freezing at the bottom.
I Systematic exploration of the parameter space is needed: Rayleigh number, buoyancy ratio,rheology, internal heating, mineral phase diagram...

Conclusions

I The heat budget of the Earth is rather well established (with known uncertainties).
I Thermal evolution models can match the observation if exotic behavior is assumed for the mantleor (preferred) a core cooling faster than the mantle.
I Implies long-term persistence of a basal magma ocean.
I Convection in the solid greatly influenced by the possibility of melting and freezing at the bottom.
I Systematic exploration of the parameter space is needed: Rayleigh number, buoyancy ratio,rheology, internal heating, mineral phase diagram...