Approaching Earth’s core conditions in high-resolution geodynamo simulations

SUMMARY

The geodynamo features a broad separation between the large scale at which Earth’s magnetic field is sustained against ohmic dissipation and the small scales of the turbulent and electrically conducting underlying fluid flow in the outer core. Here, the properties of this scale separation are analysed using high-resolution numerical simulations that approach closer to Earth’s core conditions than earlier models. The new simulations are obtained by increasing the resolution and gradually relaxing the hyperdiffusive approximation of previously published low-resolution cases. This upsizing process does not perturb the previously obtained large-scale, leading-order quasi-geostrophic (QG) and first-order magneto-Archimedes-Coriolis (MAC) force balances. As a result, upsizing causes only weak transients typically lasting a fraction of a convective overturn time, thereby demonstrating the efficiency of this approach to reach extreme conditions at reduced computational cost. As Earth’s core conditions are approached in the upsized simulations, Ohmic losses dissipate up to 97 per cent of the injected convective power. Kinetic energy spectra feature a gradually broadening self-similar, power-law spectral range extending over more than a decade in length scale. In this range, the spectral energy density profile of vorticity is shown to be approximately flat between the large scale at which the magnetic field draws its energy from convection through the QG-MAC force balance and the small scale at which this energy is dissipated. The resulting velocity and density anomaly planforms in the physical space consist in large-scale columnar sheets and plumes, respectively, co-existing with small-scale vorticity filaments and density anomaly ramifications. In contrast, magnetic field planforms keep their large-scale structure after upsizing. The small-scale vorticity filaments are aligned with the large-scale magnetic field lines, thereby minimizing the dynamical influence of the Lorentz force. The diagnostic outputs of the upsized simulations are more consistent with the asymptotic QG-MAC theory than those of the low-resolution cases that they originate from, but still feature small residual deviations that may call for further theoretical refinements to account for the structuring constraints of the magnetic field on the flow.

Core, Dynamo: theories and simulations, Numerical modelling

1 INTRODUCTION

Steadily improving numerical models of the geodynamo have appeared since the historical outset of the discipline, more than two decades ago. In retrospect, it is quite striking that a significant part of the fundamental insight on how Earth’s magnetic field may be regenerated by the underlying core convection has been gained through early models that featured a low spatial resolution. A self-sustained magnetic field with dipole reversals (Glatzmaier & Roberts 1995) and basic but plausible induction mechanisms (Olson et al. 1999) could indeed be captured with simulations resolving spatial scales down to only a thousand kilometers at Earth’s core–mantle boundary (corresponding to spherical harmonic degree 20). At first glance this may seem sufficient since geomagnetic field models using even the most recent satellite observations can resolve the main field of internal origin only up to degree 13 and its rate of change up to degree about 16, after which contributions from other sources and contamination from noise become too strong (Hulot et al. 2015). This is of course misleading because the Earth’s core is expected to be in a strongly turbulent state (e.g. Aubert et al. 2017), implying that the large-scale picture that we observe should be influenced by nonlinear interactions between much smaller, unobservable scales. This theoretical expectation has been partly confirmed by subsequent modelling efforts, where advances of the computer power allowed for calculations with a gradually increasing spatial resolution (harmonic degree in excess of 100), and more realistic physical parameters. The output from these models has indeed reached a higher level of quantitative agreement with the detailed morphology of the geomagnetic field (Christensen et al. 2010; Mound et al. 2015), and of its temporal variations, notably the westwards drift of low-latitude magnetic flux patches at the core–mantle boundary (Aubert et al. 2013; Schaeffer et al. 2017). However, this new generation of models still featured a relatively low spatial resolution and parameters quite far from Earth’s core conditions, raising the question of whether these successes were indeed obtained in physically relevant conditions, especially regarding the force balance obtained in the numerical models (e.g. Soderlund et al. 2012).

In an attempt to answer this question, it has recently been shown that the physical conditions at which these models operate and those of Earth’s core can be connected through the formulation of a unidimensional theoretical path in model parameter space (Aubert et al. 2017). The level of magnetic turbulence (as characterized by the value of the magnetic Reynolds number, see Section 2) is already realistic at path start and therefore remains constant along this path. However, a decreasing ratio between the kinematic and magnetic diffusivities along the path (decreasing magnetic Prandtl number) and increasing levels of forcing imply that hydrodynamic turbulence becomes stronger (increasing Reynolds number), and separation between the largest magnetic field scales and the smallest velocity field scales increases accordingly. Despite this increasing level of turbulence, a striking result of calculations carried out along this path is that at large-scales (below spherical harmonic degree 30), and timescales longer than the convective overturn time, the morphology of the magnetic field, of its rate-of-change and underlying core flows, as well as the leading-order force balance remain approximately invariant as we progress towards Earth’s core conditions (Aubert et al. 2017; Aubert 2018). This result rationalizes the early successes of geodynamo modelling that were obtained at low spatial resolution, but raises further questions regarding the role and importance of the broad scale separation that is expected in Earth’s core, and indeed observed in recent landmark simulations at extreme conditions and high resolution (harmonic degrees in excess of 200 and up to 1000; Sakuraba & Roberts 2009; Yadav et al. 2016; Schaeffer et al. 2017; Sheyko et al. 2018). One key point is that only a small fraction (the first 15 per cent) of this parameter space path has been covered by fully resolved direct numerical simulations (DNS), and a larger portion (up to 50 per cent) has been explored using lower resolution, large-eddy simulations (LES) employing a form of eddy viscosity that prescribes the smallest possible scales in the system. For these simulations, hyperdiffusion has been preferred over sophisticated subgrid parametrizations (e.g. Baerenzung et al. 2008, 2010; Matsui & Buffett 2013), on the basis of physical scaling arguments (Davidson 2013) that predict that the dominant force balance in the system is achieved at a relatively large scale (spherical harmonic degree about 10) that does not significantly evolve along the path. The realism of LES simulations has been successfully tested in the path region where both DNS and LES are feasible (Aubert et al. 2017). However, this region corresponds only to a moderate level of hydrodynamic turbulence, and one may question the realism of LES as regards the details of the scale separation and the associated nonlinear energy transfers when turbulence significantly increases. This calls for high-resolution geodynamo simulations at advanced positions along this path.

High-resolution geodynamo simulations in the rapidly rotating and strongly turbulent regime pertaining to Earth’s core conditions are also desirable to advance the understanding of how rotation and a self-generated magnetic field constrain turbulence in such systems. Recent cross-fertilization between theoretical investigations (Davidson 2013; Calkins et al. 2015; Aurnou & King 2017; Calkins 2018) and numerical efforts (Aubert et al. 2017; Schaeffer et al. 2017) have led to an emerging consensus that in this regime, the system is governed at large scales by a leading-order geostrophic balance between the pressure and Coriolis forces, as opposed to the earlier conjecture of a large-scale magnetostrophic balance where the Lorentz force would also reach a leading order (e.g. Hollerbach 1996). With a leading-order geostrophic balance being enforced, it has been argued that the planforms of the dynamo are in fact broadly similar to those of non-magnetic, rotating convection (Soderlund et al. 2012). This apparently weak influence of the magnetic field was however obtained in contexts where the magnetic energy was equivalent to the kinetic energy, implying a weak magnetic control on the flow. Subsequent simulations at more extreme parameters and higher magnetic to kinetic energy ratios (stronger magnetic control) have indeed highlighted different distributions of length scales in magnetic and non-magnetic rotating turbulence (Yadav et al. 2016; Sheyko et al. 2018). Still, the Lorentz force tends to be minimized as a consequence of Lenz law, leading, for instance, to the enforcement of the Taylor constraint (Taylor 1963) at the axisymmetric level (see e.g. Aubert et al. 2017). Though there exists an important body of numerical work on rotating or magnetic turbulence in Cartesian domains (see reviews in Tobias et al. 2012; Nataf & Schaeffer 2015), the configuration of forced magnetohydrodynamic (MHD) and rotating turbulence in a spherical shell has been less studied in this so-called strong regime of interaction between the velocity and magnetic fields.

In this study, the hyperdiffusive approximation of LES simulations previously published in Aubert et al. (2017) and Aubert (2018) is gradually relaxed, and their spatial resolution is increased accordingly, with all other parameters kept the same, a process that will be referred to as upsizing. Due to computational limitations, high spatial resolution and long integration times cannot be achieved simultaneously. The upsized models are therefore computed only for a short amount of time and our first goal is to assess the duration of transients caused by upsizing and the relevance of this numerical approach. Our next goals are to evaluate the successes and shortcomings of the previously employed hyperdiffusive approximations in the light of the upsized solutions and of other available extreme DNS simulations, and propose optimal approximations for forthcoming work. The upsized simulations are also useful to evaluate the robustness of the leading-order force balance, extract the next-order force balances established at smaller scales, characterize the turbulence and scale separation through their statistical properties and assess their level in Earth’s core. The manuscript is organized as follows: Section 2 presents the numerical dynamo model and methods. Results are presented in Section 3 and discussed in Section 4.

2 MODEL AND METHODS

2.1 Model set-up, dimensionless inputs and outputs

We implement the equations of Boussinesq convection, thermochemical density anomaly transport and magnetic induction in the MHD approximation within an electrically conducting and rotating spherical fluid shell of thickness D = r_o − r_i and aspect ratio r_i/r_o = 0.35 representing the Earth’s outer core. These equations can be found in Aubert et al. (2017) (from hereafter A17). We solve for the velocity field |$\boldsymbol{\mathrm{u}}$|⁠, magnetic field |$\boldsymbol{\mathrm{B}}$| and density anomaly field C. The shell is surrounded by solid inner core of radius r_i, and a solid mantle between radii r_o and 1.83r_o, both of which are electrically conducting. These two surrounding layers are electromagnetically coupled to the outer core and gravitationally coupled between each other, and present a variable axial differential rotation with respect to the fluid shell. The moments of inertia of the three regions are set in Earth-like proportions, and the constant angular momentum of the ensemble defines the planetary rotation rate vector |$\boldsymbol{\mathrm{\Omega }}$|⁠. Convection is driven from below using a homogeneous mass anomaly flux F imposed at radius r_i, and a volumetric buoyancy sink term such that there is no homogeneous mass anomaly flux at radius r_o. The mechanical, thermochemical and electromagnetic boundary conditions are respectively of the stress-free, fixed-flux and electromagnetically conducting types at both boundaries. The complete set-up of boundary conditions, lateral heterogeneities and core-mantle-inner core couplings corresponds to the Coupled Earth model described in Aubert et al. (2013); A17; Aubert (2018) and full details on this set-up can be found in these earlier studies. All models produced a self-sustained and dipole-dominated magnetic field with Earth-like geometry (Christensen et al. 2010, A17) that did not reverse polarity during their integration time.

The main model parameters are the flux-based Rayleigh, Ekman, Prandtl and magnetic Prandtl numbers

\begin{eqnarray*} Ra_{F}&=&\dfrac{g_{o}F}{4\pi \rho \Omega ^{3}D^{4}}, \end{eqnarray*}

(1)

\begin{eqnarray*} E&=&\dfrac{\nu }{\Omega D^{2}}, \end{eqnarray*}

(2)

\begin{eqnarray*} Pr&=&\dfrac{\nu }{\kappa }, \end{eqnarray*}

(3)

\begin{eqnarray*} Pm&=&\dfrac{\nu }{\eta }. \end{eqnarray*}

(4)

Here g_o, ρ, ν, κ and η are, respectively, the gravity at radius r_o, the fluid density, viscosity, thermochemical and magnetic diffusivities. We use the path theory (A17) that bridges the parameter space gap between our previous coupled Earth model (Aubert et al. 2013) and Earth’s core conditions by relating these four parameters to a single variable ϵ:

\begin{eqnarray*} Ra_{F}&=&\epsilon Ra_{F} (\mathrm{CE}), \end{eqnarray*}

(5)

\begin{eqnarray*} E&=&\epsilon E (\mathrm{CE}), \end{eqnarray*}

(6)

\begin{eqnarray*} Pr&=&1, \end{eqnarray*}

(7)

\begin{eqnarray*} Pm&=&\sqrt{\epsilon } Pm (\mathrm{CE}). \end{eqnarray*}

(8)

Here Ra_F(CE) = 2.7 × 10⁻⁵, E(CE) = 3 × 10⁻⁵ and Pm(CE) = 2.5 are the control parameters of the coupled Earth dynamo model. Through these definitions, A17 introduced the concept of a unidimensional path in parameter space, starting at ϵ = 1, and ending at ϵ = 10⁻⁷ where reasonable estimates for Earth’s core conditions are obtained (A17). The model cases presented here (Table 1) range from ϵ = 0.1 to ϵ = 3.33 × 10⁻⁴, this latter value defining the model at the middle of the path, which to date is the closest to Earth’s core conditions computed in a numerical simulation. The numerical implementation involves a decomposition of |$\boldsymbol{\mathrm{u}}$|⁠, |$\boldsymbol{\mathrm{B}}$| and C in spherical harmonics up to degree and order ℓ_max and a discretization in the radial direction on a second-order finite-differencing scheme with NR gridpoints. We use the spherical harmonics transform library SHTns (Schaeffer 2013) freely available at https://bitbucket.org/nschaeff/shtns. Full details on the numerical implementation can be found in A17.

Table 1.

Input parameters of numerical models (see the text for definitions). Cases marked with a (^∗) symbol have been previously reported in A17.

Case	ϵ	Path position (per cent)	Run length (overturns)	E	Pm	ℓmax	NR	ℓh	qh
DNS (∗)	0.1	14	164	|$3\,\times 10^{-6}$|	0.8	256	480	–	–
DNS	|$3.33\,\times 10^{-2}$|	21	59	10−6	0.45	320	720	–	–
qDNS	|$3.33\,\times 10^{-2}$|	21	8	10−6	0.45	320	720	290	1.1
qDNS1	10−2	29	33	|$3\,\times 10^{-7}$|	0.25	256	816	196	1.15
qDNS2	10−2	29	3.3	|$3\,\times 10^{-7}$|	0.25	512	1248	448	1.15
qDNS	|$3.33\,\times 10^{-4}$|	50	0.8	10−8	0.045	640	2496	512	1.45
LES (∗)	0.1	14	1172	|$3\,\times 10^{-6}$|	0.8	133	200	30	1.045
LES (∗)	|$3.33\,\times 10^{-2}$|	21	595	10−6	0.45	133	240	30	1.0575
LES1 (∗)	10−2	29	332	|$3\,\times 10^{-7}$|	0.25	133	320	30	1.07
LES2	10−2	29	22.3	|$3\,\times 10^{-7}$|	0.25	256	320	30	1.03
LES (∗)	|$3.33\,\times 10^{-4}$|	50	196	10−8	0.045	133	624	30	1.10

Along the parameter space path, the magnetic Reynolds number Rm = UD/η (where U is the rms value of |$\boldsymbol{\mathrm{u}}$| in the shell) remains approximately constant and Earth-like (Table 2), a defining property of this path. As a consequence, the level of hydrodynamic turbulence increases as we progress towards Earth’s core conditions, as witnessed by the increase of the hydrodynamic Reynolds number Re = UD/ν = Rm/Pm reaching values up to Re ≈ 24 000. A17 explored a parameter range down to ϵ = 0.1 (the first 14 per cent of the path on a logarithmic scale) with DNS simulations. Beyond this point, DNS simulations become computationally very expensive (Schaeffer et al. 2017). Therefore, in A17 the remaining portion to the middle of the path was explored with LES simulations where a hyperdiffusive approximation was used on the velocity and density anomaly fields, but not on the magnetic field that remained fully resolved. The principle is to use values of ℓ_max and NR lower than those required by DNS, and effective diffusivities ν_eff, κ_eff that depend on the spherical harmonic degree ℓ and the molecular diffusivities ν, κ according to Nataf & Schaeffer (2015)

\begin{eqnarray*} (\nu _\text{eff},\kappa _\text{eff})&=&(\nu ,\kappa )\, \mathrm{for}\, \ell \lt \ell _{\mathrm{ h}}, \end{eqnarray*}

(9)

\begin{eqnarray*} (\nu _\text{eff},\kappa _\text{eff})&=&(\nu ,\kappa )\, q_{\mathrm{ h}}^{\ell -\ell _{\mathrm{ h}}} \, \mathrm{for}\, \ell \ge \ell _{\mathrm{ h}}. \end{eqnarray*}

(10)

Here ℓ_h is the cut-off degree below which the hyperdiffusive treatment is not applied, and q_h is the strength of hyperdiffusion. In Table 1 we recall the values of q_h, ℓ_h of a few LES cases from A17 (marked with a star symbol), where the approach was to start hyperdiffusivity at a relatively large scale ℓ_h = 30 and smoothly ramp its strength up towards smaller scales. Here we increase the model resolution (ℓ_max, NR) and test three approaches to relax hyperdiffusivity. First, we introduce new LES simulations where we keep ℓ_h = 30 but use weaker hyperdiffusivity, i.e. values of q_h weaker than in A17. Second, we introduce quasi-DNS simulations (from hereafter qDNS) where we bring ℓ_h close to the smallest numerically resolved scale ℓ_max and employ larger values of q_h than those used in A17. This latter configuration attempts to maximize the range of length scales over which hyperdiffusivity is not applied, and capture as much as possible the physically relevant length scales within this range. Third, we extend the portion of the path where full DNS computations are carried out to 21 per cent, by reporting on a DNS at E = 10⁻⁶ (ϵ = 3.33 × 10⁻²). Due to the increased computational cost, qDNS and DNS simulations have been integrated over a shorter simulation time than their LES counterparts (see run length data in Table 1, presented in units of the overturn time D/U).

Table 2.

Output parameters of numerical models (see the text for definitions). Cases marked with a (^∗) symbol have been previously reported in A17.

Case	Path position (per cent)	|$\dfrac{4\pi D^{2} p}{g_{o} F}$|	|$f_\Omega$|	|$Rm=\dfrac{UD}{\eta }$|	|$\dfrac{B}{\sqrt{\rho \mu \eta \Omega }}$|	|$\mathcal {T}$|	ℓΩ	ℓMS	ℓ⊥	ℓCIA	ℓVAC
DNS (∗)	14	0.332	0.81	1092	4.52	0.17	159	71	10	76	207
DNS	21	0.335	0.89	1071	4.41	0.11	166	70	10	110	>320
qDNS	21	0.331	0.88	1050	4.33	0.11	167	78	11	105	295
qDNS1	29	0.346	0.92	1130	4.49	0.07	171	60	10	199	212
qDNS2	29	0.340	0.93	1095	4.49	0.07	170	65	12	182	455
qDNS	50	0.343	0.97	1105	3.90	0.02	196	130	13	522	522
LES (∗)	14	0.336	0.72	1046	4.64	0.17	146	62	9	61	73
LES (∗)	21	0.343	0.76	1036	4.57	0.12	150	65	10	72	79
LES1 (∗)	29	0.348	0.80	1046	4.70	0.08	152	60	10	87	83
LES2	29	0.347	0.86	1084	4.61	0.09	162	67	12	122	127
LES (∗)	50	0.360	0.84	1089	4.48	0.04	164	81	12	101	92

Time-averaged diagnostic outputs of the simulation are presented in Table 2. Listed first are the convective power density p in units of g_oF/4πD², its fraction |$f_\Omega$| that is dissipated by ohmic losses (also represented in Fig. 1 as a function of the path parameter), the magnetic Reynolds number Rm presenting the shell rms velocity U in units of η/D, the shell rms magnetic field B in Elsasser units of |$\sqrt{\rho \mu \eta \Omega }$| (where μ is the magnetic permeability, the resulting dimensionless value representing the square root of the classical Elsasser number, A17). We also report on the level of enforcement of the Taylor constraint by computing the level of cancellation |$\mathcal {T}$| of azimuthal magnetic force integrated over geostrophic cylinders as in A17. The magnetic dissipation length scale |$d_\Omega$| is defined as in A17 (where it was labelled as d_min) from the square root of the ratio of volume-averaged magnetic energy and magnetic dissipation, and reported as an equivalent harmonic degree ℓ_Ω = πD/d_Ω.

Figure 1.

Evolution of the fraction |$f_\Omega$| of convective power dissipated by Ohmic losses with the path parameter ϵ, in LES, qDNS and DNS simulations (Earth’s core conditions are towards the left of this graph). Simulation results already reported in A17 are presented in grey, and new results are reported in blue. The shaded regions represent the ±1 std. dev. of fluctuations relative to the time average.

As previously in A17, we compute scale-dependent force balance diagrams (Fig. 2) that enable the investigation of the spatial force balance structure. To represent length scales, we prefer the harmonic degree ℓ over the harmonic order because of the invariance of ℓ-dependent energy spectra against rotation of the coordinate system, and because unlike an harmonic order, a given degree ℓ in the spectral space can be associated to a minimal corresponding length scale πD/ℓ in the physical space. It has been a common practice to use curled forces rather than the forces themselves to examine their contributions in the Navier-Stokes equation (e.g. Christensen & Aubert 2006; Davidson 2013). This approach is however misleading if these curls are then analysed in the spectral space. Indeed, the contributions at different harmonic degrees are coupled by the curl operation, thereby defeating the possibility to associate a force at degree ℓ to its action at length scale πD/ℓ in the physical space. For this reason, and also because we wish to also assess the role of the pressure force, it is mandatory to analyse the non-curled forces in the spectral space (as previously done in A17). Crossover length scales are obtained in Fig. 2 and reported in Table 2 by determining the harmonic degrees at which two forces of the system become of equal amplitude. These may be thought of as scales where a given force balance is achieved, and they usually bound scale ranges where different dynamical equilibria are enforced. The temporal variability of force balance diagrams is weak (A17) and they are therefore computed from snapshots in time. The typical error on crossover length scales is then less than 10 per cent. Here we recall some key results found in A17 for the force balance structure and introduce new crossover length scales. These results found along the parameter space path are typical of the situation most commonly found in a wide region of the parameter space (Schwaiger et al. 2019), with the exceptions being found at low convective supercriticality and low magnetic Prandtl number Pm. The different forces span up to five orders in magnitude (with unprecedently low levels of viscosity being reached in the qDNS at 50 per cent of the path, Fig. 2c), involving at least three successive orders of force balances. Fig. 2 shows that at large scales (low harmonic degrees), a leading order quasi-geostrophic balance is achieved between the Coriolis and pressure forces. The crossover of the Coriolis and Lorentz forces defines the magnetostrophic harmonic degree ℓ_MS. For harmonic degrees larger than ℓ_MS, the zeroth-order balance becomes magnetostrophic, with the pressure force equilibrating the Lorentz force (which is hence essentially a magnetic pressure) and the residual Coriolis force. Considering now the first-order balance, the crossover harmonic degree of the Lorentz and buoyancy forces is defined as ℓ_⊥. As this balance is achieved under the constraint of the part of the Coriolis force not balanced at zeroth order by pressure (the ageostrophic Coriolis force), the associated length scale is usually referred to as that of the MAC (magneto-Archimedes-Coriolis) balance (Davidson 2013; Yadav et al. 2016; Aubert 2018, A17). For ℓ < ℓ_⊥, the first-order balance is of thermal wind type between the ageostrophic Coriolis and buoyancy forces, and for ℓ > ℓ_⊥ it becomes of magnetic wind type between the ageostrophic Coriolis and Lorentz force. Moving to the next-order force balances, the crossover harmonic degrees of the inertial and buoyancy forces ℓ_CIA, and of the viscous and buoyancy forces ℓ_VAC are, respectively, representative of the CIA (Coriolis-inertia-Archimedes) and VAC (viscous-Archimedes-Coriolis) force balances that operate under the constraint of the Coriolis force residual not balanced at zeroth and first orders (the amagnetostrophic Coriolis force). A summary diagram presenting the relative positions of the main length scales is presented in Fig. 3. Scale separation in this system refers to the difference between the large scale of the first-order MAC balance at degree ℓ_⊥, and the much smaller scales of the second-order force balances at degrees ℓ_CIA, ℓ_VAC and of magnetic dissipation at degree ℓ_Ω.

Figure 2.

Scale-dependent force balance diagrams for direct and quasi-direct numerical simulations (DNS and qDNS) obtained along the parameter space path. Represented are the rms amplitudes of each force as functions of the spherical harmonic degree ℓ (see A17 for technical computation details). Forces are normalized relative to the maximum value of the pressure force. The ageostrophic Coriolis force represents the residual Coriolis force after removal of the pressure force. The amagnetostrophic Coriolis force represents the residual Coriolis force after removal of the pressure and Lorentz forces. The crossover length scales ℓ_MS, ℓ_⊥, ℓ_CIA and ℓ_VAC are defined from the crossings observed in these diagrams (see the text for interpretation, values in Table 2 and a summary diagram in Fig. 3).

Figure 3.

Relative positions of key length scales in the simulations, represented as equivalent spherical harmonic degrees. The horizontal bands represent the extent of the spherical harmonics expansion (with the right end representing the maximum degree ℓ_max), and their colour delineate the ranges where hyperdiffusivity is applied (pink) or not (grey). Note that ℓ_Ω is based on the integral length scale d_Ω encompassing magnetic diffusion in all three spatial dimensions, hence it is possible that ℓ_Ω exceeds ℓ_max. The degree ℓ_⊥ is indicative of the large scale at which the first-order MAC balance is achieved. Only the smallest of the two harmonic degrees ℓ_CIA, ℓ_VAC is represented, thereby indicating the nature of the force balance achieved at second order.

2.2 Asymptotic QG-MAC theory

The results from the upsized numerical simulations will be analysed in the light of the asymptotic theory proposed by Davidson (2013), some elements of which are recalled here. We consider a leading-order, quasi-geostrophic (QG) equilibrium and a first-order magneto-Archimedes-Coriolis (MAC) force balance, both well supported by the numerical simulations (Fig. 2), together with the dynamo power budget. Due to the zeroth-order QG equilibrium and the associated Taylor-Proudman constraint, the scale |$d_{\mathbin {\!//\!}}$| of flow structures along the rotation axis |$\boldsymbol{\mathrm{\Omega }}$| remains at the system size, that is, |$d_{\mathbin {\!//\!}}\sim D$| (see Fig. 5). The first-order MAC balance operating at the large scale d_⊥ = πD/ℓ_⊥ can be written from the curl of the Navier-Stokes equation (see A17):

\begin{eqnarray*} \dfrac{\rho \Omega U}{D} \sim \dfrac{g_{o} C}{d_{\perp }} \sim \dfrac{B^{2}}{\mu d_{\perp }^{2}}, \end{eqnarray*}

(11)

where C is a typical value for the density anomaly field. In the present context where Ohmic losses account for the essential part of the dissipated energy (Fig. 1), the balance between the convective power input and magnetic dissipation (Christensen & Tilgner 2004; Christensen & Aubert 2006) may be written as

\begin{eqnarray*} \dfrac{\eta B^{2}}{\mu d_\Omega ^2} \sim p\sim \dfrac{g_{o} F}{D^{2}}, \end{eqnarray*}

(12)

where we have used the equivalence between convective power and mass anomaly flux (Christensen & Aubert 2006; Aubert et al. 2009, A17), as attested by the almost constant value of 4πD²p/g_oF in Table 2. Combining eqs (11) and (12), and using F ∼ UCD², we obtain the vorticity equivalence

\begin{eqnarray*} U/d_{\perp }= c_{\omega } ~\eta /d_\Omega ^{2}. \end{eqnarray*}

(13)

This equivalence states that the vorticities at the large scale d_⊥ where the magnetic field is sustained by drawing convective power (balance between buoyancy and Lorentz forces) and at the small scale |$d_\Omega$| where this power is ohmically dissipated are related by a proportionality constant c_ω (interpreting |$\eta /d_\Omega ^{2}$| as a vorticity implicitly assumes that the gradient and rotational parts of |$\boldsymbol{\mathrm{u}}$| are of similar magnitude at the scale of magnetic diffusion).

2.2.1 Self-similar ranges in spatial spectra

Despite the fact that power is essentially carried by magnetic (rather than hydrodynamic) nonlinearities from the large injection scale d_⊥ down to the small dissipation scale |$d_\Omega$|⁠, we do not expect spatial spectra of the magnetic energy density to feature self-similar, power-law behaviours in the range |${}[\ell _\perp ,\ell _\Omega ]$| because the level of magnetic turbulence is only moderate (Rm is moderate; Table 2). However, the hydrodynamic Reynolds number Re = Rm/Pm reaches values up to 24 000 along the explored part of the path, suggesting that such ranges are likely in spectra of the velocity field. The vorticity equivalence (13) further supports this idea. Indeed, the constant c_ω being in principle of order 1, at large values |$d_\perp /d_\Omega$| of the scale separation we expect the spectral enstrophy density ω²(ℓ) at degree ℓ (where |$\boldsymbol{\mathrm{\omega }}=\nabla \times \boldsymbol{\mathrm{u}}$| is vorticity) to remain approximately constant, that is,

\begin{eqnarray*} \omega ^{2}(\ell ) \propto \ell ^{0}\, \mathrm{for}\, \ell _{\perp }\lt \ell \lt \ell _{\Omega }. \end{eqnarray*}

(14)

With length scales varying like ℓ⁻¹, the energy density of velocity u²(ℓ) should then also present a power-law range

\begin{eqnarray*} u^{2}(\ell ) \propto \ell ^{-2}\, \mathrm{for}\, \ell _{\perp }\lt \ell \lt \ell _{\Omega }. \end{eqnarray*}

(15)

Deviations from these ideal spectra are however likely. First, length scales are indeed anisotropic in this system, with flow structures elongated along the rotation axis (⁠|$d_{\perp } \ll d_{\mathbin {\!//\!}}$|⁠), implying that the velocity energy spectrum is less steep than ℓ⁻² if the vorticity spectrum is flat. In that sense, the spectra of the vorticity component parallel to the rotation axis |$\omega _{\mathbin {\!//\!}}=\left(\nabla \times \boldsymbol{\mathrm{u}}\right)\cdot \boldsymbol{\mathrm{\Omega }}/\Omega$| and velocity component perpendicular to the rotation axis |$\boldsymbol{\mathrm{u}}_{\perp }=\boldsymbol{\mathrm{u}}-\left(\boldsymbol{\mathrm{u}}\cdot \boldsymbol{\mathrm{\Omega }}/\Omega \right)$| should better approach (14, 15) than those of the full vorticity |$\boldsymbol{\mathrm{\omega }}$| and velocity |$\boldsymbol{\mathrm{u}}$|⁠, because the former leave out the variations along the rotation axis. Second, the scale separation |$d_\perp /d_\Omega$| only reaches values up to 20 in our simulations (Fig. 3). Finally, the correspondence between length scales and harmonic degrees is not straightforward (e.g. Schaeffer et al. 2017), because a given harmonic degree can represent different length scales at different radii, and because different dynamics may exist in the regions inside and outside the axial cylinder tangent to the inner core (the tangent cylinder). We will nevertheless use (14, 15) as guidelines for the interpretation of spatial spectra.

2.2.2 Asymptotic scaling laws along the parameter space path

We now recall the derivation of the scaling laws for U, B, d_⊥ and d_Ω along the parameter space path, as functions of the path parameter ϵ. If one requires the large- and small-scale vorticities present in eq. (13) to be independent on the system rotation rate and diffusivities, then from (12) it follows that the magnetic field itself is independent on the rotation rate and diffusivities. Davidson (2013) pointed out that dimensional analysis finally yields the following scaling, corresponding to the initial proposal of Christensen & Aubert (2006):

\begin{eqnarray*} B \sim \sqrt{\rho \mu } (g_{o}F/\rho D)^{1/3}. \end{eqnarray*}

(16)

Combining eqs (11, 12, 16) and using Ra_F ∼ ϵ leads to the scalings previously proposed in A17, which are presented here in a slightly different manner to account for our choice of dimensionless quantities:

\begin{eqnarray*} \dfrac{B}{\sqrt{\rho \mu \eta \Omega }} \sim \epsilon ^{1/12}, \end{eqnarray*}

(17)

\begin{eqnarray*} Rm=\dfrac{UD}{\eta }\sim \epsilon ^{-1/18}, \end{eqnarray*}

(18)

\begin{eqnarray*} d_{\perp }/D \sim \epsilon ^{1/9}, \end{eqnarray*}

(19)

\begin{eqnarray*} \ d_{\Omega }/D \sim \epsilon ^{1/12}. \end{eqnarray*}

(20)

The powers of ϵ entering (17–20) imply that these quantities only undergo weak variations along the path (as confirmed by Table 2), from ϵ = 1 down to ϵ = 10⁻⁷. This is another evidence of the continuum that exists along the parameter space path between numerical simulations and Earth’s core conditions. In particular, the weak variations of the length scales d_⊥ and d_Ω (see Fig. 3) led A17 to propose a scale-invariant theory (the path theory, see also Starchenko & Jones 2002) tailored for the analysis of LES simulations, where these two length scales are fixed, which in turn implies constant values of Rm and |$B/\sqrt{\rho \mu \eta \Omega }$|⁠. The weak values of scaling exponents may initially obscure the prospect of checking these theories against numerical data, but A17 have shown that this is not the case owing to the broad range of accessible path parameter values ϵ and the low scatter of numerical data. However, an uncertainty analysis will be needed for the upsized runs, which have a short integration time (Table 1) and may therefore present more scatter than the LES simulations analysed in A17.

3 RESULTS

3.1 Transients

Fig. 4 presents a typical temporal evolution sequence of output diagnostic quantities during a gradual upsizing process starting from LES1 (Table 1), a large-eddy simulation at 29 per cent of the parameter space path (Ekman number |$E=3\,\times 10^{-7}$|⁠) taken from A17. At each increase of spatial resolution, the simulation is restarted using the output of the previous step. Transients associated to this process are short, typically a fraction of a convective overturn. As we shall see, the process does not alter the leading-order QG-MAC force balance (Fig. 2, see also Fig. 9 below) and large-scale simulation structure (see Figs 5–7), leaving only the smaller scales to reach their saturated energy after the upsizing step. Integral diagnostics for power, velocity and magnetic field (Figs 4a–c) are therefore largely invariant against upsizing, with the upsized simulations even following the evolution of the previous lower-resolution case for about an overturn time after the upsizing step. The magnetic dissipation length scale d_Ω (Fig. 4d) is more sensitive to the distribution of small scales and hence more influenced by upsizing, until a high enough resolution is reached in qDNS2 to fully encompass ℓ_Ω = π/d_Ω (see Fig. 3). As a corollary, the fraction f_Ω of Ohmic dissipation in the system (Fig. 1) generally increases during the upsizing procedure, and reaches a value independent of further upsizing when d_Ω also converges (see similarity of f_Ω, ℓ_Ω for qDNS1 and qDNS2 in Table 2). At 50 per cent of the path, the upsizing process thereby leads to a value f_Ω = 0.97 indicative of almost entirely Ohmic dissipation. We are therefore confident that despite being integrated for a small amount of overturn times, the upsized simulations are fully representative of the saturated dynamical state that could be reached at much greater computational cost in a regular DNS.

Figure 4.

Temporal evolution of output diagnostics during an upsizing sequence at 29 per cent of the parameter space path (⁠|$E=3\,\times 10^{-7}$|⁠). At each upsizing step, the restart file from the previous step is used. Case LES1 from A17 (Table 1, black) is first expanded into LES2 (purple) with the same cut-off harmonic degree ℓ_h = 30 but a lower value of the hyperdiffusion strength q_h. Two quasi-DNS simulations (qDNS1, qDNS2) are then performed, first by expanding the spherical harmonic truncation to ℓ_max = 256 and the hyperdiffusivity cut-off to ℓ_h = 196 (blue), and then by expanding ℓ_max to 512 and ℓ_h to 448. For clarity, only a portion of the first three sequences (black, purple, blue) is represented. Because of an incorrect specification of the numerical code options, the diagnostic outputs presented here were not available during a short portion of the last sequence (green dashes).

Figure 5.

Planforms of velocity |$\boldsymbol{\mathrm{u}}$| components in the cylindrical radial direction (a,b), the azimuthal direction (c,d) and planforms of the axial vorticity |$\omega _{\mathbin {\!//\!}}$| (e,f). Represented are the fields in the equatorial plane, in a half meridian plane, at the outer (core–mantle) boundary of the shell at radius r_o, and at radius r/D = 0.545, close to the inner (inner core) boundary. The position of the rotation axis |$\boldsymbol{\mathrm{\Omega }}$| is marked by a black line on all panels. The left-hand column (a,c,e) refers to a snapshot taken during case LES1 from A17 at 29 per cent of the path (Ekman number |$E=3\,\times 10^{-7}$|⁠, Table 1), and the right-hand column to a snapshot of the highest-resolution case qDNS2 also obtained at 29 per cent of the path after the upsizing sequence (green lines in Fig. 4). Velocity is presented in magnetic Reynolds units of η/D, and vorticity in inverse magnetic diffusion time units of η/D².

Figure 6.

Planforms of the density anomaly field C, presented in units of F/4πΩD³. The snapshots are taken at the same instants in time as in Fig. 5, for cases LES1 (a) and qDNS2 (b) obtained at 29 per cent of the path. Same visualization conventions as in Fig. 5.

Figure 7.

Planforms of the radial component of the magnetic field |$\boldsymbol{\mathrm{B}}$|⁠, presented in Elsasser units of |$\sqrt{\rho \mu \eta \Omega }$|⁠. The snapshots are taken at the same instants in time as in Fig. 5, for cases LES1 (a) and qDNS2 (b) obtained at 29 per cent of the path. Same visualization conventions as in Fig. 5.

3.2 Planforms in physical space

Planforms of the velocity, density anomaly and magnetic fields in the physical space are presented in Figs 5–7, for the initial case LES1 and final case qDNS2 of the upsizing sequence at 29 per cent of the path (Table 1, Fig. 4). The structure of |$\boldsymbol{\mathrm{u}},\boldsymbol{\mathrm{B}},C$| at large scales is remarkably preserved after upsizing, confirming the accuracy of LES, as previously advocated for in A17. These large scales consist in upwellings taking the shape of radially elongated columnar sheets (Figs 5a,b) associated with large buoyancy plumes (Fig. 6a,b), entrained in a mostly columnar and westward azimuthal flow presenting rim-like polar circulations on the tangent cylinder (Figs 5c,d). Relative to its LES1 counterpart, the qDNS2 simulation however features significantly enriched small scales and a higher level of separation between the largest and smallest observable scales. These small scales comprise columnar flow filaments coexisting with the sheet-like cylindrical radial upwellings (Fig. 5b) and ramifications of the buoyancy anomaly plumes (Fig. 6b). The small-scale content is best seen on planforms of the axial vorticity |$\omega _{\mathbin {\!//\!}}$| (Fig. 5f) where the columnar filaments are gathered at the edges of the large-scale velocity structures. As we shall see below in Figs 10–12, the large scale in |$\boldsymbol{\mathrm{u}}$|⁠, |$\boldsymbol{\mathrm{B}}$|⁠, C is d_⊥, the scale of the MAC balance at which the magnetic field draws its power from convection, and the dominant small scale seen in |$\boldsymbol{\mathrm{\omega }}$| in the qDNS2 case corresponds to d_Ω where this power is ohmically dissipated. In the LES1 case (Figs 5a,c,d and 6a), the scale separation is much weaker, and we shall see the small scale of |$\boldsymbol{\mathrm{\omega }}$| in this case correspond to the balance of hyperdiffusive viscosity and buoyancy force (the VAC balance at harmonic degree ℓ_VAC). Figs 7(a) and (b) show that the upsizing process does not cause an increase in the scale separation of the magnetic field |$\boldsymbol{\mathrm{B}}$|⁠, the content of which remains largely at scale d_⊥. Indeed, in case of LES1 the magnetic field is already close to being fully resolved, because the low value of the magnetic Prandtl number Pm warrants a modest degree of magnetic turbulence Rm at strong forcing. Analysing the magnetic field geometry of the upsized case qDNS2 (Fig. 8), it can be seen that field lines tend to avoid regions of strong velocity and also gather at the edges of upwelling sheets or density anomaly plumes. This leads to a remarkable alignment between small-scale vorticity filaments and magnetic field lines, which is commonly referred to as dynamic alignment in the theory of strong MHD turbulence (see Tobias et al. 2012). This is the process through which the magnetic field adjusts so as to minimize its shear by the flow, as a consequence of Lenz law. The Lorentz force therefore constrains the flow but ultimately becomes dynamically irrelevant, reducing its contribution at small scales to a magnetic pressure (Fig. 2). The planforms of qDNS2 are broadly comparable to those of recently published DNSs with strong scale separation (Schaeffer et al. 2017; Sheyko et al. 2018). However, less accumulation of light material within the tangent cylinder is observed in qDNS2 than in simulation S2 from Schaeffer et al. (2017), and we do not observe the generation of reverse surface magnetic flux inside this cylinder (compare Figs 6a,b to 7 in Schaeffer et al. 2017). Case qDNS2 indeed features a stronger level of forcing (stronger hydrodynamic and magnetic Reynolds numbers) than the cases from Schaeffer et al. (2017) and Sheyko et al. (2018), and this presumably improves the exchange of material through the tangent cylinder, leading to a spatially more homogeneous system.

Figure 8.

Polar view of the axial vorticity field |$\omega _{\mathbin {\!//\!}}$| from case qDNS2 at 29 per cent of the path, presented in the entire equatorial plane and at radius r/D = 0.545, with black lines delineating the cut presented in Fig. 5(f) (same snapshot and colour coding as this figure). Magnetic field lines parallel to the local direction of |$\boldsymbol{\mathrm{B}}$| are also represented (orange), with thickness proportional to the local magnetic energy density |$\boldsymbol{\mathrm{B}}^{2}$|⁠. In order to see field lines in the whole volume, the equatorial plane is made selectively transparent with opacity proportional to |$|\omega _{\mathbin {\!//\!}}|$|⁠. This reveals an organization of field lines aligned with the columnar vorticity filaments and along the rotation axis.

3.3 Force balances at leading and first order

Fig. 9 presents scale-dependent force balances in simulations with variable degrees of upsizing (LES, qDNS, DNS) at a fixed position on this path (21 per cent). This complements Fig. 2 where force balances are presented for DNS and qDNS along the parameter space path. Invariance of the zeroth-order quasi-geostrophic and magnetostrophic balances, as well as the first-order MAC balance was observed along the path in the LES of A17, and is confirmed here in the upsized qDNS and DNS. The LES are furthermore efficient at accurately capturing these balances (compare Figs 9a and c). The value of the crossover length scale to magnetostrophy ℓ_MS is also broadly invariant in our sets of runs, except for the qDNS at 50 per cent of the parameter space path (Fig. 2c), where the force balance is perturbed by the choice of upsizing strategy (see below). Invariance of ℓ_MS is consistent with the analysis of Aurnou & King (2017) where this scale is conjectured to depend solely on the Elsasser number B²/ρμηΩ and magnetic Reynolds number Rm, both of which are about constant along the path (Table 2). The MAC length scale ℓ_⊥ is also almost invariant along the parameter space and against upsizing (see also Fig. 3). In summary, the QG-MAC force balance previously advocated for in Davidson (2013), Yadav et al. (2016), A17, Calkins (2018) is therefore confirmed here as the leading-order balance asymptotically holding in rapidly rotating regimes, and it is shown to be remarkably stable even in the high-resolution, upsized simulations.

Figure 9.

Scale-dependent force balances for LES (a), qDNS (b) and DNS (c) simulations at 21 per cent of the parameter space path (Ekman number E = 10⁻⁶). See Fig. 2 and the text for the definitions of crossover harmonic degrees, and A17 for technical computation details.

The projection of these leading and first-order force balances on axisymmetric scales yields the Taylor constraint, stating that the integral of azimuthal Lorentz force over imaginary cylinders co-axial to the rotation axis vanishes to the level of the residual inertia and viscosity (see A17). Since the magnetic field structure is largely preserved by upsizing (Figs 7a,b) the level of enforcement |$\mathcal {T}$| of the Taylor constraint is also preserved (Table 2) and increases along the parameter space path. We again emphasize that adherence to the Taylor constraint should not be misinterpreted as a sign of leading-order magnetostrophy, as it follows here from leading-order quasi-geostrophy.

3.4 Second-order force balance and optimal approximations

Differences between upsizing strategies emerge at the next orders in the force balance diagrams. DNS simulations (Figs 2a and 9c, see also A17) show that the second-order force balance is of the CIA nature. It has been argued that the VAC balance could influence the dynamics of dynamo simulations at moderate parameter values (King & Buffett 2013) but this is clearly not the case here as this balance only comes at the third order in amplitude. Quasi-DNS simulations constructed such that the hyperdiffusive cut-off ℓ_h largely exceeds ℓ_CIA (such as qDNS at 21 per cent and qDNS2 at 29 per cent of the path; Fig. 3) have a force balance structure that is undistinguishable from that of the full DNS down to the second order in amplitude (compare Figs 9b and c, see also Fig. 2b). At third order, the VAC balance length scale ℓ_VAC is improperly resolved in these qDNS but this has no consequence on the partitioning of dissipation between Ohmic and viscous losses (compare the f_Ω values in qDNS and DNS at 21 per cent of the path in Fig. 1). This type of qDNS may therefore be safely used to fully replace a DNS, but it still involves a sizeable computational cost because ℓ_CIA becomes large along the path (Fig. 3). As an acceptable tradeoff between cost and accuracy, a qDNS simulation where the cut-off is close to, or smaller than the true value of ℓ_CIA, but still larger than ℓ_Ω may be performed (respectively, cases qDNS1 at 29 per cent of the path and case qDNS at 50 per cent of the path, Fig. 3). This type of qDNS still captures the Ohmic dissipation fraction f_Ω well, and in any case better than LES (compare LES, qDNS1 and qDNS2 at 29 per cent of the path in Fig. 1), but the stronger approximation implies a perturbed force balance structure at high harmonic degree (Fig. 2c) and an inaccurate determination of ℓ_CIA (compare qDNS1 and qDNS2 in Fig. 3). Finally, LES simulations (Fig. 9a) tend to render a mix of CIA and VAC force balances at second order, with ℓ_CIA ≈ ℓ_VAC and both values being significantly decreased relative to DNS (Table 2, Fig. 3). At path positions of 29 per cent and beyond the ordering of these two length scales is even reversed compared to the DNS ordering (i.e. ℓ_VAC < ℓ_CIA, see violet symbols in Fig. 3), implying that VAC has replaced CIA as the second-order force balance. This implies that the determination of f_Ω is inaccurate (Fig. 1), as already reported in A17.

3.5 Spatial energy spectra and power-law ranges

Differences between upsizing strategies can be further assessed by examining spatial spectra of the kinetic and magnetic energy densities (Fig. 10). Here again, a representation using the spherical harmonic degree ℓ is used, but similar results can be obtained by using the spherical harmonic order instead. All simulations have broadly similar magnetic energy spectra throughout the spherical harmonic degree range, and also similar kinetic energy spectra up to spherical harmonic degree ℓ ≈ 60, underlining the relevance of LES at such large scales and the weak impact of upsizing on these scales. In particular, the peak of kinetic and magnetic energy spectra for ℓ > 1 occurs at ℓ = ℓ_⊥, confirming that ℓ_⊥ is the dominant scale of velocity, density anomaly and magnetic field planforms (Figs 5–7). There is no evidence of a self-similar, power-law range in spectra of the magnetic energy density. This is expected given the modest level Rm ≈ 1000 of magnetic turbulence (Table 2). In kinetic energy spectra of the LES simulations, a power-law decay range cannot be observed because it is obscured by hyperdiffusion kicking in at ℓ_h = 30 (red and pink curves in Figs 10a,b), but it appears gradually in qDNS and DNS simulations where the hyperdiffusive cut-off ℓ_h is increased close to ℓ_max or removed altogether. It is again advisable to use a sufficiently resolved qDNS simulation (dark green curves in Figs 10a,b) for best accuracy (compare the dark green and black curves in Fig. 10a). In particular, choosing ℓ_h above ℓ_Ω, but close to or below ℓ_CIA causes energy stacking effects where the tail of the spectrum overshoots the true value (compare light green and dark green curves in Fig. 10b). This in turn causes the buoyancy and Coriolis force lines to overshoot their true value in the force balance diagrams (Fig. 2c).

Figure 10.

Instantanous kinetic (KE) and magnetic (ME) energy density spectra, as functions of the spherical harmonic degree ℓ, in DNS, qDNS and LES simulations performed at 21 per cent (a) and 29 per cent (b) of the parameter space path (Ekman numbers E = 10⁻⁶ and |$E=3\,\times 10^{-7}$|⁠). Spectra are normalized such that the total magnetic energy is one, illustrating the increasing dominance of the magnetic energy as we progress along the parameter space path.

Plotting together uncompensated and compensated kinetic energy spectra (Figs 11a,b) of our best-resolved qDNS and DNS simulations along the path further reveals that the simulations are supportive of a gradually broadening range with power-law decay ℓ^−3/2 starting from the large scale ℓ_⊥. Less steep decay ranges have been previously observed in strongly scale-separated dynamo simulations (Schaeffer et al. 2017; Sheyko et al. 2018), with spectral indices from −1 to −4/3, and evidence of dependence of spectral index with the simulation forcing has also been reported. The steeper decay observed here occurs in a context of even stronger forcing (stronger Re and Rm) than in these previous simulations and approaches closer to the index −2 predicted from vorticity equivalence (Section 2.2.1). In our three simulations where resolution is high enough to avoid energy stacking at the spectrum tail (black, red and green curves in Fig. 11), the extent of the power-law range is measured through the harmonic degree at which the compensated spectra of Fig. 11(b) fall below 85 per cent of the plateau energy. We then find that the corresponding harmonic degrees match the lowest of the two scales |$\ell _\mathrm{CIA}, \ell _\Omega$| (see their relative position in Fig. 3). This is theoretically expected because vorticity equivalence should be enforced down to the scale of magnetic dissipation, unless the MAC balance is perturbed at small scales by the second-order CIA balance. In our qDNS simulation at 50 per cent of the parameter space path (blue curve in Fig. 11) the same analysis cannot be carried out because of energy stacking at the spectrum tail, but we note that stacking again starts at ℓ_Ω which in this case is significantly lower than ℓ_CIA (Figs 2c and 3). Power transfers from injection at ℓ_⊥ down to Ohmic dissipation at ℓ_Ω are well controlled in this simulation (hence the flat compensated spectral range), but the transfer to even smaller scales of the residual power (3 per cent in this case, Fig. 1) for viscous dissipation via inertial forces is improperly rendered (because the true value of ℓ_CIA is not resolved), thereby causing the stacking of energy. Though this stacking appears significant in the compensated spectra (Fig. 11b), it should not be overemphasized as it represents a small fraction of the kinetic energy (Fig. 11a). Stacking is also commonly observed in fully resolved DNS (see the tails of the black and red curves in Figs 11a,b), and in previously published simulations with moderate resolution (see fig. 3b of Sheyko et al. 2018) without bearing too much consequence on the simulated dynamics.

Figure 11.

(a) Instantaneous kinetic energy density spectra in DNS and best-resolved qDNS simulations along the parameter space path, as functions of the spherical harmonic degree ℓ. Spectra are normalized in each case relatively to the total kinetic energy. (b) kinetic energy density spectra compensated with ℓ^3/2. The horizontal dashed line in (b) represents 85 per cent of the flat level of compensated kinetic energy density.

In order to further establish the link between Davidson’s vorticity equivalence (eq. 13) and the power-law range observed in Figs 10 and 11, spatial spectra of enstrophy ω² are plotted together with kinetic energy spectra in Fig. 12 for case qDNS2 at 29 per cent of the parameter space path. At scales between the MAC scale ℓ_⊥ and the magnetic dissipation scale ℓ_Ω, the enstrophy ω² presents a weakly increasing power-law range ω²∝ℓ^1/3. Given that scale separation covers a decade in this example, this implies that the amplitude of vorticity varies by less than a factor 2 between the large and small scales, consistent with the idea of a flat spectrum underlain by vorticity equivalence (eq. 14). The vorticity spectrum peaks at ℓ = ℓ_Ω, showing that the dominant scale of the vorticity planforms observed in Fig. 5(f) is indeed the scale of magnetic dissipation (the same analysis carried out in the case of the LES case of Fig. 5(e) reveals that ℓ_VAC is the dominant scale). Note that a length scale defined as δ(ℓ) = u(ℓ)/ω(ℓ) ∼ Dℓ^−11/12 does not exactly vary like ℓ⁻¹, underlining weak but measurable effects of the anisotropic flow structure (Figs 5b,d). Computing the spectra of axial enstrophy |$\omega _{\mathbin {\!//\!}}^{2}$| removes the influence of the long length scale parallel to the rotation axis, and reveals a better respected vorticity equivalence |$\omega _{\mathbin {\!//\!}}^{2} \sim \ell ^{1/4}$| in the range [ℓ_⊥, ℓ_Ω]. Likewise, the kinetic energy spectrum of flow perpendicular to the rotation axis follows |$u_{\perp }^{2} \propto \ell ^{-7/4}$|⁠, better approaching the spectral index −2 (eq. 15), and the associated perpendicular length scale |$\delta _{\perp }(\ell )=u_\perp (\ell )/\omega _{\mathbin {\!//\!}}(\ell )$| now indeed varies like ℓ⁻¹, confirming the relevance of spherical harmonic degrees to analyse length scales.

Figure 12.

Instantaneous energy density spectra of the enstrophy ω², kinetic energy u², axial enstrophy |$\omega _{\mathbin {\!//\!}}^{2}$| and kinetic energy of flow perpendicular to the rotation axis |$u_{\perp }^{2}$|⁠, in case qDNS2 at 29 per cent of the parameter space path (Ekman number |$E=3\,\times 10^{-7}$|⁠). Spectra are normalized relatively to the total kinetic energy.

3.6 Scaling properties

The outputs from our upsized qDNS and DNS simulations at highest resolution are expected to approach the asymptotic QG-MAC scaling theory of Davidson (2013), and to depart from the scale-independent path theory outlined in A17 and verified on LES simulations (see Section 2.2.2). Fig. 3 indeed shows that ℓ_⊥ in qDNS and DNS simulations tends to be slightly higher than in LES simulations, and to slightly increase along the path, both results being consistent with these expectations. However these results are not fully systematic, presumably because of our evaluation of ℓ_⊥ on snapshots and the associated uncertainty, which while being small (about 10 per cent) remains on the order of the observed variations. Furthermore, from eq. (19) a doubling of ℓ_⊥ = π/d_⊥ is predicted between 14 and 50 per cent of the path, but clearly not observed. The departure of our results for ℓ_⊥ from the QG-MAC prediction simply reflects the fact that the scale separation d_⊥/d_Ω in our simulations actually increases along the path, while from eqs (19 and 20), it is (rather counter-intuitively) expected to decrease in the QG-MAC theory.

To refine this analysis, we further examine in Fig. 13 the evolution along the path of the dimensionless velocity Rm, magnetic field |$B/\sqrt{\rho \mu \eta \Omega }$| (corrected as in A17 with the inverse square root of |$f_\Omega$| to account for incompletely Ohmic dissipation at the start of path), and magnetic dissipation length scale d_Ω/D. To evaluate the uncertainties that are due to the shortness of DNS and qDNS runs, our longest LES simulation at 14 per cent of the path (Table 1) is used as a reference. Error bars for other cases are then computed by evaluating in this reference case the maximum deviation of outputs averaged over the short duration of these cases versus the long time average of the reference. Power laws of DNS and qDNS are then obtained from the numerical data by using least-squares weighted with the resulting uncertainties. The magnetic Reynolds number of upsized simulations (Fig. 13a) remains about constant along the path, with a scaling Rm ∼ ϵ^{0.004 ± 0.04}. This result is close to the prediction from the path theory (constant Rm, A17), and deviates from the asymptotic QG-MAC prediction Rm ∼ ϵ^−1/18 ∼ ϵ^−0.056 (eq. 18, yellow line in Fig. 13a) in a significant manner even with uncertainties taken into account. The magnetic field of upsized simulations (Fig. 13b) follows the f_Ω-corrected scaling |$B/\sqrt{f_{\Omega }\rho \mu \eta \Omega } \sim \epsilon ^{0.04\pm 0.06}$| and uncorrected scaling |$B/\sqrt{\rho \mu \eta \Omega } \sim \epsilon ^{0.02\pm 0.06}$|⁠, not fundamentally deviating from the constant prediction of the path theory (A17) but this time also marginally compatible with the QG-MAC prediction |$B/\sqrt{\rho \mu \eta \Omega } \sim \epsilon ^{1/12}{\sim \epsilon ^{0.083}}$| (eq. 17). Similar results are found for the magnetic dissipation length scale (Fig. 13c) that follows d_Ω/D ∼ ϵ^{0.04 ± 0.04}, marginally compatible with the constant prediction from the path theory and the QG-MAC prediction d_Ω/D ∼ ϵ^1/12 (eq. 20). In summary, output diagnostics from the upsized qDNS and DNS simulations indeed better approach the QG-MAC theory of Davidson (2013) than their LES counterparts, but fail to fully adhere to this theory, with the most important deviations being observed on the diagnostic for flow velocity.

Figure 13.

Evolution of the magnetic Reynolds number Rm (a), magnetic field amplitude |$B/\sqrt{\rho \mu \eta \Omega }$| corrected as in A17 with the inverse square root of the ohmic fraction |$f_\Omega$| (b), and magnetic dissipation length scale d_Ω (c) with the path parameter ϵ (Earth’s core conditions are towards the left on these graphs). Presented are the set of LES from A17 (black circles), the set of DNS from A17 (black stars) and the best-resolved qDNS and DNS from this work (blue diamonds), with the blue shaded regions representing the uncertainty range associated to short-term time averaging of these runs (see the text). Power laws obtained by least-squares fitting of the A17 results are also reported (orange and green lines for DNS and LES, respectively). The power-law fit obtained from DNS simulations (orange line) closely follows the predictions from the QG-MAC theory (Davidson 2013, A17). For the new DNS and qDNS simulations, best-fitting power laws (purple lines) are obtained by using weighted least squares taking the uncertainty range into account.

4 DISCUSSION

By increasing the resolution and relaxing the hyperdiffusive approximation of the large-eddy simulations presented in A17, we have presented high-resolution numerical simulations of the geodynamo that approach closer to the rapidly rotating and turbulent conditions of Earth’s core than earlier extreme simulations (Yadav et al. 2016; Schaeffer et al. 2017; Sheyko et al. 2018). The high-resolution models could be integrated only for a short time, but they are still relevant as they reach dynamical saturation after transients lasting only a fraction of an overturn time (Fig. 4). This rapid equilibration owes to the fact that the large-eddy simulations that they originate from are already in the correct large-scale force balance (Fig. 9). In that sense, large-eddy simulation followed by upsizing is an highly efficient strategy to reach extreme conditions at a fraction of the computer cost that would be involved in equilibrating the simulation from scratch. By reaching extremely low levels of the fluid viscosity in much of the harmonic degree range (up to five orders of magnitude below those of the leading forces, Fig. 2), the upsized simulations also reach a state where the injected convective power is almost entirely dissipated by Ohmic losses (Fig. 1). Combined with the previous results from A17, the new high-resolution simulations provide insight into force balance, scale separation and turbulence in Earth’s core.

4.1 QG-MAC force balance in Earth’s core

A quasi-geostrophic equilibrium between Coriolis and pressure forces (Fig. 2) has been confirmed as the leading-order force balance enforced at large scales as we progress towards Earth’s core conditions with high-resolution simulations. Though the magnetic force also reaches leading order at smaller scales (as previously advocated by Aurnou & King 2017), its contribution is essentially a magnetic pressure because of dynamic alignment of vorticity filaments and magnetic field lines (Fig. 8), implying that from a dynamical standpoint the whole length scale range in fact remains close to quasi-geostrophic. These results further rule out the possibility that Earth’s core is in a system-scale magnetostrophic state, and should encourage further theoretical developments centred on quasi-geostrophy (Calkins et al. 2015; Calkins 2018). At the next order, the magneto-Archimedes-Coriolis force balance already found in Yadav et al. (2016) and Aubert et al. (2017) has also been confirmed in high-resolution simulations.

4.2 Scale separation in Earth’s core

In the vicinity of the large scale ℓ_⊥ of the MAC balance, the upsized simulations preserve the large-scale features and axial invariance of the large-eddy simulations that they originate from. They however reveal an important scale separation in their velocity, vorticity and density anomaly fields (Figs 5 and 6) while the magnetic field remains at large scales (Fig. 7). The most interesting features revealed by upsizing are vorticity filaments that are most prominent at the small scale ℓ_Ω of magnetic dissipation (Fig. 12). These exist at the edges of the large-scale, sheet-like convective upwelling and induce ramifications of density anomaly plumes. The filaments align with the large-scale magnetic field lines (Fig. 8), thereby minimizing the non-pressure part of the associated Lorentz force at small scales (Fig. 2). The distribution of velocity and vorticity across scales has been analysed in the spectral space by using a decomposition of energy into spherical harmonic degrees. This approach is justified because at the strong forcing of our simulations, the effects of spatial heterogeneity and anisotropy are both reduced (Figs 5, 6 and 12) relative to simulations at weaker forcing (Schaeffer et al. 2017; Sheyko et al. 2018). The spectral analysis (Figs 11 and 12) has partly confirmed the principle of vorticity equivalence from the QG-MAC theory (Davidson 2013), with approximately flat enstrophy spectra between degree ℓ_⊥ (large scale d_⊥) at which the magnetic field draws its power from convection and degree ℓ_Ω (small scale d_Ω) where this power is dissipated. High-resolution simulations feature a ℓ^−3/2 power-law decay in the range [ℓ_⊥, ℓ_Ω], approaching closer to the spectral index −2 predicted from vorticity equivalence than earlier simulations at lower forcing (Schaeffer et al. 2017; Sheyko et al. 2018). On a side note, Schaeffer et al. (2017) mention that with a kinetic energy density power-law decay not much steeper than ℓ⁻¹, integral length scales built by weighting each harmonic degree ℓ with this density (Christensen & Aubert 2006; King & Buffett 2013) become ill-posed. The present high-resolution simulations confirm this view and suggest that such measures should be abandoned altogether as they do not match any of the physically relevant crossover length scales identified here.

From the evolution of key length scales along the parameter space path (Figs 3 and 13), at the end of the parameter space ϵ = 10⁻⁷ (i.e. at Earth’s core conditions) we predict large-scale features at scale ℓ_⊥ ≈ 10 and small-scale features at ℓ_Ω ≈ 300 (respectively |$d_{\perp } =\pi D / \ell _{\perp } \approx 700 \mathrm{km}$| and |$d_{\Omega }=\pi D / \ell _{\Omega } \approx 20 \mathrm{km}$|⁠), hence a significant but not extremely large level of scale separation. This implies that the large-scale structures currently imaged by core flow modelling (see Holme 2015 for a review) are representative of the large scales of the geodynamo, and that numerical dynamo simulations can provide robust a priori constraints on the spectral decay of the flow in this inverse problem (e.g. Aubert 2013, 2014, 2015; Fournier et al. 2015).

4.3 Optimal strategies for large-scale approximation of Earth’s core dynamics

In our system, the convective power p is essentially carried by magnetic nonlinearities associated to the Lorentz force from the injection scale ℓ_⊥ to the Ohmic dissipation scale ℓ_Ω (Section 2.2). To further render the transfer of the residual power (1 − f_Ω)p to even smaller viscous dissipation scales through hydrodynamic nonlinearities, one needs to resolve at least the Coriolis-Inertia-Archimedes scale ℓ_CIA (Figs 1 and 9). From the simulation cases where ℓ_CIA is correctly resolved (Fig. 3) we can infer the dependency ℓ_CIA ∼ ϵ^−0.4, leading at Earth’s core conditions to ℓ_CIA ≈ 20 000, or the equivalent length scale |$\pi D/ \ell _\mathrm{CIA} \approx 400 \mathrm{m}$|⁠. Though the need to resolve such a small scale obviously complexifies the numerical simulation problem, we note that this residual part of power (1 − f_Ω)p becomes vanishingly small in our simulations and therefore at Earth’s core conditions (Fig. 1), thereby motivating further research towards approximated simulations of moderate maximal resolution ℓ_max ≈ ℓ_Ω ≪ ℓ_CIA. Introducing quasi-direct numerical simulations where the hypderdiffusive cut-off ℓ_h has been set to a value larger than ℓ_Ω, but smaller than ℓ_CIA, we have for instance obtained an acceptable trade-off between computational cost and accuracy, with an accurate value of the Ohmic dissipation fraction f_Ω (Fig. 1) but some stacking towards the kinetic energy spectrum tail (Figs 10 and 11) and some degradation of the small-scale force balance (Fig. 2). The upsized results have also confirmed the quality of earlier large-eddy simulations from A17, which capture the large-scale planforms (Figs 5–7), the leading order and first-order force balances (Fig. 9) and dynamical equilibria such as the enforcement of the Taylor constraint (Table 2) to an excellent level of accuracy, but fail to render the partition of dissipation between Ohmic and viscous losses (Fig. 1) because they introduce a viscous-Archimedes-Coriolis force balance at second order.

4.4 Rotating magnetohydrodynamic turbulence in the spherical shell geometry

It is interesting to compare the turbulence properties observed here to the non-magnetic, but rotating case. In classical 2-D and quasi-geostrophic turbulence (e.g. Nataf & Schaeffer 2015), the enstrophy cascades to smaller scales through vortex stretching. In the present high-resolution simulations, vortex stretching is strongly constrained by dynamic alignment of vortex filaments with magnetic field lines (Fig. 8). This may explain why the velocity spectrum decay ℓ^−3/2 (Figs 11 and 12) is less steep than the decay ℓ^−5/3 previously reported in a strongly forced simulation of nonmagnetic and rotating convection (Sheyko et al. 2018), and much less steep than the decay m⁻⁵ obtained along spherical harmonic orders in simulations of quasi-geostrophic Rossby wave turbulence (Schaeffer & Cardin 2005).

Examining our results within the framework of non-rotating, MHD turbulence, we first note that our simulations are in the strong turbulence regime where magnetic field lines are bent by velocity fluctuations. Cartesian simulations of strong MHD turbulence commonly feature energy density spectra decaying like |$k_{\perp }^{-3/2}$|⁠, where k_⊥ is the field-perpendicular wavenumber, and scale-dependent dynamic alignment of vorticity structures along field lines is essential to reconcile this spectral decay with existing theories (see e.g. Tobias et al. 2012). It is interesting to note that dynamic alignment in the physical space is also clearly observed here, and that Alfvén waves, the energy carriers of strong MHD turbulence, have been highlighted as important energy carriers in our simulations, both the axisymmetric and non-axisymmetric levels (Aubert 2018). Interpreting the common spectral decay index −3/2 of our results and the strong MHD turbulence theory is however subject to caution, as the relationship between the cartesian wavenumber k_⊥ and the spherical harmonic degree ℓ in volume-averaged spectra is not straightforward (Schaeffer et al. 2017).

4.5 Towards an ultimate asymptotic scaling theory

The outputs from the best-resolved upsized simulations approach the scaling predictions of the QG-MAC theory (Davidson 2013) better than their large-eddy counterparts, but do not fully adhere to this theory. The deviations are clearest for the length scale d_⊥ (Fig. 3) and the flow velocity or magnetic Reynolds number Rm (Fig. 13). They already appear in a DNS carried out at 21 per cent of the path, suggesting that this is not the result of residual hyperdiffusivity effects. It remains possible that the magnetic equilibration of the upsized simulations takes place on a timescale that is still longer than the short integration times that could be achieved here, in which case the error bars specified of Fig. 13 would be underestimated. If such is the case however, the DNS at 21 per cent of the path remains puzzling as its magnetic field amplitude appears to almost match the QG-MAC prediction while the velocity field amplitude strongly deviates from this theory. Deviations from the QG-MAC theory may be needed somehow, because we have seen that this theory fails to predict an increase of scale separation along the path as Earth’s core conditions are approached, in apparent contradiction with the present results (Fig. 11) and previous high-resolution numerical simulations (Sakuraba & Roberts 2009; Yadav et al. 2016; Schaeffer et al. 2017; Sheyko et al. 2018). It has also been shown in A17 that the predictions from the scale-invariant path theory match the amplitude of Earth’s core magnetic and velocity fields strikingly well, and better than the QG-MAC predictions. In this respect, it is possible that the vorticity equivalence (spectral index −2 in the kinetic energy spectrum) is incompatible with dynamic alignment of vorticity and magnetic field lines (spectral index −3/2). When stacked at advancing positions along the parameter space path (Fig. 11), kinetic energy spectra do indeed appear to converge towards index −3/2 and there is no evidence that the index could evolve towards −2 upon further progress along the path. If index −3/2 is the asymptotic value, then the vorticity equivalence constant c_ω in eq. (13) may have some residual dependence on the path parameter ϵ, while still remaining of order one at the end of path. This indicates that some refinements may still be needed for an ultimate scaling theory of Earth’s core dynamo. Further work should account for the gradual enforcement of structuring constraints exerted by the magnetic field on the flow as we progress towards Earth’s core conditions.

ACKNOWLEDGEMENTS

JA wishes to thank Thomas Gastine for insightful discussions and suggestions, and acknowledges support from the Foundation Simone et Cino Del Duca of Institut de France (2017 Research Grant). Numerical computations were performed at S-CAPAD, IPGP and using high-performance computing resources from Grand Equipement National de Calcul Intensitf (GENCI)-IDRIS, GENCI-TGCC and GENCI-CINES (Grants A0020402122 and A0040402122). This is IPGP contribution 4022.