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

TL;DR

This study uses high-resolution simulations to explore Earth's core conditions, revealing how scale separation influences magnetic and flow structures, and improving the alignment with theoretical models of the geodynamo.

Contribution

It advances geodynamo modeling by increasing resolution and relaxing hyperdiffusive approximations, approaching more realistic Earth's core conditions and validating quasi-geostrophic and MAC force balances.

Findings

Kinetic energy spectra broaden with increased resolution.

Large-scale magnetic structures remain stable after upsizing.

Simulations align better with QG-MAC theory, with some residual deviations.

Abstract

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 Earth's core conditions are approached in the upsized simulations, 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 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.

Tables

Table 1. Table 1: Input parameters of numerical models (see text for definitions). Cases marked with a (*) symbol have been previously reported in A17.

Case	$ϵ$	Path position (percent)	run length (overturns)	$E$	$Pm$	${\ell }_{\max }$	NR	${\ell }_h$	$q_h$
DNS (*)	0.1	14	164	$3{10}^{-6}$	0.8	256	480	-	-
DNS	$3.33{10}^{-2}$	21	59	${10}^{-6}$	0.45	320	720	-	-
qDNS	$3.33{10}^{-2}$	21	8	${10}^{-6}$	0.45	320	720	290	1.1
qDNS1	${10}^{-2}$	29	33	$3{10}^{-7}$	0.25	256	816	196	1.15
qDNS2	${10}^{-2}$	29	3.3	$3{10}^{-7}$	0.25	512	1248	448	1.15
qDNS	$3.33{10}^{-4}$	50	0.8	${10}^{-8}$	0.045	640	2496	512	1.45
LES (*)	0.1	14	1172	$3{10}^{-6}$	0.8	133	200	30	1.045
LES (*)	$3.33{10}^{-2}$	21	595	${10}^{-6}$	0.45	133	240	30	1.0575
LES1 (*)	${10}^{-2}$	29	332	$3{10}^{-7}$	0.25	133	320	30	1.07
LES2	${10}^{-2}$	29	22.3	$3{10}^{-7}$	0.25	256	320	30	1.03
LES (*)	$3.33{10}^{-4}$	50	196	${10}^{-8}$	0.045	133	624	30	1.10

Table 2. Table 2: Output parameters of numerical models (see text for definitions). Cases marked with a (*) symbol have been previously reported in A17.

Case	Path position (percent)	$\frac{4\pi D^{2}p}{g_{o}F}$	$f_{\Omega }$	$Rm=\frac{UD}{\eta }$	$\frac{B}{\sqrt{\rho \mu \eta \Omega }}$	$𝒯$	${\ell }_{\Omega }$	${\ell }_{\mathrm{MS}}$	${\ell }_{⟂}$	${\ell }_{\mathrm{CIA}}$	${\ell }_{\mathrm{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

Equations

$$R{a}_F$$

$$E$$

$$Pr$$

$$Pm$$

$$R{a}_F$$

$$E$$

$$Pr$$

$$Pm$$

$$({\nu }_{\text{eff}},{\kappa }_{\text{eff}})$$

$$\frac{\rho \Omega U}{D}\sim \frac{g_{o}C}{d_{\perp }}\sim \frac{B^2}{\mu d_{\perp }^2},$$

$$\frac{\eta B^2}{\mu d_{\Omega }^2}\sim p\sim \frac{g_{o}F}{D^2},$$

$$U/d_{\perp }=c_{\omega }\eta /d_{\Omega }^2.$$

$${\omega }^2(\ell )\propto {\ell }^0\mathrm{for}{\ell }_{\perp }<\ell <{\ell }_{\Omega }.$$

$$u^2(\ell )\propto {\ell }^{-2}\mathrm{for}{\ell }_{\perp }<\ell <{\ell }_{\Omega }.$$

$$B\sim \sqrt{\rho \mu }(g_{o}F/\rho D{)}^{1/3}.$$

$$\frac{B}{\sqrt{\rho \mu \eta \Omega }}\sim {ϵ}^{1/12},$$

$$Rm=\frac{UD}{\eta }\sim {ϵ}^{-1/18},$$

$$d_{\perp }/D\sim {ϵ}^{1/9},$$

$$d_{\Omega }/D\sim {ϵ}^{1/12}.$$

Full text

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 non-linear 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 westward 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 characterised 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 time scales 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 rationalises 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 percent) of this parameter space path has been covered by fully resolved direct numerical simulations (DNS), and a larger portion (up to 50 percent) 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 parameterisations (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 non-linear 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-fertilisation 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 minimised 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 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 the present study, the hyperdiffusive approximation of LES simulations previously published in Aubert et al. (2017); 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, characterise the turbulence and scale separation through their statistical properties and assess their level in Earth’s core. The manuscript is organised 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 magnetohydrodynamic 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

[TABLE]

Here $g_{o}$, $\rho$, $\nu$, $\kappa$ and $\eta$ are respectively the gravity at radius $r_{o}$, the fluid density, viscosity, thermo-chemical 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 $\epsilon$:

[TABLE]

Here $Ra_{F}(CE)=2.7~{}10^{-5}$, $E(CE)=3~{}10^{-5}$ 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 $\epsilon=1$, and ending at $\epsilon=10^{-7}$ where reasonable estimates for Earth’s core conditions are obtained (A17). The model cases presented here (Table 1) range from $\epsilon=0.1$ to $\epsilon=3.33~{}10^{-4}$, 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 $\ell_{\mathrm{max}}$ and a discretisation in the radial direction on a second-order finite-differencing scheme with $NR$ grid points. 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.

Along the parameter space path, the magnetic Reynolds number $Rm=UD/\eta$ (where $U$ is the r.m.s 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/\nu=Rm/Pm$ reaching values up to $Re\approx 24000$. A17 explored a parameter range down to $\epsilon=0.1$ (the first 14 percent 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 which remained fully resolved. The principle is to use values of $\ell_{\mathrm{max}}$ and $NR$ lower than those required by DNS, and effective diffusivities $\nu_{\text{eff}},\kappa_{\text{eff}}$ that depend on the spherical harmonic degree $\ell$ and the molecular diffusivities $\nu,\kappa$ according to (Nataf & Schaeffer, 2015)

[TABLE]

Here $\ell_{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},\ell_{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 $\ell_{h}=30$ and smoothly ramp its strength up towards smaller scales. Here we increase the model resolution $(\ell_{\mathrm{max}},NR)$ and test three approaches to relax hyperdiffusivity. First, we introduce new LES simulations where we keep $\ell_{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 $\ell_{h}$ close to the smallest numerically resolved scale $\ell_{\mathrm{max}}$ and employ larger values of $q_{h}$ than those used in A17. This latter configuration attempts to maximise 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 percent, by reporting on a DNS at $E=10^{-6}$ ($\epsilon=3.33~{}10^{-2}$). 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$).

Time-averaged diagnostic outputs of the simulation are presented in Table 2. Listed first are the convective power density $p$ in units of $g_{o}F/4\pi D^{2}$, 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 r.m.s. velocity $U$ in units of $\eta/D$, the shell r.m.s magnetic field $B$ in Elsasser units of $\sqrt{\rho\mu\eta\Omega}$ (where $\mu$ 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_{\mathrm{min}}$) from the square root of the ratio of volume-averaged magnetic energy and magnetic dissipation, and reported as an equivalent harmonic degree $\ell_{\Omega}=\pi D/d_{\Omega}$.

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 $\ell$ over the harmonic order because of the invariance of $\ell$-dependent energy spectra against rotation of the coordinate system, and because unlike an harmonic order, a given degree $\ell$ in the spectral space can be associated to a minimal corresponding length scale $\pi D/\ell$ in the physical space. It has been 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 $\ell$ to its action at length scale $\pi D/\ell$ 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 percent. 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., submitted), 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 percent 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 $\ell_{\mathrm{MS}}$. For harmonic degrees larger than $\ell_{\mathrm{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 $\ell_{\perp}$. 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 $\ell<\ell_{\perp}$, the first order balance is of thermal wind type between the ageostrophic Coriolis and buoyancy forces, and for $\ell>\ell_{\perp}$ 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 $\ell_{\mathrm{CIA}}$, and of the viscous and buoyancy forces $\ell_{\mathrm{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 on Fig. 3. Scale separation in this system refers to the difference between the large scale of the first-order MAC balance at degree $\ell_{\perp}$, and the much smaller scales of the second-order force balances at degrees $\ell_{\mathrm{CIA}}$, $\ell_{\mathrm{VAC}}$ and of magnetic dissipation at degree $\ell_{\Omega}$.

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{\!/\mkern-4.0mu/\!}}$ of flow structures along the rotation axis $\boldsymbol{\mathrm{\Omega}}$ remains at the system size i. e. $d_{\mathbin{\!/\mkern-4.0mu/\!}}\sim D$ (see Fig. 5). The first-order MAC balance operating at the large scale $d_{\perp}=\pi D/\ell_{\perp}$ can be written from the curl of the Navier-Stokes equation (see A17):

[TABLE]

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

[TABLE]

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\pi D^{2}p/g_{o}F$ in Table 2. Combining equations (11) and (12), and using $F\sim UCD^{2}$, we obtain the vorticity equivalence

[TABLE]

This equivalence states that the vorticities at the large scale $d_{\perp}$ 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_{\omega}$ (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) non-linearities from the large injection scale $d_{\perp}$ 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 24000 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_{\omega}$ being in principle of order 1, at large values $d_{\perp}/d_{\Omega}$ of the scale separation we expect the spectral enstrophy density $\omega^{2}(\ell)$ at degree $\ell$ (where $\boldsymbol{\mathrm{\omega}}=\nabla\times\boldsymbol{\mathrm{u}}$ is vorticity) to remain approximately constant i.e.

[TABLE]

With length scales varying like $\ell^{-1}$, the energy density of velocity $u^{2}(\ell)$ should then also present a power-law range

[TABLE]

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{\!/\mkern-4.0mu/\!}}$), implying that the velocity energy spectrum is less steep than $\ell^{-2}$ if the vorticity spectrum is flat. In that sense, the spectra of the vorticity component parallel to the rotation axis $\omega_{\mathbin{\!/\mkern-4.0mu/\!}}=\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_{\perp}$, and $d_{\Omega}$ along the parameter space path, as functions of the path parameter $\epsilon$. If one requires the large- and small-scale vorticities present in equation (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):

[TABLE]

Combining equations (11,12,16) and using $Ra_{F}\sim\epsilon$ 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:

[TABLE]

The powers of $\epsilon$ entering (17-20) imply that these quantities only undergo weak variations along the path (as confirmed by Table 2), from $\epsilon=1$ down to $\epsilon=10^{-7}$. 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_{\perp}$ and $d_{\Omega}$ (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 $\epsilon$ 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 percent of the parameter space path (Ekman number $E=3~{}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 (Fig. 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_{\Omega}$ (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 $\ell_{\Omega}=\pi/d_{\Omega}$ (see Fig. 3). As a corollary, the fraction $f_{\Omega}$ of Ohmic dissipation in the system (Fig. 1) generally increases during the upsizing procedure, and reaches a value independent on further upsizing when $d_{\Omega}$ also converges (see similarity of $f_{\Omega}$, $\ell_{\Omega}$ for qDNS1 and qDNS2 in Table 2). At 50 percent of the path, the upsizing process thereby leads to a value $f_{\Omega}=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.

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 percent 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 (Fig. 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 (Fig. 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{\!/\mkern-4.0mu/\!}}$ (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,11,12, the large scale in $\boldsymbol{\mathrm{u}}$, $\boldsymbol{\mathrm{B}}$, $C$ is $d_{\perp}$, 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_{\Omega}$ where this power is ohmically dissipated. In the LES1 case (Figs. 5a,c,d, 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 $\ell_{\mathrm{VAC}}$). Fig. 7a,b shows 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_{\perp}$. Indeed, in case 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 magnetohydrodynamic turbulence (see Tobias et al., 2012). This is the process through which the magnetic field adjusts so as to minimise 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. 6b, 7b to Fig. 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); Sheyko et al. (2018) and this presumably improves the exchange of material through the tangent cylinder, leading to a spatially more homogeneous system.

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 percent). 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 Fig. 9a and 9c). The value of the crossover length scale to magnetostrophy $\ell_{\mathrm{MS}}$ is also broadly invariant in our sets of runs, except for the qDNS at 50 percent of the parameter space path (Fig. 2c), where the force balance is perturbed by the choice of upsizing strategy (see below). Invariance of $\ell_{\mathrm{MS}}$ is consistent with the analysis of Aurnou & King (2017) where this scale is conjectured to depend solely on the Elsasser number $B^{2}/\rho\mu\eta\Omega$ and magnetic Reynolds number $Rm$, both of which are about constant along the path (Table 2). The MAC length scale $\ell_{\perp}$ 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.

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 (Fig. 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, 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 $\ell_{h}$ largely exceeds $\ell_{\mathrm{CIA}}$ (such as qDNS at 21 percent and qDNS2 at 29 percent 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 Fig. 9b and 9c, see also Fig. 2b). At third order, the VAC balance length scale $\ell_{\mathrm{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_{\Omega}$ values in qDNS and DNS at 21 percent 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 $\ell_{\mathrm{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 $\ell_{\mathrm{CIA}}$, but still larger than $\ell_{\Omega}$ may be performed (respectively, cases qDNS1 at 29 percent of the path and case qDNS at 50 percent of the path, Fig. 3). This type of qDNS still captures the Ohmic dissipation fraction $f_{\Omega}$ well, and in any case better than LES (compare LES, qDNS1 and qDNS2 at 29 percent 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 $\ell_{\mathrm{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 $\ell_{\mathrm{CIA}}\approx\ell_{\mathrm{VAC}}$ and both values being significantly decreased relative to DNS (Table 2, Fig. 3). At path positions of 29 percent and beyond the ordering of these two length scales is even reversed compared to the DNS ordering (i.e. $\ell_{\mathrm{VAC}}<\ell_{\mathrm{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_{\Omega}$ 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 $\ell$ 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 $\ell\approx 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 $\ell>1$ occurs at $\ell=\ell_{\perp}$, confirming that $\ell_{\perp}$ 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\approx 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 $\ell_{h}=30$ (red and pink curves in Fig. 10a,b), but it appears gradually in qDNS and DNS simulations where the hyperdiffusive cut-off $\ell_{h}$ is increased close to $\ell_{\mathrm{max}}$ or removed altogether. It is again advisable to use a sufficiently resolved qDNS simulation (dark green curves in Fig. 10a,b) for best accuracy (compare the dark green and black curves in Fig. 10a). In particular, choosing $\ell_{h}$ above $\ell_{\Omega}$, but close to or below $\ell_{\mathrm{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).

Plotting together uncompensated and compensated kinetic energy spectra (Fig. 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 $\ell^{-3/2}$ starting from the large scale $\ell_{\perp}$. 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. 11b fall below 85 percent 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 percent 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 $\ell_{\Omega}$ which in this case is significantly lower than $\ell_{\mathrm{CIA}}$ (Figs. 2c, 3). Power transfers from injection at $\ell_{\perp}$ down to Ohmic dissipation at $\ell_{\Omega}$ are well controlled in this simulation (hence the flat compensated spectral range), but the transfer to even smaller scales of the residual power (3 percent in this case, Fig. 1) for viscous dissipation via inertial forces is improperly rendered (because the true value of $\ell_{\mathrm{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 Fig. 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.

In order to further establish the link between Davidson’s vorticity equivalence (equation 13) and the power-law range observed in Figs. 10,11, spatial spectra of enstrophy $\omega^{2}$ are plotted together with kinetic energy spectra in Fig. 12 for case qDNS2 at 29 percent of the parameter space path. At scales between the MAC scale $\ell_{\perp}$ and the magnetic dissipation scale $\ell_{\Omega}$, the enstrophy $\omega^{2}$ presents a weakly increasing power-law range $\omega^{2}\propto\ell^{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 (equation 14). The vorticity spectrum peaks at $\ell=\ell_{\Omega}$, showing that the dominant scale of the vorticity planforms observed in Fig. 5f is indeed the scale of magnetic dissipation (the same analysis carried out in the case of the LES case of Fig. 5e reveals that $\ell_{\mathrm{VAC}}$ is the dominant scale). Note that a length scale defined as $\delta(\ell)=u(\ell)/\omega(\ell)\sim D\ell^{-11/12}$ does not exactly vary like $\ell^{-1}$, underlining weak but measurable effects of the anisotropic flow structure (Fig. 5b,d). Computing the spectra of axial enstrophy $\omega_{\mathbin{\!/\mkern-4.0mu/\!}}^{2}$ removes the influence of the long length scale parallel to the rotation axis, and reveals a better respected vorticity equivalence $\omega_{\mathbin{\!/\mkern-4.0mu/\!}}^{2}\sim\ell^{1/4}$ in the range $[\ell_{\perp},\ell_{\Omega}]$. 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$ (equation 15), and the associated perpendicular length scale $\delta_{\perp}(\ell)=u_{\perp}(\ell)/\omega_{\mathbin{\!/\mkern-4.0mu/\!}}(\ell)$ now indeed varies like $\ell^{-1}$, confirming the relevance of spherical harmonic degrees to analyse length scales.

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 $\ell_{\perp}$ 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 $\ell_{\perp}$ on snapshots and the associated uncertainty, which while being small (about 10 percent) remains on the order of the observed variations. Furthermore, from equation (19) a doubling of $\ell_{\perp}=\pi/d_{\perp}$ is predicted between 14 percent and 50 percent of the path, but clearly not observed. The departure of our results for $\ell_{\perp}$ from the QG-MAC prediction simply reflects the fact that the scale separation $d_{\perp}/d_{\Omega}$ in our simulations actually increases along the path, while from equations (19,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_{\Omega}/D$. To evaluate the uncertainties that are due to the shortness of DNS and qDNS runs, our longest LES simulation at 14 percent 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\sim\epsilon^{0.004\pm 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\sim\epsilon^{-1/18}{\sim\epsilon^{-0.056}}$ (equation 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_{\Omega}$-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}}$ (equation 17). Similar results are found for the magnetic dissipation length scale (Fig. 13c) that follows $d_{\Omega}/D\sim\epsilon^{0.04\pm 0.04}$, marginally compatible with the constant prediction from the path theory and the QG-MAC prediction $d_{\Omega}/D\sim\epsilon^{1/12}$ (equation 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.

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 centered 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); 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 $\ell_{\perp}$ 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,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 $\ell_{\Omega}$ 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 minimising 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,12) relative to simulations at weaker forcing (Schaeffer et al., 2017; Sheyko et al., 2018). The spectral analysis (Figs. 11,12) has partly confirmed the principle of vorticity equivalence from the QG-MAC theory (Davidson, 2013), with approximately flat enstrophy spectra between degree $\ell_{\perp}$ (large scale $d_{\perp}$) at which the magnetic field draws its power from convection and degree $\ell_{\Omega}$ (small scale $d_{\Omega}$) where this power is dissipated. High-resolution simulations feature a $\ell^{-3/2}$ power-law decay in the range $[\ell_{\perp},\ell_{\Omega}]$, 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 $\ell^{-1}$, integral length scales built by weighting each harmonic degree $\ell$ 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,13), at the end of the parameter space $\epsilon=10^{-7}$ (i.e. at Earth’s core conditions) we predict large-scale features at scale $\ell_{\perp}\approx 10$ and small scale features at $\ell_{\Omega}\approx 300$ (respectively $d_{\perp}=\pi D/\ell_{\perp}\approx 700\hskip 2.84526pt\mathrm{km}$ and $d_{\Omega}=\pi D/\ell_{\Omega}\approx 20\hskip 2.84526pt\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 $\ell_{\perp}$ to the Ohmic dissipation scale $\ell_{\Omega}$ (section 2.2). To further render the transfer of the residual power $(1-f_{\Omega})p$ to even smaller viscous dissipation scales through hydrodynamic nonlinearities, one needs to resolve at least the Coriolis-Inertia-Archimedes scale $\ell_{\mathrm{CIA}}$ (Figs. 1, 9). From the simulation cases where $\ell_{\mathrm{CIA}}$ is correctly resolved (Fig. 3) we can infer the dependency $\ell_{\mathrm{CIA}}\sim\epsilon^{-0.4}$, leading at Earth’s core conditions to $\ell_{\mathrm{CIA}}\approx 20000$, or the equivalent length scale $\pi D/\ell_{\mathrm{CIA}}\approx 400\hskip 2.84526pt\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_{\Omega})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 $\ell_{\mathrm{max}}\approx\ell_{\Omega}\ll\ell_{\mathrm{CIA}}$. Introducing quasi-direct numerical simulations where the hypderdiffusive cut-off $\ell_{h}$ has been set to a value larger than $\ell_{\Omega}$, but smaller than $\ell_{\mathrm{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_{\Omega}$ (Fig. 1) but some stacking towards the kinetic energy spectrum tail (Figs. 10,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 two-dimensional 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 $\ell^{-3/2}$ (Figs. 11, 12) is less steep than the decay $\ell^{-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^{-5}$ 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, magnetohydrodynamic 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_{\perp}$ 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 magnetohydrodynamic 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_{\perp}$ and the spherical harmonic degree $\ell$ 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_{\perp}$ (Fig. 3) and the flow velocity or magnetic Reynolds number $Rm$ (Fig. 13). They already appear in a DNS carried out at 21 percent 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 time scale 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 percent 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_{\omega}$ in equation (13) may have some residual dependence on the path parameter $\epsilon$, 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 Fondation Simone et Cino Del Duca of Institut de France (2017 Research Grant). Numerical computations were performed at S-CAPAD, IPGP and using HPC resources from GENCI-IDRIS, GENCI-TGCC and GENCI-CINES (Grants A0020402122 and A0040402122). This is IPGP contribution 4022.