On the Evolution of Thermally Stratified Layers at the top of Earth’s Core

Stable stratification at the top of the Earth’s outer core has been suggested based upon seismic and geomagnetic observations, however, the origin of the layer is still unknown. In this paper we focus on a thermal origin for the layer and conduct a systematic study on the thermal evolution of the core. We develop a new numerical code to model the growth of thermally stable layers beneath the CMB, integrated into a thermodynamic model for the long term evolution of the core. We conduct a systematic study on plausible thermal histories using a range of core properties and, combining thickness and stratification strength constraints, investigate the limits upon the present day structure of the thermal layer. We find that whilst there are a number of scenarios for the history of the CMB heat flow, Qc, that give rise to thermal stratification, many of them are inconsistent with previously published exponential trends in Qc from mantle evolution models. Layers formed due to an exponentially decaying Qc are limited to 250-400 km thick and have maximum present-day Brunt-Väisälä periods, TBV = 8 − 24 hrs. When entrainment of the lowermost region of the layer is included in our model, the upper limit of the layer size is reduced and can fully inhibit the growth of any layer if our non-dimensional measure of entrainment, E > 0.2. The period TBV is insensitive to the evolution and so our estimates remain distinct from estimates arising from a chemical origin. Therefore, TBV should be able to discern between thermal and chemical mechanisms as improved seismic constraints are obtained.


1
The Earth's large scale magnetic field is generated within the liquid iron outer core by 2 the geodynamo process, which converts the mechanical energy of fluid motion into magnetic 3 energy. Spatial and temporal variations of the field observed at Earth's surface reflect 4 processes at the top of the core and so establishing the structure and dynamics of this region 5 is of particular importance. Much debate has focused on the presence of stable stratification period. A dynamo powered solely by thermal convection cannot be sustained if the CMB to the evolution timescale of the core (Braginsky and Roberts, 1995;Gubbins et al., 2003;Nimmo, 2015). In the convecting core lateral density fluctuations are thought to be much smaller than the radial density variation (Stevenson, 1987) and are assumed to average out. 136 This assumption is also applied to the stable region, which essentially ignores effects arising where an overbar denotes a mole fraction, A x is the atomic mass of element x, A is the 153 mean atomic mass of the mixture, and the superscript denotes liquid or solid phase. Core 154 temperature and transport properties are calculated self-consistently for each composition. 155 All parameter values are listed in Tables 1 and 2. 156 Global conservation of energy through the core requires that where k(r) is thermal conductivity, ρ(r) is the density, C p the specific heat at constant 158 pressure, ψ(r) the gravitational potential referred to zero potential at the CMB, α l x the 159 expansion coefficient for element x in the liquid phase, L = T ∆S Fe the latent heat coefficient 160 with ∆S Fe the entropy of melting for pure iron, V the volume of the whole core, and S the 161 surface of the core with outward normal n. Subscripts i, c, rs and conv denote quantities 162 evaluated at r i , r c , r s and over the convecting core respectively. Equation (2) states that the 163 heat Q c leaving the core across the CMB is balanced by the heat sources within the core: 164 the sensible heat Q s , gravitational energy Q g released as light elements left in the liquid 165 at the ICB mix the core, and latent heat Q L released on freezing at the ICB. In the Q g . 168 We have also neglected radiogenic heating due to 40 K since recent calculations suggest that 169 only minor amounts of potassium will partition into the core (Xiong et al., 2018). 170 The global energy balance can be divided into contributions from the stable layer and 171 the remainder of the core. All of the latent heat released at the ICB passes through the 172 CMB (Davies and Gubbins, 2011). We follow Lister and Buffett (1998) by assuming that 173 any gravitational energy change due to rearrangement of mass within the stable layer is 174 small enough to neglect. With these assumptions Q L and Q g are apportioned to the energy 175 balance of the well-mixed core and the global energy balance can be written where Q rs = − k(r s )∇T (r s )n · dS is the heat leaving the well-mixed region. The first 177 integral in equation (3) is evaluated using the temperature profile from the stable layer 178 while Q rs is evaluated from the parameterisation of the well-mixed region.

179
The energy budget does not contain any information about the magnetic field and therefore cannot predict if a dynamo may be sustained. Whilst a magnetic field is generated through the induction process, electric currents in the core give rise to resistive heating. This energy loss from ohmic dissipation is transferred as heat throughout the core and so does not represent any energy transfer in/out of the core. To evaluate the potential for the geodynamo to operate an entropy balance can be constructed where the ohmic dissipation does enter the equation due to being a non-reversible process. The entropy change within the core is where T i is the ICB temperature, T c is the CMB temperature, i 2 is the square of the mass 180 flux vector, and α D x is the barodiffusion coefficient for element x given by where D x and µ x are the molecular diffusivity and chemical potential for element x. can be evaluated using information from the convecting region and the CMB temperature.

196
The terms E k and E α both contain contributions from stable and well-mixed regions. The

197
ohmic dissipation E J is calculated as the remainder of equation (4) where T cen is the temperature at the center of the core, γ is the Grüneisen parameter, φ is 207 the seismic parameter and g is gravity. The total adiabatic heat flow at the CMB is which, along with Q c determines the onset of thermal stratification. The exponential in 209 equation (6)  The contributions from the well-mixed region to all terms on the right-hands side of 215 equations (2) and (4) can be expressed in terms of the cooling rate at the centre, dT cen /dt.

216
The rate of change of the inner core radius is given by (Gubbins et al., 2003) 217 where T m is the melting temperature of the core alloy. This equation defines the quantity 218 C r , which relates the core cooling rate to the inner core growth rate. The rate of change of 219 light element x in the liquid is obtained from conservation of mass and is (Gubbins et al., where M conv is the mass of the convecting core.

222
With the above definitions the energy balance for the well-mixed region can be written where V s (t) is the volume of the core below r s (t). If no stable layer exists, Q rs = Q c and 224 r s = r c . Q rs is either known based on the temperature profile at the base of the stable layer ∆S Fe (P ) = ∆S 0 + ∆S 1 P + ∆S 2 P 2 + . . . ∆S N P N .
For ρ the polynomial coefficients are all assumed constant in time with the exception of 232 ρ o 0 which is adjusted to ensure mass is conserved as the inner core radius changes. g(r) 233 and ψ(r) are found by integrating the density polynomials where g(0) = 0 and ψ(r c ) = 0.

234
The pressure P (r) is found by numerically integrating the hydrostatic pressure gradient 235 dP/dr = −ρg, subject to a specified CMB pressure of 135 GPa.

236
The melting temperature T m of the core alloy is written as where ∆T x is the depression of the melting point by impurity x and we have assumed that wherec s x is the mole fraction of element x in the solid. Relatingc s x andc l x requires knowledge 241 of how light elements partition between the liquid and solid as the inner core grows. We where µ

248
The adiabatic temperature profile is calculated at each timestep and its gradient dT a /dr 249 is used to calculate the stable layer evolution. If no stable layer is present, E J is directly 250 calculated at this stage by equation (4). Inner core nucleation occurs when T a (r = 0) = 251 T m (r = 0) and r i is thereafter defined as the radius where T a (r) = T m (r). We assume that the 252 core solidifies from the inside out and hence the radial gradient in the melting temperature 253 is necessarily steeper than the adiabat. Within the stable layer we assume that heat transport is governed by thermal conduction: where ρ s and T s are the density and temperature in the stable layer and the thermal con-256 ductivity k is allowed to vary with radius. Composition is assumed to have a uniform value, At the time-dependent stable layer interface, r s (t), the situation is more complicated. Dy- where E is the entrainment coefficient. Both upper and lower boundary conditions are 272 therefore of the Neumann type.
where ∆t is the timestep. If this condition is not satisfied by the current ∆t then a smaller   The density in (20) is derived from the temperature in the stable layer at the previous 287 iteration as where ρ and T a are respectively the PREM density and adiabatic temperature extrapolated 289 through the stable layer from the convecting region.

290
At this point the adiabatic and stable layer temperatures at the new time, T a (r, t + ∆t) 291 and T s (r, t + ∆t), will in general be discontinuous at r s (t), which will no longer be the point If fluid at any radius within the layer satisfies equation (25) or is more dense than a fluid 296 parcel from the adiabatic region would be when raised to its level (ρ s (r, t+∆t) > ρ(r, t+∆t)) 297 then the unstable fluid is assumed to mix into the bulk; the layer thickness decreases and 298 r s (t + ∆t) is moved to the point of neutral stability, ∂T a (r, t + ∆t)/∂r = ∂T s (r, t + ∆t)/∂r. If 299 the entire stable layer satisfies equation (25) then the stable region thickens and r s (t + ∆t) 300 is set as the radius where T a (r, t + ∆t) = T s (r, t + ∆t) ( Figure 2).

301
To obtain the temperature between r s (t) and r s (t + ∆t) we linearly interpolate between 302 T a (r, t + ∆t) and T s (r, t + ∆t). Consequently the temperature profile across the core at 303 the end of each iteration will be continuous, but the temperature gradient will only be 304 piecewise continuous at r s (t + ∆t). Since the individual layers generally cool by only a 305 fraction of a degree over a timestep of 1 million years the discontinuity in ∂T /∂r is orders of 306 magnitude smaller than the absolute temperature gradient. We have investigated different 307 interpolation schemes that allow continuity of T and ∂T /∂r at r s , however these higher 308 order schemes generally permit unphysical behaviour such as unstable gradients in the stable and the solution for a fixed temperature gradient at r = a is (Crank, 1979, equation 6.45 ) where α n are defined by the n th root of aα n cot(aα n ) = 1.

324
Numerical solutions were run in a spherical shell with 0.001 ≤ r ≤ a = 1 to avoid the 325 singularity at the origin, which we found to adequately represent the full-sphere geometry where I(r 1 ) and I(r 2 ) are the values of the integral I(r) at r 1 and r 2 given by and κ varies in radius such that  (Table 1). Figure   348 4 shows how the layer quickly grows and then converges to the radii at which Q rs = Q c .

349
The temperature profile in the layer is elevated above the adiabat until it merges with the 350 adiabat at r s as expected.

351
Finally, we reproduce the results of Labrosse et al. (1997). We parameterise their CMB 352 heat flow in the form where q 0 =75 mW m −2 and β = -3.5 W m −2 s −1 . The thermal conductivity of the core is  for T a , T m,Fe and k for the three cases are given in Table 1.
where A and B are the present day CMB heat flow and the linear decrease in Q c over time.

405
The best fit linear decrease in Q c over the last 0.7 Gyrs for the histories shown in Figure 6  at the present day. We find that the required adjustment to the initial temperature is very 422 small, typically less than 20 K, and so we do not expect any significant dependence of our 423 results on the initial core temperature. achieve a positive E J just prior to inner core nucleation. 466 We calculated stable layer properties for the 3 sets of core properties in Table 1 . 490 We assume that the linear time-dependence of Q c used to obtain the results in Figure 8

513
We take models that satisfy this constraint as being compatible with the published 514 models in Figure 6, limiting the selection to those models that give Q i c < 70 TW, with 515 maximum layer thicknesses for a range of ∆ρ and E values shown in Table 3. When E = 0, 516 the maximum layer thickness is ∼250-300 km for ∆ρ = 600 and 800 kg m −3 , and ∼400 km 517 for ∆ρ = 1000 kg m −3 . Increasing E quickly lowers this upper limit since thicker layers are 518 only found in regions of the parameter space that give progressively higher values for Q i c .

519
When E = 0.1, the maximum layer thickness is just <60 km for ∆ρ = 600 and 800 kg m −3 , 520 and ∼200 km for ∆ρ = 1000 kg m −3 . Only models with ∆ρ = 1000 kg m −3 produce a stable 521 layer when E = 0.2, at a maximum of just 12 km, and no models at E = 0.3 produce a 522 layer, given the constraint upon Q i c . 523 Figure 10 shows the peak Brunt-Väisälä period for all models. The maximum thermal 524 anomaly always occurs at the present-day directly below the CMB (e.g. Figure 7) and so the gravitational energy release, though its effect on transport properties has not been calculated.

541
Umemoto and Hirose (2020) suggest that hydrogen becomes relevant if the ICB temperature 542 is in the range 4800 − 5400 K, which is low compared to the range 5300 − 5900 K considered 543 here. Naively we might expect the temperature drop from 5300 K to 4800 K to produce could also address our assumption that all gravitational energy is released in the adiabatic 599 region of the core, though we do not expect this to bear strongly on our conclusions since 600 the stable layer thickness remains relatively thin.

601
The main result from this work is that thermally stable layers in Earth's core driven by     The ICB is at the radius r i , the stable layer interface at r s , and the CMB at r c . The adiabatic region is defined as 0 ≤ r ≤ r s and the stable layer at r s ≤ r ≤ r c .       Grey indicates that no stable layer forms. Black contours indicate the value for Q c at t = 500 Myr assuming that the present day rate of change in Q c were due to an exponential decay in Q c over the last 3.5 Gyrs (see text for details).