Scaling Laws for Regional Stratification at the Top of Earth's Core

Seismic and geomagnetic observations have been used to argue both for and against a global stratified layer at the top of Earth's outer core. Recently, we used numerical models of turbulent thermal convection to show that imposed lateral variations in core‐mantle boundary (CMB) heat flow can give rise to regional lenses of stratified fluid at the top of the core while the bulk of the core remains actively convecting. Here, we develop theoretical scaling laws to extrapolate the properties of regional stratified lenses measured in simulations to the conditions of Earth's core. We estimate that regional stratified lenses in Earth's core have thicknesses of up to a few hundred kilometres and Brunt‐Väisälä frequencies of hours, consistent with independent observational constraints. The location, thickness, and strength of the stratified regions would change over geological time scales in response to the slowly evolving CMB heat flux heterogeneity imposed by mantle convection.

Convection in the core is controlled by heat flow across the CMB. Compared to the dynamics of the relatively low-viscosity core, solid-state convection in the overlying mantle is associated with long time scales and large temperature variations, such that the core is subjected to large lateral variations in CMB heat flux (Nakagawa & Tackley, 2008;Olson et al., 2015;Stackhouse et al., 2015;Zhang & Zhong, 2011). This CMB heat flux heterogeneity would interact with, and potentially disrupt, any inherent core stratification (e.g., Christensen, 2018;Cox et al., 2019;Gibbons & Gubbins, 2000;Gubbins et al., 2015;Lister, 2004;Olson et al., 2017) and can have a significant influence on the pattern of core convection and hence the geomagnetic field (e.g., Davies et al., 2008;Glatzmaier et al., 1999;Gubbins & Gibbons, 2004;Gubbins et al., 2007;Olson et al., 2015;Olson & Christensen, 2002).
An alternative view of core stratification has recently been suggested from numerical modeling in which stratification is caused, rather than opposed, by lateral CMB heat flow variations; furthermore, the resultant stratification is found to be confined into regional lenses, rather than a global layer (Mound et al., 2019). In some cases, 1D averaging over strong and laterally extensive regional inversion lenses can produce an apparent global stratification despite there being radial motion throughout the core including its outermost regions. Regional inversion lenses are ubiquitous in our simulations; however, estimation of their expected thickness L and Brunt-Väisälä frequency N in the Earth requires extrapolation from the computationally accessible parameter regime to that characteristic of Earth's core.
Three nondimensional parameters control the dynamic behavior in our numerical model of rotating nonmagnetic convection in a spherical shell . The Prandtl number Pr = ν/κ is the ratio of ©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. the fluid's kinematic viscosity ν and its thermal diffusivity κ. The strength of convective driving is described by the Rayleigh number f Ra¼αg o β=2Ωκ, where α is the thermal expansivity of the fluid, g o is the gravitational acceleration on the outer boundary (r = r o ), Ω is the planetary rotation rate, and β¼r 2 o q ave =k, where q ave is the average heat flux across the outer boundary and k = κρC P is the thermal conductivity of the fluid, with ρ and C P the fluid density and specific heat, respectively. The importance of the fluid viscosity relative to rotation is described by the Ekman number E = ν/2Ωh 2 , where h = r o − r i is the shell thickness. We describe the amplitude of heat flux heterogeneity at the CMB using q ⋆ = (q max − q min )/q ave , where q max and q min are the maximum and minimum heat flux, respectively (with outward heat flux defined to be positive). Scaling analysis requires simulations spanning a sufficient region of parameter space; focusing on nonmagnetic convection to reduce computational cost enabled us to explore the influence of E, f Ra, and q ⋆ (in all cases, we hold Pr = 1). Lorentz forces will also play a role in the Earth's core, a point that we will return to in the discussion.
Our previous work considered two patterns of CMB heat flux heterogeneity, one derived from seismic tomography (Masters et al., 1996) and an east-west hemispheric pattern, described by a spherical harmonic pattern of degree and order 1 that could be relevant when considering the mantle flow associated with supercontinent assembly (Zhong et al., 2007). Values of q ⋆ = 2.3 or 5.0 were considered for the amplitude of the CMB heat flux heterogeneity. We produced a suite of nonmagnetic rotating convection simulations covering E = {10 −4 ,10 −5 ,10 −6 }, f Ra up to several hundred times the critical value for the onset of convection.
Although our simulations approach the limit of what is computationally feasible, they remain far from the parameter regime for the Earth's core. In particular, estimates of the relevant parameters suggest that f Ra may be far larger and E far smaller in the Earth than in our simulations (Mound et al., 2019). The value of q ⋆ is uncertain in the Earth as it requires knowledge of both the temperature structure and thermal conductivity of the lowermost mantle and the total superadiabatic CMB heat flow; nevertheless, its value in the Earth may be an order of magnitude larger than in our simulations (Mound et al., 2019). In this work, we first establish the theory relating L and N to the underlying physical parameters of the convecting system. We then show that our simulations match this theoretical expectation, enabling us to extrapolate to parameter values plausibly representative of the Earth's core.

Scaling Theory
The dynamics of convection falls into qualitatively different regimes depending on what combination of forces are important and which play a subdominant or inconsequential role. Scaling laws relating emergent behaviors to the imposed control parameters differ between dynamic regimes; so, care must be taken when extrapolating simulation results to planetary conditions (e.g., Aubert et al., 2017;Gastine et al., 2016;Jones, 2015;King et al., 2013). We focus on the regime of turbulent rotating convection where Inertial, Archimedean buoyancy, and Coriolis forces all play an important role in the dynamics (the IAC balance), which holds in 34 of our previously presented simulations (Long et al., 2020;Mound & Davies, 2017).
In a fluid where density decreases with increasing radius, a fluid parcel displaced radially with be returned to its original depth by buoyancy forces with a characteristic Brunt-Väisälä frequency given by where ρ 0 is a reference density. For radial density variations arising from purely thermal effects For our Boussinesq models with fixed-flux thermal boundary conditions, the strength of thermal stratification is approximately set by the temperature gradient associated with the value of q min imposed at the CMB. We note that along some radial profiles, the maximum temperature gradient occurs some distance below the outer boundary; nevertheless, we will use ∂T/∂r ≈ −q min /k to estimate the maximum value of N expected in our simulations. For a simple pattern of CMB heat flux variation and our definition of q ⋆ , we expect and hence Using 2Ω as our frequency scaling leads to For the Earth, it will be the average superadiabatic heat flux q þ ave that controls the vigor of convection; so, it is useful to recast our expression for N as where q ptp = q max − q min is the peak-to-peak variation in CMB heat flux. For comparison to our simulations, it is useful to combine the physical parameters into the relevant control parameters giving an expected scaling of Only for sufficiently strong heat flux heterogeneity will there be regions of the CMB beneath which convection is entirely suppressed. In the Earth, this requires regions of sufficiently hot lowermost mantle such that the imposed temperature gradient is subadiabatic. This requirement enters the equations above via the need for q ⋆ > 2 or, equivalently, q ptp > 2q þ ave in order to ensure N is a positive real number.
The vigor of convection within the core is ultimately controlled by the imposed heat flux boundary conditions. The average CMB heat flux in our simulations is sufficiently supercritical with respect to the onset of convection that the fluid interior would be effectively well mixed throughout the bulk of the shell if q ⋆ = 0. As a result, the radial transport of heat at depth will be primarily associated with the advective term of the temperature equation (u · ∇T). Sufficiently close to the boundary, the imposed pattern of heat flux heterogeneity will directly influence the dynamics. Under regions where the heat extracted from the core is higher than average, stronger convective motion will be driven at the top of the core, with correspondingly larger radial heat advection. Beneath regions with progressively lower CMB heat flux, the vigor of convection at the top of the core will be correspondingly reduced and eventually suppressed entirely.
In our simulations, q ⋆ is sufficiently large that a thermally stratified regional inversion lens must exist immediately beneath some of the CMB. Within a lens, heat transport is radially inward and dominated by the diffusive term of the temperature equation (κ∇ 2 T). The more anomalous the amplitude of q min in our simulations, the greater the depth over which convection is suppressed by the conductive temperature profile arising due to the CMB heat flux minimum. We define the thickness of the regional inversion lenses by finding the depth of neutral stability (time-averaged dT/dr = 0 in our thermally driven Boussinesq model) and expect that this point is determined via competition between the temperature profiles associated with the advective heat transport in the interior and the conductive heat transport imposed at the CMB. This suggests a scaling of the form where U is the characteristic velocity of convection, ℓ is the characteristic length scale of convection, T′ is the characteristic convective temperature fluctuation, ΔT′ is the total temperature anomaly across the thickness of the lens, and L is the characteristic lens thickness. Multiplying each side by 1/h 2 and rearranging gives 10.1029/2020GL087715

Geophysical Research Letters
The average advective heat flux in the interior of our models will be determined by the average imposed heat flux at the CMB; so, we expect ρC P UT′ ∼ q ave . As noted in the discussion of the scaling for the Brunt-Väisälä frequency, the strongest inverted temperature gradient (and hence conductive heat transport) can be associated with the minimum imposed CMB heat flux such that kΔT′/h ∼ q min . Making use of these associations and Equations 3 and 4, we can rewrite Equation 9 as The expected scaling of ℓ will depend on which force balance describes the convective dynamics. For our simulations that sit within the IAC regime, the convective length scale is expected (Aubert et al., 2001) to scale as where the flux Rayleigh number Ra F ¼ f Ra=E (Mound & Davies, 2017). Therefore, we expect that the thickness of the regional inversion lenses should scale as In terms of the underlying physical parameters, this scaling becomes The scaling laws for N and L depend on a number of physical parameters, values of which are listed in Table 1. To compare our Boussinesq model with the Earth, it is the superadiabatic heat flow across the CMB that should be used to determine the relevant thermal forcing β = Q + /4πk, a quantity that is poorly constrained with even the sign of Q + uncertain (Jones, 2015;Olson, 2015). Regional inversion lenses will occur only if the combination of Q + and q ptp result in a CMB heat flux pattern that has both superadiabatic and subadiabatic regions. Based on a scaling argument for core velocity, Jones (2011) estimated Q + ≈ 0.6 TW; whereas a comparison of core adiabatic heat flow estimates  and total CMB heat flow estimates (Nimmo, 2015) suggest values as large as Q + ≈ 3 TW are possible. We will use both of these estimates to bound our extrapolations. Similarly, the lateral variation in heat flux, q ptp should be considered relative to the average superadiabatic flux, q þ ave ¼Q þ =4πr 2 o when determining q ⋆ ; here we adopt q ptp = 0.14 W/m 2 (Stackhouse et al., 2015) leading to q ⋆ ≈ 10 or 35 for our chosen values of Q + .

Scaling Results
We have determined the thickness and maximum Brunt-Väisälä frequency for the regional inversion lenses in our simulations at the centre of the lens for the hemispheric boundary forcing and at two locations beneath the CMB for our tomographic boundary forcing. For the tomographic case, we will refer to the locations as African (0°N, 0°E) and Pacific (0°N, 180°E). We first establish that our simulations obey the expected scaling by restricting ourselves to the subset of simulations with a hemispheric pattern of CMB heat flux heterogeneity. For this set of simulations, q min is located where we measure N and L for the regional inversion lens and the pattern of heat flux heterogeneity obeys the adopted geometric assumptions (Equations 3 and 4); therefore, we expect this subset of lenses should best conform to the derived scalings. In Figures 1a and  Dziewonski and Anderson (1981). b Aoki et al. (1982). c Gubbins et al. (2003). d Pozzo et al. (2012).

Geophysical Research Letters
MOUND AND DAVIES 1d, we plot N and L for the hemispheric simulations against the predicted scaling; the agreement between simulations and theory is indeed excellent.
For the tomographic pattern of CMB heterogeneity, q min is located beneath the south Pacific; however, a regional inversion lens may still form beneath Africa provided q ptp is sufficiently large. When considering the African (Figures 1b and 1e) and Pacific (Figures 1c and 1f) regional inversion lenses in our tomographic simulations, the developed scaling laws do not fit the measurements of N and L as well as they do when considering the hemispheric pattern. For the Pacific regional inversion lens characteristics, there is a systematic offset between the q ⋆ = 2.3 and the q ⋆ = 5.0 simulations in their scalings for N and L, perhaps because we have not evaluated the lens properties directly beneath q min and our geometric assumptions about q (Equations 3 and 4) do not hold as well in this case. Nevertheless, each combination of CMB heat flux pattern and lens location does follow the expected scaling and the best-fit prefactors for each lens location do a reasonable job of explaining L and N across all of our simulations falling in the IAC regime.
The scaling fits to our simulations allow us to extrapolate to the conditions relevant to the Earth's core for our two choices of Q + (gray symbols in Figure 1, values in Table 2). For the lower value of superadiabatic CMB heat flow (Q þ low ¼0:6 TW, gray stars) the extrapolated L and N are somewhat larger than for Q þ high ¼3 TW (gray squares). For the chosen tomographic boundary condition, the heat flux low under the Pacific is deeper than the low under Africa. As a result, the predicted thickness and Brunt-Väisälä frequency of the Pacific lens are always larger than the African lens, with scaling predictions of N/2Ω = O(1) and L = O(100) km in all cases.

Discussion
The developed theoretical scalings do a good job of fitting our simulations of regional inversion lenses and allow us to extrapolate to Earth's core conditions, predicting values of L and N that are geophysically plausible. Observational constraints on L and N for a global stratification at the top of the core have been derived from both seismic and geomagnetic observations. Seismic evidence allows a layer of anomalously slow P wave speed up to 450 km thick (Kaneshima, 2018); when combined with a model for chemical enrichment (Helffrich & Kaneshima, 2010), the Brunt-Väisälä frequency is inferred to be N/2Ω ≈ 3.5 − 7.35. Magnetic-Archimedean-Coriolis (MAC) waves in a stable layer 130-140 km thick with N/2Ω ≈ 0.37 − 0.42 have been suggested to explain certain periodic variations of the magnetic field (Buffett et al., 2016). The fundamental difference in our scenario is that stratified lenses arise only under regions of anomalously low CMB heat flux and are absent where the CMB heat flux is superadiabatic. The lateral temperature difference between the stratified and unstratified regions at the top of the core will drive thermal winds within the fluid core. Such flows are seen within our simulations and, when extrapolated to the Earth, suggest a vertical velocity gradient (∂u/∂z) on the order of 10 −5 (m/s)/m (Mound et al., 2019). At depths on the order of 100 m, this corresponds to flow speeds on the order of 1 mm/s, compatible with those inferred from observations of magnetic secular variation for flow at the top of the free stream (e.g., Lesur et al., 2015).
The smaller the total superadiabatic heat flow across the CMB, the broader, stronger, and thicker the regional inversion lenses are expected to be. Ascertaining the existence and extent of regional inversion lenses at the top of the core would, therefore, provide constraints on the values of Q + and q ptp and hence the thermal state of the lowermost mantle. As the pattern and strength of CMB heat flux heterogeneity evolves over geological time in response to ongoing mantle convection, the predicted locations, thicknesses, and strengths of regional inversion lenses would similarly evolve in response.
If a mechanism of chemical stratification is also active within the core, then this additional stabilizing contribution might enable regional inversion lenses to form for smaller values of q ⋆ or promote thicker and stronger lenses. A global stratified layer might even result despite the CMB heat flux heterogeneity, provided that the process of light element enrichment at the top of the core is sufficiently strong to shut down convection beneath regions of superadiabatic CMB heat flow. Determining whether any such mechanism is active in Earth's core and its influence on regional inversion requires additional numerical and theoretical investigation beyond the scope of this study. Extrapolation of regional inversion lens thickness to the Earth requires understanding of the appropriate force balance, both for the Earth and the given suite of simulations. For a given convecting system, the dominant force balance will depend on its physical properties and boundary conditions and can differ between boundary layers and the interior and with the length-scale considered (e.g., Aubert, 2019;Aurnou & King, 2017;Gastine et al., 2016;Grossmann & Lohse, 2000;Schwaiger et al., 2019). The goodness of our fits ( Figure 1) and previous analysis (Long et al., 2020;Mound & Davies, 2017) indicate that the IAC balance holds in the bulk of the fluid interior for the selected simulations. We have run additional simulations (see, Mound & Davies, 2017) that do not fall in the IAC regime. For insufficiently supercritical f Ra , the simulations sit within the weakly nonlinear regime with relatively sluggish convection that may not fill the spherical shell. For sufficiently large forcing, the rotational constraint on convection starts to break down and the simulations sit within the transitional regime between rapidly rotating and nonrotating convection. Although thick regional inversion lenses do form in simulations outside the IAC regime, they do not follow the thickness scaling developed here (Equation 13), which suggests that the global force balance does play an important role in setting L.
Our simulations reach Reynolds numbers (up to O(10 3 )) and Rossby numbers (down to O(10 −4 )) that are characteristic of turbulent rotating convection, as expected in Earth's core.
The global force balance enters into the scaling for L by determining the small length-scale associated with convection. It is possible to consider another force balance, such as one incorporating the flux heterogeneity. The dashed lines and R 2 values are based on the best-fit prefactors for each case and the exponents predicted by our theory. These fits are used to extrapolate to the Earth for Q þ low ¼0:6 TW (star) or Q þ high ¼3 TW (square). influence of the magnetic field on the dynamics within the Earth's core. For example, Davidson (2013) derived a scaling for ℓ assuming a MAC balance holds; substitution of his equation (12) into our Equation 10 results in As with the scaling based on the IAC balance, this scaling depends only weakly on E, Pr, and f Ra and hence the associated physical parameters and has the same expected scaling between L/h and q ⋆ . Therefore, we expect that if this force balance holds, a sufficiently large amplitude of CMB heat flux heterogeneity would result in the presence of thick (O(100) km) regional inversion lenses at the top of the core.

Data Availability Statement
The simulation output data used to produce Figures 1 and 2 can be found in the supplementary information and can be accessed via DOI: 10.17605/OSF.IO/ES47H.