Peter Olson
email: olson@jhu.edu

Maylis Landeau

Evan Reynolds

Dynamo tests for stratification below the core-mantle boundary

Keywords: Partially stratified dynamos, outer core stratification, core-mantle boundary, core heat flux

Evidence from seismology, mineral physics, and core dynamics suggests a layer with an overall stable stratification in the Earth's outer core, possibly thermal in origin, extending below the core-mantle boundary (CMB) for several hundred kilometers. Yet vigorous deep mantle convection with locally elevated heat flux implies unstable thermal stratification in some regions below the CMB, consistent with interpretations of non-dipole geomagnetic field behavior that favor upwelling flows in places below the CMB. To resolve this apparent inconsistency, we investigate the structure of convection and magnetic fields in the core using numerical dynamos with laterally heterogeneous boundary heat flux. Strongly heterogeneous boundary heat flux generates localized convection beneath the CMB that coexists with an overall stable stratification there.

Our partially stratified dynamos are distinguished by their time average magnetic field structures. Without stratification or with stratification confined to a thin layer, the octupole component is small and the CMB magnetic field structure includes polar intensity minima. With more extensive stratification, the octupole component is large and the magnetic field structure includes intense patches or high intensity lobes in the polar regions. Comparisons with the time-averaged geomagnetic field are generally favorable for partial stratification in a thin (< 400 km) layer but unfavorable for stratification in a thick (∼1000 km) layer beneath the CMB.

The possibility of a distinct layer below Earth's core-mantle boundary (CMB), [START_REF] Braginsky | Mac-oscillations of the hidden ocean of the core[END_REF] so-called hidden ocean, has been the subject of numerous studies using a variety of seismic, geomagnetic, and mineral physics approaches, with the twin objectives of resolving the properties of this layer and understanding its dynamical significance. The majority of these studies conclude that the layer extends one to two hundred kilometers below the CMB [START_REF] Whaler | Does the whole of the Earth's core convect?[END_REF][START_REF] Lay | The stably-stratified outermost core revisited[END_REF][START_REF] Gubbins | Geomagnetic constraints on stratification at the top of Earth's core[END_REF][START_REF] Tanaka | Possibility of a low p-wave velocity layer in the outermost core from global smks waveforms[END_REF][START_REF] Buffett | Geomagnetic fluctuations reveal stable stratification at the top of the Earth's core[END_REF] although some claim it extends to far greater depths, perhaps three hundred kilometers [START_REF] Helffrich | Outer-core compositional stratification from observed core wave speed profiles[END_REF][START_REF] Gomi | The high conductivity of iron and thermal evolution of the Earth's core[END_REF][START_REF] Tang | Seismological evidence for a non-monotonic velocity gradient in the topmost outer core[END_REF] or more [START_REF] Kaneshima | Array analyses of SmKS waves and the stratification of Earth's outermost core[END_REF]. Interpretations include stable (subadiabatic) thermal stratification [START_REF] Gomi | The high conductivity of iron and thermal evolution of the Earth's core[END_REF][START_REF] Buffett | Geomagnetic fluctuations reveal stable stratification at the top of the Earth's core[END_REF] as well as stable compositional stratification due to excess light element concentrations in the layer [START_REF] Helffrich | Causes and consequences of outer core stratification[END_REF][START_REF] Gubbins | The stratified layer at the core-mantle boundary caused by barodifusion of oxygen, sulfur and silicon[END_REF].

It is important to point out, however, that not every study supports the existence of such a layer, or at least, there are reasons to doubt that the region below the CMB is uniformly stable to all convective motions. Interpretations of the geomagnetic secular variation are most easily accommodated by core flows including upwelling and downwelling motions that extend to within 100 km below the CMB or shallower [START_REF] Gubbins | Geomagnetic constraints on stratification at the top of Earth's core[END_REF][START_REF] Amit | Can downwelling at the top of the Earth's core be detected in the geomagnetic secular variation?[END_REF][START_REF] Lesur | Are geomagnetic data consistent with stably stratified flow at the core-mantle boundary?[END_REF]). Likewise, the proliferation and rapid evolution of reverse flux spots in the present-day geomagnetic field on the CMB [START_REF] Olsen | The chaos-4 geomagnetic feld model[END_REF] argue for flux expulsion related to upwellings and downwellings [START_REF] Olson | Changes in Earth's dipole[END_REF].

In addition, it is necessary to consider the effects of the mantle circulation on the geodynamo. Mantle global circulation models (Zhang and Rudolph, 2015;[START_REF] Nakagawa | Implications of high core thermal conductivity on Earth's coupled mantle and core evolution[END_REF]2015) predict vigorous deep mantle convection with locally elevated heat flux that is large enough to sustain unstable thermal stratification in some regions beneath the CMB (Olson et al., 2015), even if recent estimates of high thermal conductivity in the core (Ohata et al., 2016) apply. Alternatively, with lower thermal conductivity [START_REF] Konopkova | Direct measurement of thermal 583 conductivity in solid iron at planetary core conditions[END_REF], thermal conditions may be unstable everywhere, but in that case a small accumulation of light elements at the top of the core [START_REF] Buffett | Stratification of the top of the core due to chemical interactions with the mantle[END_REF] could provide the stable stratification.

These issues raise the question of whether it is possible to reconcile seemingly contradictory inferences: a layer providing overall stable stratification on the one hand, with radial motions in the fluid outer core below the CMB on the other. In this paper we address this apparent incongruity using numerical dynamos with a particular suite of boundary conditions that (1) model inner core boundary buoyancy release as the source of the main convection,

(2) produce an overall (i.e., spherically-averaged) stable thermal stratification below the outer boundary, and (3) generate strong lateral heterogeneity in the stratification, including localized convection. We find that this combination of boundary conditions implies a style of convection in the outer core that dynamically maintains stably stratified conditions in limited regions below the CMB, yet allows for radial motions in places as well as generating a dipole-dominant magnetic field. We call these partially stratified dynamos. In addition, we demonstrate that partially stratified dynamos have distinctive high latitude magnetic field structures, allowing the strength of the stratification below the CMB to be inferred remotely, using the geomagnetic field on the CMB.

Our study is limited to the types of stratification that are produced when the destabilizing effects of inner core boundary buoyancy release are comparable in the stratified region to the stabilizing effects of subadiabatic CMB heat flux. This regime has been explored previously using numerical dynamos with homogeneous outer boundary conditions [START_REF] Christensen | Models of magnetic field generation in partly stable planetary cores: Applications to Mercury and Saturn[END_REF][START_REF] Nakagawa | 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[END_REF][START_REF] Nakagawa | An implication for the origin of stratification below the core-mantle boundary region in numerical dynamo simulations in a rotating spherical shell[END_REF] and magnetoconvection models [START_REF] Takehiro | Penetration of steady fluid motions into an outer stable layer excited by MHD thermal convection in rotating spherical shells[END_REF]. It has been shown that stable stratification tends to filter the non-axisymmetric non-dipolar fields, and if the stratified layer is thick, also reduces the intensity of the axial dipole field [START_REF] Christensen | Models of magnetic field generation in partly stable planetary cores: Applications to Mercury and Saturn[END_REF][START_REF] Nakagawa | 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[END_REF]2015). [START_REF] Christensen | Geodynamo models with a thick stable layer and heterogeneous CMB heat flow[END_REF] used combinations of subadiabatic mean boundary heat flux plus lateral boundary variations to produce dynamos with stratification extending below the CMB to 20-40% of the outer core depth. Under these conditions he finds thin horizontal circulations that mediate the boundary heat flux heterogeneity, but little mixing of the stratification. In contrast, stratified magnetoconvection calculations by [START_REF] Takehiro | Penetration of steady fluid motions into an outer stable layer excited by MHD thermal convection in rotating spherical shells[END_REF] produce strong flows capable of penetrating most of the stable region.

The stratification analyzed in this study refers to radial density gradients that deviate from adiabatic (i.e., uniform entropy) thermal conditions. Temperature gradients resulting from self-compression of the fluid are therefore excluded from our dynamo calculations, and must be factored in before comparing our results to Earth's core. Dynamos that explicitly include adiabatic density and thermal gradients due to compressibility of the fluid also show stratification effects, particularly when the adiabatic density variation is large across the fluid [START_REF] Jones | Anelastic convection-driven dynamo benchmarks[END_REF]Gastine et al., 2012[START_REF] Yadav | Consistent scaling laws in anelastic spherical shell dynamos[END_REF]. Because the density scale height of the outer core is greater than its depth, the direct effects of compressibility are not included in our dynamos. In addition, our study does not consider the situation in which the stabilizing effects of stratification vastly outweigh the destabilizing effects of inner core growth, as would be the case for strong, pre-existing compositional stratification [START_REF] Landeau | Core merging and stratification after giant impacts[END_REF] or compositional stratification that develops over time through rapid chemical diffusion [START_REF] Nakagawa | On the thermo-chemical origin of the stratified region at the top of the Earth's core[END_REF]. With such strong stratification, lateral variations in heat flux at the CMB would likely play a more limited role in determining the global structure of the outer core and its overall dynamical behavior.

Partially stratified dynamos

To model dynamo action with thermal and compositional buoyancy originating at the inner core boundary (ICB) due to inner core growth plus dynamically-regulated thermal stratification below the CMB, all within the context of the Boussinesq approximation, we follow standard procedures [START_REF] Jones | Thermal and compositional convection in the core[END_REF], defining the codensity C in the outer core as the sum of densities due to temperature and light element concentration:

C = ρ o (αT + βχ) , (1)
where ρ o is outer core average density, T is temperature relative to the core adiabat, χ is the outer core light element concentration relative to its mean, and α and β are volumetric expansion coefficients for T and χ, respectively. We let χo and Ṫo denote the time rate-ofchange of the mean (volume-averaged) light element concentration and temperature of the outer core, each assumed to be constant over the time span of a single dynamo simulation, so that Ċo = ρ o (α Ṫo +β χo ) is the volume-averaged rate-of-change of codensity over a simulation, also assumed constant. Further, let Ω denote angular velocity of Earth's rotation, g gravity at the CMB, D = r cmb -r icb the depth of the outer core fluid, r cmb and r icb being the radii of the CMB and the ICB, respectively, and let ν and κ denote outer core kinematic viscosity and codensity diffusivity, respectively.

With these definitions, the Boussinesq equations for conservation of momentum including rotation, conservation of mass, and codensity transport in a rotating spherical shell (see

E = ν ΩD 2 ; P r = ν κ ; Ra = βgD 5 χo ν 2 κ . ( 2
)
Here E is the Ekman number, P r is the Prandtl number, and Ra is the Rayleigh number, and the factors D, D 2 /ν and D 2 ρβ χo /ν scale length, time, and codensity variations, respectively.

Two additional control parameters appearing in the magnetic induction equation and the codensity equation are the magnetic Prandtl number P m and the dimensionless volumetric codensity source/sink that quantifies the rate of buoyancy absorbed in the outer core from the mixing of light elements, secular cooling of the outer core, curvature of the core adiabat, and radioactive heat sources. Our dynamos are driven by the combination of light element release at r icb and secular cooling, so that

= -(1 + α Ṫo β χo ).
(3)

The magnetic Prandtl number P m is defined by

P m = ν η (4)
where η is the magnetic diffusivity of the outer core. Magnetic fields are scaled by ρ o Ω/σ, where σ is electrical conductivity.

Additional control parameters arise in defining the boundary conditions. At the ICB we assume no-slip velocity conditions and a uniform codensity, C icb . At the CMB also assume no-slip velocity conditions, zero compositional flux, and we specify the heat flux there to be the sum of a spherical mean part (denoted by an overbar) and a deviation from the spherical mean (denoted by a prime):

q = q + q (φ, θ) (5)
where φ and θ are longitude and co-latitude, respectively, and q is measured relative to the heat flux down the core adiabat, with q > 0 being superadiabatic heat flux and q < 0 being subadiabatic heat flux. The variable q in (5) specifies the pattern and the amplitude of the CMB heat flux heterogeneity. In the same way we can write the codensity as the sum of a spherical mean part C and a laterally varying part C . Then using (1) and ( 5) and assuming Fourier's law for conduction, the spherical mean and laterally varying codensity gradients on the CMB can be written as

∂ C ∂r = - αρ 0 k q (6)
and

∂C ∂r = - αρ 0 k q (7)
where r is the radial coordinate and k is the thermal conductivity of the outer core, assumed constant. In dimensionless form ( 6) and ( 7) become

∂ C * ∂r * = Ra q Ra = S, (8)
where

S = - αν Dk χo q (9)
is the dimensionless stratification parameter we use in this study, and

∂C * ∂r * = Ra q Ra L(θ, φ) = S L. (10)
In ( 8)

Ra q = - αgD 4 kνκ q (11)
is the Rayleigh number based on the spherical mean CMB heat flux, and in ( 10)

Ra q = αgD 4 kνκ ∆q (12)
is the Rayleigh number based on the peak-to-peak variation ∆q of the laterally varying CMB heat flux, L(θ, φ) being the pattern (or planform) of its lateral variation in terms of colatitude θ and longitude φ, S is the boundary heterogeneity counterpart to S, and the asterisks denote dimensionless variables. Note that both S and Ra q are defined to be positive when q is negative, i.e., when the mean CMB heat flux is stabilizing.

In this study we restrict consideration to dynamos with P r = 1, P m = 6, and = -0.8, the latter appropriate for dominantly compositional convection but with some secular cooling. In most cases we choose Ekman numbers E = 1 × 10 -4 , with a few cases at E = 3 × 10 -5 . We fix the aspect ratio of the fluid outer core to be r icb /r cmb = 0.351.

The solid region r ≤ r icb representing the inner core is assumed to have the same electrical conductivity σ as the fluid, and the solid region r ≥ r cmb representing the mantle is assumed to be electrically insulating.

By fixing the codensity on the inner boundary, we allow the codensity flux there to freely adjust, which permits it to respond to the heterogeneous CMB heat flux in a dynamically consistent way. For modeling convective interaction between the inner core and the mantle, this boundary condition has advantages over other simplified boundary conditions, such as fixing the ICB codensity flux a-priori, because it naturally yields lateral variations in the ICB codensity flux, which is a dynamo model proxy for lateral variations in the rate of inner core growth. Its main disadvantage in the context of our study is that it requires two global Rayleigh numbers with different definitions to characterize these dynamos: Ra q for the CMB forcing that depends inversely on viscosity and Ra for the ICB forcing that depends inversely on viscosity squared. Their ratio, the stratification parameter S, therefore depends linearly on the fluid viscosity, which unfortunately is poorly constrained in the outer core.

Similarly, by setting the sink function to a constant we let the dynamics determine the thickness and the gravitational stability of the stratified region, through interactions between the control parameters Ra, S, S . An alternative approach is to prescribe the thickness of the stratification by tuning the sink function, as in [START_REF] Nakagawa | An implication for the origin of stratification below the core-mantle boundary region in numerical dynamo simulations in a rotating spherical shell[END_REF]. Later we show that these two approaches yield broadly consistent results in terms of the time average dynamo magnetic field structure.

We examine dynamo action over ranges of the control parameters Ra and S, for three different CMB heat flux planforms L. The first planform, denoted by L0, corresponds to uniform CMB heat flux, with C * = 0 everywhere on the outer boundary. These uniform boundary cases serve as references for the heterogeneous boundary cases. The second planform, denoted by L2, is defined by C * on the outer boundary consisting of spherical harmonic degree two components, and is shown in Figure 1a. The third planform, denoted by L4, consists of C * components up to and including spherical harmonic degree four, is shown in Figure 1b.

The L4 planform in Figure 1b was obtained from the present-day CMB heat flux pattern produced by the mantle GCM (mantle global circulation model) HF5 of [START_REF] Rudolph | History and dynamics of net rotation of the mantle and lithosphere[END_REF] and [START_REF] Zhong | On the temporal evolution of long-wavelength mantle structure of the Earth's mantle since the Early Paleozoic[END_REF], after truncating that CMB heat flux pattern at spherical harmonic degree and order four. The mantle GCM HF5 includes variable mantle viscosity, compositionally dense material at the mantle base, plus plate motion surface velocity constraints starting around 400 Ma and continuing to present-day. This particular mantle GCM has been used by Olson et al. (2015) to model core evolution. The mean and standard deviation of its present-day CMB heat flux are 79 ± 24 mW/m 2 , and its maximum and minimum are 131 and 36 mW/m 2 , respectively. The L2 planform in Figure 1a is a modified version of L4 planform using only spherical harmonic degree =2 components at orders m=0 and m=2, adjusted to create a pattern with bilateral (i.e., 2-fold) azimuthal symmetry.

The L2 pattern in Figure 1a corresponds to the largest scale, lower mantle heterogeneity structure advocated by some seismologists [START_REF] Dziewonski | Mantle anchor structure: an argument for bottom up tectonics[END_REF]. It has the same mean value as the L4 pattern, whereas its extreme values are 137 and 31 mW/m 2 , respectively. It is essentially the same planform used by [START_REF] Olson | Mantle superplumes induce geomagnetic superchrons[END_REF] in their study of the effects of dense basal mantle piles on magnetic polarity reversal behavior.

The stratification parameter S in our dynamos are determined according to (8) and ( 9).

The amplitudes of the boundary heterogeneity S of the L2 and L4 patterns in our dynamos are scaled using the following ratio:

Γ = |S L| ∂ C * ∂r * icb -1 , (13)
where || denotes peak-to-peak variation. The Γ ratio can be estimated by combining the results of our mantle GCM with a model of the thermodynamic state of the core. We use a core state model with the CMB heat flux obtained from the mantle GCM (which includes an adiabatic part, assumed to be 96 mW/m 2 ) to estimate the implied CMB and ICB codensity fluxes and hence the ratio Γ. We find that the mean cmb heat flux from mantle GCMs with plate motion constraints [START_REF] Zhong | On the temporal evolution of long-wavelength mantle structure of the Earth's mantle since the Early Paleozoic[END_REF] is typically comparable to the adiabatic core heat flux, such that Γ (r icb /r cmb ) 2 0.123, approximately. We calculate the spherical mean ICB codensity flux from our numerical dynamo at a given Ra with Ra q = Ra q =0, and then adjust |S L| in (13) so that Γ=0.123 for that dynamo. The above steps yield a constant value of |S L|=0.58 for all the L2 dynamo cases. A similar procedure is followed for the L4 dynamo cases. Amplitude coefficients of the spherical harmonics that generate these L2 and L4 heat flux planforms for our dynamos are given in Table A1.

We have carried out a systematic parameter sweep of stratified thermochemical convection and dynamo action with the control parameters just described, using the MAGIC code [START_REF] Wicht | Inner-core conductivity in numerical dynamo simulations[END_REF], varying the parameters Ra and S for the three outer boundary heat flux heterogeneity patterns L0, L2 and L4. Table A2 summarizes the control parameter ranges of our sweep. Most of the calculations were done at E = 1 × 10 -4 and run for at least 5 viscous diffusion times in order that the run averages approximate true time averages. For these calculations we used a numerical grid with (n r , n θ , n φ ) = (81, 128, 256) in the fluid shell and spherical harmonic truncation ( , m) max =85. We also ran one case at Ra = 2 × 10 7 with S=1.0, which was strongly subcritical for convection. The majority of our calculations used the L2 boundary heterogeneity, although uniform L0 as well as some L4 cases were included for comparison purposes. Overall, we found that the L4 cases mostly added shorter wavelength heterogeneity to the L2 case results, somewhat complicating their interpretation without changing their behavior in a fundamental way. We also computed three S = 0.1

cases at E = 3 × 10 -5 with L0, L2, and L4 boundary heterogeneity, respectively, which are given in in Table A2. For the smaller E calculations we used a numerical grid with (n r , n θ , n φ )

= (121,192,384) in the fluid shell and spherical harmonic truncation ( , m) max =128. tangent cylinder in each hemisphere, plus conspicuous polar intensity minima in both north and south hemispheres. In contrast, the high latitude radial magnetic field in the S = 0.2 dynamo consists of patches of high intensity field that terminate short of the pole and no polar intensity minima, while the high latitude radial magnetic field in the S = 0.3 dynamo consists of lobes of high intensity field that extend all the way to the poles, forming polar intensity maxima.

Dynamo results

The differences in these field structures can be directly attributed to differences in the high latitude pattern of radial velocity in the two cases. As seen in Figure 5b, there are strong polar upwellings in the S = 0 case, whereas Figure 5e shows much weaker polar upwellings in the S=0.2 case, and Figure 5h shows polar downwellings in the more strongly stabilized S = 0.3 case. Accordingly, the magnetic field becomes concentrated near the pole in the most stratified case, the two high intensity flux lobes being located along longitudes that approximately correspond to the bands of radial downward flow, whereas in the neutral S=0 case the field is mostly concentrated by the circular downwellings that occur along the tangent cylinder. The intermediate S=0.2 case produces a field structure that results from a mixture of the flows seen in the two more extreme cases.

Patterns of the azimuthally and time averaged internal structure are shown in Figures 5c,f, and i for each of the three dynamos. These images reveal the underlying dynamics that produce distinct radial velocity structures, which in turn produce the polar magnetic intensity minima, patch, or lobe radial magnetic field structures. In the S=0 case the tangent cylinder region is supercritical for convection and has the familiar combination of polar upwellings plus tangent cylinder downwellings in its azimuthal averaged flow structure [START_REF] Olson | A polar vortex in the Earth's core[END_REF][START_REF] Sreenivasan | Structure and dynamics of the polar vortex in the earth's core[END_REF]. Together these up-and downflows create a circulation that tends to expel poloidal magnetic field from inside the tangent cylinder region. In contrast, in the strongly stabilized S=0.3 case the tangent cylinder is (or appears to be) subcritical for free convection and the flow directions are reversed there relative to the S=0 case. Magnetic flux concentration, rather than flux expulsion, is active inside the tangent cylinder in this case. In the intermediate S=0.2 case there is a weak polar upwelling, slightly concentrating the field inside the tangent cylinder. Consequently, the radial field becomes more concentrated in the polar regions with stronger stratification, eliminating the polar intensity minimum and replacing it with high intensity field that, with the L2 boundary heterogeneity, generates high intensity patches or lobes. The high intensity flux lobes or patches are located approximately at the longitudes that correspond to the bands of radial downward flow (Figures 5e andh). We note that the dynamics of the flows that underlie the structural changes in the field here are similar to what has been found in dynamos when the ratio of inner boundary-to-outer boundary radius is changed, as is expected to be the situation before versus after the inner core nucleates [START_REF] Landeau | The signature of inner-core nucleation on the geodynamo[END_REF].

As with the magnetic field structure, major differences in the structure of the internal azimuthal flows in Figures 5c,f, and i lie at polar latitudes, where the neutral S=0 case includes a strong retrograde (westward) polar vortex that extends to the outer boundary, whereas in the S=0.3 case the flow in the same region is strongly attenuated in amplitude and weakly prograde (eastward). Equatorial westward jets are present in all three cases, but they are stronger with stratification, and unlike the polar flows, penetrate to near the outer boundary, even in the most stabilized S=0.3 case. Qualitatively, this behavior is in accord with previous findings by [START_REF] Nakagawa | 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[END_REF] and [START_REF] Takehiro | Penetration of steady fluid motions into an outer stable layer excited by MHD thermal convection in rotating spherical shells[END_REF].

In order to systematize how these morphologic differences in the magnetic field structure vary with our dynamo control parameters, we show in Figure 6 and that these structures are particularly sensitive to the stratification parameter S, although some weaker dependence on the Rayleigh number Ra is evident from the locations of the transitional cases. In particular, the patches P-regime appears to attenuate and possibly disappear at strongly supercritical Rayleigh numbers. With more modest supercriticality, the sequence with increasing S consists of M for unstable, neutral, or very weakly stratified cases, transitioning to P with slightly stronger stratification, then to L and finally to N when the boundary-induced stratification is strong enough to kill dynamo action. With reference to Table A2, our three S = 0.1 cases at E = 3 × 10 -5 and Ra = 2 × 10 8 are all in the M-regime, consistent with the results shown in Figure 6.

The distinctions between M, P, and L field structures involves arbitrary considerations in some of the transitional cases, but in most situations, identification can be made on the basis of contours of B r on the dynamo outer boundary in the polar regions. Figure 7 shows contours of the time-averaged radial magnetic field intensity on the outer boundary for the three dynamos shown in Figure 5. Black curves are contours that enclose the high field intensity regions in each case. If two contours are needed to enclose the B r -intensity maxima and both contours include the pole, the structure is classified as M-type, as in Figure 7a. Note that the two bounding contours should not be very different in radius, so as to exclude from consideration very small, inconsequential field intensity minima. Alternatively, if two bounding contours are needed and neither contour includes the pole, the structure is classified as P-type, as in Figure 7b. Finally, if a single bounding contour encloses the intense field and also includes the pole, the structure is classified as L-type, as in Figure 7c.

The three magnetic field structures in Figure 7 are distinct because of the combined action of the stratification measured by S or alternatively Ra q, plus the L2 lateral boundary heterogeneity measured by S or alternatively Ra q . Other parameter combinations could lead to different results. For example, eliminating the lateral boundary heterogeneity yields axisymmetric time averaged magnetic structures, so that the M-and L-structures transform to axisymmetric polar minima and axisymmetric polar maxima, respectively, while the regime with clearly defined P-structures probably gets lost. The radial magnetic field structures in figures 3 and 5 of [START_REF] Nakagawa | An implication for the origin of stratification below the core-mantle boundary region in numerical dynamo simulations in a rotating spherical shell[END_REF], calculated using homogeneous outer boundary conditions in the same range of Ekman and magnetic Prandtl numbers, support this interpretation. This behavior, along with the inference from Figure 6 that the P-regime becomes attenuated at highly supercritical Ra suggests that the high intensity patch regime may be less robust than polar intensity minima or maxima. In addition, reducing the symmetry of the boundary heterogeneity can produce different results, by obscuring the clear-cut differences in the high latitude structures found with L2 symmetry. For example, we find that the loss of bilateral symmetry with L4 boundary heterogeneity yields a nominally P-style field morphology, but with a single rather than two high field intensity patches. However, in spite of these complications, our results suggest it is possible to probe Earth's outer core for stratification using time average geomagnetic field morphology in the polar regions.

Figure 8 shows additional stratification diagnostics that involve ratios G 0 of the axial (m=0) Gauss coefficients of the magnetic field, obtained by time-averaging our partially stratified dynamos. Figure 8a shows axially symmetric (m=0) Gauss coefficient ratios G 0 = g( , 0)/g(1, 0) versus spherical harmonic degree from time-averaged dynamos with Ra = 6 × 10 7 and L2 boundary heterogeneity, for various stratifications S. Figure 8b shows Gauss coefficient ratio G30 (axial octupole over axial dipole) versus stratification S for the same dynamos. M,P,L denote polar minima, patches, and lobes, respectively, in the high latitude magnetic field structure. Stable/unstable refers to the sign of the gravitational stability below the outer boundary. In Figure 8a the only appreciable (i.e., observable) G-ratio is G30. The G50 ratios change sign with stratification but are too small to be observed. In

Stratification scaling

Before applying our results to thermal stratification in Earth's outer core, a necessary first step is to derive scaling laws that summarize the strength and extent of the stratification in our numerical dynamos as a function of the control parameters. The two control parameters that we varied substantially are the Rayleigh number Ra and the stratification parameter S, or alternatively, the boundary Rayleigh number Ra q. Furthermore, most of our dynamos used the L2 boundary heat flux heterogeneity with the other control parameters fixed (that is, E = 1 × 10 -4 , P r = 1, P m = 6). Accordingly, our fits are biased toward dynamos of this type. In addition, we focus attention on two parameterizations of the stratification: (1) the spherically averaged thickness of the stratified region beneath the outer boundary and (2) the spherically averaged gravitational stability of the layer, both given in Table A2 for each case.

We define the spherically averaged dimensionless thickness of the stratified region to be

δ * = r cmb -r min r cmb (14)
where r cmb is, as before, the outer boundary radius and r min is the radius where C * reaches its local minimum value below the outer boundary, C * min , as seen in Figure 3d. Similarly, we define the characteristic gravitational stability of this region in terms of the following parameter:

N * 2 = δC * δ * , (15)
where δC * = C * cmb -C * min is the dimensionless codensity increase across the stratified region. The choice of the notation here is motivated by the connection between N * and the buoyancy frequency in stratified fluids [START_REF] Turner | Buoyancy effects in fluids[END_REF], the asterisk emphasizing that this is a nondimensional parameter.

Figure 9 shows fits of the dimensionless stratified layer thickness to power laws in Ra, Ra q, and S. In Figure 9a we assume a power law of the form

δ * = aRa b qRa c ( 16
)
where a is a (constant) coefficient and b and c are (constant) exponents. The fit was obtained by minimizing the function

F = δ * 2 σ 2 δ (log δ * -log a -b log Ra q -c log Ra) 2 , ( 17
)
where σ δ , the characteristic uncertainty on δ * , was assumed to be the same for all cases. The symbols in Figure 9a have the same meaning as in Figure 2 The only systematic discrepancies between stratified layer thickness and our scaling occur at low S values, for cases with heterogeneous heat flux. This can be explained as follows.

The lateral heat flux variations produce a mix of regions that are convectively unstable below the outer boundary and regions that are convectively stable. Because the thermal boundary layer in convecting regions is thinner than the stratified layer in stable regions, the volumeaveraged codensity field is biased toward the properties of the stable region. Accordingly, the globally averaged stratified thickness is non-zero with strong heterogeneous boundary conditions, even in cases with S=0.

The gravitational stability immediately below the CMB is measured by the squared buoyancy frequency defined in terms of CMB heat fluxes:

N 2 cmb = αg k (q ad -qtotal ), (19)
where q ad and qtotal refer to the adiabatic and mean CMB heat fluxes, the latter including the adiabatic contribution. The property values in Table A3 yield N 2 cmb = 2.4 × 10 -8 (rad/s) 2 , so that the ratio of the squared buoyancy frequency to the square of the Coriolis parameter is approximately (N cmb /2Ω) 2 1.1. However, this measure of the gravitational stability applies only in a very limited depth range beneath the CMB, essentially within a thin boundary layer region where the codensity gradient is conductive. In particular, it does not apply over the entire depth range of the stratification, because the codensity gradient there is affected by convection.

A better measure of the overall gravitational stability is the average of the buoyancy frequency over the stratified region, given by (15). Figure 10 shows the dependence of N * 2 defined by (15) on Ra, Ra q, and S for all cases with S <0.4 The symbols and color scheme are the same as in Figure 9. Figure 10a shows that the best-fitting constants in a power-law relationship with Ra and Ra q similar to ( 16) are (a, b, c)=(1.1,1.05,-1.07). Within uncertainties, the condition c = -b again holds in this case, indicating that the gravitational stability also depends primarily on S. The latter result implies that the mean stratification across the layer is proportional to the stabilizing codensity gradient imposed at the CMB .

Refitting the stability data to

N * 2 = a N S b (20)
with b = 1 yields a N = 0.72, the fit shown by the dashed line in Figure 10b. The deviations from linear scaling at S ≥0.4 result from weak or non-existent convection from the inner core buoyancy source.

Extrapolation to the outer core

Table A3 gives values of the physical properties needed to calculate the stratification parameter in the outer core according to the definition of S in (9). We use the core state model described in Olson et al. (2015) with the adiabatic thermal gradient as modified by [START_REF] Labrosse | Thermal evolution of the core with a high thermal conductivity[END_REF], which for a thermal conductivity k=100 W/m/K gives and adiabatic heat flux at the CMB of q ad = 96 mW/m 2 . We select qtotal = 79 mW/m 2 , representative of the mean CMB heat flux from mantle GCMs [START_REF] Nakagawa | Influence of combined primordial layering and recycled MORB on the coupled thermal evolution of Earth's mantle and core[END_REF][START_REF] Zhong | On the temporal evolution of long-wavelength mantle structure of the Earth's mantle since the Early Paleozoic[END_REF], which includes the adiabatic conductive contribution. Together these imply a subadiabatic heat flux at the CMB, with q = qtotal -q ad = -17 mW/m 2 . The core state model then yields Ċo =6.8×10 -16 kg/s for the rate of codensity change in the outer core due to the combined effects of cooling and light element increase.

Substituting Ċo , q, and the other parameter values from Table A3 into (9) yields a low value of S=0.018 for an outer core viscosity of ν = 2 × 10 -6 m 2 /s, an intermediate value of S=0.064 for ν = 3.5 × 10 -6 m 2 /s and a high value of S=0.2 for an outer core viscosity of ν = 2.2 × 10 -5 m 2 /s. Here we have adjusted the diffusivity so that κ = 10ν in all cases.

Application of our scaling law (18) with a δ =1.82 and b=1.2 yields stratified layer thickness predictions of δ = 52 km, 230 km, and 920 km, respectively, for these three viscosity choices.

The stratified layer thickness also depends sensitively on qtotal . For example, reducing qtotal to 50 mW/m 2 and assuming an outer core kinematic viscosity of ν = 2×10 -6 or ν = 7×10 -6 m 2 /s along with the other properties in Table A3, the predicted stratified layer thickness would be 255 km or 1150 km, respectively.

The dependence of stratified layer thickness on viscosity in our scaling is a consequence of the definition of S, which is the ratio of two Rayleigh numbers, one based on the outer boundary codensity flux, the other based on the rate of codensity change in the outer core as a whole. It could be argued that Ohmic rather than viscous dissipation should control the layer thickness in the Earth's core, and the dependence of layer thickness on viscosity is therefore a model artifact. Nevertheless, the above range encompasses most previous estimates of the thickness of the outer core stratified region, and as we demonstrate below, the gravitational stability of the stratified layer predicted by our scaling is nearly independent of viscosity.

To apply our scalings to the gravitational stability of a stratified layer in the outer core, we combine ( 9) and ( 20) to get, in terms of dimensional properties,

N 2 ave = a N αg k (q ad -qtotal ). ( 21
)
The property values listed in Table A3 give N 2 ave = 1.7 × 10 -8 (rad/s) 2 . Note that the layer average gravitational stability ( 21) is reduced from its value immediately below the CMB (20) by the factor a N =0.72, but otherwise its dependence on the properties of the outer core remains the same. Also, (21) implies that, unlike the thickness of the stratified layer, the gravitational stability does not depend explicitly on the inner core buoyancy source or the outer core viscosity. This is a consequence of our ignoring the (very weak) dependence of N * 2 on Ra in Figure 10. Finally, there appears to little dependence of the stratification parameters on the Ekman number, especially the buoyancy frequency. For example, N * 2 = 0.0649 for the L2, S = 0.1 E = 3 × 10 -5 case in Table A2, compared to an average value of N * 2 = 0.0678 all of the L2 S = 0.1 cases at E = 1 × 10 -4 , just a 4% change in gravitational stability accompanying three-fold change in rotation rate.

6 Testing for outer core stratification

Our results show that the combination of subadiabatic average heat flux plus large amplitude lateral variations in heat flux generates an overall stable stratification below the CMB, yet allows for localized radial motions where the CMB heat flux is particularly large. This partial stratification contrasts with the usual assumption of homogeneously stratified conditions, in which radial motions would be suppressed uniformly.

Our partially stratified dynamos show clearly identifiable transitions in the structure of the time-averaged radial magnetic field in polar regions that correspond to increasing strength of stable stratification, starting from polar minima without stratification (or with weak stratification), to patches of intense field for moderate stratification, finally to lobes of intense field for stronger stratification. Dynamos with homogeneous outer boundaries lack permanent patch or lobe structures, and switch from polar minima to polar maxima with increasing stratification [START_REF] Nakagawa | An implication for the origin of stratification below the core-mantle boundary region in numerical dynamo simulations in a rotating spherical shell[END_REF].

The structural transitions we find with heterogeneous outer boundary conditions offer the possibility of inferring the stratification below the CMB based on the morphology of the time-averaged geomagnetic field on the CMB. In order to realize this possibility, however, it is necessary to image the high latitude geomagnetic field with sufficient resolution and over a sufficiently long period of time, in order to image its nonaxisymmetric structure. It is not obvious that current models of the time average geomagnetic field satisfy these requirements.

The present-day geomagnetic field on the CMB has been imaged up to a maximum spherical harmonic degree of approximately max =14, as shown in Figure 12a. At this resolution, polar minima are clearly evident, and there is even reversed magnetic flux in the polar regions of both hemispheres. The northern hemisphere structure has previously been interpreted as evidence of a convective upwelling associated with a polar vortex [START_REF] Olson | A polar vortex in the Earth's core[END_REF].

When the geomagnetic field is time-averaged over the past 400 yr, as in field model gufm1 by [START_REF] Jonkers | Four centuries of geomagnetic secular variation from historical records[END_REF], both polar minima are strongly attenuated compared to the presentday, and the time-averaged field structure is more patch-like at high latitudes. This trend continues in field reconstructions with increasingly long time averages. In archeomagnetic field reconstructions such as CALS10k [START_REF] Korte | Reconstructing the Holocene Geomagnetic Field[END_REF], time averaging over thousands of years results in high latitude field structures that are distinctly more lobe-like, compared to the present-day, with no expression of polar minima. Finally, paleomagnetic field reconstructions that average the geomagnetic field over the past five million years (0-5Ma) typically only show broad, reduced intensity lobe structures in the polar regions [START_REF] Johnson | The time-averaged geomagnetic field as recorded by lava flows over the last 5 Ma[END_REF]. But even these modest deviations from axial symmetry have been questioned, as it remains unclear that the paleomagnetic data absolutely requires their existence [START_REF] Johnson | The Time-Averaged Field and Paleosecular Variation[END_REF].

Superficially, the tendency for geomagnetic and paleomagnetic field reconstructions to exhibit broad, high latitude lobes with increasingly long time averaging would implicate moderate or strong stratification beneath the CMB, characterized by a stratification parameter S > 0.2, according to our results. However, there is an alternative interpretation.

Disappearance of the polar intensity minima and the transition to broad lobe structures when averaged are taken over increasingly long times might simply be a consequence of a lack of adequate spatial coverage in the geomagnetic and paleomagnetic data, especially at high latitudes.

To illustrate this effect, Figure 11 shows time average radial magnetic fields on the outer boundary from the Table A2 dynamo with Ra = 6 × 10 7 , L2 boundary heterogeneity and neutral (S=0) stratification parameter. Panel a shows the field structure with truncation at spherical harmonic degree and order ( , m) max = 24; panel b is the same field with truncation at ( , m) max =12, whereas panel c is the same field truncated at ( , m) max =6. This figure illustrates the hazards of using limited-resolution geomagnetic field representations to interpret high latitude structure on the CMB in terms of outer core stratification. The fully-resolved field structure is squarely in the M (polar minima) regime, as are the weakly filtered images in Figure 11a andb. However, the strongly filtered image Figure 11c removes the polar minima, making the high latitude field appear more like the L (lobe) regime.

Similar trends emerge when the present-day geomagnetic field is subjected to increasingly severe spherical harmonic truncation. Figure 12 Because full resolution of the ancient core field is problematic, it is useful to consider simpler diagnostics that would indicate stratification. Axially symmetric departures from an axial dipole have been extensively investigated over the entire paleomagnetic record [START_REF] Evans | Test of the dipolar nature of the geomagnetic field throughout Phanerozoic time[END_REF]McFadden and Reid, 1982;[START_REF] Veikkolainen | On the low-inclination bias of the Precambrian geomagnetic field[END_REF][START_REF] Johnson | The Time-Averaged Field and Paleosecular Variation[END_REF] The general consensus is that the ratio of axial quadrupole to axial dipole Gauss coefficients, denoted by G20, is nonzero at several epochs, whereas the ratio of axial octupole to axial dipole Gauss coefficients, denoted by G30, is quite small when averaged over the past 0-5 Ma, less than 0.05 in absolute value, with considerable debate concerning its sign. This is generally consistent with the more recent history of G30, which is approximately -0.045 when averaged over the past decade and even smaller when averaged over the past ten millennia, -0.006 according to archeomagnetic field reconstruction CALS10k [START_REF] Korte | Reconstructing the Holocene Geomagnetic Field[END_REF]. When averaged over separate polarity chrons, paleomagnetic inclination data indicate G30 0.01 during the Brunhes chron, but somewhat larger G30 0.05 during the Matuyama chron [START_REF] Aubert | Observations and models of the long-term evolution of Earth's magnetic field[END_REF]. There is some evidence that G30 may have been far larger in the Precambrian [START_REF] Evans | Testing the GAD throughout geological time[END_REF][START_REF] Veikkolainen | On the low-inclination bias of the Precambrian geomagnetic field[END_REF] but back then the influence of inner core buoyancy release may have been less or non-existent. Indeed, numerical dynamos with a passive inner core (Heimpel and Evans, 2013) or no inner core at all [START_REF] Landeau | The signature of inner-core nucleation on the geodynamo[END_REF] generally produce larger and more positive values of G30 compared to similar dynamos with an active inner core.

The key question here is whether this evidence of a modest-sized G30 in the time-averaged field is large enough to validate stable stratification beneath the CMB. According to Figure 8b, the observation that G30 ≤ 0.05 when the geomagnetic field is averaged over the past few million years suggests a stratification parameter beneath the CMB of S <0.1. According to (18), the size of G30 constrains the thickness of such a layer to 400 km or less, approximately, and would seem to preclude very thick thermal stratification, such as the 1000 km layer proposed by [START_REF] Gomi | The high conductivity of iron and thermal evolution of the Earth's core[END_REF], for which S 0.2 and G30 0.15. However, this conclusion is based on a limited set of dynamo calculations, and further exploration of partially stratified dynamos with different boundary conditions and control parameters is needed for confirmation.

Core flow inversions offer another way to probe for outer core stratification in localized regions beneath the CMB, including polar regions. Typically, flow core inversions based on the assumption of frozen magnetic flux find westward zonal flow in the polar regions (e.g. [START_REF] Eymin | On core surface flows inferred from satellite magnetic data[END_REF][START_REF] Olson | Changes in Earth's dipole[END_REF], as in our low S cases, but unlike our higher S cases. By this measure, our neutral S = 0 or stratified S = 0.1 dynamos seem more Earthlike in terms of core flow, compared to the larger S dynamos. Is this strong enough evidence to say that the outermost core is at most weakly stratified? A recent inversion for core flow concluded that a purely horizontal (i.e., toroidal) flow fails to explain the global geomagnetic basically consistent with weak stratification. Similarly, an even more localized frozen flux study by [START_REF] Chulliat | Magnetic flux expulsion from the core as a possible cause of the unusually large acceleration of the north magnetic pole during the 1990s[END_REF] confined to the north polar region argued for radial magnetic diffusion supported by an underlying polar upwelling, i.e., at most weak stratification, in basic agreement with [START_REF] Olson | A polar vortex in the Earth's core[END_REF].

We can also compare our results with inferences of outer core stratification derived from idealized analytical models of thermochemical convection and from dynamical interpretations of the geomagnetic secular variation. According to the convection model of [START_REF] Lister | Stratification of the outer core at the core-mantle boundary[END_REF], the present-day core is predicted to have δ 190 km of thermal stratification for the adiabatic and total heat fluxes listed in Table 1. For these same heat fluxes their model predicts a maximum value of (N/2Ω) 2 0.5 at the cmb, compared to our maximum value of 1.1 calculated using (19). [START_REF] Buffett | Geomagnetic fluctuations reveal stable stratification at the top of the Earth's core[END_REF] or more in thickness in order to produce L-type lobes in the time average high latitude geomagnetic field, according to our results. For the octupole field component, Figure 8 predicts a value of G30 0.08 for S=0.04, somewhat higher than typically inferred from the time-averaged paleomagnetic field over the past two polarity chrons, but probably within the range of the combined dynamo model and observational uncertainties.

We conclude that a deep, thermally stratified layer beneath the CMB (nominally 1000 km thick) would be detectable with our methods but is at variance with the observed structure of the geomagnetic field. A thin thermally stratified layer beneath the CMB (a few hundred km thick or less) is more consistent with observations but would be difficult to confirm using our methods, because of uncertainties in long time average geomagnetic field reconstructions.

Detection of this amount of stratification may well be possible in the future, with better resolution of the time-averaged geomagnetic field combined with more extensive modeling of partially stratified dynamos.

Appendix Governing equations

Using the notation from Section 2, we define the codensity inthe outer core as the sum of a spatially uniform backgrond state (C o + C icb ) with a time-dependent part denoted by subscript o, plus a deviation from that state C, defined so that

C o + C = ρ o α(T o + T ) + β(χ o + χ) , (22)
where T and χ are the outer core variations in temperature and light element concentration, respectively. We take the background temperature to be adiabatic, with its mean value

E du * dt * -∇ 2 u * + 2ẑ × u * + ∇P * = EP r -1 Ra r r cmb C * + P m -1 (∇ × B * × B * ), (24)
∂B * ∂t * = ∇ × (u * × B * ) + P m -1 ∇ 2 B * , (25)
and

dC * dt * = P r -1 ∇ 2 C * + , (26) where
= Ra ad Ra -1 (27) in which Ra = βgD 5 χo ν 2 κ , (28)
as in Section 2, and

Ra ad = - αgD 5 Ṫad ν 2 κ (29)
is the Rayleigh number measuring the cooling of the core. The other dimensionless control parameters are Figure 9: a: Dimensionless thickness of the stratified layer as a function of its best fit in terms of Ra q and Ra as defined in the main text. The layer thickness is normalized by the core radius such that its maximum value is 0.6491. b: Dimensionless thickness of the stratified layer as a function of S = Ra q/Ra alone. Dotted lines shows the least square best fits. Circles, squares and crosses denote L0, L4 and L2 boundary conditions, respectively. The symbol color varies progressively with Ra from light grey for Ra = 10 7 to black for Ra = 10 8 . Simulations with S > 0.4 have been excluded for these fits since the stratified layer reaches its maximal value near S = 0.4 (where it invades the whole outer core).

E = ν ΩD 2 ; P r = ν κ ; P m = ν η . (30)
Figure 10: a: Dimensionless stability parameter as a function of its best fit in terms of Ra q and Ra as defined in the main text. b: Same parameter as a function of S = Ra q/Ra alone. Dotted lines shows the least square best fits. Circles, squares and crosses denote L0, L4 and L2 boundary conditions, respectively. Colors denote the different regimes for the time-averaged magnetic field structures: red for minima (regime M), green for patches (regime P), blue for lobes (regime L) and black for no dynamo (N). The symbol color varies progressively with Ra from light grey for Ra = 10 7 to black for Ra = 10 8 . Simulations with S ≥ 0.4, in which the layer invades the whole outer core have been excluded from the fits.

Figure 2

 2 Figure 2 shows E = 1 × 10 -4 cases in terms of the control parameters we varied. Axes are the Rayleigh number Ra and the stratification parameter S, defined as positive for stable boundary stratification and negative for unstable boundary stratification. The dashed line marks neutral (i.e., adiabatic) boundary flux cases. L0 (squares), L2 (crosses) and L4 (circles) denote spherical harmonic representation of the outer boundary heat flux heterogeneity as described in the previous section. The combined effects on the fluid motions of stable stratification and lateral boundary heterogeneity can be seen by comparing Figures 3 and 4, which show equatorial plane and global surface views of the time average flow structure at E = 1×10 -4 and Ra = 2×10 7 , with L0 (uniform) and L2 boundary heterogeneity, respectively, for two stratification parameters.With a uniform outer boundary, the time average codensity is spherically symmetric for both stabilizing (Figure3a, b, c) and destabilizing (Figure3d, e, f) stratification parameter. The only difference between the two cases is the presence of stable stratification in the S=0.2 case, extending downward from the outer boundary for a distance equal to about one quarter of the fluid shell depth. This stable stratification has major influences on the fluid velocities, as can be seen in the radial (u r ) and azimuthal (u φ ) images in Figure3. In both cases the radial velocities are highest near the inner boundary, where the buoyancy release is greatest, but in the S=0.2 case the radial velocities are truncated at the depth where stratification begins, and fail to penetrate to the outer boundary. Comparable differences are seen in the azimuthal velocity near the outer boundary, where in the S=-0.1 case the highest velocities occur at high latitudes, within and near the inner core tangent cylinder, whereas in the S=0.2 case the highest velocities occur in a retrograde (westward flowing) equatorial jet.

Figure 4

 4 Figure4has the same control parameters as Figure3but with L2 boundary heterogeneity added. In the destabilizing S=-0.1 case the stratification is highly destabilizing below the regions of higher than average boundary heat flux, and beneath these regions the time average radial velocity is strongly negative (downward flow). Conversely, below the regions of lower than average boundary heat flux the stratification is weakly stable and the time average radial velocity is positive (upward flow), from the inner boundary all the way to the outer. In contrast, in the S=0.2 case the thermal stratification is stabilizing at all longitudes immediately below the outer boundary and the radial velocity pattern is layered. Downwellings located below high boundary heat flux regions penetrate only about one quarter of the fluid depth before terminating, and the laterally broad upwellings below regions with low outer boundary heat flux show layering, with the strongest upwelling layers occurring immediately below the outer boundary and also just above the inner boundary. The azimuthal velocity patterns beneath the outer boundary for both S=-0.1 and S=0.2 cases include, at low latitudes, thin lenses of eastward (prograde) flow immediately beneath the high boundary heat flux regions but displaced slightly in the downstream direction of the heterogeneity, plus generally westward (retrograde) flow beneath the low boundary heat flux regions, again displaced slightly downstream of the center of the boundary heterogeneity. The other major difference between cases with versus without stabilizing boundary heat flux, which we later show is crucial for stratification detection, is the pattern and strength of the circulation inside the tangent cylinder of the inner boundary. The cases with destabilizing outer boundary fluxes in both Figures 3 and 4 have stronger azimuthal flows inside the tangent cylinder compared to their counterparts with stabilizing outer boundary fluxes.This difference indicates that the pattern of convection at depth inside the tangent cylinder is different with versus without boundary stabilization, and as we demonstrate next, this has observable effects on the magnetic fields generated in the tangent cylinder region.

Figure 5

 5 Figure5shows global views of the time average dynamo structure on and just below the outer boundary and in cross section, at E = 1 × 10 -4 and Ra = 6 × 10 7 , for stratification parameters S= 0 (neutral), S=0.2, and S=0.3, respectively, all with L2 boundary heterogeneity. The radial magnetic field on the outer boundary in the neutral case Figure5ahas higher intensity and shorter length scales in its structure, compared to its counterparts with stabilizing boundary heat flux, because without stratification, the short length scale fluctuations are more intense and require longer run times to average out. More importantly, the

  a regime diagram of the magnetic field in the polar region for the E = 1 × 10 -4 and L2 boundary heterogeneity cases. Axes are Rayleigh number Ra and outer boundary stratification parameter S, once again positive for stable stratification, negative for unstable. Other symbols denote the structure of the time average radial magnetic field structure on the outer (core-mantle) boundary at high latitudes. M denotes existence of polar intensity minima, P denotes high field intensity patches without polar minima, L denotes high field intensity lobes without polar minima, and N (no dynamo) denotes subcritical for dynamo action due to excessively strong boundary-induced stratification. Dual symbols in Figure6denote high latitude magnetic field structures that appear to be transitional, M/P for polar minima transitioning to patches, and P/M for the reverse. This figure demonstrates there are several (3 or more) distinctive magnetic field structures that occur systematically as the control parameters vary,

Figure 8b ,

 8b Figure 8b, G30 is positive for all S-values considered, although it might become negative for very strongly unstable stratification. More significantly, G30 increases with S, perhaps saturating around 0.15 near S=0.2.

  and the color scheme indicates the different magnetic field structures. Cases with δ * =0.6491 are saturated, meaning that the stratification spans the entire fluid layer from r cmb to r icb . Excluding saturated cases yields (a, b, c)=(1.26,1.2,-1.18) for the best-fitting constants in (16).

Figure

  Figure9areveals that the thickness of the stratified region increases strongly with Ra q and decreases about equally strongly with Ra, as anticipated on physical grounds. Significantly, the exponents in Figure9anearly satisfy the relationship c = -b, which implies that the thickness of the stratified layer is solely a function of the stratification parameter S. We explain the latter result by considering mass conservation in the stratified layer : assuming zero mass anomaly flux at the base of the stratified layer (i.e. zero codensity gradient), one obtains from mass balance that the volume of the layer evolves as -Raq/(Ra ) =

  shows the radial component of the geomagnetic field on the CMB at epoch 2010 from the CHAOS field model[START_REF] Olsen | The chaos-4 geomagnetic feld model[END_REF]. Panel a is the complete field model, with representation to spherical harmonic degree and order ( , m) max = 14. Intensity minima are evident at both poles at this resolution. Panels b and c show the same field model truncated at spherical harmonic degree and order ( , m) max = 12 and 6, respectively. The polar minima are strongly attenuated with increasingly severe truncation, such that the southern polar minima is completely removed at ( , m) max = 6 in panel c, replaced by two broad, high intensity lobes. Since the data coverage necessarily degrades in going from present-day to historical to archeomagnetic to paleomagnetic time scales, a plausible interpretation of the trends shown in Figures11 and 12is that the longer time span geomagnetic reconstructions lack the resolution necessary to resolve the true high latitude field structure, would lead to an over-estimation of the stratification parameter S, and therefore would overestimate the amount of stratification presently below the CMB.

  has proposed that a portion of the time-dependent zonal flow in the outer core as well as the time-dependent deviations from uniform decrease of the historical geomagnetic dipole moment are expressions of MAC oscillations -waves in the outer core governed by magnetic, buoyancy (i.e., Archimedes) and Coriolis forces that propagate in the north-south direction within a layer beneath the CMB, their propagation characteristics depending on the stratification in the layer. The best-fitting model for MAC wave propagation determined by Buffett et al. (2016) consists of an approximately 140 km thick layer with peak stratification immediately below the CMB given by (N/2Ω) 2 0.21 and a layer average gravitational stability of (N/2Ω) 2 0.1, approximately. Their MAC model stability is about 20% of our value and their stratified layer average value is about 13% of ours, based on (21). In terms of our scaling law (18) that relates the stratified layer thickness to the stratification parameter, the Buffett et al. (2016) δ 140 km corresponds to δ * 0.062. Inverting (18) with a=1.82 and b =1.2 yields S 0.042. Referring to Figure 6, this value of the stratification parameter is expected produce a high latitude, time average radial magnetic field structure of the M-type with polar minima, much like the present-day geomagnetic field. For reference, such a stratified layer would have to increase to 600 km

  denoted by T o = T ad . The rate-of-change of the background codensity (assumed constant in each calculation) is then Ċo = ρ o α Ṫad + β χo . (23) Scaling length, time, codensity, and magnetic field with D, D 2 /ν, D 2 ρβ χo /ν and ρ o Ω/σ as in Section 2, the dimensionless Boussinesq equations of motion are

  × 10 7 0.1 2 0.824921 0.815779 0.1254 0.0729 M 4.0 × 10 7 0.2 2 0.896273 0.864884 0.2326 0.1349 P 4.0 × 10 7 0.3 2 0.999349 0.913303 0.3866 0.2225 L 4.0 × 10 7 0.4

b

  [START_REF] Masters | On the resolution of density within the Earth[END_REF]); c Ancellini et al. (2013); d Poirier (1990); e Vocadlo et al. (2003); f Hirose et al. (2013); g Perriallt et al. (2010).

Figure 1 :

 1 Figure 1: Core-mantle boundary heat flux patterns L2 (a) and L4 (b) used in this study.

Figure 2 :

 2 Figure 2: Regime diagram showing numerical dynamos used in this study as a function of the primary control parameters Ra and S. Insert shows symbols that refer to the type of outer boundary heterogeneity.

Figure 3 :

 3 Figure 3: Equatorial plane and global Aitoff projection views of the time average flow structure at E = 1 × 10 -4 and Ra = 2 × 10 7 with L0 boundary heterogeneity, for boundary stratification parameters S=-0.1 (left column) and S=0.2 (right column). From top to bottom images show equatorial plane codensity (a,d), radial velocity (b,e) and azimuthal velocity at 0.95r cmb (c,f). Velocity scales are in dimensionless Reynolds number uD/η units, u being the appropriate dimensional velocity component. Longitude increases anti-clockwise starting from the right edge (3 pm) in each equatorial image and to the right of the centerline in the global projections.

Figure 4 :

 4 Figure 4: Equatorial plane and global Aitoff projections of the time average flow structure at E = 1 × 10 -4 and Ra = 2 × 10 7 with L2 boundary heterogeneity, for boundary stratification parameters S=-0.1 (left column) and S=0.2 (right column). From top to bottom images show equatorial plane codensity (a,d), radial velocity (b,e) and azimuthal velocity at 0.95r cmb (c,f). Velocity scales are in Reynolds number units, as in Figure 3. Longitude increases anticlockwise starting from the right edge (3 pm) in each equatorial image and to the right of the centerline in the global projections. Plus (+) and minus (-) signs indicate the longitudes of maximum and minimum boundary heat flux, respectively.

Figure 5 :

 5 Figure 5: Global Aitoff projections and zonal averages of the time average dynamo structure at E = 1 × 10 -4 and Ra = 6 × 10 7 with L2 boundary heterogeneity, for boundary stratification parameters S=0 (left column), S=0.2 (middle column) and S=0.3 (right column). Top row images (a,d,g) show radial magnetic field at r cmb . Middle row images (b,e,h) show radial velocity at 0.95r cmb . Bottom row images (c, f, i) show zonal averages of codensity (C), meridional streamlines over azimuthal velocity contours (U), and poloidal magnetic field lines over azimuthal field contours, , with (blue,red) and (dashed, solid) contours indicating (negative, positive) values, respectively. Magnetic field scale bars are in dimensionless Elsasser number σ/ρ o ΩB r units, where σ is electrical conductivity; likewise the velocity scale bars are in dimensionless Reynolds number uD/η units, u being the appropriate dimensional velocity component. Longitude increases to the right of the centerline in the global projections.

Figure 6 :

 6 Figure 6: Regime diagram showing E = 1 × 10 -4 , L2 dynamo results for the time average radial magnetic field structure on the outer (core-mantle) boundary at high latitudes.

Figure 7 :

 7 Figure 7: Polar views of time-averaged radial magnetic fields on the outer boundary from the dynamos in Figure 5. Solid black lines indicate bounding contours used for labeling field structures. a: Polar minimum M-structure with S=0 stratification parameter; b: Patch Pstructure with S=0.2; c: Lobe L-structure with S=0.3. Magnetic intensity scale bars are in dimensionless Elsasser number units; red crosses mark the time-averaged geomagnetic pole.

Figure 8

 8 Figure8: a: Axially symmetric (m=0) Gauss coefficient ratios versus spherical harmonic degree from time averaged dynamos with Ra = 6 × 10 7 and L2 boundary heterogeneity, for various stratifications. b: Gauss coefficient ratio G30 versus stratification parameter S for the same dynamos. M,P,L denote polar minima, patches, and lobes, respectively, in the high latitude magnetic field structure. Stable/unstable refers to the stratification below the outer boundary.

Figure 11 :

 11 Figure 11: Time average radial magnetic fields on the outer boundary from the dynamo case with Ra = 6 × 10 7 and L2 boundary heterogeneity and neutral (S=0) stratification parameter. a: truncation at spherical harmonic degree and order ( , m) max = 24; b: truncation at ( , m) max =12; c: truncation at ( , m) max =6.

Figure 12 :

 12 Figure 12: Radial magnetic fields on the CMB from the CHAOS4 2010 geomagnetic field model (Olsen et al., 2014) a: truncation at spherical harmonic degree and order ( , m) max = 14; b: truncation at ( , m) max =12; c: truncation at ( , m) max =6.

Table A2 :

 A2 Dynamo Cases Ra S L CMB C * Min. C * δ *

					N * 2	Structure
	1.0 × 10 7 0.1 4 0.704284 0.697211 0.1104 0.0640 M/P
	1.0 × 10 7 0.1 2 0.693749 0.68662 0.1033 0.0690 M
	1.0 × 10 7 -0.1 2 0.517744 0.517744 0	NA	M
	1.0 × 10 7 -0.1 4 0.521391 0.521391 0	NA	M
	1.0 × 10 7 0	2 0.600992 0.600992 0	NA	M
	1.0 × 10 7 0	4 0.604899 0.604899 0	NA	M
	2.0 × 10 7 0.1 2 0.764032 0.757908 0.0922 0.0663 M
	2.0 × 10 7 0.2 2 0.862358 0.832471 0.2209 0.1352 P
	2.0 × 10 7 0.2 0 0.859579 0.822195 0.2417 0.1546 -
	2.0 × 10 7 0.3 2 0.978527 0.893147 0.3755 0.2273 L
	2.0 × 10 7 0.4 2 1.23117	1	0.6491 0.3561 N
	2.0 × 10 7 0.5 2 1.51669	1	0.6491 0.7960 N
	2.0 × 10 7 1	4 2.93827	1	0.6491 2.986 N
	2.0 × 10 7

Table A2

 A2 Cases with E = 3 × 10 -5 ; all others with E = 1 × 10 -4 .

	Table A3: Core Properties
	: Continued Input Properties	Notation Value
	Ra ICB radius	S	L CMB C * Min. C * δ * r icb	N * 2 1220 km a Structure
	6.0 × 10 7 CMB radius 6.0 × 10 7 Density at core center 0 2 0.812963 0.812782 0.0571 0.0032 M r cmb 3480 km a 0 ρ c 12500 kg.m -3 a 0 0.815215 0.81507 0.0331 0.0044 M 6.0 × 10 7 Density at zero pressure ρ 0 7500 kg.m -3 0.1 2 0.852811 0.844548 0.1273 0.0648 M 6.0 × 10 7 Mean core density ρ o 11040 kg.m -3 a 0.1 0 0.856437 0.846691 0.1383 0.0704 M 6.0 × 10 7 Gravity at the CMB g 10.68 m.s -2 a 0.2 2 0.913068 0.882692 0.2495 0.1217 P 6.0 × 10 7 Compositional density jump at the ICB ∆ρ 500 kg.m -3 b 0.2 0 0.919266 0.884893 0.2657 0.1293 -6.0 × 10 7 0.3 2 1.01299 0.9276 Incompressibility at zero pressure K 0 4.75 × 10 11 Pa 0.3788 0.2254 L 6.0 × 10 7 0.3 0 1.02607 Melting temperature at the ICB T melt 5800 K c 0.93166 0.4041 0.2336 -6.0 × 10 7 0.4 2 1.23071 1 Entropy of melting ∆S 120 J.kg -1 .K -1 d 0.6491 0.3554 N 6.0 × 10 7 0.4 0 1.24977 1 Grüneisen parameter γ 1.5 e 0.6491 0.3848 N 6.0 × 10 7 0.5 2 1.51621 1 Specific heat C p 850 J.kg -1 .K -1 e 0.6491 0.7953 N 6.0 × 10 7 0.5 0 1.53475 1 Thermal expansion coefficient α 1.3 × 10 -5 K -1 e 0.6491 0.8238 N 6.0 × 10 7 -0.1 2 0.777253 0.777253 0 NA Compositional expansion coefficient β 1 M 6.0 × 10 7 -0.1 0 0.77527 0.77527 0 NA Thermal conductivity k 100 W.m -1 .K -1 M 7.0 × 10 7 0 Density length scale r ρ 6600 km a 2 0.822773 0.822546 0.0506 0.0045 M 7.0 × 10 7 Temperature length scale r T 6040 km c 0.1 2 0.865818 0.854998 0.1455 0.0743 M 7.0 × 10 7 Outer core light elements χ 9.8 wt.% f 0.2 2 0.922185 0.889587 0.2741 0.1189 P/M 7.0 × 10 7 0.3 2 1.0149 Mass of the core M c 1.95 × 10 24 kg 0.929291 0.3905 0.2192 L 8.0 × 10 7 Outer core kinematic viscosity ν (2,7,22)×10 -6 m 2 .s -1 g 0.1 2 0.873372 0.862076 0.1624 0.0695 M 8.0 × 10 7 Outer core diffusivity κ 10ν 0.2 2 0.927032 0.89635 0.2534 0.1211 M 8.0 × 10 7 0.3 2 1.01797 Mean CMB heat flux qtotal 79 mW.m -2 0.932767 0.3820 0.2230 L 9.0 × 10 7 0 2 0.84325 0.842954 0.0760 0.0039 M Adiabatic CMB heat flux q ad 96 mW.m -2
	9.0 × 10 7 Output Parameters 0.1 2 0.877763 0.868198 0.1539 0.0621 M Notation Value
	9.0 × 10 7 9.0 × 10 7 Light element concentration change rate χo 0.2 2 0.929826 0.897254 0.2534 0.1285 M 0.3 2 1.01873 3.2 × 10 -16 s -1 0.933266 0.3801 0.2248 L 1.0 × 10 8 Cooling rate Ṫo -3.7 × 10 -15 K.s -1 0.05 2 0.865793 0.861533 0.1117 0.0381 M 1.0 × 10 8 Codensity change rate Ċo 6.8 × 10 -16 kg s -1 0.1 2 0.886055 0.874903 0.1754 0.0636 M 1.0 × 10 8 0.15 2 0.90753 Codensity sink -0.8 0.888956 0.1903 0.0976 M 2.0 × 10 8 a 0.1 0 0.85105 Rayleigh number Ra (9 × 10 29 , 2 × 10 28 , 6 × 10 26 ) 0.84252 0.120 0.0710 M 2.0 × 10 8 a 0.1 2 0.85196 Stratification parameter S (0.018, 0.064, 0.2) 0.84456 0.114 0.0649 M 2.0 × 10 8 a 0.1 4 0.85411 Stratified layer thickness δ (52, 230, 920) km 0.84590 0.127 0.0646 M Stratified layer stability (CMB, ave.) N 2 (2.4, 1.7) × 10 -8 rad 2 .s -2
	a Dziewonski and Anderson (1981);

a

Acknowledgments

This research was supported by Frontiers in Earth System Dynamics grant EAR-1135382 from the National Science Foundation. The dynamo calculations were made at the Maryland Advanced Research Computer Center (MARCC). We benefited from insightful comments and suggestions by Hagay Amit, and from thorough reviews by Phil Livermore and Bruce Buffett.

Author contributions

PO and ML defined the numerical dynamo parameter sweep; ER did the dynamo calculations, collected and organized the numerical data; PO and ML did the scaling analysis, interpretation, and wrote the paper.