Outer Core Stratification From the High Latitude Structure of the Geomagnetic Field

The presence of stable stratiﬁcation has broad implications for the thermal and compositional state of the outer core, the evolution of Earth’s deep interior, and the energetics of the geodynamo. Yet the origin, strength, and depth extent of stratiﬁcation in the region below the core-mantle boundary remain open questions. Here we compare magnetic ﬁelds produced by numerical dynamos that include heterogeneous stable thermal stratiﬁcation below their outer boundary with models of the geomagnetic ﬁeld on the core-mantle boundary, focusing on high latitude structures. We demonstrate that the combination of high magnetic ﬁeld intensity regions and reversed magnetic ﬂux spots, especially at high latitudes, constrains outer core stratiﬁcation below the core-mantle boundary. In particular, we ﬁnd that the negative contribution to the axial dipole from reversed ﬂux spots is a strong inverse function of the stratiﬁcation. Comparison of our numerical dynamo results to the structure of the historical geomagnetic ﬁeld suggests up to 400 km of permeable, laterally heterogeneous thermal stratiﬁcation below the core-mantle boundary.

This raises multiple questions for the dynamics of the core. First, is the outer core stratification inferred by recent seismic studies compatible with the geomagnetic field and its secular variation? Core flow inversions based on the geomagnetic secular variation are best accommodated by including upwelling and downwelling motions extending very close to the core-mantle boundary (Gubbins, 2007;Amit, 2014;Lesur et al., 2015;Huguet et al., 2016). For example, Gubbins (2007) argued that the production of reversed flux spots on the coremantle boundary, which are rapidly evolving in the present-day geomagnetic field Olsen et al., 2014;Terra-Nova et al., 2015;Metman et al., 2018), limits the depth extent of the stratification to less than 150 km, assuming no radial motion in that layer and that the reversed flux spots on the coremantle boundary result from the expulsion of magnetic flux from the outer core.
Second, can numerical dynamos provide independent constraints on the strength and depth extent of the stratification? There are relatively few systematic investigations of the geodynamo in the presence of stratification (Sreenivasan and Gubbins, 2008;Nakagawa, 2011Nakagawa, , 2015Olson et al., 2017;Christensen, 2018). However, stratification effects have been extensively studied in the context of the solar dynamo (e.g., Browning et al., 2006Browning et al., , 2007Käpylä et al., 2008;Tobias et al., 2008;Brummell et al., 2010;Masada et al., 2013), Jupiter (Zhang and Schubert, 2000), Saturn (Christensen and Wicht, 2008;Stanley, 2010), and also Mercury (Christensen, 2006;Manglik et al., 2010). All these investigations found that the presence of a stratified layer affects the morphology of the magnetic field. In particular, a stratified layer below a convective region is key to generating a large-scale magnetic field in solar dynamo simulations (Browning et al., 2006(Browning et al., , 2007Käpylä et al., 2008), where strong zonal flows in the stratified layer stretch the poloidal magnetic field in the convective region into a large-scale toroidal magnetic field through an ω-effect. Other investigations have reported the generation of strong azimuthal flows within a stratified layer adjacent to a convective region (Zhang and Schubert, 2000;Takehiro and Lister, 2002;Couston et al., 2018), which attenuate high-frequency, nonaxisymmetric magnetic field components in the stratified layer (Christensen, 2006;Christensen and Wicht, 2008;Stanley, 2010).
Because stratification affects the magnetic field structure, dynamo simulations are useful in constraining the stratification in Earth's core. In a previous paper (Olson et al., 2017) we conducted a systematic investigation of the flow and the time average magnetic field in the presence of thermal stratification. We showed that the high latitude structures of the time average magnetic fields in numerical dynamos are sensitive to the strength and depth extent of thermal stratification below the dynamo upper boundary. This sensitivity offers the means to infer the properties of stratification below the core-mantle boundary (CMB) in terms of the time average structure of the geomagnetic field. In this paper we quantitatively compare the high latitude CMB structure of the COV-OBS geomagnetic field model (Gillet et al., 2013) to a suite of thermally stratified numerical dynamos. Extending the analysis in Olson et al. (2017), we compute the correlation of the high latitude structures of the time average magnetic field in the COV-OBS model and in our numerical dynamos. In addition, we analyze the time varying field, focusing on the effects of reversed flux spots on the axial dipole. These comparisons favor the existence of stratification below the CMB but also indicate that substantial radial motions are present there, implying that the stratification is rather weak and permeable to outer core convection.

NUMERICAL DYNAMOS WITH THERMAL STRATIFICATION BELOW THE OUTER BOUNDARY
The stratification analyzed in this study is due to thermal gradients that deviate from adiabatic (i.e., uniform entropy) conditions and are maintained by the heat flux imposed at the outer boundary. We include lateral variations of the boundary heat flux, following the results of mantle global circulation models Tackley, 2013, 2015;Zhong and Rudolph, 2015) that yield vigorous deep mantle convection with locally variable heat flux on the core-mantle boundary that is large enough in some places to sustain unstable thermal stratification , even if the thermal conductivity of the outer core is high (Ohta et al., 2016).
We model stratified thermochemical convection in the outer core with heterogeneous heat flux at the CMB using the formulation in Olson et al. (2017). Outer core density variations are expressed in terms of the codensity, i.e., density variations due to the combination of temperature and light element concentration variations: where ρ o is fluid mean density, T is temperature relative to the adiabat with mean T o , χ is the fluid light element concentration with mean χ o , and α and β are volumetric expansion coefficients for T and χ, respectively. In terms of these, the governing equations for thermochemical convection and dynamo action in a rotating spherical shell (with the Boussinesq approximation) include the following dimensionless control parameters: Here E is the Ekman number, Pr is the Prandtl number, Pm is the magnetic Prandtl number, and ǫ is the volumetric codensity source. In (2), denotes angular velocity of rotation, D = r o − r i is the depth of the fluid shell, r o and r i , the radii of the inner and outer fluid boundaries, with ν, η, and κ denoting kinematic viscosity, magnetic diffusivity, and codensity diffusivity, respectively.
At the inner boundary r i we assume no-slip velocity conditions and a uniform codensity C i . At the outer boundary we also assume no-slip velocity conditions, zero light element flux, and we specify the heat flux q to be the sum of a spherical mean part (denoted by an overbar) and a deviation from the spherical mean (denoted by a prime): where φ and θ are longitude and colatitude, respectively, andq is measured relative to the heat flux down the adiabat, withq > 0 being superadiabatic heat flux andq < 0 being subadiabatic heat flux. This formulation yields three additional dimensionless parameters that control the convection: a Rayleigh number based on the rate of increase of light element concentration in the fluid a second Rayleigh number based on the spherical mean heat flux at the outer boundary and a third Rayleigh number based on the peak-to-peak variation q ′ of the laterally varying boundary heat flux In (4-6), g is gravity at the outer boundary and k is thermal conductivity. In the numerical dynamos, the factors D 2 ρ o βχ o /ν and √ ρ o /σ (where σ is electrical conductivity) non-dimensionalize codensity variations and magnetic field intensity, respectively, and ν/D non-dimensionalizes the fluid velocity. In what follows, we retain these scalings for codensity and magnetic field, but we use η/D to scale the fluid velocity. With these factors, the scaling for velocity and magnetic field intensities are referred to as magnetic Reynolds number and Elsasser number units, respectively.
In Olson et al. (2017) we introduced a parameter describing the spherical mean stratification: defined to be positive whenq is negative, i.e., when the spherical mean boundary heat flux is stabilizing. There is also a related stratification parameter describing the effects of the boundary heat flux heterogeneity: We analyze dynamos with E = 10 −4 , Pr = 1, Pm = 6, and ǫ = −0.8, the latter appropriate for dominantly compositional convection but with some secular cooling included. The aspect ratio of the fluid shell is fixed at r i /r o = 0.351. The solid region r ≤ r i representing the inner core is assumed to have the same electrical conductivity as the fluid, and the solid region r ≥ r o representing the mantle is assumed to be electrically insulating. The boundary heat flux pattern is defined by a spherical mean part plus a heterogeneous part consisting mostly of spherical harmonic degree ℓ = 2 components at orders m = 0 and m = 2, adjusted so as to produce a pattern with nearly bilateral (i.e., 2-fold) azimuthal symmetry. The resulting boundary heat flux pattern is shown in Figure 1A and corresponds to the largest scale of lower mantle heterogeneity structure determined by Dziewonski et al. (2010). It is basically the same planform used by Olson and Amit (2015) in their study of the influences of lower mantle piles on magnetic polarity reversal behavior. Dynamo calculations are made at Ra = 6 × 10 7 and Ra = 9 × 10 7 for stratification parameters S ranging from −0.1, corresponding to superadiabatic CMB heat flux, to +0.3, corresponding to strongly subadiabatic CMB heat flux, using the MAGIC dynamo code (Wicht, 2002). We assume a constant value of S ′ = 0.58 for all cases. We find that by preserving S ′ , key properties of these dynamos such as the r.m.s. dipole axis tilt are nearly independent of S, while other properties such as the contribution of reversed flux to the axial dipole are relatively insensitive to Ra. We use a numerical grid with (n r , n θ , n φ ) = (81, 128, 256) in the fluid shell and spherical harmonic truncation (ℓ, m) max = 85. All the calculations were run for at least one magnetic diffusion time, in order that the run averages approximate true time averages.

STRATIFICATION DIAGNOSTICS
For comparison with the geomagnetic field, we focus on properties of the dynamo magnetic field structures, particularly at high latitudes. Previously, Olson et al. (2017) found that the high latitude dynamo magnetic fields are especially sensitive to stratification beneath the outer boundary, and the effects of stratification produce distinct and readily identifiable structures, both inside the tangent cylinder of the inner core and beyond, down to latitudes of approximately 45 o . In contrast, some dynamo magnetic field structures at low latitudes are not so easily related to stratification. Accordingly, most of our comparisons between numerical dynamos and the geomagnetic field are based on the variable B r cos(θ ), where θ is colatitude, which is the kernel of the axial dipole moment density on the CMB . For our applications, B r cos(θ ) is superior to the radial component of the magnetic field B r because the cos(θ ) factor adds weight to the high latitude field structure.
We characterize our numerical dynamos in terms of the structure of B r cos(θ ) on the outer boundary, along with the magnitude of the stratification and the upwelling below the outer boundary. One important diagnostic is the ratio of reversed to normal polarity flux on the outer boundary. The two individual contributors to this ratio are given by where the superscripts N and R denote the signs of B r (positive or negative) that define the dominant and the subordinate components, respectively, of the axial dipole moment, and A is the outer boundary surface area. The ratio of these two fluxes is defined so that 0 ≤ F * ≤ 1, the lower limit indicating zero contribution to the axial dipole moment from reversed flux, the upper limit corresponding to a vanishingly small axial dipole. Another magnetic diagnostic we use is the distribution of high latitude, high intensity B r cos(θ )-structures. Our previous study (Olson et al., 2017) documented that the morphology of high latitude, high intensity B r -structures in time average dynamo magnetic fields can be used to constrain the stratification parameter S. In the next sections we demonstrate that B r cos(θ ) is even more sensitive to S, both in snapshots and in time averages. We measure the stratification in our numerical dynamos using the spherical mean thickness of the stratified region and its gravitational stability. The dimensionless spherical mean thickness of the stratified region of the dynamo is defined as where r min is the radius where the dimensionless spherical mean codensityC * reaches its local minimum value below the outer boundary. Likewise, we define the gravitational stability of the stratified layer in terms of the dimensionless buoyancy frequency squared: where δC * =C * o −C * min is the dimensionless codensity increase across the stratified region. In Olson et al. (2017) we derived the following scaling laws for these quantities: in which (a δ , b δ ) = (1.82, 1.2), plus in which (a N , b N ) = (0.72, 1). Lastly, the r.m.s. upwelling strength below the outer boundary is used to characterize the effects of stratification on the flow. We define the dimensionless outer boundary (or CMB) upwelling strength as where ∇ H and u are the dimensionless horizontal divergence and the fluid velocity, respectively, and || denotes r.m.s. average over the spherical shell at 0.95r o . Figure 1 shows the pattern of heat flux applied to the outer dynamo boundary and the resulting radial velocity pattern from two dynamos with Ra = 9 × 10 7 but different amounts of stratification. Figure 1A shows the nearly bilaterally symmetric boundary heat flux in dimensionless form, with a great circle of elevated heat flux that includes both polar regions, separating two large, low latitude regions with reduced heat flux. The reduced heat flux regions correspond approximately to the Large Low Shear Velocity Provinces (LLSVPs) imaged by seismic tomography (Garnero and McNamara, 2008)  In the snapshots, the effects of stratification are most evident in the difference in magnitude of the radial velocities. In the S = 0 dynamo, dimensionless radial velocities in Figure 1B exceed 400 in places, with an r.m.s. at this depth of approximately 180. In the S = 0.3 dynamo, in contrast, dimensionless radial velocities in Figure 1D nowhere exceed 30, and the r.m.s. at this depth is approximately 9. Clearly, the stabilizing effects of the boundary heat flux suppress the radial velocity below the outer boundary, reducing the r.m.s. strength of upwellings and downwellings there by a nearly a factor of 20 between the two cases.

DYNAMOS WITH STRATIFICATION
The strong reduction in radial velocity caused by stratification that is seen in the snapshots is less extreme in the time averages in Figures 1C,E. Overall, the patterns of radial velocity are more similar in these time averages compared to their corresponding snapshots, because the boundary heat flux heterogeneity plays a relatively greater role in structuring the time average velocities. The greatest differences between the two dynamos in terms of their time average radial velocities are found at high latitudes. In the S = 0 dynamo, there are strong polar upwellings and strong downwellings along the inner core tangent cylinder in both hemispheres (Figure 1C), structures that are missing from the strongly stratified S = 0.3 dynamo (Figure 1E). Figure 2 shows snapshots and time averages of the radial magnetic field intensity on the outer boundary at Ra = 9 × 10 7 for stratification parameter S varying between 0 and 0.3. Unlike the radial velocity, for which the amplitude of the upwellings and downwellings show the strongest influence of stratification, the magnetic field on the outer boundary mainly responds to the stratification through changes in its structure, rather than its amplitude. For example, in the snapshot field structures in Figure 2 there is a progressive reduction in the number and the intensity of reversed flux spots with increasing S, such that the S = 0.3 dynamo snapshot ( Figure 2G) is entirely lacking in reversed flux at high latitudes in both hemispheres, yet the overall magnetic field intensity barely changes with S. The other expression of structural change at high latitudes is seen in the time average field structures. In the S = 0 and S = 0.1 dynamos (Figures 2B,D) the high latitude structure consists of rings of high intensity field located near the tangent cylinder surrounding deep intensity minima, with localized reversed flux at the poles. In the S = 0.2 dynamo ( Figure 2F) the polar minima are gone and the high intensity field is localized in patches, two in each hemisphere. Lastly, in the strongly stratified S = 0.3 case ( Figure 2H) the two patches in each hemisphere have merged into a single high intensity lobe, positioned such that there is a field intensity maximum located at each pole.
The trends in the time average magnetic field structure in Figure 2 can be explained in terms of the changes in the internal dynamo structure with increasing stratification. Figure 3 compares the azimuthally averaged structure of Ra = 9 × 10 7 dynamos with S = 0 and S = 0.3, respectively. The internal structure of the S = 0 dynamo (Figures 3A-C) includes an adverse (i.e., destabilizing) codensity gradient, strong thermal wind circulations with meridional overturning inside the tangent cylinder in both hemispheres, and low magnetic field intensity near the outer boundary inside the tangent cylinder, locally reversed at each pole. The polar reversed flux, the low field intensity inside the tangent cylinder, and the high intensity field along the tangent cylinder, can be explained in this dynamo in terms of incomplete flux expulsion by the meridional circulations inside each tangent cylinder region. This circulation advects the poloidal magnetic field away from the poles and concentrates it along the tangent cylinder, producing the high latitude pattern seen in Figure 2B. In contrast, the azimuthally averaged internal structure of the S = 0.3 dynamo (Figures 3D-F) includes stable stratification below the outer boundary at all latitudes, a two-layer meridional circulation pattern at low and middle latitudes, and reversed circulations inside the tangent cylinder that exchange fluid with the meridional circulations outside. The meridional circulations inside the tangent cylinder region include polar downwellings that produce horizontal convergence beneath the outer boundary. These circulations concentrate poloidal magnetic flux close to the pole, producing polar intensity maxima in both hemispheres, as seen in Figure 2H. Figure 4 shows B r cos(θ ) on the core-mantle boundary from the COV-OBS geomagnetic field model (Gillet et al., 2013; http://www.spacecenter.dk/files/magnetic-models/COV-OBS/ Northern and Southern hemispheres, respectively. Data sources for this geomagnetic field model include space-borne magnetometer measurements during low altitude satellite orbits plus annual means from ground-based observatories. The COV-OBS core field is represented at epochs spaced 2 years apart, and is complete to spherical harmonic degree and order 14. We treat Figures 4A,B as snapshots of the present-day core field, for comparison with our dynamo snapshots. The maps of B r cos(θ ) in Figures 4C,D are averages over 88 epochs, but their 174 year time span is far shorter than the averaging times in our dynamos, which are of the order of a few hundred thousand years. Nevertheless, in what follows we treat the geomagnetic field average as a true time average for purposes of comparison with the dynamo averages. Figure 5 shows snapshots and time averages of B r cos(θ ) from numerical dynamos at Ra = 6 × 10 7 for stratification parameters ranging from S = −0.1, corresponding to a superadiabatic thermal gradient at the CMB, to S = 0.3, corresponding to a strongly subadiabatic thermal gradient at the CMB. The top row of maps in Figure 5 are Northern hemisphere B r cos(θ ) snapshots, the middle row are Southern hemisphere snapshots at the same times, and the bottom row are Northern hemisphere time averages. Southern hemisphere time averages differ insignificantly from their northern counterparts and are not shown.

COMPARISONS WITH THE GEOMAGNETIC FIELD AT THE CMB
The top and middle rows in Figure 5 show the same qualitative trends as in Figure 2 in terms of the disappearance of reversed flux with increasing stratification parameter. To demonstrate this quantitatively, Figure 6 shows F * , the ratio of reversed to normal flux defined by Equations (9) and (10) vs. stratification parameter S, for the Ra = 6 × 10 7 dynamos in Figure 5 and the Ra = 9 × 10 7 dynamos in Figure 2. The error bars indicate the standard deviation of F * based on six to eight snapshots from each dynamo. Although there is some dependence on the Rayleigh number at S = 0 and S = −0.1, the reversed to normal flux ratios at both Rayleigh numbers decrease strongly with increasing S, rapidly converging toward zero at larger S. Reversed flux patches are generally non-axisymmetric structures. Therefore, this decrease in F * with increasing S agrees with previous studies that found that stratification removes not only reversed flux (Sreenivasan and Gubbins, 2008), but also other non-axisymmetric components of the magnetic field (Christensen, 2006;Christensen and Wicht, 2008;Stanley, 2010). We also show in Figure 6 the reversed to normal flux ratio on the CMB from the Gillet et al. (2013) COV-OBS geomagnetic field model at epoch 2014. Dynamos with S = 0.1 best match the present-day geomagnetic field structure in terms of the relative contribution of reversed flux to the axial dipole.
There are several important differences between the numerical dynamos and the core field model that need to be factored out in order to make the comparison in Figure 6 more direct. First, the core field model is truncated at spherical harmonic degree 14, whereas the numerical dynamos used for F * in Figure 6 represent the field to spherical harmonic degree 85. Second, ambiguities arise in the calculation of F * that depend on the choice of the geographic equator vs. the magnetic equator. All of the values of F * in Figure 6 are based on the geographic equator, whereas the standard methods for calculating reversed flux on the CMB make use of the magnetic equator (Terra-Nova et al., 2015; Metman  , 2018). The most obvious consequence of the choice of equator is the contribution to reversed flux from the tilt of the dipole axis. Dipole axis tilt contributes to the inventory of reversed flux when using the geographic equator, but it need not when using the magnetic equator. Third, the value of F * changes with time in the core field model, being generally smaller in the past, whereas the averaging of widely spaced snapshots removes most (or all) of the secular drift in F * from the numerical dynamos.
For these reasons, we show in Figure 7 comparisons between numerical dynamos, the COV-OBS core field model, and two other core field models, based on a modified reversed to normal flux ratio, F * C . For the core field model COV-OBS, F * C is just F * with the equatorial dipole terms removed. Removing the equatorial dipole represents the lowest order correction to the magnetic equator. F * C from COV-OBS is shown at epochs 2014 and 1964, to illustrate the magnitude of the drift in this parameter with time. MLM in Figure 7 corresponds to the mean value of F * C calculated by Metman et al. (2018) for epoch 2015 using their definition of magnetic equator on core field model COV-OBS.x1 (Gillet et al., 2015). TN in Figure 7 corresponds to the value of F * C calculated by Terra-Nova et al. (2015) using their definition of magnetic equator on the present-day (zero age) limit of archeomagnetic field model CALSk.4b (Korte and Constable, 2011). For the numerical dynamos, F * C in Figure 7 is F * with the equatorial dipole terms removed and with a crustal filter applied, such that the magnetic field amplitude decreases by a factor of e with each spherical harmonic degree above 14. We note that the effects of removing the equatorial dipole from the numerical dynamos and the modern core field models are comparable, because the r.m.s. dipole axis tilt of the numerical dynamos (10 degrees at Ra = 6 × 10 7 and 12 degrees at Ra = 6 × 10 7 ) are comparable to the time average dipole axis tilt in the historical geomagnetic field. Finally, we calculate F * C for the dynamos and for field model COV-OBS using the same 1.5 x 1.5 degree grid.
The effects of crustal filtering and correction to the magnetic equator are to reduce F * C relative to F * , for the core field models as well as the numerical dynamos. Yet the same trends evident in Figure 6 are seen in Figure 7, with perhaps greater clarity. The numerical dynamos with S = 0.1 are compatible with all three core field models, in spite of the differences in processing that went into calculating reversed and normal flux in each case. There is some suggestion in Figure 7 that neutrally stratified dynamos with S = 0 may also be compatible, although this comparison is less convincing. And, just like Figure 6, this comparison argues against the more strongly stratified dynamos Frontiers in Earth Science | www.frontiersin.org 8 October 2018 | Volume 6 | Article 140 with S = 0.2 and greater. In short, Figures 6, 7 imply that strong thermal stratification below the CMB, characterized by S ≥ 0.2, as well as strongly superadiabatic conditions below the CMB characterized by S ≤ −0.1, are incompatible with the present-day structure of the geomagnetic field insofar as the amount of reversed flux is concerned, whereas on this same basis, the present-day field is compatible with weak stratification characterized by S = 0.1 or perhaps a bit less. Disappearance of reversed flux with increasingly strong stratification is a direct consequence of the reduction in strength FIGURE 6 | Ratios of reversed to normal flux contributions to the axial dipole F * vs. stratification parameter S from numerical dynamos (symbols), compared to F * from the COV-OBS geomagnetic field model on the core-mantle boundary at epoch 2014 (dashed line). Error bars denote one standard deviation of dynamo snapshots.
FIGURE 7 | Modified ratios F * C of reversed to normal flux contributions to the axial dipole vs. stratification parameter S from numerical dynamos (symbols), compared to F * C from geomagnetic field models on the core-mantle boundary. Geomagnetic field F * C values labeled COV-OBS, MLM, and TN are explained in the text. Symbol error bars denote one standard deviation of dynamo snapshots.
of the radial velocity below the outer boundary. Figure 8 shows F * C vs. the CMB upwelling strength W * defined by (15). The CMB upwelling is given in dimensionless form, in units of η/D 2 . The color and symbol schemes in Figure 8 are the same as in Figure 6, and only the snapshot averaged values of W * are plotted because the variation between snapshots is no larger than the symbols. Figure 8 shows a strong, positive and approximately linear correlation between the dynamo reversed to normal flux ratio and CMB upwelling. CMB upwelling less than a few hundred hardly produces any reversed flux, whereas for CMB upwelling above W * ≃ 1700, reversed flux reduces the axial dipole by 10% or more. Figure 8 also shows the range in F * C from the core field models in Figure 7. The best matching S = 0.1 and S = 0 dynamos intersect the dynamo trend at dimensionless FIGURE 8 | Modified ratios F * C of reversed to normal flux contributions to the axial dipole vs. dimensionless r.m.s. upwelling W * below the core-mantle boundary from numerical dynamos (symbols), compared to the geomagnetic field models on the core-mantle boundary described in the text. CMB upwelling strengths of W * = 800-1500, with W * ≃ 1000 being a representative value. In addition to reversed flux in snapshots, the polar structure of the time average geomagnetic field is also sensitive to core stratification. Based on visual comparison of the time averages in Figures 4, 5, the S = 0.1 dynamo best replicates the polar field structure of COV-OBS geomagnetic field model. The high latitude structure of that dynamo in Figure 5 includes two partially isolated high field intensity patches enclosing a polar intensity minimum, much like the high latitude geomagnetic field structures in Figure 4. In contrast, the dynamos with S ≤ 0 in Figure 5 have ring-shaped high intensity field regions, while the dynamos with S ≥ 0.2 lack polar intensity minima or in the extreme case, have polar intensity maxima.
A quantitative test of this visual interpretation can be made using the cross correlation between a time average dynamo magnetic field and a 174 year geomagnetic field average. Figure 9 shows global cross correlations of time average B r cos(θ ) between the Ra = 6 × 10 7 dynamos and the COV-OBS geomagnetic field model vs. longitude shift in degrees, with positive and negative denoting westward and eastward shifts, respectively, of the dynamo relative to the geomagnetic field model. It is helpful to include longitude pattern shifts in this analysis, since the longitudes of the high field intensity patches vary with the dynamo control parameters. Allowance for some longitude pattern shift mitigates the bias from this variation. The spectra of the time average dynamo fields on the outer boundary contain little power above spherical harmonic degree 14, so crustal filtering is not necessary here. The cross correlations were preconditioned for weak field suppression by masking boundary regions with field intensity below 20% of the maximum intensity, in order to add weight to the high field intensity regions. Figure 9 indicates there is some dependence of the correlation on longitude shift, but for shifts of 20 o or less the effect is relatively minor. More significantly, there is a substantial difference in this correlation between unstratified and weakly stratified dynamos vs. the strongly stratified dynamos, with the former group correlating above 0.5 and the latter group below 0.5. Interestingly, the best correlation is found for the unstratified S = 0 dynamo and the second best is the S = −0.1 dynamo, although their correlations differ very little from the S = 0.1 dynamo overall.

IMPLICATIONS FOR OUTER CORE STRATIFICATION
Our comparisons between numerical dynamos and the geomagnetic field on the CMB favor the existence of outer core stratification with stratification parameter S close to 0.1. Equally significant, these same comparisons argue against stronger outer core stratification, as would be characterized by S ≥ 0.2, say. Although our study does not consider situations in which the stabilizing effects of stratification vastly outweigh the destabilizing effects of inner core growth, as would be the case for strong compositional stratification (Landeau et al., 2016;Nakagawa, 2017;Christensen, 2018), the fact that we can exclude thermally stratified dynamos with large S suggests our results might also be applicable for constraining outer core compositional stratification.
Assuming that S = 0.1 in the region below the CMB, our previously-derived dynamo scaling laws yield estimates of the thickness of the stratified layer and its gravitational stability. In dimensional terms, our scaling laws for stratified layer thickness (13) and squared buoyancy frequency (14) are, for the outer core plus and from the definition (7) of S, the subadiabatic heat flux on the CMB is Using the core property values in Table 1 with S = 0.1, (16) gives δ ≃ 400 km, (18) gives q ad −q cmb ≃ 17 mW.m −2 , and (17) gives N 2 ≃ 1.7×10 −8 rad 2 .s −2 .
A 400 km layer may seem excessively thick for a thermal stratification, but it is important to note that this value refers to the full spherical mean thickness of the layer, from the CMB to the depth where the spherically averaged codensity profile has a local minimum. Furthermore, although our results favor S = 0.1  Perrillat et al. (2010).
Frontiers in Earth Science | www.frontiersin.org stratification, they are also marginally consistent with somewhat weaker stratification, S = 0.05 for example. In that case, the stratified layer would be substantially thinner, with δ ≃ 170 km. In dynamical terms, such a layer would not prevent upward radial motions reaching close to the CMB, as evidenced by our finding that the r.m.s. CMB upwelling strength W * ≃ 1000. For the geodynamo, in dimensional units, W ≃ 1000η/(r cmb − r icb ) 2 , where the subscripts cmb and icb denote outer and inner core radii, respectively. In terms of the values of core properties in Table 1, this corresponds to W ≃ 0.5/century for the r.m.s upwelling below the CMB, within the range of the estimates of the r.m.s. CMB upwelling obtained from frozen flux inversions of the geomagnetic secular variation, which vary between 0.1/century and 4/century r.m.s. Amit and Pais, 2013). Even with S=0.1 stratification, superadiabatic thermal conditions may be present beneath approximately 5% of the CMB, according to the boundary heat flux pattern in Figure 1. If so, thermal instabilities originating at the CMB can penetrate the layer in these regions, making the thermal stratification somewhat permeable to outer core convection and allowing the formation of reversed flux spots as observed in the geomagnetic core field.
Permeable stratification distributed over several hundred kilometers beneath the CMB is consistent with other fluid dynamical effects, in particular, the upward penetration of convection through a weakly stratified layer (Takehiro and Lister, 2002;Rogers and Glatzmaier, 2005). Two scalings for the penetration distance have been proposed; because it is unclear which applies best to the core, we consider both. The first, by Takehiro and Lister (2002), predicts that convection penetrates a distance given by δ p ∼ 2 λ/N, where λ is the horizontal flow length scale. Using (17), we estimate 2 /N ∼ 1 for thermal stratification in the Earth's core. This implies a weak stratification, where the effects of stable stratification below the CMB are only about as strong as Coriolis effects from rotation. With this stratification, the Takehiro and Lister (2002) scaling predicts that convective eddies wider than about 400 km will penetrate to the CMB. The second scaling is derived from numerical models of solar convection (Hurlburt et al., 1994;Rogers and Glatzmaier, 2005). These studies find that the penetration distance scales with the ratio of the unstable to the stable stratification, i.e., δ p ∼ DS −1 in our notation. This scaling also predicts that convective motions easily penetrate a 400 km layer with S = 0.1. The Rayleigh number Ra, the Ekman number E and the magnetic Prandtl number Pm in our numerical dynamos are orders of magnitude away from Earth's core values. This raises a standard question for dynamo modelers: How sensitive are our conclusions to our parameter choices? Assuming reversed flux spots originate from toroidal flux expulsion (Gubbins, 2007), we expect the flux ratio at the CMB (either F * or F * C ) to scale as the flux ratio measured in the underlying convective region modulated by the radial velocity in the stratified region relative to that in the convective region. For dipole-dominated dynamos, the relative strength of the dipole varies only marginally with Ra, E, and Pm (Aubert et al., 2009). We hypothesize that the flux ratios in the convective region are only weakly sensitive to these parameters. In addition, the radial velocity in the stratified region relative to that in the convective region depends only on the ratio of the stratified layer thickness δ to the penetration distance of the convection δ p . Using the scalings discussed above for δ p and relation (16), we infer that δ/δ p , and therefore F * and F * C depend only on S and possibly /N. And, in contrast to Ra, E, and Pm, the values of /N and S in our dynamos are in the range expected for thermal stratification at the top of Earth's core (Takehiro and Lister, 2002;Buffett et al., 2016). Provided these expectations are met, our conclusions about stratification are applicable to the core. This can be tested by extending our analysis to stratified dynamos with more realistic values of Ra, E, and Pm.

AUTHOR CONTRIBUTIONS
PO and ML designed the study, did the analysis, and wrote the paper. ER managed the dynamo calculations and produced the dynamo data products.