# Convective Erosion of a Primordial Stratification Atop Earth's Core

M. Bouffard, M. Landeau, A. Goument

<figure><img src="image_1.png" /></figure>

To cite this version:

M. Bouffard, M. Landeau, A. Goument. Convective Erosion of a Primordial Stratification Atop Earth's Core. Geophysical Research Letters, 2020, 47 (14), 10.1029/2020GL087109. http://hal-03202947

HAL Id: hal-03202947

https://hal.science/hal-03202947

Submitted on 11 Dec 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Geophysical Research Letters

<figure><img src="image_3.png" /></figure>

## RESEARCH LETTER

10.1029/2020GL087109

**Key Points:**

* We perform the first simulations on the erosion by thermal convection of a chemically stratified layer atop Earth's core

* Our numerical scaling suggests that a primordial chemical stratification in Earth's core can persist up to the present day

* The stratification at the top of Earth's core, suggested by seismic and geomagnetic observations, could be a vestige of core formation

**Supporting Information:**

* Supporting Information S1

* Movie S1

* Movie S2

* Movie S3

* Movie S4

* Movie S5

* Movie S6

**Correspondence to:**

M. Bouffard,
bouffard@irphe.univ-mrs.fr

**Citation:**

Bouffard, M., Landeau, M., & Goument, A. (2020). Convective erosion of a primordial stratification atop Earth's core. *Geophysical Research Letters*, 47, e2020GL087109. https://doi.org/10.1029/2020GL087109

Received 22 JAN 2020
Accepted 30 APR 2020
Accepted article online 3 MAY 2020

## Convective Erosion of a Primordial Stratification Atop Earth's Core

M. Bouffard<fnref n="1" />,² ID, M. Landeau³ ID, and A. Goument<fnref n="4" /> ID

**Abstract** Seismic and geomagnetic observations suggest the presence of a stratified layer atop Earth's core. Previous laboratory experiments showed that this layer could be primordial, produced by a collision between the primitive Earth and a giant impactor. However, paleomagnetic data require turbulent flow motions in the core for the last 3.4 Ga. Such flows can erode an existing stratification. It is therefore unclear whether a primordial stratification still exists nowadays. Here, we use numerical simulations to investigate the erosion by thermal convection of a chemical layer atop Earth's core. Our scaling law predicts that a primordial layer thicker than 1 km with a density anomaly above 0.01% can survive 4.5 Ga of convective erosion. We conclude that the observed present-day stratification could be a vestige of core formation. We also observe strong double-diffusive flows in the layer. These might reconcile the existence of a stratification with the present-day structure of the geomagnetic field.

### 1. Introduction

Seismic and geomagnetic observations have long suggested the presence of a stratified layer at the top of Earth's core (Buffett, 2014; Buffett et al., 2016; Braginsky, 1994; Lay & Young, 1990; Souriau & Poupinet, 1991; Tanaka, 2007; Whaler, 1980). Although some found no evidence for stratification (Alexandrakis & Eaton, 2010; Irving et al., 2018), recent seismic investigations indicate low P-waves velocities within a 100 to 500 km thick layer below the core-mantle boundary (CMB) (Helffrich & Kaneshima, 2010; Kaneshima, 2018; Kaneshima & Helffrich, 2013; Kaneshima & Matsuzawa, 2015; Tang et al., 2015). Previous investigations have shown that such a stratified layer strongly affects the morphology of the geomagnetic field (Christensen & Wicht, 2008; Nakagawa, 2011; Olson et al., 2017) and the long-term evolution of the core (Labrosse, 2015).

The origin of a stable layer at the top of the core remains enigmatic. The outermost core may be thermally stratified if the CMB heat flow is less than that lost by conduction (Labrosse, 2015). Alternatively, the layer could be chemically stratified and produced by accumulation of light elements emanating from the inner core boundary (Braginsky, 1994; Bouffard et al., 2019; Moffatt & Loper, 1994), barodiffusion (Braginsky, 2006; Gubbins & Davies, 2013; Fearn & Loper, 1981), or chemical interactions with the mantle (Buffett & Seagle, 2010). More recently, Landeau et al. (2016) proposed that a chemically stratified layer could be produced early in the Earth's history, as the result of the incomplete turbulent mixing between the terrestrial core and that of a giant impactor. This scenario can produce thick stratified layers, up to several hundreds of kilometers, in agreement with recent seismic investigations (Kaneshima, 2018; Tang et al., 2015). In addition, a giant impact can produce a layer with both low density and low wave speed anomalies. Brodholt and Badro (2017) showed that these constraints are met when the impactor core is enriched in one light element (e.g., C, O, or S) but depleted in another (e.g., Si), compared to the protocore. In contrast, barodiffusion or light element accumulation produce a layer enriched in all light elements, which cannot explain the observed low wave speed anomalies (Brodholt & Badro, 2017).

Primordial layers are also supported by models of Earth's accretion. These predict that the metal, progressively added to the core, varies in composition through time (Rubie et al., 2011, 2015). Injecting scalings from fluid dynamical experiments (Landeau et al., 2016) into models of core formation (Jacobson et al., 2017), we find that chemically stratified layers are almost unavoidable in the primitive core, even when we account for mixing by impacts (supporting information, Text S1).

<footnote n="1">¹Max Planck Institute for Solar System Research, Göttingen, Germany, ²Aix Marseille Univ, CNRS, Centrale Marseille, IRPHE, Marseille, France, ³Institut de physique du globe de Paris, CNRS, Université de Paris, Paris, France,</footnote>

<footnote n="4">⁴Département Sciences de la Terre et de l'Univers, Université de Nantes, Nantes, France</footnote>

<footnote>©2020. The Authors.
This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.</footnote>
However, paleomagnetic data indicate that the Earth already possessed a magnetic field 3.4 Ga ago (Biggin et al., 2011; Tarduno et al., 2010). This implies that, to sustain the geodynamo, vigorous flow motions must have taken place in the fluid core. Such turbulent flow motions can gradually erode a stratified layer. This is a well-known process in oceanic, atmospheric, and astrophysical sciences (Cristini et al., 2017; Fernando, 1991; Lareau & Horel, 2015). Whether a primordial chemically stratified layer at the top of the core could have survived up to the present day is an open question.

Most investigations on the erosion of stratification by turbulent motions focused on planar or cylindrical geometries, in the absence of rotation (Cristini et al., 2017; Meakin & Arnett, 2007; Pratt et al., 2017; Singh et al., 1995). Only a few experimental studies have included rotation (Fleury et al., 1991; Levy & Fernando, 2002; Masuda, 1983; Maxworthy, 1986) but never in the spherical geometry of the outer core. Other authors have modeled the effect of stable layers on convection in rotating spherical shells (Christensen, 2006, 2018; Manglik et al., 2010; Masada et al., 2013; Mound et al., 2019; Nakagawa, 2015; Olson et al., 2018; Stanley & Mohammadi, 2008). However, these investigations focus on a quasi-steady stratified layer, with a constant thickness. Hence, they do not capture the time-dependent erosion of a stratified layer. Furthermore these investigations either neglect or underestimate the diffusivity difference between temperature and composition.

Here, we conduct the first investigation on the erosion of a chemically stratified layer by thermal convection in a rotating spherical shell, accounting for the large diffusivity contrast between temperature and composition. We obtain a scaling law for the erosion rate as a function of a Richardson number, which measures the strength of stratification relative to that of convection. Our scaling predicts that core convection likely eroded less than 1 km of primordial stratification over the age of the Earth.

## 2. Numerical Model

When the core formed and acquired its primordial structure, it was likely entirely liquid (Nakajima & Stevenson, 2015; Tonks & Melosh, 1993). Core evolution models suggest that it remained fully liquid for most of Earth's history, for inner core likely nucleated during the last 0.3 to 1 Ga (Davies, 2015; Labrosse, 2015; Nimmo, 2007). We therefore consider stratification in a fully liquid core, prior to inner core nucleation. We model a primordial chemically stratified layer by imposing a linear initial compositional profile at the top of the core:

$$ \xi_0 = \xi(t=0) = \begin{cases} 0 & \text{if } r < r_o - h_0 \\ \frac{\Delta\xi_0}{h_0}(r - r_o + h_0) & \text{if } r \ge r_o - h_0, \end{cases} \quad (1) $$

where r is the radius, $r_o$ is the CMB radius, $\xi$ is the mass fraction of light elements, $\Delta\xi_0$ is the initial compositional difference across the layer, and $h_0$ is the initial thickness of the layer.

To model the erosion of this chemical layer by thermal convection, we use the code PARODY (Aubert et al., 2008a), which solves equations for the conservation of mass, momentum, and energy in a 3D spherical shell for a Newtonian rotating fluid in the Boussinesq approximation. We neglect magnetic effects. For numerical reasons, we retain a small inner core with radius $r_i$ so that $r_i/r_o = 0.1$. Such a small inner core has a negligible impact on the solution (Landeau & Aubert, 2011). We impose no-slip boundary conditions and neglect any chemical interaction with the inner core and mantle, so that no light elements are introduced nor removed at the boundaries. To model thermal convection, driven by the secular cooling of the Earth, we assume a constant and spatially uniform internal heating, while imposing a fixed heat flow $Q_o$ at the CMB and a null heat flow at the inner boundary (Aubert et al., 2009; Dietrich & Wicht, 2013; Landeau & Aubert, 2011). The initial overshoot in kinetic energy causes massive entrainment of the layer. To avoid this, we first run simulations without the stable layer until we reach a statistically steady state for pure thermal convection. Then, we introduce the chemical layer. In all simulations but two (see Table S1), we set the layer thickness to 20% of the shell thickness $D = (r_o - r_i)$.

We scale times and distances using the rotation frequency $1/\Omega$ and the thickness of the shell $D$, respectively. Composition is scaled by the initial compositional difference $\Delta\xi_0$ across the layer. We use $T^* = Q_o/(4\pi\rho C_p\Omega D^3)$ as the temperature scale, with $\rho$ the core density, $C_p$ the heat capacity. Five independent dimensionless numbers control the dynamics. The Ekman number, $Ek = v/(\Omega D^2)$, where v is the kinematic viscosity, measures the ratio of viscous to Coriolis forces. We vary this number within the range $10^{-5}-3 \times 10^{-4}$.
The thermal Rayleigh number, $Ra_T = \alpha g_0 Q_0 / (4\pi \rho C_p \Omega^3 D^4)$, characterizes the vigor of thermal convection below the stratified layer, with $\alpha$ the thermal expansion coefficient and $g_0$ the gravitational acceleration at the CMB. We consider thermal Rayleigh numbers that are 11 to 270 times supercritical, using the critical Rayleigh numbers in Landeau and Aubert (2011). The strength of chemical stratification is controlled by the compositional Rayleigh number, $Ra_\xi = g_0 \beta \Delta\xi_0 / (\Omega^2 D)$, where $\beta$ is the compositional expansion coefficient. We explore a range of $Ra_\xi$ corresponding to ratios $Ra_\xi/Ra_T$ comprised between $6.67 \times 10^{-3}$ and $1.67$. The Prandtl number, $Pr = v/\kappa^T$, is the ratio of momentum to thermal molecular diffusivity. We fix $Pr = 1$. The ratio of thermal to compositional molecular diffusivities is the Lewis number, $Le = \kappa^T/\kappa^\xi$. Since this number is large in the Earth's core (Bouffard et al., 2019), we explore the end member $Le = +\infty$. This implies that we totally neglect compositional molecular diffusion and focus on purely advective erosion mechanisms. To achieve infinite Lewis number numerically, we use a particle-in-cell method recently implemented into the PARODY code to describe the compositional field (Bouffard et al., 2017). We ran a total of 32 simulations summarized in Table S1.

### 3. Flow Regimes

Before we introduce the chemical layer, the whole core convects by secular cooling. Cold plumes form at the outer boundary and descend at velocity $U$. When we introduce a stable chemical layer, the cold plumes emerging from the CMB become enriched in light elements, with a compositional excess $\Delta\xi_0$. This causes a stabilizing buoyancy force, which counteracts the thermal buoyancy and slows down the plumes in a time $t_c \propto U/(\Delta\xi_0 \beta g_0)$. We measure these competing effects by the bulk Richardson number $Ri_0$, which is the ratio of the plumes advective time $t_a = D/U$ to the plumes braking time $t_c$:

$$Ri_0 = t_a/t_c = \frac{\Delta\xi_0 \beta g_0 D}{U^2}. \qquad (2)$$

As shown in Figure 1, the bulk Richardson number controls the dynamical regime in our simulations. Using the scalings from Aubert et al. (2009) (Figure S8), $Ri_0$ can be expressed as a function of the input parameters:

$$Ri_0 \approx 1.94 Ra_\xi Ra_T^{-0.84}. \qquad (3)$$

When the bulk Richardson is less than a critical value, $Ri_{0,c} \sim 90$, stratification is too weak to resist thermal convection. In this regime, the preexisting thermal plumes are little affected by composition and continue their path across the chemical layer (Figure 2 top row and Movies S1–S2). The radial kinetic energy has a comparable magnitude inside and outside the chemical layer. The entire layer is rapidly entrained and mixed with the rest of the core within a few advection times. Following Davaille et al. (2003), we refer to this case as the whole-layer convection regime. The critical Richardson number, $Ri_{0,c} \approx 90$, is well constrained by our data at $Ek = 3 \times 10^{-4}$. However, the uncertainty is larger for lower values of the Ekman number. The value of $Ri_{0,c}$ is then in the range 70–300 (Figure S2). Yet we observe no resolvable dependency of $Ri_{0,c}$ on $Ek$.

For $Ri_0 > Ri_{0,c}$, the compositional stratification is sufficiently strong to stop preexisting thermal plumes and suppress thermal convection in the layer. In this regime, the layer is stable for many advective times. However, convective motions near the base of the stratification entrain portions of the layer downward, converting kinetic energy into gravitational potential. This process progressively erodes the layer from underneath, a regime we name erosion (Figure 2 middle row and Movies S3–S6). We observe that the rate of erosion depends on latitude and on the strength of stratification (Figure S3 and Movies S3–S6). We observe no or little erosion in polar regions compared to the rest of the core. At midlatitudes, columnar vortices suck in the base of the layer, forming spikes that are stretched along the rotation axis (Figure S3a) and analogous to the structures observed in previous laboratory experiments (Fleury et al., 1991; Levy & Fernando, 2002; Masuda, 1983). In the equatorial region, we most often observe entrainment by counter rotating eddies, which sweep and stretch thin filaments from the bottom of the layer (Figures S3b–S3d). Less frequently, we observe shear instabilities (Figure S3e, see also Mory, 1991) or wave breaking (Carruthers & Hunt, 1986; Strang & Fernando, 2001). The nature and efficiency of the entrainment mechanisms vary with latitude, creating a lateral topography at the base of the layer (Figures 2, S3a, S4b, and Movies S3–S4).

In some cases of the erosion regime, we observe waves that form within the layer and eventually break into convective motions (Movies S5–S6). These flow motions are of similar strength to the thermal convection
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Regime diagram as a function of the thermal Rayleigh number $Ra_T$, normalized by its value $Ra_{T,c}$ at the onset of convection, and the bulk Richardson number $Ri_0$. Each symbol corresponds to one simulation. Blue diamonds, orange circles, and orange squares denote the whole-layer convection, erosion, and erosion with double-diffusive convection regimes, respectively. The erosion with double-diffusive convection regime is defined by strong convective radial motions that develop in more than 50% of the volume occupied by the layer. The vertical dashed line separates the whole-layer convection and erosion regimes.</figcaption></figure>

underneath, but they involve much smaller length scales (Figure 2 bottom row and Movies S5–S6). Since the layer is chemically stable and thermally unstable, it is prone to the “diffusive” regime of double-diffusive convection (Radko, 2013; Turner, 1979). We therefore interpret these small-scale convective motions as diffusive convection and call the corresponding cases *erosion with double-diffusive convection*. Double-diffusive convection extends throughout the entire layer only for $Ri_0 \le 10^3$ and for Ekman numbers $Ek \le 3 \times 10^{-5}$ (Table S1). We propose an explanation for this in Text S2.

In the early Earth’s core, the density anomaly was likely larger than 0.01% (Jacobson et al., 2017; Landeau et al., 2016). Assuming convective velocities close to $10^{-3}$ m s$^{-1}$, similar to those in the present-day core (Christensen & Aubert, 2006; Finlay & Amit, 2011), we obtain a bulk Richardson number $Ri_0$ larger than $10^{10}$. These are orders of magnitude larger than the critical value $Ri_{0,c} \approx 90$. Our results therefore suggest that the early core was in the erosion regime (Figure 1), where the layer is progressively eroded from underneath. We therefore focus solely on this regime in what follows.

## 4. Erosion Scaling

Since erosion proceeds from below, we assume that the gradient of composition inside the layer does not change over time. This approximation is common in laboratory experiments (Deardorff et al., 1980; Kato & Phillips, 1969; Levy & Fernando, 2002; Turner, 1986) and is well verified in our simulations, though to a lesser extent in the presence of double-diffusive motions (Figure S4). With this assumption, the mass balance in the convective region writes

$$ \frac{4}{3} \pi \bar{\xi}_{\text{conv}}(t) [(r_o - h(t))^3 - r_i^3] = \int_{r_o-h_0}^{r_o-h(t)} 4\pi r^2 \xi_0(r) dr, \quad (4) $$

where $\bar{\xi}_{\text{conv}}$ is the mean composition in the convective region. We measure the value of $\bar{\xi}_{\text{conv}}(t)$ in our simulations, and we then deduce the layer thickness $h(t)$, using equation (4). In the erosion regime, we observe that the layer thickness $h$ always decreases with time (Figure 3a), indicating that the stratified layer is gradually eroded.
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Various flow regimes. The top row corresponds to the whole-layer convection regime (Case 28 in Table S1) after 3.1 advection times, the middle row to the erosion regime (Case 30 in Table S1) after 9.8 advection times, and the bottom row to the erosion regime with double-diffusive convection (Case 29 in Table S1) after 6.3 advection times. Left column: azimuthal section of the compositional field. Second column: compositional field in the equatorial plane. Third column: azimuthal section of the radial velocity. Fourth column: equatorial section of the radial velocity. Color bars apply to both projections of the velocity, which is given in viscous timescale.</figcaption></figure>

As erosion proceeds and removes materials from the bottom of the layer, a jump in composition, and hence in density, builds up at the interface between the convective and stratified regions (Figure S4b). Previous investigations found that the erosion rate primarily depends on a local Richardson number, which measures the buoyancy associated with the density jump versus the strength of the flow (e.g., Fleury et al., 1991; Kato & Phillips, 1969; Kantha et al., 1977; Strang & Fernando, 2001). Following these findings, we define a local Richardson number $Ri$, based on the density jump $\Delta\rho$ at the interface with the convective region, rather than the total density jump $\Delta\rho_0$ across the entire layer. Accordingly, $Ri = \Delta\xi\beta gD/U^2$, where $\Delta\xi$ is the mean jump in composition at the base of the layer.

Because composition is larger at the top of the layer than at its base, the jump in composition $\Delta\xi$, and the Richardson number $Ri$, increase with time, as the layer gets thinner. This makes it harder for convection to entrain fluid from the base of the layer, which explains why the erosion speed $-\partial h/\partial t$ decreases with time in our simulations (Figure 3a).

To quantify this, we write the jump in composition as a function of the layer thickness, its initial value, and the initial difference in composition across the layer, $\Delta\xi(t) = \Delta\xi_0(1 - h(t)/h_0)$. This relation is valid when the gradient of composition within the layer is constant over time. Injecting this into the definition of the local Richardson number, we obtain $Ri = R_{i0}(1 - h/h_0)$.

Following previous investigations on mixing in a stratified environment (e.g., Fleury et al., 1991; Kato & Phillips, 1969; Kantha et al., 1977; Strang & Fernando, 2001), we define the dimensionless erosion rate $E = -(\partial h/\partial t)/U$, where $U$ is the typical convective velocity, measured as the time average of the root mean square velocity in the shell. For each simulation, we find a smooth best-fit curve for $h(t)$ from which we extract the erosion rate $\partial h/\partial t$ at a given time, when the convective velocity is quasi-stationary (detailed in
<figure><img src="image_6.png" /><figcaption>**Figure 3.** (a) Yellow curve: layer thickness as a function of time normalized by the advective time D/U, for a given simulation (Case 8, Ek = 3 × 10⁻⁴, Ra_T = 9 × 10⁻⁴, Ra_ε = 3 × 10⁻¹). Solid black curve: least-square best-fitting curve for h(t) using the model detailed in Text S3. From the best-fit curve, we extract the erosion rate (slope of the black curve) at a time where the convective velocity is quasi-stationary (Text S4). (b) Erosion rate E = -(dh/dt)/U as a function of the local Richardson number Ri = Δεβg₀D/U² = Ri₀(1 - h/h₀), as defined in the main text. The solid curve represents the least-square best-fitting power law E = cRi⁻ⁿ, where c = 60 ± 10 and n = 1.7 ± 0.3.</figcaption></figure>

Text S4 and illustrated in Figure 3a). From these, we deduce the local erosion rate E and the local Richardson number Ri.

Figure 3b shows that the local Richardson number Ri controls the erosion rate in our simulations. We obtain the least-square best-fitting power law:

$$E = cRi^{-n}, \qquad (5)$$

where c = 60 ± 10 and n = 1.7 ± 0.3. The positive value of n indicates that the stronger the convection, the faster it erodes the stratified layer.

The dispersion of our data in Figure 3b, measured as the relative standard deviation to the power law (5) in logarithmic scale, is 15%. Such a large dispersion is common in previous investigations (see for instance scaling laws by Deardorff et al., 1980; Kato & Phillips, 1969; Strang & Fernando, 2001). This may originate from flows inside the layer, which modify the initial compositional profile, or from bursts in the convective velocity, which induce sudden variations in the erosion rate. Both effects cause uncertainties in the measurement of the mean value of h, ∂h/∂t, and Ri. Despite this dispersion, our best-fit exponent n = 1.7 ± 0.3 is within the range reported in the literature (Fernando, 1991). Some investigations have reported a secondary dependency of the erosion rate on rotation (Fleury et al., 1991; Levy & Fernando, 2002). Yet careful analysis of our data reveals no residual dependency on the thermal Rayleigh number, the chemical Rayleigh number or the Ekman number (Text S5 and Figures S5–S7). This suggests that rotation and viscosity little affect the erosion rate in our simulations. We also find that the erosion scaling (5) does not depend on the initial layer thickness h₀ (Figures S5 and S6).

## 5. Implications for Core Stratification

The convective velocity U in the core depends on the power generated by buoyancy forces (Aubert et al., 2009) that we estimate from the energetics of the core (Lister, 2003; Labrosse, 2015). We therefore convert the scaling law of Figure 3b into a diagram, shown in Figure 4, for the mean erosion rate, as a function of the power available for convection and the initial chemical density difference across the layer, Δρ₀/ρ.

Core evolution models predict a convective power of 1 – 3 TW at present day (Davies et al., 2015; Labrosse, 2015), a range compatible with the intensity of the present-day geomagnetic field (Christensen & Aubert, 2006). Prior to inner core nucleation, the convective power takes values as low as 0.01 TW, when convection is driven by secular cooling alone (Davies et al., 2015; Labrosse, 2015), and up to about 3 TW, when we allow for exsolution of magnesium or silicon oxides (Badro et al., 2018). Such convective power values produce a
<figure><img src="image_7.png" /></figure>

**Figure 4.** Mean erosion rate, equal to the initial layer thickness divided by the total erosion time, as a function of the initial density anomaly $\Delta\rho_0/\rho$ across the layer, for various values of the power available for convection (in TW). From equation (5), we compute the mean erosion rate given by equation (13) in Text S3. We relate the convective power to the convective velocity using the power law from Aubert et al. (2009) and assuming $\Omega \sim 7.3 \times 10^{-5} \text{ s}^{-1}$, $D = 3.4 \times 10^6 \text{ m}$, and $g = 10.68 \text{ m}^2\text{s}^{-1}$. Yellow region: plausible stratification after core formation (Jacobson et al., 2017; Landeau et al., 2016). Green region: stratified layer at the top of the present-day core, as inferred from seismic observations and mineral physics (Brodholt & Badro, 2017). Blue region: stratified layer inferred from geomagnetic observations (Buffett et al., 2016; Olson et al., 2018), which suggest a layer of 140–400 km with a maximum buoyancy frequency in the range 0.84–2Ω.

magnetic field compatible with paleointensity data for the last 3.4 Ga (Aubert et al., 2009). In Figure 4, we therefore assume that a convective power in the range 0.01–3 TW is representative of that explored by the core, from its formation to the present day.

Geochemical models of core formation can estimate the primordial stratification (Jacobson et al., 2017; Rubie et al., 2011, 2015). During Earth accretion, planetary impacts brought new metal to the growing Earth. Before being added to the core, this metal equilibrates with the mantle silicates. Its abundances in light elements depend on the oxidation state and on the mantle temperature and pressure (Fischer et al., 2015; Siebert et al., 2012), which all vary during accretion (Rubie et al., 2011, 2015). Thus, the newly added metal is likely enriched or depleted in light elements compared to the protocore. This results in the formation of compositionally stratified layers up to a few hundreds kilometer thick at the top of the early core, with a density anomaly larger than 0.01% and up to 2%–3% (Jacobson et al., 2017; Landeau et al., 2016). Such primordial layers correspond to the yellow region in Figure 4, for which we predict an erosion rate of less than 100 m in 4.5 Ga. Even when accounting for the uncertainty on the power law exponent in Figure 3b, convection erodes less than about 5 km in 4.5 Ga (Figure S9). We conclude that primordial layers, inherited from core formation, can easily persist up to the present day.

Do we observe such primordial layers in the present-day core? Seismic observations suggest a stratified layer with a density anomaly of up to 1% atop Earth's core (Brodholt & Badro, 2017). The green region in Figure 4 locates such a layer, for the present-day convective power. The blue region corresponds to the layer inferred from geomagnetic observations (Buffett et al., 2016; Olson et al., 2018). The blue and green regions overlap at erosion rates always smaller than 1 km in 4.5 Ga (Figure 4). This is negligible compared to the 100 km or more of the layer, implying that the present-day stratification could easily have survived to 4.5 Ga of core convection. The stratified layer at the top of the core could therefore be a vestige of core formation.

A recent seismic investigation by Irving et al. (2018) suggests that a stratified layer atop Earth's core is not required to explain observations of normal modes and body waves travel times. Therefore, a presently unstratified outermost core is still possible. However, our results indicate that this scenario is not easy to
reconcile with a stratified primordial core, which is predicted by accretion models (Jacobson et al., 2017; Landeau et al., 2016). Figure 4 shows that convection cannot efficiently mix the stratification inherited from core formation. Hence, if the present-day outermost core is not stratified, an additional and presently unknown mixing mechanism is needed.

Dynamo simulations recently brought additional constraints on core stratification. Gastine et al. (2019) and Yan and Stanley (2018) showed that a thick stratified layer at the top of the core is hardly compatible with the structure of the present-day geomagnetic field, unless a more complex dynamics takes place in the layer. Manglik et al. (2010) showed that fingering convection within a thermally stratified layer enhances the poloidal magnetic field and produces a stronger magnetic field outside the core. This suggests that a double-diffusive layer above a dynamo region may affect the magnetic field in a way that is more complex than a pure skin effect. In some simulations, we observe strong double-diffusive radial flows inside the chemically stratified layer, similar to those observed in Earth's oceans and proposed for gas giants and stars (Leconte & Chabrier, 2012; Merryfield, 1995; Rosenblum et al., 2011; Woodgate et al., 2007). Such dynamics had never been reported in a rotating shell. Despite the absence of a magnetic field in our simulations, we suggest that these flows might reconcile the presence of a thick stratified layer below the CMB with the morphology of the geomagnetic field (Gubbins, 2007; Gastine et al., 2020; Olson et al., 2017). In the future, it is therefore crucial to investigate double-diffusive flows in the core and their role in models of the geodynamo.

Seismic observations also indicate a compositional stratification, the F-layer, at the bottom of the present-day core, with density anomalies up to 6% (Gubbins et al., 2008; Souriau & Poupinet, 1991; Wong et al., 2018; Zou et al., 2008). The F-layer could be produced on long timescales by inner core translation (Alboussière et al., 2010) or by multiphase flows at the inner core boundary (Deguen et al., 2007; Fearn et al., 1981; Wong et al., 2018). Our scaling law (5) has been established for a chemical layer at the top of the core. A chemical layer at the bottom of the core may have different dynamics due to the continuous release of light elements caused by inner core crystallization. Therefore, erosion of a chemical layer at the bottom of the core deserves a full investigation. Yet many of the entrainment mechanisms listed in section 3, such as the suction exerted by cyclonic vortices or the entrainment of thin filaments, are also expected for a layer at the base of the outer core. If our scaling law applies to such a layer, the results in Figure 4 suggest that the F-layer would likely survive convective erosion on long timescales.

Our simulations neglect the effect of compositional molecular diffusion. Using previous scalings on diffusion-limited erosion (Noh & Fernando, 1993), we find that accounting for molecular diffusion in the core does not affect the above conclusions (Text S6). The ratio of inertial to viscous forces, called the Reynolds number, varies in the range 31–739 in our simulations (Table S1). In this range, we observe no dependency of the erosion rate on the Reynolds number (Text S5 and Figures S5 and S6). Yet our values of the Reynolds number are lower than in the core. More extreme simulations are required to investigate the nature and efficiency of the erosion process in turbulent regimes. Furthermore, our simulations perform at unrealistically high Ekman numbers. Previous investigations find that increasing the rotation rate at fixed Richardson number decreases the erosion rate (Fleury et al., 1991; Levy & Fernando, 2002). Thus, decreasing the Ekman number toward Earth-like values might lower the erosion rate. With even lower erosion rates, our main conclusions from Figure 4 remain unchanged. Earth's magnetic field is another feature absent in our simulations. It may affect flow motions in the core and therefore the erosion of a primordial stratification. This effect should be investigated in the future.

## Acknowledgments

We are grateful to two anonymous reviewers for their valuable comments that significantly improved the manuscript. We also thank P. Cardin, U.R. Christensen, J. Wicht, W. Dietrich, S. Sanchez, A. Fournier and T. Gastine for fruitful discussions and suggestions on that topic. Numerical simulations were performed on the Hydra and Cobra clusters of the MPCDF computing facility (Garching, Germany). This project has received funding from the European Union's Horizon 2020 research and innovation program under the Marie-Sklodowska-Curie Grant Agreement 703767 and through Grant No. 681835-FLUDYCO-ERC-2015-CoG. Data displayed on Figure 3b can be accessed at https://zenodo.org/badge/ latestdoi/235595334. For the PARODY-PIC code and simulations data, please ask the corresponding author.

## References

* Alboussière, T., Deguen, R., & Melzani, M. (2010). Melting-induced stratification above the Earth's inner core due to convective translation. *Nature*, **466**, 744–747.

* Alexandrakis, C., & Eaton, D. W. (2010). Precise seismic-wave velocity atop Earth's core: No evidence for outer-core stratification. *Physics of the Earth and Planetary Interiors*, **180**(1-2), 59–65.

* Aubert, J., Amit, H., Hulot, G., & Olson, P. (2008a). Thermochemical flows couple the Earth's inner core growth to mantle heterogeneity. *Nature*, **758–761**.

* Aubert, J., Labrosse, S., & Poitou, C. (2009). Modelling the palaeo-evolution of the geodynamo. *Geophysical Journal International*, **179**(3), 1414–1428.

* Badro, J., Aubert, J., Hirose, K., Nomura, R., Blanchard, I., Borensztajn, S., & Siebert, J. (2018). Magnesium partitioning between Earth's mantle and core and its potential to drive an early exsolution geodynamo. *Geophysical Research Letters*, **45**, 13–240. https://doi.org/10.1029/2018GL080405
Biggin, A. J., de Wit, M. J., Langereis, C. G., Zegers, T. E., Voute, S., Dekkers, M. J., & Drost, K. (2011). Palaeomagnetism of Archaean rocks of the Onverwacht Group, Barberton Greenstone Belt (southern Africa): Evidence for a stable and potentially reversing geomagnetic field at ca. 3.5 Ga. *Earth and Planetary Science Letters*, **302**(3-4), 314–328.

Bouffard, M., Choblet, G., Labrosse, S., & Wicht, J. (2019). Chemical convection and stratification in the Earth's outer core. *Frontiers in Earth Science*, **7**, 99. https://doi.org/10.3389/feart.2019.00099

Bouffard, M., Labrosse, S., Choblet, G., Fournier, A., Aubert, J., & Tackley, P. J. (2017). A particle-in-cell method for studying double-diffusive convection in the liquid layers of planetary interiors. *Journal of Computational Physics*, **346**, 552–571.

Braginsky, S. (1994). MAC-oscillations of the hidden ocean of the core. *Journal of Geomagnetism and Geoelectricity*, **45**(11-12), 1517–1538.

Braginsky, S. I. (2006). Formation of the stratified ocean of the core. *Earth and Planetary Science Letters*, **243**(3), 650–656.

Brodholt, J., & Badro, J. (2017). Composition of the low seismic velocity E' layer at the top of the Earth's core. *Geophysical Research Letters*, **44**, 8303–8310. https://doi.org/10.1002/2017GL074261

Buffett, B. (2014). Geomagnetic fluctuations reveal stable stratification at the top of the Earth's core. *Nature*, **507**(7493), 484. https://doi.org/10.1038/nature13122

Buffett, B., Knezek, N., & Holme, R. (2016). Evidence for MAC waves at the top of Earth's core and implications for variations in length of day. *Geophysical Journal International*, **204**(3), 1789–1800.

Buffett, B. A., & Seagle, C. T. (2010). Stratification of the top of the core due to chemical interactions with the mantle. *Journal of Geophysical Research: Solid Earth*, **115**, B04407. https://doi.org/10.1029/2009JB006751

Carruthers, D., & Hunt, J. (1986). Velocity fluctuations near an interface between a turbulent region and a stably stratified layer. *Journal of Fluid Mechanics*, **165**, 475–501.

Christensen, U. R. (2006). A deep dynamo generating Mercury's magnetic field. *Nature*, **444**(7122), 1056–1058.

Christensen, U. R. (2018). Geodynamo models with a stable layer and heterogeneous heat flow at the top of the core. *Geophysical Journal International*, **215**(2), 1338–1351.

Christensen, U., & Aubert, J. (2006). Scaling properties of convection-driven dynamos in rotating spherical shells and application to planetary magnetic fields. *Geophysical Journal International*, **117**, 97–114.

Christensen, U. R., & Wicht, J. (2008). Models of magnetic field generation in partly stable planetary cores: Applications to Mercury and Saturn. *Icarus*, **196**(1), 16–34.

Cristini, A., Meakin, C., Hirschi, R., Arnett, D., Georgy, C., Viallet, M., & Walkington, I. (2017). 3D hydrodynamic simulations of carbon burning in massive stars. *Monthly Notices on the Royal Astronomical Society*, **471**(1), 279–300.

Davaille, A., Le Bars, M., & Carbonne, C. (2003). Thermal convection in a heterogenous mantle. *Comptes Rendus Geoscience*, **335**(1), 141–156.

Davies, C. (2015). Cooling history of Earth's core with high thermal conductivity. *Physical Earth Planetary International*, **247**, 65–79.

Davies, C., Pozzo, M., Gubbins, D., & Alfè, D. (2015). Constraints from material properties on the dynamics and evolution of Earth's core. *Nature Geoscience*, **8**(9), 678–685.

Deardorff, J., Willis, G., & Stockton, B. (1980). Laboratory studies of the entrainment zone of a convectively mixed layer. *Journal of Fluid Mechanics*, **100**, 41–64.

Deguen, R., Alboussiere, T., & Brito, D. (2007). On the existence and structure of a mush at the inner core boundary of the Earth. *Physics of the Earth and Planetary Interiors*, **164**(1), 36–49.

Dietrich, W., & Wicht, J. (2013). A hemispherical dynamo model: Implications for the Martian crustal magnetization. *Physics of the Earth and Planetary Interiors*, **217**, 10–21.

Fearn, D. R., & Loper, D. E. (1981). Compositional convection and stratification of Earth's core. *Nature*, **289**(5796), 393–394.

Fearn, D., Loper, P., & Roberts, D. (1981). Structure of the Earth's inner core. *Nature*, **292**(5820), 232–233.

Fernando, HarindraJ. S. (1991). Turbulent mixing in stratified fluids. *Annual Review of Fluid Mechanics*, **23**(1), 455–493.

Finlay, C. C., & Amit, H. (2011). On flow magnitude and field-flow alignment at Earth's core surface. *Geophysical Journal International*, **186**(1), 175–192.

Fischer, R., Nakajima, Y., Campbell, A., Frost, D., Harries, D., Langenhorst, F., et al. (2015). High pressure metal-silicate partitioning of Ni, Co, V, Cr, Si, and O. *Geochemistry, Cosmochimica Acta*, **167**, 177–194.

Fleury, M., Mory, M., Hopfinger, E., & Auchere, D. (1991). Effects of rotation on turbulent mixing across a density interface. *Journal of Fluid Mechanics*, **223**, 165–191.

Gastine, T., Aubert, J., & Fournier, A. (2019). Dynamo-based limit to the extent of a stable layer atop Earth's core. *Geophysical Journal International*. https://doi.org/10.1093/gji/ggaa250

Gubbins, D. (2007). Geomagnetic constraints on stratification at the top of Earth's core. *Earth, Planets and Space*, **59**(7), 661–664.

Gubbins, D., & Davies, C. (2013). The stratified layer at the core–mantle boundary caused by barodiffusion of oxygen, sulphur and silicon. *Physics of the Earth and Planetary Interiors*, **215**, 21–28.

Gubbins, D., Masters, G., & Nimmo, F. (2008). A thermochemical boundary layer at the base of Earth's outer core and independent estimate of core heat flux. *Geophysical Journal International*, **174**(3), 1007–1018.

Helffrich, G., & Kaneshima, S. (2010). Outer-core compositional stratification from observed core wave speed profiles. *Nature*, **468**, 807–810.

Irving, J. C., Cottaar, S., & Lekić, V. (2018). Seismically determined elastic parameters for Earth's outer core. *Science Advances*, **4**(6), eaa2538.

Jacobson, S. A., Rubie, D. C., Hernlund, J., Morbidelli, A., & Nakajima, M. (2017). Formation, stratification, and mixing of the cores of Earth and Venus. *Earth and Planetary Science Letters*, **474**, 375–386.

Kaneshima, S. (2018). Array analyses of SmKS waves and the stratification of Earth's outermost core. *Physics of the Earth and Planetary Interiors*, **276**, 234–246.

Kaneshima, S., & Helffrich, G. (2013). Vp structure of the outermost core derived from analysing large-scale array data of SmKS waves. *Geophysical Journal International*, **193**(3), 1537–1555.

Kaneshima, S., & Matsuzawa, T. (2015). Stratification of Earth's outermost core inferred from SmKS array data. *Progress in Earth and Planetary Science*, **2**(1), 15.

Kantha, L., Phillips, O., & Azad, R. (1977). On turbulent entrainment at a stable density interface. *Journal of Fluid Mechanics*, **79**, 753–768.

Kato, H., & Phillips, O. (1969). On the penetration of a turbulent layer into stratified fluid. *Journal of Fluid Mechanics*, **37**, 643–655.

Labrosse, S. (2015). Thermal evolution of the core with a high thermal conductivity. *Physics of the Earth and Planetary Interiors*, **247**, 36–55.

Landeau, M., & Aubert, J. (2011). Equatorially asymmetric convection inducing a hemispherical magnetic field in rotating spheres and implications for the past Martian dynamo. *Physical Earth Planetary International*, **185**(3-4), 61–73.

Landeau, M., Olson, P., Deguen, R., & Hirsh, B. H. (2016). Core merging and stratification following giant impact. *Nature Geoscience*, **9**(10), 786–789.
Lareau, N. P., & Horel, J. D. (2015). Turbulent erosion of persistent cold-air pools: Numerical simulations. *Journal of the Atmospheric Sciences*, 72(4), 1409–1427.

Lay, T., & Young, C. J. (1990). The stably-stratified outermost core revisited. *Geophysical Research Letters*, 17(11), 2001–2004.

Leconte, J., & Chabrier, G. (2012). A new vision of giant planet interiors: Impact of double diffusive convection. *Astronomy & Astrophysics*, 540, A20.

Levy, M., & Fernando, H. (2002). Turbulent thermal convection in a rotating stratified fluid. *Journal of Fluid Mechanics*, 467, 19–40.

Lister, J. R. (2003). Expressions for the dissipation driven by convection in the Earth's core. *Physical Earth Planetary International*, 140(1-3), 145–158.

Manglik, A., Wicht, J., & Christensen, U. R. (2010). A dynamo model with double diffusive convection for Mercury's core. *Earth and Planetary Science Letters*, 289(3-4), 619–628.

Masada, Y., Yamada, K., & Kageyama, A. (2013). Effects of penetrative convection on solar dynamo. *The Astrophysical Journal*, 778(1), 11.

Masuda, A. (1983). Entrainment by grid-generated turbulence in rotating stratified fluids. *Kyushu University Research Institute Applied Mechanics Reports*, 31, 29–45.

Maxworthy, T. (1986). On turbulent mixing across a density interface in the presence of rotation. *Journal of Physical Oceanography*, 16(6), 1136–1137.

Meakin, C. A., & Arnett, D. (2007). Turbulent convection in stellar interiors. I. Hydrodynamic simulation. *The Astrophysical Journal*, 667(1), 448.

Merryfield, W. J. (1995). Hydrodynamics of semiconvection. *The Astrophysical Journal*, 444, 318–337.

Moffatt, H., & Loper, D. (1994). The magnetostrophic rise of a buoyant parcel in the Earth's core. *Geophysical Journal International*, 117(2), 394–402.

Mory, M. (1991). A model of turbulent mixing across a density interface including the effect of rotation. *Journal of Fluid Mechanics*, 223, 193–207.

Mound, J., Davies, C., Rost, S., & Aurnou, J. (2019). Regional stratification at the top of Earth's core due to core-mantle boundary heat flux variations. *Nature Geoscience*, 12(7), 575–580.

Nakagawa, T. (2011). Effect of a stably stratified layer near the outer boundary in numerical simulations of a magnetohydrodynamic dynamo in a rotating spherical shell and its implications for Earth's core. *Physics of the Earth and Planetary Interiors*, 187(3), 342–352.

Nakagawa, T. (2015). An implication for the origin of stratification below the core–mantle boundary region in numerical dynamo simulations in a rotating spherical shell. *Physics of the Earth and Planetary Interiors*, 247, 94–104.

Nakajima, M., & Stevenson, D. J. (2015). Melting and mixing states of the Earth's mantle after the Moon-forming impact. *Earth and Planetary Science Letters*, 427, 286–295.

Nimmo, F. (2007). *Treatise on geophysics*, vol. 8: Elsevier Amsterdam.

Noh, Y., & Fernando, H. J. (1993). The role of molecular diffusion in the deepening of the mixed layer. *Dynamics of Atmospheres and Oceans*, 17(2-3), 187–215.

Olson, P., Landeau, M., & Reynolds, E. (2017). Dynamo tests for stratification below the core-mantle boundary. *Physics of the Earth and Planetary Interiors*, 271, 1–18.

Olson, P., Landeau, M., & Reynolds, E. (2018). Outer core stratification from the high latitude structure of the geomagnetic field. *Frontiers Earth Science*, 6, 140.

Pratt, J., Baraffe, I., Goffrey, T., Constantino, T., Viallet, M., Popov, M., et al. (2017). Extreme value statistics for two-dimensional convective penetration in a pre-main sequence star. *Astronomy & Astrophysics*, 604, A125.

Radko, T. (2013). *Double-diffusive convection*. Cambridge: Cambridge University Press.

Rosenblum, E., Garaud, P., Traxler, A., & Stellmach, S. (2011). Turbulent mixing and layer formation in double-diffusive convection: Three-dimensional numerical simulations and theory. *The Astrophysical Journal*, 731(1), 66.

Rubie, D. C., Frost, D. J., Mann, U., Asahara, Y., Nimmo, F., Tsuno, K., et al. (2011). Heterogeneous accretion, composition and core–mantle differentiation of the Earth. *Earth and Planetary Science Letters*, 301(1-2), 31–42.

Rubie, D. C., Jacobson, S. A., Morbidelli, A., O'Brien, D. P., Young, E. D. D., de Vries, J., et al. (2015). Accretion and differentiation of the terrestrial planets with implications for the compositions of early-formed solar system bodies and accretion of water. *Icarus*, 248, 89–108.

Siebert, J., Badro, J., Antonangeli, D., & Ryerson, F. (2012). Metal–silicate partitioning of Ni and Co in a deep magma ocean. *Earth and Planetary Science Letters*, 321-322, 189–197.

Singh, H. P., Roxburgh, I. W., & Chan, K. L. (1995). Three-dimensional simulation of penetrative convection: Penetration below a convection zone. *Astronomy & Astrophysics*, 295, 703.

Souriau, A., & Poupinet, G. (1991). The velocity profile at the base of the liquid core from a PKP(BC+Cdiff) data: An argument in favor of radial inhomogeneity. *Geophysical Research Letters*, 18, 2023–2026.

Stanley, S., & Mohammadi, A. (2008). Effects of an outer thin stably stratified layer on planetary dynamos. *Physics of the Earth and Planetary Interiors*, 168, 179–190.

Strang, E., & Fernando, H. (2001). Entrainment and mixing in stratified shear flows. *Journal of Fluid Mechanics*, 428, 349–386.

Tanaka, S. (2007). Possibility of a low P-wave velocity layer in the outermost core from global SmKS waveforms. *Earth and Planetary Science Letters*, 259(3), 486–499.

Tang, V., Zhao, L., & Hung, S.-H. (2015). Seismological evidence for a non-monotonic velocity gradient in the topmost outer core. *Scientific reports*, 5, 8613.

Tarduno, J., Cottrell, R., Davis, W., Nimmo, F., & Bono, R. (2010). A Hadean to Paleoarchean geodynamo recorded by single zircon crystals. *Science*, 349(6247), 521–524.

Tonks, W. B., & Melosh, H. J. (1993). Magma ocean formation due to giant impacts. *Journal of Geophysical Research*, 98(E3), 5319–5333.

Turner, J. S. (1979). *Buoyancy effects in fluids*. Cambridge: Cambridge University Press.

Turner, J. S. (1986). Turbulent entrainment: The development of the entrainment assumption, and its application to geophysical flows. *Journal Fluids Mechanics*, 173, 431–471.

Whaler, K. (1980). Does the whole of the Earth's core convect? *Nature*, 287(5782), 528.

Wong, J., Davies, C., & Jones, C. (2018). A Boussinesq slurry model of the F-layer at the base of Earth's outer core. *Geophysical Journal International*, 214(3), 2236–2249.

Woodgate, R. A., Aagaard, K., Swift, J. H., Smethie Jr, W. M., & Falkner, K. K. (2007). Atlantic water circulation over the Mendeleev Ridge and Chukchi Borderland from thermohaline intrusions and water mass properties. *Journal of Geophysical Research*, 112, C02005. https://doi.org/10.1029/2005JC003416

Yan, C., & Stanley, S. (2018). Sensitivity of the geomagnetic octupole to a stably stratified layer in the Earth's core. *Geophysical Research Letters*, 45, 11,005–11,011. https://doi.org/10.1029/2018GL078975
Zou, Z., Koper, K., & Cormier, V. (2008). The structure of the base of the outer core inferred from seismic waves diffracted around the inner core. *Journal of Geophysical Research, 113*, B05314. https://doi.org/10.1029/2007JB005316