Convective Erosion of a Primordial Stratification Atop Earth's Core

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.

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;, 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(Rubie et al., , 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).
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(Christensen, , 2018Manglik 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.

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: where r is the radius, r o is the CMB radius, is the mass fraction of light elements, Δ 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∕ and the thickness of the shell D, respectively. Composition is scaled by the initial compositional difference Δ 0 across the layer. We use T * = Q o ∕(4 C p D 3 ) as the temperature scale, with the core density, C p the heat capacity. Five independent dimensionless numbers control the dynamics. The Ekman number, Ek = ∕( D 2 ), where is the kinematic viscosity, measures the ratio of viscous to Coriolis forces. We vary this number within the range 10 −5 −3×10 −4 .
The thermal Rayleigh number, Ra T = g o Q o ∕(4 C p 3 D 4 ), characterizes the vigor of thermal convection below the stratified layer, with the thermal expansion coefficient and g o 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, We explore a range of Ra corresponding to ratios Ra ∕Ra T comprised between 6.67 × 10 −3 and 1.67. The Prandtl number, Pr = ∕ 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 = T ∕ . Since this number is large in the Earth's core (Bouffard et al., 2019), we explore the end member Le = +∞. 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.

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 Δ 0 . This causes a stabilizing buoyancy force, which counteracts the thermal buoyancy and slows down the plumes in a time t c ∝ U∕(Δ 0 g o ). 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 : 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: When the bulk Richardson is less than a critical value, Ri 0,c ∼ 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 ≃ 90, is well constrained by our data at Ek = 3 × 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 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.
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 ≲ 10 3 and for Ekman numbers Ek ≤ 3 × 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 ≃ 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.

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 4 3̄c onv (t) wherēc onv is the mean composition in the convective region. We measure the value of̄c onv (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.  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.
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 Δ at the interface with the convective region, rather than the total density jump Δ 0 across the entire layer. Accordingly, Ri = Δ gD∕U 2 , where Δ 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 Δ , 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 − h∕ 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, Δ (t) = Δ 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 = Ri 0 (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 = −( h∕ 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 h∕ t at a given time, when the convective velocity is quasi-stationary (detailed in 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: 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 disperion 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 0 ( Figures S5 and S6).

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, Δ 0 ∕ .
Core evolution models predict a convective power of 1 − 3 TW at present day 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 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 4. Mean erosion rate, equal to the initial layer thickness divided by the total erosion time, as a function of the initial density anomaly Δ 0 ∕ 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 ∼ 7.3 × 10 −5 s −1 , D = 3.4 × 10 6 m, and g = 10.68 m 2 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., 2011Rubie et al., , 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(Rubie et al., , 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 10.1029/2020GL087109 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;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. 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.