- Split View
-
Views
-
Cite
Cite
N. Schaeffer, D. Jault, H.-C. Nataf, A. Fournier, Turbulent geodynamo simulations: a leap towards Earth’s core, Geophysical Journal International, Volume 211, Issue 1, October 2017, Pages 1–29, https://doi.org/10.1093/gji/ggx265
- Share Icon Share
Summary
We present an attempt to reach realistic turbulent regime in direct numerical simulations of the geodynamo. We rely on a sequence of three convection-driven simulations in a rapidly rotating spherical shell. The most extreme case reaches towards the Earth’s core regime by lowering viscosity (magnetic Prandtl number Pm = 0.1) while maintaining vigorous convection (magnetic Reynolds number Rm > 500) and rapid rotation (Ekman number E = 10−7) at the limit of what is feasible on today’s supercomputers. A detailed and comprehensive analysis highlights several key features matching geomagnetic observations or dynamo theory predictions—all present together in the same simulation—but it also unveils interesting insights relevant for Earth’s core dynamics. In this strong-field, dipole-dominated dynamo simulation, the magnetic energy is one order of magnitude larger than the kinetic energy. The spatial distribution of magnetic intensity is highly heterogeneous, and a stark dynamical contrast exists between the interior and the exterior of the tangent cylinder (the cylinder parallel to the axis of rotation that circumscribes the inner core). In the interior, the magnetic field is strongest, and is associated with a vigorous twisted polar vortex, whose dynamics may occasionally lead to the formation of a reverse polar flux patch at the surface of the shell. Furthermore, the strong magnetic field also allows accumulation of light material within the tangent cylinder, leading to stable stratification there. Torsional Alfvén waves are frequently triggered in the vicinity of the tangent cylinder and propagate towards the equator. Outside the tangent cylinder, the magnetic field inhibits the growth of zonal winds and the kinetic energy is mostly non-zonal. Spatio-temporal analysis indicates that the low-frequency, non-zonal flow is quite geostrophic (columnar) and predominantly large-scale: an m = 1 eddy spontaneously emerges in our most extreme simulations, without any heterogeneous boundary forcing. Our spatio-temporal analysis further reveals that (i) the low-frequency, large-scale flow is governed by a balance between Coriolis and buoyancy forces—magnetic field and flow tend to align, minimizing the Lorentz force; (ii) the high-frequency flow obeys a balance between magnetic and Coriolis forces; (iii) the convective plumes mostly live at an intermediate scale, whose dynamics is driven by a three-term MAC balance—involving Coriolis, Lorentz and buoyancy forces. However, small-scale (≃E1/3) quasi-geostrophic convection is still observed in the regions of low magnetic intensity.
1 INTRODUCTION
Earth’s magnetic field is generated by a turbulent flow of liquid metal in the core. The pioneering work of Glatzmaier & Roberts (1995) launched the first generation of numerical simulations of the geodynamo. For the first time, a self-sustained magnetic field was produced by a self-consistent convective geodynamo model. It exhibited several key features of the Earth’s magnetic field: a mostly dipolar field aligned with the rotation axis, and polarity reversals. These findings were somewhat surprising, considering the huge gap between the parameters used in the simulations and their expected values in the Earth’s core.
Second-generation numerical simulations of the geodynamo explored the parameter space and derived scaling laws (Christensen & Aubert 2006). As the mean properties of the simulations could be cast in laws that excluded diffusion properties (viscosity, thermal and magnetic diffusivities), it was argued that the simulations had reached an asymptotic regime, which allowed extrapolation to the Earth (Christensen et al. 2010), planets (Christensen 2010), and some stars (Christensen et al. 2009).
Rapidly, several teams questioned that these simulations really reached an asymptotic regime, and that it was the right asymptotic regime for the Earth. Indeed, Soderlund et al. (2012) highlighted the important role of viscosity in their simulation set, and that magnetic forces did not seem to play a major role, contrary to what is expected for the Earth’s core. Reanalysing the large suite of simulations of Christensen & Aubert (2006) and King & Buffett (2013) confirmed that viscous dissipation was far from negligible in these simulations; Cheng & Aurnou (2016) showed that the diffusionless scaling laws were hiding an actual dependency upon viscosity; Oruba & Dormy (2014) derived alternative scaling laws in which the magnetic field intensity depends upon viscosity and rotation rate.
In the meantime, analyses of geomagnetic observations suggested a quasi-geostrophic (invariant along the rotation axis) description of the flow (Gillet et al. 2011). Using this assumption, the large-scale flow at the top of the core inferred from geomagnetic observations can be continued throughout the core, leading to a root-mean-square velocity of about 10 km yr−1, that is |$\bar{U}_{\text{rms}} \simeq 3 \times 10^{-4}$| m s−1, encompassing large-scale flow, and consistent with previous estimates of the flow speed at the top of the core. In addition, several remarkable features were found (Gillet et al. 2015; Pais et al. 2015): the analysis of Pais & Jault (2008) yielded a large-scale off-centred quasi-geostrophic gyre around the inner core (also see: Schaeffer & Pais 2011); Finlay & Jackson (2003) showed evidence for equatorial magnetic waves propagating at decadal timescales; Gillet et al. (2010) detected torsional oscillations and, from their propagation velocity, inferred the magnetic field inside the core. Indeed, the propagation speed of torsional waves scale with the intensity Brms of the magnetic field perpendicular to the rotation axis. Assuming isotropy, Gillet et al. (2010, 2015) found Brms = 4 mT, a value that includes the contribution of all length-scales.
These discoveries prompted modellers to look for these signatures in their simulations, but also highlighted that magnetic energy Em in the Earth’s core is much larger than kinetic energy Ek. Indeed, estimating the former using Brms and the latter assuming only the large-scale flow |$\bar{U}_{\rm {rms}}$| leads to Em/Ek ≃ 104. Even allowing for hidden and energetic small scales of about |$U_{\text{rms}} \simeq 3\,\bar{U}_{\text{rms}}$| leads to Em/Ek ≃ 1000.
This strong magnetic field revived the search for inviscid magnetostrophic dynamos—where Pressure, Coriolis and Lorentz forces form the main balance (Livermore et al. 2013; Wu & Roberts 2015). Taking an opposite approach, Dormy (2016) obtained strong-field dynamos for large magnetic Prandtl numbers Pm—the ratio of kinematic viscosity to magnetic diffusivity. Increased computational power also made it possible to push the parameters towards more realistic values (Kageyama et al. 2008; Sakuraba & Roberts 2009; Sheyko 2014; Yadav et al. 2016b). Fig. 1 displays the current parameter-space coverage that has been achieved by several groups, in terms of Ekman number E (ratio of viscous to Coriolis forces), magnetic Prandtl number Pm, magnetic Reynolds number Rm (ratio of induction over ohmic dissipation), and Alfvén number A (ratio of fluid velocity over Alfvén wave velocity, or equivalently the square root of kinetic over magnetic energies). It shows that having a large Rm is difficult to achieve at low viscosity (low E and low Pm, as in the Earth’s core), because the flow becomes turbulent, implying that smaller and smaller scales need to be resolved, leading to tremendous computing costs. Similarly, few simulations lie in the low Pm, low A region where the Earth’s core is found. In contrast to direct simulations that are committed to resolve all dynamical scales, Aubert et al. (2017) have used a form of eddy viscosity to suppress the smallest scales while keeping the largest ones unaffected. This allowed them to reach very low values of the large-scale viscosity and argue for a continuous path connecting today’s simulations with the Earth’s core.
Here, we present an attempt to reach turbulent regime in direct numerical simulations of the geodynamo. By not artificially damping the smallest flow scales, we ensure unbiased dynamics at all scales. We have been able to reach strongly forced dynamos at low viscosity, where the magnetic energy becomes one order of magnitude larger than the kinetic energy. In this way, all obvious timescales have the same ordering than what is inferred from observations. From the shortest to the longest time: the rotation period 2π/Ω, the Alfvén time |$\tau _b=D\sqrt{\mu _0\rho }/B$| (typical time for an Alfvén wave to propagate across the core), the turnover time τu = D/U (time for a fluid parcel to travel across the core), the magnetic diffusion time (D2/η), and the viscous damping time (D2/ν). Fast and slow variations of the magnetic field can be observed, while the flow exhibits tall and thin structures under the effect of the strong global rotation, with a wide range of excited scales. We also observe torsional waves, which appear to be forced near the tangent cylinder. Using intensive post-processing, we analyse the fields and forces as a function of time- and length-scale in the two dynamically distinct regions separated by the imaginary cylinder aligned with rotation axis and tangent to the inner-core.
The rest of this paper is organized as follows. In the next section, we describe our numerical model, as well as the simulation targets. Then we describe and analyse our set of simulations, focusing on both averaged fields and rapid dynamics. We conclude the paper with a discussion.
2 MODEL EQUATIONS, TARGET, AND NUMERICAL METHODS
2.1 Model
In order to simulate the liquid core of the Earth, confined between the solid inner core and mantle, we consider a spherical shell with an inner boundary at radius r = ri and an outer boundary at r = ro. The aspect ratio is fixed to ri/ro = 0.35, as for the present Earth. The solid boundaries are at rest in the reference frame rotating at a constant rate Ω along the z-axis. The fluid has an electrical conductivity σ and a kinematic viscosity ν, both being constants (in space and time). Its magnetic permeability μ0 is that of empty space.
The conductive codensity profile C0(r) is obtained using the thermochemical model of Aubert et al. (2009), as detailed in the Appendix.
At both inner and outer boundaries, the codensity deviation from C0 is imposed to have a zero gradient ∂rC = 0 (effectively fixing the constant uniform heat flux set by C0), the velocity is zero (no-slip), and the magnetic field is matched to a potential field outside the fluid domain (the inner-core and mantle are both modelled as electrically insulating). Note that fixing the codensity flux also enforces uniform crystallization of the inner-core.
2.2 Target features for our simulations
The Earth’s core operates at low viscosity (as measured by E ≃ 10−15 and Pm ≃ 10−6), so we should try to lower E and Pm as much as possible. Fig. 1 (left) illustrates the difficulty to lower E and Pm: the magnetic Reynolds number Rm (represented by the size of the circles) decreases and at some point it is not enough to maintain a magnetic field against ohmic dissipation.
Another outstanding feature of Earth’s dynamo is the small ratio of kinetic to magnetic energy, which is the squared Alfvén number A2 ≃ 10−4. It is apparent in Fig. 1 (right) that low A are difficult to obtain for Pm < 1. Low Alfvén number dynamos are readily obtained with Pm ≥ 1 (e.g. Kageyama et al. 2008; Dormy 2016) and just above the convection threshold (e.g. Takahashi & Shimizu 2012). Sakuraba & Roberts (2009) argued that imposing the heat flux rather than the temperature at the boundaries also helps to get a strong magnetic field (i.e. a low A). However, the convective power was not the same for their fixed flux and fixed temperature cases. Indeed, for fixed flux the convective power is proportional to Ra, whereas for fixed temperature the convective power is proportional to (Nu − 1)Ra, leading to weaker driving power in the fixed flux case for the same value of super-criticality Ra/Rac (Aubert et al. 2017, Fig. 1). According to the scaling laws of Christensen & Aubert (2006), the Alfvén number A increases when the convective power increases, and the trend is in broad agreement with the field strength obtained in our simulations.
Our goal is to obtain vigorous convection (far above its onset) and low Pm dynamo with a strong magnetic field (A < 1).
When trying to compute geodynamo models as close as possible to the parameters of the Earth’s core, the computation cost increases not only because of the higher and higher spatial resolution required, but also because the time-step is smaller and smaller compared to the magnetic diffusion time. Hence, in order to reach a statistically stationary dynamo regime, the time needed for a simulation to run increases prohibitively.1 To reach extreme parameters in our simulations, we prepare the initial conditions close to the expected statistical equilibrium state, in order to reduce the duration of transients. Those initial conditions are obtained by applying previously established scaling laws (Christensen & Aubert 2006) to the output of a lower resolution simulation at parameters further from the Earth’s core. This procedure can be repeated to achieve simulations that are closer and closer to the conditions of the Earth’s core.
2.3 Numerical implementation
We use the XSHELLS code, available as free software.2 It passes the dynamo benchmark of Christensen et al. (2001) as reported by Matsui et al. (2016). It uses second order finite differences in radius and pseudo-spectral spherical harmonic expansion. The time-stepping scheme is second order in time, and treats the diffusive terms implicitly, while the nonlinear and Coriolis terms are handled explicitly. We have carefully optimized the code for speed. The SHTns library (Schaeffer 2013) performs all the spherical harmonic transforms and, thanks to its low memory requirement and high performance, we can reach high resolutions (up to harmonic degree ℓmax = 1000 in this study). A domain decomposition in the radial direction allows efficient parallel execution using multiple processes (using the MPI standard). In addition, within each process, multiple threads (using the OpenMP standard) are used for an added level of parallelism. The most demanding simulation presented below has been run routinely on up to 8192 cores with a good scaling of the run time with the number of cores. More details about performance can be found in Matsui et al. (2016) and in the user manual.2
3 THE SIMULATIONS
3.1 Overview
We have run three dynamo simulations (S0, S1 and S2) that have decreasing Ekman numbers, E = 10−5, 10−6 and 10−7 respectively. They all have a magnetic Reynolds number Rm between 500 and 700. This is achieved by lowering the magnetic Prandtl number Pm and increasing the Rayleigh number to preserve high supercriticality (Ra/Rac > 4870 is almost constant, where Rac is the Rayleigh number for the onset of non-magnetic convection). In order to avoid transients, initial fields for simulations S1 and S2 are produced from S0 and S1 respectively, by multiplying each field by a scaling factor to reach the energy levels expected from the scaling laws of Christensen & Aubert (2006). Note that the S0 simulation is a continuation of the simulation presented by Fournier et al. (2012), and that parts of both S0 and S1 were used by Bouligand et al. (2016) with the same name.
In order to quantify the influence of the magnetic field, we have run a fourth simulation, S1*, which has the same parameters as S1, but without magnetic field. S1* is started from velocity and codensity fields of S1 at a given time, without any rescaling.
Table 1 summarizes the input parameters of our simulations, and compares it to the simulations of Kageyama et al. (2008), Sakuraba & Roberts (2009), Sheyko (2014), Yadav et al. (2016b) and Aubert et al. (2017). The table also includes some output parameters computed from averaged quantities. No magnetic polarity reversal occurred in our simulations.
. | K08 . | S09 . | S14 . | Y16 . | A17 . | S0 . | S1 . | S2 . | Earth . |
---|---|---|---|---|---|---|---|---|---|
Nr | 511 | 160 | 528 | 181 | 624 | 256 | 512 | 1280 | |
Nθ | 1024 | 384 | 384 | 640 | 200 | 448 | 720 | 1504 | |
Nϕ | 2048 | 768 | 768 | 1280 | 400 | 720 | 1440 | 2688 | |
time | 0.016 | 0.027 | 0.05 | 0.12 | 2.2 | 0.51 | 0.052 | ||
E | 9.4 × 10−7 | 2.4 × 10−6 | 5.9 × 10−7 | 10−6 | 10−8 | 10−5 | 10−6 | 10−7 | 3 × 10−15 |
Ra/Rac | ∼300 | 144 | 470 | 400* | 4879 | 5770 | 6310 | 106 ? | |
Pm | 1 | 0.2 | 0.05 | 0.4 | 0.045 | 0.4 | 0.2 | 0.1 | 2 × 10−6 |
Pr | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 0.1–10 |
Rm | 700 | 120 | 274 | 346 | 1082 | 671 | 546 | 514 | 2000 |
A | 0.42 | 0.46 | 0.35 | 0.3 | 0.11 | 1.2 | 0.43 | 0.27 | 0.01 |
Re | 700 | 600 | 5480 | 865 | 24000 | 1680 | 2730 | 5140 | 109 |
Ro | 7.6 × 10−4 | 1.4 × 10−3 | 3.2 × 10−3 | 8.7 × 10−4 | 2.4 × 10−4 | 0.017 | 2.7 × 10−3 | 5.1 × 10−4 | 3 × 10−6 |
Le | 1.8 × 10−3 | 2.2 × 10−3 | 9.2 × 10−3 | 2.9 × 10−3 | 2.1 × 10−3 | 0.014 | 6.4 × 10−3 | 1.9 × 10−3 | 10−4 |
Λ | 3 | 0.8 | 3.6 | 3.4 | 20 | 8.0 | 8.2 | 3.7 | ≳10 |
. | K08 . | S09 . | S14 . | Y16 . | A17 . | S0 . | S1 . | S2 . | Earth . |
---|---|---|---|---|---|---|---|---|---|
Nr | 511 | 160 | 528 | 181 | 624 | 256 | 512 | 1280 | |
Nθ | 1024 | 384 | 384 | 640 | 200 | 448 | 720 | 1504 | |
Nϕ | 2048 | 768 | 768 | 1280 | 400 | 720 | 1440 | 2688 | |
time | 0.016 | 0.027 | 0.05 | 0.12 | 2.2 | 0.51 | 0.052 | ||
E | 9.4 × 10−7 | 2.4 × 10−6 | 5.9 × 10−7 | 10−6 | 10−8 | 10−5 | 10−6 | 10−7 | 3 × 10−15 |
Ra/Rac | ∼300 | 144 | 470 | 400* | 4879 | 5770 | 6310 | 106 ? | |
Pm | 1 | 0.2 | 0.05 | 0.4 | 0.045 | 0.4 | 0.2 | 0.1 | 2 × 10−6 |
Pr | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 0.1–10 |
Rm | 700 | 120 | 274 | 346 | 1082 | 671 | 546 | 514 | 2000 |
A | 0.42 | 0.46 | 0.35 | 0.3 | 0.11 | 1.2 | 0.43 | 0.27 | 0.01 |
Re | 700 | 600 | 5480 | 865 | 24000 | 1680 | 2730 | 5140 | 109 |
Ro | 7.6 × 10−4 | 1.4 × 10−3 | 3.2 × 10−3 | 8.7 × 10−4 | 2.4 × 10−4 | 0.017 | 2.7 × 10−3 | 5.1 × 10−4 | 3 × 10−6 |
Le | 1.8 × 10−3 | 2.2 × 10−3 | 9.2 × 10−3 | 2.9 × 10−3 | 2.1 × 10−3 | 0.014 | 6.4 × 10−3 | 1.9 × 10−3 | 10−4 |
Λ | 3 | 0.8 | 3.6 | 3.4 | 20 | 8.0 | 8.2 | 3.7 | ≳10 |
. | K08 . | S09 . | S14 . | Y16 . | A17 . | S0 . | S1 . | S2 . | Earth . |
---|---|---|---|---|---|---|---|---|---|
Nr | 511 | 160 | 528 | 181 | 624 | 256 | 512 | 1280 | |
Nθ | 1024 | 384 | 384 | 640 | 200 | 448 | 720 | 1504 | |
Nϕ | 2048 | 768 | 768 | 1280 | 400 | 720 | 1440 | 2688 | |
time | 0.016 | 0.027 | 0.05 | 0.12 | 2.2 | 0.51 | 0.052 | ||
E | 9.4 × 10−7 | 2.4 × 10−6 | 5.9 × 10−7 | 10−6 | 10−8 | 10−5 | 10−6 | 10−7 | 3 × 10−15 |
Ra/Rac | ∼300 | 144 | 470 | 400* | 4879 | 5770 | 6310 | 106 ? | |
Pm | 1 | 0.2 | 0.05 | 0.4 | 0.045 | 0.4 | 0.2 | 0.1 | 2 × 10−6 |
Pr | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 0.1–10 |
Rm | 700 | 120 | 274 | 346 | 1082 | 671 | 546 | 514 | 2000 |
A | 0.42 | 0.46 | 0.35 | 0.3 | 0.11 | 1.2 | 0.43 | 0.27 | 0.01 |
Re | 700 | 600 | 5480 | 865 | 24000 | 1680 | 2730 | 5140 | 109 |
Ro | 7.6 × 10−4 | 1.4 × 10−3 | 3.2 × 10−3 | 8.7 × 10−4 | 2.4 × 10−4 | 0.017 | 2.7 × 10−3 | 5.1 × 10−4 | 3 × 10−6 |
Le | 1.8 × 10−3 | 2.2 × 10−3 | 9.2 × 10−3 | 2.9 × 10−3 | 2.1 × 10−3 | 0.014 | 6.4 × 10−3 | 1.9 × 10−3 | 10−4 |
Λ | 3 | 0.8 | 3.6 | 3.4 | 20 | 8.0 | 8.2 | 3.7 | ≳10 |
. | K08 . | S09 . | S14 . | Y16 . | A17 . | S0 . | S1 . | S2 . | Earth . |
---|---|---|---|---|---|---|---|---|---|
Nr | 511 | 160 | 528 | 181 | 624 | 256 | 512 | 1280 | |
Nθ | 1024 | 384 | 384 | 640 | 200 | 448 | 720 | 1504 | |
Nϕ | 2048 | 768 | 768 | 1280 | 400 | 720 | 1440 | 2688 | |
time | 0.016 | 0.027 | 0.05 | 0.12 | 2.2 | 0.51 | 0.052 | ||
E | 9.4 × 10−7 | 2.4 × 10−6 | 5.9 × 10−7 | 10−6 | 10−8 | 10−5 | 10−6 | 10−7 | 3 × 10−15 |
Ra/Rac | ∼300 | 144 | 470 | 400* | 4879 | 5770 | 6310 | 106 ? | |
Pm | 1 | 0.2 | 0.05 | 0.4 | 0.045 | 0.4 | 0.2 | 0.1 | 2 × 10−6 |
Pr | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 0.1–10 |
Rm | 700 | 120 | 274 | 346 | 1082 | 671 | 546 | 514 | 2000 |
A | 0.42 | 0.46 | 0.35 | 0.3 | 0.11 | 1.2 | 0.43 | 0.27 | 0.01 |
Re | 700 | 600 | 5480 | 865 | 24000 | 1680 | 2730 | 5140 | 109 |
Ro | 7.6 × 10−4 | 1.4 × 10−3 | 3.2 × 10−3 | 8.7 × 10−4 | 2.4 × 10−4 | 0.017 | 2.7 × 10−3 | 5.1 × 10−4 | 3 × 10−6 |
Le | 1.8 × 10−3 | 2.2 × 10−3 | 9.2 × 10−3 | 2.9 × 10−3 | 2.1 × 10−3 | 0.014 | 6.4 × 10−3 | 1.9 × 10−3 | 10−4 |
Λ | 3 | 0.8 | 3.6 | 3.4 | 20 | 8.0 | 8.2 | 3.7 | ≳10 |
Fig. 2 shows both kinetic and magnetic energies as a function of time, for our four simulations. It is clear that the ratio of magnetic over kinetic energy increases significantly from S0 to S2. In addition, the kinetic energy of S1* is about 10 times larger than that of S1, showing the strong influence of the magnetic field. In S1 and S2, the magnetic energy has fluctuations of larger amplitudes but lower frequencies than the kinetic energy. Furthermore, all the dynamo simulations exhibit much larger variations of their energy levels compared to the non-magnetic case S1*, with much longer correlation times. This seems to be an inherent property of turbulent dynamo simulations. We note that the magnetic energy in S1 and S2 increases above the level expected by the scaling laws of Christensen & Aubert (2006), but not by more than a factor 3, which corresponds also to the fluctuations observed in S0. However, the magnetic energy of all three simulations correspond to the trend exhibited by the Christensen & Aubert (2006) data set, which also displays a scatter by a factor 3 to 10 (not shown).
We have checked the spatial convergence of S2 by multiplying the number of radial shells by 1.5, and the maximum degree of spherical harmonics by 1.12. No notable difference was found in the energy levels. Furthermore, the highest degrees of the spherical harmonic spectra have two orders of magnitude lower energy than the most energetic degree, at each radius. As a supplementary convergence check, the time-averaged total energy dissipation rate is equal to the time-averaged power of the buoyancy force within less than 1 per cent in all our dynamo simulations (values listed in Table 2).
. | S0 . | S1 . | S1* . | S2 . |
---|---|---|---|---|
E | 10−5 | 10−6 | 10−6 | 10−7 |
Ra | 6.34 × 109 | 1.27 × 1011 | 1.27 × 1011 | 2.54 × 1012 |
Ra/Rac | 4879 | 5770 | 5770 | 6310 |
Pm | 0.4 | 0.2 | 0 | 0.1 |
Rm | 671 | 546 | 0 | 514 |
A | 1.2 | 0.43 | 0.27 | |
Re | 1680 | 2730 | 8760 | 5140 |
Ro | 0.017 | 2.7 × 10−3 | 8.8 × 10−3 | 5.1 × 10−4 |
Le | 0.014 | 6.4 × 10−3 | 0 | 1.9 × 10−3 |
Λ | 8.0 | 8.2 | 0 | 3.7 |
Nu | 31 | 45 | 42 | 59 |
Pconv | 2.643 × 1011 | 5.579 × 1012 | 5.54 × 1012 | 10.3 × 1013 |
Dη | 1.540 × 1011 | 4.451 × 1012 | 0 | 8.71 × 1013 |
Dν | 1.085 × 1011 | 1.113 × 1012 | 5.48 × 1012 | 1.44 × 1013 |
fohm = Dη/(Dν + Dη) | 0.59 per cent | 0.80 per cent | 0.00 per cent | 0.86 per cent |
Lu | 0.020 | 0.011 | 0.014 | 5.2 × 10−3 |
Lb | 0.022 | 0.025 | 0.023 | |
mc | 15 | 32 | 32 | 67 |
|$\bar{\ell }$| | 26 | 46 | 24 | 86 |
Ek(m = 0)/Ek | 0.09 | 0.12 | 0.50 | 0.06 |
fdip | 0.49 | 0.72 | 0.67 | |
Bsurf(ℓ = 1)/Brms | 0.11 | 0.16 | 0.10 | |
|$\langle |\mathcal {T}(s,t)|\rangle _{s,t}$| | 0.447 | 0.161 | 0.119 | |
|$\langle |\langle \mathcal {T}(s,t)\rangle _t|\rangle _s$| | 0.257 | 0.0247 | 0.0219 | |
Pconv/(Dν + Dη) | 1.007 | 1.003 | 1.011 | 1.010 |
Nr | 256 | 512 | 512 | 1280 |
ℓmax | 297 | 479 | 479 | 1000 |
mmax | 238 | 479 | 479 | 893 |
million core.hours | ≃0.4 | ≃1.4 | ≃0.4 | ≃10 |
. | S0 . | S1 . | S1* . | S2 . |
---|---|---|---|---|
E | 10−5 | 10−6 | 10−6 | 10−7 |
Ra | 6.34 × 109 | 1.27 × 1011 | 1.27 × 1011 | 2.54 × 1012 |
Ra/Rac | 4879 | 5770 | 5770 | 6310 |
Pm | 0.4 | 0.2 | 0 | 0.1 |
Rm | 671 | 546 | 0 | 514 |
A | 1.2 | 0.43 | 0.27 | |
Re | 1680 | 2730 | 8760 | 5140 |
Ro | 0.017 | 2.7 × 10−3 | 8.8 × 10−3 | 5.1 × 10−4 |
Le | 0.014 | 6.4 × 10−3 | 0 | 1.9 × 10−3 |
Λ | 8.0 | 8.2 | 0 | 3.7 |
Nu | 31 | 45 | 42 | 59 |
Pconv | 2.643 × 1011 | 5.579 × 1012 | 5.54 × 1012 | 10.3 × 1013 |
Dη | 1.540 × 1011 | 4.451 × 1012 | 0 | 8.71 × 1013 |
Dν | 1.085 × 1011 | 1.113 × 1012 | 5.48 × 1012 | 1.44 × 1013 |
fohm = Dη/(Dν + Dη) | 0.59 per cent | 0.80 per cent | 0.00 per cent | 0.86 per cent |
Lu | 0.020 | 0.011 | 0.014 | 5.2 × 10−3 |
Lb | 0.022 | 0.025 | 0.023 | |
mc | 15 | 32 | 32 | 67 |
|$\bar{\ell }$| | 26 | 46 | 24 | 86 |
Ek(m = 0)/Ek | 0.09 | 0.12 | 0.50 | 0.06 |
fdip | 0.49 | 0.72 | 0.67 | |
Bsurf(ℓ = 1)/Brms | 0.11 | 0.16 | 0.10 | |
|$\langle |\mathcal {T}(s,t)|\rangle _{s,t}$| | 0.447 | 0.161 | 0.119 | |
|$\langle |\langle \mathcal {T}(s,t)\rangle _t|\rangle _s$| | 0.257 | 0.0247 | 0.0219 | |
Pconv/(Dν + Dη) | 1.007 | 1.003 | 1.011 | 1.010 |
Nr | 256 | 512 | 512 | 1280 |
ℓmax | 297 | 479 | 479 | 1000 |
mmax | 238 | 479 | 479 | 893 |
million core.hours | ≃0.4 | ≃1.4 | ≃0.4 | ≃10 |
. | S0 . | S1 . | S1* . | S2 . |
---|---|---|---|---|
E | 10−5 | 10−6 | 10−6 | 10−7 |
Ra | 6.34 × 109 | 1.27 × 1011 | 1.27 × 1011 | 2.54 × 1012 |
Ra/Rac | 4879 | 5770 | 5770 | 6310 |
Pm | 0.4 | 0.2 | 0 | 0.1 |
Rm | 671 | 546 | 0 | 514 |
A | 1.2 | 0.43 | 0.27 | |
Re | 1680 | 2730 | 8760 | 5140 |
Ro | 0.017 | 2.7 × 10−3 | 8.8 × 10−3 | 5.1 × 10−4 |
Le | 0.014 | 6.4 × 10−3 | 0 | 1.9 × 10−3 |
Λ | 8.0 | 8.2 | 0 | 3.7 |
Nu | 31 | 45 | 42 | 59 |
Pconv | 2.643 × 1011 | 5.579 × 1012 | 5.54 × 1012 | 10.3 × 1013 |
Dη | 1.540 × 1011 | 4.451 × 1012 | 0 | 8.71 × 1013 |
Dν | 1.085 × 1011 | 1.113 × 1012 | 5.48 × 1012 | 1.44 × 1013 |
fohm = Dη/(Dν + Dη) | 0.59 per cent | 0.80 per cent | 0.00 per cent | 0.86 per cent |
Lu | 0.020 | 0.011 | 0.014 | 5.2 × 10−3 |
Lb | 0.022 | 0.025 | 0.023 | |
mc | 15 | 32 | 32 | 67 |
|$\bar{\ell }$| | 26 | 46 | 24 | 86 |
Ek(m = 0)/Ek | 0.09 | 0.12 | 0.50 | 0.06 |
fdip | 0.49 | 0.72 | 0.67 | |
Bsurf(ℓ = 1)/Brms | 0.11 | 0.16 | 0.10 | |
|$\langle |\mathcal {T}(s,t)|\rangle _{s,t}$| | 0.447 | 0.161 | 0.119 | |
|$\langle |\langle \mathcal {T}(s,t)\rangle _t|\rangle _s$| | 0.257 | 0.0247 | 0.0219 | |
Pconv/(Dν + Dη) | 1.007 | 1.003 | 1.011 | 1.010 |
Nr | 256 | 512 | 512 | 1280 |
ℓmax | 297 | 479 | 479 | 1000 |
mmax | 238 | 479 | 479 | 893 |
million core.hours | ≃0.4 | ≃1.4 | ≃0.4 | ≃10 |
. | S0 . | S1 . | S1* . | S2 . |
---|---|---|---|---|
E | 10−5 | 10−6 | 10−6 | 10−7 |
Ra | 6.34 × 109 | 1.27 × 1011 | 1.27 × 1011 | 2.54 × 1012 |
Ra/Rac | 4879 | 5770 | 5770 | 6310 |
Pm | 0.4 | 0.2 | 0 | 0.1 |
Rm | 671 | 546 | 0 | 514 |
A | 1.2 | 0.43 | 0.27 | |
Re | 1680 | 2730 | 8760 | 5140 |
Ro | 0.017 | 2.7 × 10−3 | 8.8 × 10−3 | 5.1 × 10−4 |
Le | 0.014 | 6.4 × 10−3 | 0 | 1.9 × 10−3 |
Λ | 8.0 | 8.2 | 0 | 3.7 |
Nu | 31 | 45 | 42 | 59 |
Pconv | 2.643 × 1011 | 5.579 × 1012 | 5.54 × 1012 | 10.3 × 1013 |
Dη | 1.540 × 1011 | 4.451 × 1012 | 0 | 8.71 × 1013 |
Dν | 1.085 × 1011 | 1.113 × 1012 | 5.48 × 1012 | 1.44 × 1013 |
fohm = Dη/(Dν + Dη) | 0.59 per cent | 0.80 per cent | 0.00 per cent | 0.86 per cent |
Lu | 0.020 | 0.011 | 0.014 | 5.2 × 10−3 |
Lb | 0.022 | 0.025 | 0.023 | |
mc | 15 | 32 | 32 | 67 |
|$\bar{\ell }$| | 26 | 46 | 24 | 86 |
Ek(m = 0)/Ek | 0.09 | 0.12 | 0.50 | 0.06 |
fdip | 0.49 | 0.72 | 0.67 | |
Bsurf(ℓ = 1)/Brms | 0.11 | 0.16 | 0.10 | |
|$\langle |\mathcal {T}(s,t)|\rangle _{s,t}$| | 0.447 | 0.161 | 0.119 | |
|$\langle |\langle \mathcal {T}(s,t)\rangle _t|\rangle _s$| | 0.257 | 0.0247 | 0.0219 | |
Pconv/(Dν + Dη) | 1.007 | 1.003 | 1.011 | 1.010 |
Nr | 256 | 512 | 512 | 1280 |
ℓmax | 297 | 479 | 479 | 1000 |
mmax | 238 | 479 | 479 | 893 |
million core.hours | ≃0.4 | ≃1.4 | ≃0.4 | ≃10 |
3.2 Mean fields
Time and longitude averaged fields are represented in Fig. 3, for the three dynamo simulations S0, S1 and S2, and for the non-magnetic convection simulation S1*.
The root-mean-square (rms) poloidal magnetic field is comparable to the toroidal one, with a poloidal over toroidal ratio of 1.73, 0.95 and 0.83 for simulations S0 to S2. In all three cases, the amplitude of meridional velocity is weak, about 4–6 per cent of the zonal one.
A north–south asymmetry remains in S0 and S2 after averaging over the time Tavail when the required data are available (in terms of turnover time τu = D/U: Tavail/τu = 1390 for S0, 190 for S1, and 8 for S2). Note that a similar asymmetry is present in snapshots of S1 (not shown).
The imaginary cylinder tangent to the inner-core and aligned with the rotation axis (hereafter named tangent cylinder or TC) separates two regions of different dynamics. The lower the Ekman number, the sharper the transition between these two regions. This dichotomy is visible on both velocity components as well as on the toroidal magnetic field, for S1 (E = 10−6) and S2 (E = 10−7). A sharp shear layer, associated with a meridional circulation materializes the tangent cylinder, and is reminiscent of Stewartson layers (Stewartson 1966). The zonal mean flow inside the tangent cylinder has amplitudes about 10 times larger than outside. This contrasts with the findings of Aubert (2005) who argued for isorotation along magnetic field lines and had comparable flow amplitudes inside and outside the tangent cylinder. However, we checked that the scaling of the zonal flow proposed by Aubert (2005) seems to hold.
Interestingly, the total codensity shows that lighter fluid is trapped within the tangent cylinder (see also Fig. 7), with strong lateral gradients. This feature is much less pronounced in the non-magnetic case S1*. A more quantitative view is given by the codensity profiles of Fig. A2 in the Appendix. It is noteworthy that all dynamo simulations present an adverse codensity gradient within the tangent cylinder. This emergence of a stable stratification within TC is most likely an effect of the strong magnetic field, as the stratification remains unstable everywhere in the non-magnetic S1*. Note that in S0, light material is not well contained and appears to leak (especially in the northern hemisphere, see Fig. 3 top). In S1, the leak is barely visible and the confinement seems even better in S2.
3.2.1 Outside the tangent cylinder
Outside the tangent cylinder, the averaged flow is weak in both S1 and S2. The zonal wind amplitude is Ro ≃ 5 × 10−5 in S2, about 10 times smaller than the overall rms averaged Rossby number (see Table 2). In particular, the zonal wind is 10 times smaller outside than inside the tangent cylinder, and 30 times smaller in S1 than in the corresponding non-magnetic thermal convection of S1* (see Fig. 3), showing the strong effect of the magnetic field on the mean flow. Note also that the weak mean flow next to the equator of the core–mantle boundary (CMB) is rather eastward (prograde) in S1, but westward (retrograde) in S2. The mean toroidal magnetic field is also weak in S1 (corresponding Elsasser number Λ ≃ 0.05) and even more so in S2 (Λ ≃ 0.01), but the poloidal field is rather strong, especially near the inner-core (where the associated Elsasser number reaches 0.25 in S2, see Fig. 3 bottom). The much lower values of Elsasser numbers here compared to the global root-mean-square Elsasser given in Table 1 emphasize the weak average toroidal magnetic field outside TC.
3.2.2 Inside the tangent cylinder
Inside the tangent cylinder, a strong azimuthal flow takes place, in the prograde direction close to the inner-core, and in the retrograde direction close to the mantle, with an associated meridional circulation (one cell). The zonal flow amplitude is Ro ≃ 5 × 10−4 in S2, similar to the overall rms Rossby number (see Table 2). They form what has been coined as the polar vortex (Olson & Aurnou 1999; Aubert 2005). This twisted polar vortex is anticyclonic next to the CMB. It has been associated with a low influence of inertia and a strong magnetic field (Sreenivasan & Jones 2006). This flow resembles the Taylor vortices described in rotating thermal convection (e.g. Grooms et al. 2010), or the Von-Karman flow generated by two impellers co-rotating at different speeds (in order to obtain one meridional circulation cell). In this respect, the VKS experiment (Monchaux et al. 2009) may be not so far from the flow inside one hemisphere of the tangent cylinder. This contrasts with the non-magnetic case S1*, where the zonal flow is almost invariant along the rotation axis. Conversely, codensity shows little variations along the rotation axis in the dynamo simulations, but important ones in S1*. The twisted zonal flow within the tangent cylinder also reaches higher speeds in S1 than in S1*, suggesting that the strong toroidal field there allows the flow to break the Taylor–Proudman constraint imposed by the global rotation.
The toroidal magnetic field is concentrated here, suggesting that an omega effect is associated with the strong twisted vortices that dominate the mean flow. The Elsasser number associated with this strong toroidal field is close to unity in S0, S1 and S2. At Elsasser close to one, convection onsets more easily and its length-scale increases dramatically from order E1/3 to order 1 (Chandrasekhar 1961; Busse & Simitev 2011; Aujogue et al. 2015). All the conditions for such an effect are met within the tangent cylinder of S1 and S2. Note also that the toroidal magnetic field in S2 seems to have switched to another topology (only one sign in each hemisphere of S2 instead of a change of sign in S1), with minor impact on the mean poloidal field (but see Fig. 8).
3.2.3 Non-zonal mean flows
Non-zonal mean flows appear to be a prominent feature of S2. Fig. 4 shows velocity field averaged in time and along the rotation axis (z-axis—note that the average along z spans only one hemisphere within the tangent cylinder). In S1, the z-averaged mean flow is dominated by a prograde jet at the tangent cylinder and a retrograde circulation within. A weak large-scale non-zonal flow is still visible outside the TC in S1. When averaging over a shorter time-span (24 turnover times instead of 190), the non-zonal flow dominates. The least viscous and strongly magnetized simulation S2 is dominated by a large non-zonal gyre outside the TC. Although arguably limited, the time span available for averaging the S2 flow is significant and would translate to about thousand years if the turnover time is used for scaling to Earth values. Similarly, the non-zonal z-invariant flow seems to last about 24 turnover times in S1, corresponding to more than 5000–10 000 yr. One outstanding difference between S1 and S2 is that the anticyclonic gyre leads to strong westward velocities at the equator in S2, whereas it is weak and mostly eastward in S1.
3.3 Instantaneous fields
We now turn to instant snapshots of the fields, which are represented in the equatorial plane for simulations S1 and S1* in Fig. 5, contrasting the differences between dynamo (S1) and non-magnetic convection (S1*). The codensity field C exhibits very small scales near the inner-core where the plumes originate in both S1 and S1*. Further away from the inner-core, the codensity field exhibits much larger structures in S1 than in S1*. The plumes also reach further out in S1 whereas they seem to be stopped by the zonal winds in S1* (see Ur in Fig. 5). However, the overall state is better mixed in S1* (lower contrasts in variations of C). This illustrates the effect of the magnetic field on the convection.
Fig. 6 shows similar views for S2 at two different times characterized by moderate and strong magnetic fields (respectively marked by circle and square in Fig. 2). 3-D renderings of the strongest field situation are shown in Fig. 7, displaying also meridional cuts.3 Very small-scale buoyant plumes originate near the inner-core for the moderate and strong field snapshots, but further away the scales appear larger where the strongest field reigns. Radial velocities are also weaker in these regions. There are notable velocity field patterns with azimuthal wavenumbers much smaller than m = 67—the critical wavenumber at onset—especially in the strong field regions. These large scales coexist with smaller scales close to the inner-core, but also in regions of weaker magnetic field.
This enlargement of convection scale is in broad agreement with the observations made by Matsui et al. (2014) and Yadav et al. (2016b) at higher Ekman numbers and lower forcing (see also Hori et al. 2012). The mechanism proposed by Matsui et al. (2014) does also fit our simulations. Namely, the presence of important variations of codensity over large regions (Fig. 6 top) produces non-axisymmetric thermal winds that convert the poloidal field into azimuthal field. The field inhibits motions, preventing the codensity anomalies to mix, thus sustaining the phenomena. Large scale motions induce less shear and are thus favoured compared to small-scale convection.
We would like to emphasize that the large scale codensity anomalies can build up only because the zonal flow is suppressed by the magnetic field. Otherwise codensity anomalies in the azimuthal direction are quickly sheared away by the zonal jets (which effectively improve lateral mixing), as seen in the non-magnetic convection of S1* (Fig. 5). This is also supported by the codensity profiles (see Fig. A2) that are better mixed (flatter) in S1 than in S1*.
In our simulations S1 and S2, magnetic field intensity is largely inhomogeneous (see |B| in Figs 5–7). This does not seem to affect the surface field which is already an order of magnitude smaller than the bulk average. Because the strength of the magnetic field is expected to control the length-scale of convection, then convection in planetary core would span a wide range of length-scales, possibly from the viscous scale E1/3 to the planetary scale. It is important to take this into account when conceptualizing turbulence in planetary cores (e.g. Nataf & Schaeffer 2015).
We remark also that because the kinetic energy spectra are not steep, using the length-scale diagnostic introduced by Christensen & Aubert (2006, and widely used afterwards) fails to capture accurately the change in length-scale, especially since the viscous length-scale may still be present in the system when the Rayleigh number is highly super-critical as in our case, in contrast to the work of Matsui et al. (2014). This failure may also explain why the viscous length-scale was found to play a role in previous dynamo studies (King & Buffett 2013) despite large Elsasser numbers.
3.4 Magnetic field at the core surface
It is important to look at the magnetic field at the surface of the core, as it can be readily compared with the geomagnetic field. Following (Christensen & Aubert 2006), the fraction of axial dipole in the observable spectrum (up to ℓ = 13) is given by fdip in Table 2, and is in reasonable agreement with the Earth’s value (fdip ≃ 0.68). Fig. 8 shows the magnetic field at the surface of S1 and S2 at a given time (the same time as in Figs 5 and 7). Interestingly, the surface magnetic field displays similarities with the surface codensity (see Fig. 7), in particular sharp gradients near the tangent cylinder. We note that the snapshot of S2 in Fig. 8 displays a reversed flux patch at high latitude (within the tangent cylinder) in the northern hemisphere. This feature is not persistent, as can be seen in the movie (see caption of Fig. 8). We remark that magnetic field models derived from satellite measurements show an inverse flux patch (albeit smaller) within TC (see e.g. Hulot et al. 2015, fig. 25).
An important feature of the Earth’s magnetic field is its secular variation, and in particular the westward drift of magnetic flux patches near the equator of the Atlantic hemisphere. The evolution in time of the field at the equator is represented in Fig. 9 for S1 and S2. Whereas there is mostly eastward drift in S1, S2 displays an important region of westward drift (roughly between 30° and 180° longitude) together with a region with no net drift (around 300°). These plots also highlight the larger longitudinal extent of patches in S2 than in S1 and their longer life-span. There is a one-to-one correspondence between the westward drifting region in S2 seen in Fig. 9 and the location of the strong westward circulation in Fig. 4. This demonstrates that large-scale westward drift arises naturally in our less viscous simulation. We also note the appearance of alternating polarity magnetic field with a well determined wave-length drifting west in the upper left part of the bottom left panel of Fig. 9, akin to the equatorial slow waves of Finlay & Jackson (2003).
3.5 Fluctuations and helicity
Velocity fluctuations are large in the polar regions just below the CMB and near the inner-core boundary (ICB), but rather weak in between. This suggests that the velocity fluctuations inside TC are associated to variability of the twisted polar vortices. Outside TC, velocity fluctuations are always small near the equator of the CMB, and increasing gradually towards the tangent cylinder, where the mean poloidal field is also concentrated.
To link these fluctuations with a possible poloidal magnetic field generation, we turn to helicity (Fig. 10 bottom and Fig. 11), which is often associated with alpha effect (e.g. Moffatt 1978; Jones 2008) whereby poloidal magnetic field is produced from toroidal field. Helicity maps exhibit a gradual change. In S0 helicity is maximum near the boundaries and extends towards the equator just outside the TC. In S2 the maximum helicity in the outside region is located near the inner-core. S1 shows both boundary bound helicity together with local maximum in the bulk (Fig. 10 bottom). Outside the tangent cylinder, relative helicity is larger in S2 than in S1. This suggests a transition from helicity mostly due to recirculation associated with the strong zonal jets just outside TC in S0 and S1* (see Busse 1975) to bulk helicity in S1 and S2, which might be a hint for the mechanism proposed by Sreenivasan & Jones (2011, relying on magnetic pumping), or the one proposed by Davidson (2014, relying on inertial waves). Inside TC, it may be simply the contribution of the helical convective vortices (see e.g. Aurnou et al. 2003; Grooms et al. 2010). In any case, we can confidently say that the S2 dynamo does not rely on helicity induced by Ekman pumping. Similarly, Ekman pumping cannot be an important source of helicity in Earth’s core (e.g. Schaeffer & Cardin 2006). In Fig. 11, strong helicity is found at the boundaries. The outer boundary is the Ekman layer, and scales as E1/2 (barely visible in Fig. 11). The inner boundary is much larger and we could not find any scaling that convincingly collapse the three dynamos. The profiles are significantly different inside and outside TC. The effect of the magnetic field on these profiles seems weak outside TC, but seizable inside TC (compare S1 and S1* in Fig. 11).
3.6 Spectra
Fig. 12 shows kinetic and magnetic energy spectra for all our simulations. These time- and radius-averaged spectra are misleading in several ways. First, they encompass regions with vastly different dynamics and length-scale (inside and outside the tangent cylinder). Second, even in these regions, fluctuations are largely inhomogeneous (see Fig. 10). Finally, a given harmonic degree ℓ corresponds to different length-scale depending on the radius.
We remark however that the non-magnetic simulation S1* displays a significant range of kinetic energy spectrum obeying a ℓ−3/2 law. This contrasts with the kinetic energy spectra in the dynamo simulations S1 and S2, which are much less steep, being almost flat at large scales, and then displaying more or less a ∼ℓ−1 range. Note that the exponents quoted here are purely empirical, and that no clear relation to theoretical Cartesian spectra can be made because of the spherical geometry and the anisotropy of such rotating flows. Christensen & Aubert (2006) introduced an average harmonic degree|$\bar{\ell } = \sum _\ell {\ell E_\ell (u)}/\sum _\ell {E_\ell (u)}$|. With a kinetic energy spectrum Eℓ(u) ∼ ℓ−1, |$\bar{\ell }$| will be determined by the smallest scales, although larger scales are present (see Fig. 4 and 12). This may explain why King & Buffett (2013) found that |$\bar{\ell }$| seems to scale as the inverse of the viscous scale in a large database of dynamo simulations.
Fig. 13 shows the energy spectra at 4 different radii in the S2 simulation. Near the inner-core (r = 0.588), the spectra are dominated by the fields inside of the tangent cylinder. The strength of both velocity and magnetic fields at large scales (ℓ ≤ 5) is much larger than at larger radii.
In the middle of the shell (r = 1.04), the magnetic spectrum is not clearly dominated by the dipole component but is rather flat up to degree 20 to 30, and then decays quickly beyond ℓ ∼ 50. The kinetic energy spectrum is about 20 to 30 times smaller than the magnetic spectrum, and also peaks at degree ℓ ≃ 20.
In the outer part (r = 1.4) the magnetic spectrum is not yet dominated by the dipole component ℓ = 1, but peaks for degree 5 ≤ ℓ ≤ 10, just as the velocity spectrum.
At the surface (r = 1.535—just below the Ekman layer), the magnetic spectrum is clearly dominated by the dipole. However, magnetic energy is much weaker than deeper in the shell. It is actually barely larger than the kinetic energy, unlike in the Earth. If we exclude the dipole (ℓ = 1) and octupole (ℓ = 3) parts, the magnetic spectrum is rather flat up to degree ℓ = 10. Beyond ℓ = 10, both velocity and magnetic spectra decay slower than at lower radii.
Finally, note that the peaks in the spectra are always located at wavenumber significantly smaller than those expected at the onset of thermal convection (mc = 67).
3.7 Spatio-temporal Fourier analysis
Because of the vastly different dynamics occurring inside and outside the tangent cylinder, we need to analyse these regions separately. This cannot be done with simple harmonic degree spectra. Furthermore, the temporal fluctuations visible in Fig. 2 also prompt us to analyse the fields and spectra beyond their time-averages.
We have recorded snapshots of the fields at regular time intervals. A full snapshot of S2 needs 32 gigabytes of disk space, and we could not afford to save many such large snapshots. For this reason, the many recorded snapshots are truncated at a spherical harmonic degree ℓtr that is smaller than the maximum ℓmax resolved by the simulation: ℓtr = 399 for S1 and ℓtr = 299 for S2. To save even more disk space, we also use single-precision to store the snapshots. We apply a Fourier transform in time and longitude (the two homogeneous directions) to the whole snapshots series. Unfortunately, some snapshots have been lost. For S1, we could effectively use a series of 523 regularly spaced snapshots sampled at period 25/Ω (every ∼4 rotation periods or 8000 samples per magnetic diffusion time). For S2, only a series of 102 snapshots sampled at period 80/Ω (every 12.7 rotation periods or 12 500 samples per magnetic diffusion time). We then study separately five different regions, represented in Fig. 14, inspired by the behaviour of mean fields described above: an inner and outer boundary layer region, located near the inner-core and near the mantle respectively, with a radial extent of 10 per cent of the outer shell radius; a tangent cylinder region, excluding the previous boundary layers and spanning cylindrical radii s from 0.3ro to 0.4ro; an inner region for s < 0.3 and an outer region for s > 0.4, both also excluding the two boundary layers. This allows us to describe quantitatively the temporal behaviour of the fields in our simulations. For each snapshot, we also compute the different terms in the evolution equations and perform the same Fourier analysis, allowing to describe the dynamical balances in our simulations in terms of spatial and temporal scales (Nataf & Schaeffer 2015).
3.7.1 Velocity and magnetic fields
The critical azimuthal order for the onset of convection is mc = 15, 32 and 67 for S0, S1 and S2, respectively. These are exact values computed by SINGE (Vidal & Schaeffer 2015) and are consistent with mc = 0.31 E−1/3. Figs 15 and 16 show that the azimuthal order of convection outside the tangent cylinder is around m = 10 for S1 and m = 20 for S2. These correspond to length-scales three times larger than the ones for the onset of convection, although they are still dependent on the viscosity, scaling roughly as m ≃ 0.1 E−1/3. However, in addition to this small scale convection, a non-zonal mean flow is produced in S2, and low frequency, low m flows are produced outside the tangent cylinder in S1 and S2. It is remarkable that in S2 the mean flow is dominated by m = 1 and 2 rather than m = 0. This mean flow is represented in Fig. 4.
As already noted, the dynamics is very different within the tangent cylinder. There, large scale motions dominate (the twisted polar vortices), and the time-variability peaks at m = 1, meaning that this vortex keeps a rather constant amplitude but deviates from the global rotation axis. Axisymmetric terms (m = 0) also dominate the magnetic field spectra. However, small-scale convection is also present where the magnetic field is weaker (as visible near the inner-core on Fig. 7, bottom).
The effect of the magnetic field can be seen by comparing Figs 16 and 17. Without magnetic field, the mean zonal flow dominates both inside and outside of the tangent cylinder, while more variability in both space and time is observed when the magnetic field is present, with important fluctuations of the m = 1 flow inside and outside TC.
3.7.2 Dynamical balances
We apply a Fourier transform, in both azimuthal and time directions, to the different terms in the vorticity equation (the curl of eq. 1). Note that examining the terms in the vorticity equation effectively extracts the second order balance, after the main geostrophic balance in which most of the Coriolis and pressure force have cancelled each other. Fields were truncated to ℓmax = 299 before computing the terms, so we lose the contribution of the smallest scales to the nonlinear terms. The spatial spectra shown in Fig. 12 suggest that this is not a problem within the bulk, as the spectra start to decay well before ℓ = 100, but we might lose significant contributions in the boundary layers. However, we have checked that truncating further at ℓmax = 100 does not alter the broad picture, although it lowers the relative contribution from nonlinear terms (Lorentz and advection).
The resulting 2-D spectra are represented in Fig. 18 for each region and each term of simulation S2. Clearly, the different regions have very different dynamics, and averaging over all of them will definitely blur the analysis. Because the analysed series spans only 1300 rotation periods for S2 (2081 for S1), we focus on the so-called rapid dynamics.
The broad picture, shown in Fig. 18, consists of a clear dominance of Coriolis, Buoyancy and Lorentz terms in the interior regions (far from boundaries). This is consistent with the so-called MAC (Magneto–Archimedes–Coriolis) balance.
Close to the boundaries the balance involves mostly viscosity and Lorentz forces perturbed by the Coriolis term, which is compatible with the balance in a Hartmann layer. We note that all forces are smaller for m = 0 than for m > 0, except for the time-average, which stands clearly out as a thermal wind balance, even outside the tangent cylinder. Remember however that the zonal flow is very much affected by the magnetic field (compare S1 and S1* in Fig. 3) even though the Lorentz force is absent from the main thermal wind balance. This analysis also illustrates the complexity of the spatio-temporal balances, which is completely overlooked when averaging in time. Note that for S1*, the balances do not depend on frequency (not shown).
With our spatio-temporal decomposition, we are able to refine the main MAC balance of the two interior regions. Fig. 19 shows the relative importance of the Lorentz and Buoyancy terms. Specifically, the Coriolis term is mostly balanced by the Lorentz term in purple regions, and by the buoyancy term in orange regions.
Inside the tangent cylinder, the Lorentz term is somewhat sub-dominant at large scales, where a balance between Coriolis and buoyancy terms describes well the largest scales and their slow evolution. Through alignment of velocity and magnetic field, the Lorentz force shapes the flow while keeping itself weak. At shorter time- and length-scales, the Lorentz term does eventually overcome buoyancy.
Outside the tangent cylinder, although the very weak mean flow is also controlled by a balance between Coriolis and buoyancy terms, the Lorentz force is more important. In fact, buoyancy plays a primary role at intermediate length-scales only (around m = 20). Note that this azimuthal wavenumber is much smaller than that of the onset of convection (mc = 67). A magnetostrophic balance, involving mostly Coriolis and Lorentz terms, controls the rapid dynamics of the flow outside the tangent cylinder.
3.7.3 Flow invariance along the rotation axis
Fig. 20 represents Gm within the tangent cylinder and outside the tangent cylinder. First we see that the large-scale flow in S1* have large values of Gm, indicating a high z-invariance of the flow, even within the tangent cylinder. This invariance is lost at smaller scales, when the local Rossby number Ro(m) = Ro m/δ > 0.1—where δ = ro for the outer region and ri for the inner region. When influenced by a strong magnetic field, the z-invariance is completely lost inside the tangent cylinder, while it is only lowered outside, all the more as the Lehnert number is decreased (compare results for S1 and S2, and see also Fig. 4), except for the axisymmetric mean flow. The quasi-geostrophic description (where the flow dynamics is treated as invariant along the rotation axis, e.g. Schaeffer & Cardin 2006; Gillet et al. 2011; Labbé et al. 2015) is thus relevant outside the tangent cylinder, even for timescales which would correspond to several decades or even centuries (when rescaled using Alfvén or advective timescales).
3.8 Torsional Alfvén waves
The high G values obtained for zonal flow at various frequencies together with a magnetic energy much larger than the kinetic energy suggest the presence of torsional waves. Fig. 21 shows a space–time diagram of z-averaged zonal flow. Outside the tangent cylinder, the propagation of waves from the tangent cylinder (s = 0.54) towards the equator (s = 1.54) is striking. There are also a few events that propagate in the other direction, which seem to be spawned in the fluid interior. Similar waves are seen in S1, although the signal is less clear (not shown). In S0, however, no such wave pattern appears. Interestingly, because the mean zonal flow is weak in our S1 and S2 simulations, we observe the torsional wave propagation outside the tangent cylinder without even removing the mean or employing any filtering. This contrasts with previous studies, where the mean flow was subtracted or filtered out (Wicht & Christensen 2010; Teed et al. 2014). We also note, that the frequency spectra of axisymmetric (m = 0) motions do not display significant peaks that would correspond to torsional modes. This absence of distinct peaks was also noted by Gillet et al. (2015) for core flows inverted from geomagnetic data. However, the measure of geostrophy Gm peaks at several frequencies that are compatible with such modes (see Fig. 20 bottom right).
The torsional waves seem to originate mostly near the tangent cylinder (s = 0.54). Referring to Fig. 18, we can see that for m = 0, the Lorentz force dominates the region close to the inner-core (named bli), whereas it stays weak in the vicinity of the tangent cylinder away from the boundary (region named tc). This hints to the Lorentz force close to the inner-core as the main source for torsional waves, which is in line with recent findings (Teed et al. 2015).
Inside the tangent cylinder, the signal is dominated by the variability of the strong twisted polar vortex system, and no waves are seen there. Note that with a calm TC, torsional waves are also expected there (as observed by Teed et al. 2015). Reflection of the waves at the equator depends on both the magnetic Prandtl number Pm (Schaeffer et al. 2012) and the mantle conductivity (Schaeffer & Jault 2016). For an insulating mantle, only weak reflection is expected even for Pm as low as 0.1. Indeed no obvious reflection is observed in our simulations.
3.9 Taylor constraint
Fig. 23 shows different views of this quantity |$\mathcal {T}$|. The space–time diagram is dominated by torsional waves outside the tangent cylinder, while inside the tangent cylinder |$\mathcal {T}$| is smaller. The time-average of the absolute value of |$\mathcal {T}$| shows this trend as well, and also shows that |$\langle |\mathcal {T}|\rangle$| decreases significantly from simulation S1 to S2. This may be linked to lower viscosity and stronger magnetic fields in S2. However, the maximum close to the equator, which scales as 1/h(s) (where |$h(s) = \sqrt{r_o^2 - s^2}$| is the height of the geostrophic cylinder) does not seem to change between S1 and S2. It is in fact not surprising that |$\mathcal {T}$| reaches 1 near the equator. There, the height of geostrophic cylinders becomes too small for the magnetic torque to have sufficient variations in z to cancel. Hence, the measure of |$\mathcal {T}$| is spoiled by geometric effects close to the equator, reaching values close to 1 independently of viscosity or inertia. Fortunately, when averaging over the volume, it has less effect. We find |$\mathcal {T}_{S1}=0.16$| and |$\mathcal {T}_{S2}=0.12$|. These numbers are significantly larger and decrease slower than those reported by Aubert et al. (2017) for stress-free boundary conditions. Considering that the Reynolds stress and viscosity in the bulk are negligible in S2, but that viscous forces balance magnetic forces near the outer boundary (see Fig. 18out and blo regions for m = 0), it suggests that viscous torques in the boundary layer (Ekman friction) dominates the torque balance on geostrophic cylinders.
The quantity |$\langle |\mathcal {T}|\rangle$| includes contribution from torsional waves and other time-dependent fluctuations, and as such is not representative of the long-term equilibrium. We checked that the time-averages of |$\mathcal {T}$| (sign included) for S1 and S2 do not tend to zero but rather to very similar profiles when averaged over various time spans. These averages |$\langle \mathcal {T}\rangle$| are shown for S0, S1 and S2 in Fig. 23 (bottom). The amplitude is significantly smaller than that of |$\langle |\mathcal {T}|\rangle$|, and of the same level inside and outside the tangent cylinder. The prominent peaks attached to the tangent cylinder scale as E1/3 (or E1/4, which are difficult to distinguish with one decade of E), pointing to an important role of viscosity there, compatible with a Stewartson layer.
4 DISCUSSION
We have reached the E = 10−7 milestone for a geodynamo simulation with vigorous convection. Our simulations do not significantly deviate from the scaling laws based on convective power input proposed by Christensen & Aubert (2006). However, the low viscosity associated with detailed analyses including force balances at different time- and length-scales brings a renewed picture of planetary dynamos.
The internal dynamics observed in the less viscous simulation S2 is summarized in Fig. 24. Importantly, global scale, columnar, non-zonal circulations emerge on intermediate timescales and are discussed in Section 4.3.2. The zonal mean flow is strongly suppressed outside the TC, as seen in Fig. 3, where it is replaced by torsional Alfvén waves (Fig. 21) triggered by Lorentz torque fluctuations at the ICB (Fig. 18). Inside the TC, the zonal flow is enhanced by the magnetic field and takes the form of a polar twisted vortex (Fig. 3), with important fluctuations in time (Fig. 15). Note also that the flow is by no means invariant along the rotation axis inside the TC, contrarily to the non-magnetic case. Torsional Alfvén waves are not visible in this very active inner region.
4.1 Inhomogeneities (in space and in time)
Our simulations point to strong spatial inhomogeneities (see Fig. 7) and to large temporal fluctuations of convective spherical dynamos (see Fig. 2). The TC separates regions with different dynamics. The inside of the TC is characterized by strong zonal flows and toroidal fields (see Fig. 3). Albeit occupying a small portion of the fluid volume, it bears a significant fraction of the total kinetic (about 1/10) and magnetic energies (about 1/4).
Our study also points to large heterogeneities outside the tangent cylinder, where regions of intense magnetic fields are observed next to regions of almost zero field (see Figs 6 and 7). This is linked to the magnetic field outside the tangent cylinder being mostly non-zonal (70 per cent of the magnetic energy) in contrast with the field at the core surface, which is mostly an axisymmetric dipole (Fig. 8). The heterogeneous magnetic field strength translates into heterogeneous length-scales for the convection. To quantify this effect, we compare some local quantities associated with the radial flow in the equatorial plane to the local magnetic field. We find that the most convincing correlations occur with the azimuthal field Bϕ, the component mostly perpendicular to the plume sheets. In Fig. 25, it is clear that the amplitude of the radial flow decreases by a factor of about 10 where the azimuthal field is strong, hence suppressing convection (Fig. 25 left). Meanwhile, the length-scale δϕ increases up to 3 times from the viscous scale E1/3 (Fig. 25 right). We therefore conclude that the primary effect of the magnetic field outside the TC is to damp convection. As a secondary effect, the length-scale increases. Note that in S1 these correlations are less spectacular (not shown).
Lastly, the fields are also heterogeneous in time. This aspect has been largely overlooked as databases for previous calculations often record time-averaged values. These large temporal fluctuations are even more intriguing as they are completely absent from our non-magnetic convection case (compare Fig. 17 with Fig. 16; see also Fig. 2). It appears that the dynamo regime is accompanied with global long-term correlation, otherwise absent from the hydrodynamics.
4.2 Comparison with the state of the art
At this stage, it is worthwhile to replace our findings in the context of the recent studies conducted by Yadav et al. (2016b), Sheyko et al. (2016) and Aubert et al. (2017), whose properties are recalled in Table 1.
4.2.1 Force balance
Yadav et al. (2016b) performed a systematic parameter survey down to an Ekman number equal to 10−6 and demonstrated convincingly that in their dynamo simulations, the first order force balance (once the zeroth-order geostrophic balance had been removed) took place between the Lorentz force, the buoyancy force and the part of the Coriolis force that is uncompensated by the pressure gradient. Yadav et al. (2016b) stressed that, in the bulk of the spherical shell, inertial forces and viscous forces represented less than 5 per cent of these three dominant actors of the dynamics (outside the boundary layers). This statement was based on an integrated analysis, that is, on the average rms values of the various forces at work in the simulations. Our simulation S2 refines this finding at more extreme conditions. Not surprisingly, our spatial-temporal analysis of the vorticity equation reveals that this three-term balance occurs in the bulk of the domain (‘in’, ‘tc’, and ‘out’ regions in Fig. 14), while the dynamics of the boundary layers is primarily governed by viscous and magnetic forces.
Inspection of the scale-dependent balance (space–time spectra of Fig. 18) shows that we must distinguish the regions inside and outside the TC. Inside the TC, at large scale (wavenumbers m ≲ 10) and low frequency (representative of the most energetic flow, see Fig. 15) the bulk three-term MAC balance is mostly a two-term balance between buoyancy and the uncompensated Coriolis force, representative of the thermal wind balance discussed above. At smaller scale or higher frequency, magnetic effects overcome buoyancy. To be precise, since we are dealing with the vorticity equation in our analysis, the curl of the Lorentz force becomes of comparable magnitude as the uncompensated Coriolis and larger than the buoyancy force (Fig. 19). Outside the TC, the picture is different. For the energetic small-scale convection plumes (20 ≲ m ≲ 70, see Fig. 15), a MAC balance is observed where Lorentz and buoyancy terms combine to balance the Coriolis term (Fig. 19). However, the Lorentz force does not contribute much to the large-scale, low-frequency force budget, which is again governed by a thermal wind balance. We also remark that outside TC, inertia (advection) is not completely negligible. In both regions, at frequencies higher than the Alfvén frequency (about 2 × 10−3 for S2), the buoyancy becomes insignificant and the Lorentz forces alone balance the Coriolis terms at all length-scales.
The analysis of Aubert et al. (2017) separates neither the regions inside and outside the TC, nor the contributions at various timescales. Instead their decomposition in harmonic degree ℓ mixes contributions from the two regions. To ease comparison with their work, we have transformed our dynamical balance analysis into a more global one to show averaged spatial or temporal dependence in Fig. 26. The global force balance analysis of (Aubert et al. 2017, their fig. 2b) at E = 3 × 10−6 shows that buoyancy does not enter the main balance at small scales (ℓ > 10), a behaviour not seen in our analysis in azimuthal wavenumber m, but maybe related to the frequency dependence in S2 (Fig. 26 bottom). In this integrated picture, the three main terms (Coriolis, Buoyancy and Lorentz) clearly dominate for 1 ≤ m ≲ 70, matching the MAC balance paradigm inside the TC. However, outside the TC, the stacking of Fig. 26 brings the terms closer together, suggesting Coriolis and Lorentz dominate, with significant contribution of both buoyancy and advection. It is thus important to separate time and space contributions to uncover the finer balance described in the previous paragraph.
4.2.2 Flow length-scale
In previous studies, a characteristic length-scale of convection is often defined by an average over the whole domain. Figs 6 and 7 show the wide variety of length-scales present in simulation S2. There are tiny plumes that form at the ICB (with measured width close to the E1/3 viscous length-scale) and progressively widen while traveling across the shell. These small plumes coexist with much larger ones (Fig. 6), and even with global-scale eddies (Fig. 4). This is most likely due to the heterogeneous magnetic field, with large regions of low field next to regions of stronger field (see also Fig. 25). Hence, even though the small viscous length-scale E1/3 is present in our simulations, the kinetic energy spectrum peaks at much larger scales (Fig. 12). In fact, the kinetic energy spectra peak at ℓ ∼ 10 for all our dynamo simulations, much like those of Aubert et al. (2017). A similar peak at m ≳ 10 appears in the spatial-temporal kinetic energy spectra outside the TC (Figs 15 and 16). Clearly, in this inhomogeneous, anisotropic system, reducing the convective motion to a single length-scale would be misleading.
Yadav et al. (2016b) pointed out the fact that when decreasing the Ekman number while maintaining a significant level of supercriticality (as large as allowed by computing resources), simulations of dynamo action driven by an imposed temperature drop across the shell were characterized by large-scale flows in the bulk of the shell. They ascribed this property to the emergence of the MAC balance (hence a stronger feedback of the Lorentz force on the flow). This result is in contrast with the previous work of Sakuraba & Roberts (2009), who found (at lower levels of supercriticality) that the emergence of large scales in the bulk of the domain was only possible if fixed-flux conditions were applied. Here we also rely on fixed-flux conditions; therefore, we cannot entirely corroborate the conclusions of Yadav et al. (2016b). That being said, our spatio-temporal analysis clearly demonstrates that the flow becomes weaker and larger scale when and where the magnetic field is strong (recall Figs 6 and 25).
4.2.3 Reversals
The simulations of Sheyko et al. (2016) are also in the rapidly rotating regime, with nominal values of the Ekman number as small as 2.4 × 10−6 when adopting our definition of E. The ratio of kinetic to magnetic energy in these calculations is always larger than 1: it varies between 3 and 30 depending on the parameters. Accordingly, the kinetic energy spectrum is always above the magnetic spectrum (their fig. 4a), and no sign of scale separation between the velocity and magnetic fields is to be seen, despite the small value of the magnetic Prandtl number (0.04−0.05). Because the magnetic energy is small, these dynamos can operate in a reversing mode driven by Parker waves (see also Busse & Simitev 2006; Schaeffer & Cardin 2006; Tobias & Cattaneo 2013; Schrinner et al. 2014). This cyclic behaviour is akin to that of the solar dynamo, whose polarity changes every 11 years. The field being weak, the emergence of a large scale flow in the bulk of the domain, as highlighted here when the field is strong and by Yadav et al. (2016b), does not seem to be happening. In addition, these simulations do not exhibit a strong dynamical contrast between the interior of the tangent cylinder and its exterior (our ‘in’ and ’out’ regions). Furthermore, we emphasize that, to our knowledge, no simulation with strong field (i.e. magnetic energy larger than kinetic energy) has ever exhibited polarity reversals. It is therefore questionable whether existing reversing simulations are relying on the same mechanism as the one governing the Earth’s core dynamo. Finally, we note that it would be interesting to study the role of the highly active polar regions for the dynamo action and for polarity reversals.
4.3 From numerical simulations to Earth’s core regime
Present-day numerical simulations, including the ones presented in this article, remain far from matching the physical parameters expected for the Earth core. However, some encouraging comparisons can be made that also suggest simulations can shed some light on the Earth’s core dynamics.
4.3.1 Invariance along the rotation axis
Assuming z-invariance of motions within the core has proven to be a fecund way of analysing the secular variation of the Earth’s magnetic field (e.g. Gillet et al. 2010, 2011). If this is characteristic of the rapidly rotating regime in which the Earth’s core operates, it is important to see if simulations reproduce this trend. Fig. 20 shows that the magnetic field reduces the z-invariance of the flow. However, outside of the tangent cylinder, the most energetic motions (compare Fig. 20 with Figs 16 and 15) are still mostly z-invariant. Shortest time- and length-scales, which are also less energetic, lose their z-invariance, as expected (Nataf & Schaeffer 2015). We note that the flow in S2 is more z-invariant than in S1, due to a lower influence of the magnetic field as measured by both the Elsasser and Lehnert numbers. Because the Earth’s core has a larger Elsasser but a smaller Lehnert number than S2, it is difficult to conclude whether the Earth’s core is globally more or less z-invariant than S2, but our simulations suggest that the large scale, slow motions can be well approximated by columnar flows outside the TC.
4.3.2 Global circulation and westward drift
A prominent feature of the geomagnetic field is the westward drift of the magnetic flux patches, especially near the equator and in the Atlantic hemisphere (e.g. Finlay & Jackson 2003). This westward drift has been associated with a large-scale, z-invariant, eccentric gyre (e.g. Pais & Jault 2008; Gillet et al. 2015; Pais et al. 2015).
Similarly, a long-term, large-scale, z-invariant circulation forms spontaneously in S2 (Figs 4 and 24) and is responsible for the westward drift of the magnetic field in one hemisphere (see Fig. 9). We emphasize that our setup and its boundary conditions have all a strict spherical symmetry (i.e. no heterogeneous thermal boundary conditions at the mantle or the inner-core). In contrast, the Coupled Earth simulations initiated by Aubert et al. (2013) produce these features using a heterogeneous heat flux at the boundaries together with a gravitational coupling between inner-core and mantle that must exceed the direct coupling at the CMB (see detailed analysis by Pichon et al. 2016). We remark however that even though the intense jet localized near the tangent cylinder in S1 and S2 (see Fig. 4) looks similar to the ones inferred from geomagnetic data (e.g. Pais & Jault 2008; Schaeffer & Pais 2011; Gillet et al. 2015; Pais et al. 2015; Livermore et al. 2017), it flows in the opposite direction. Nevertheless, it is interesting that a simple homogeneous model like S2 can reproduce features like the westward drift and a non-zonal global circulation.
4.3.3 Influence of the magnetic field
Our simulations indicate that the magnetic over kinetic energy ratio increases as we get closer to Earth’s parameters, in agreement with the scaling laws of Christensen & Aubert (2006). This is reassuring, but we get a maximum ratio of about 12 while it approaches 104 in the core.
Although the force balance may suggest a mild effect of the magnetic field on the mean flow, this is clearly not the case. As shown in Fig. 2, when the magnetic field is turned off (everything else being kept the same), the kinetic energy quickly increases by a factor 10. One also observes that the non-magnetic convection in S1* is more evenly distributed in the whole shell, with the formation of a strong zonal jet outside the tangent cylinder, instead of a strong twisted polar vortex in the dynamo case (see Figs 3 and 5). In the presence of magnetic fields, the flows are organized to limit induction, by aligning magnetic and velocity field wherever possible. This effect can be understood as a generalization of Ferraro’s law of isorotation (Aubert 2005). In this way, the strong magnetic field obtained in S1 and S2 completely changes the flow. Not only does it completely suppress the geostrophic jets, but the convective plumes extend further into the shell (because they are not blocked by the zonal jets anymore). The plumes are also of larger scale in the bulk, although they often emerge at a very small scale from the ICB, even smaller than without the magnetic field.
Inside the TC, the mean magnetic field intensity reaches a value that corresponds to the Elsasser number Λ > 1 (Fig. 3), and the mean field and flow are largely aligned, limiting the amount of induction. Of course, substantial induction is still occurring to sustain the strong field there (locally up to Λ ≳ 20, see Fig. 7), but much less than if the fields were unaligned. Indeed, focusing on the forces acting on the mean flow that are symmetric about the rotation axis, we find a dominant balance between the Coriolis and buoyancy forces (see Figs 18 and 19). This indicates that the zonal flows within the TC are very likely thermal winds driven by the density contrast between a light fluid inside the TC and a heavier one outside, which are nonetheless ‘shaped’ by the strong magnetic field. The polar vortices inferred within the Earth’s core (Olson & Aurnou 1999) suggest that this regime found inside the TC in our dynamo simulations is indeed realistic.
4.3.4 Density segregation
Interestingly, we observe that the establishment of a strong lateral density gradient requires intense magnetic fields, even though the latter do not directly enter the thermal wind balance (see Section 4.3.3). Furthermore, we remark that a stably stratified region appears within the TC, but only in the dynamo simulations (Figs 3 and A2). Analysis of the zonal force balance within the TC shows that the Lorentz force slightly helps the Coriolis force to balance the buoyancy force. It may thus be worth studying further the impact of magnetic forces on density separation and how they affect the ‘thermal’ wind balance (eq. 12) in planetary cores.
4.3.5 The role of viscosity
Despite reaching very low viscosities, our simulations are not devoid of viscous effects. They are of course important in the dynamics of the boundary layers, as seen in Fig. 18. This may result in torques that contribute to deviations from the Taylor state (see Section 3.9). But viscosity also appears to control two other aspects. First, the region of transition at the tangent cylinder seems to be controlled by viscosity, scaling with the Ekman number as a Stewartson layer (Fig. 23). Second, where the magnetic field is weak, convective plumes occur at the viscous scale E1/3 (see also Fig. 25). Whether these features are significant for the Earth’s core are open questions.
4.4 Perspectives and future work
Our direct simulations may serve as benchmarks for developing large eddy simulations schemes (see e.g. Buffett 2003; Nataf & Schaeffer 2015; Aubert et al. 2017) and for validating asymptotic models (e.g. Calkins et al. 2016), which are both needed to further progress towards realistic dynamos. Several terabytes of data will be held available for a few years to allow further analysis. It would be also instructive to characterize the turbulence and the dynamo mechanism. Calkins et al. (2016) have recently emphasized the significance of the ordinary Prandtl number in the selection of the dynamo mechanism in a rotating layer. We also have to investigate whether our results hold for Pr ≠ 1.
Because of the apparently important role played by the TC, it would be of great interest to revisit the role of a conducting inner-core in this parameter range. Indeed, it might enable the strong toroidal field to reach the boundary and alter the scale of the plumes originating there. More importantly, if the dynamics within the TC influence significantly the magnetic field generation, we should compare dynamos with and without inner-core, the latter having relevance for the early Earth or other planets. It may also be interesting to isolate the region within the tangent cylinder to study its own dynamics decoupled from the outer region.
Acknowledgments
The XSHELLS, SHTns and SINGE codes used in this paper are freely available at https://bitbucket.org/nschaeff. We thank Julien Aubert and the geodynamo team for fruitful discussions and comments. We also thank two anonymous referees for their constructive comments that helped improve this paper. We acknowledge PRACE and GENCI for awarding us access to resource Curie based in France at TGCC (grant pa1413 and t2014047258); GENCI awarded further access to Occigen (CINES) and Turing (IDRIS) under grants x2015047382 and x2016047382. Part of the computations were also performed on the Froggy platform of CIMENT (https://ciment.ujf-grenoble.fr), supported by the Rhône-Alpes region (CPER07_13 CIRA), OSUG@2020 LabEx (ANR10 LABX56) and Equip@Meso (ANR10 EQPX-29-01). Heavy data post-processing was partly performed on the S-CAPAD platform, IPGP, France. This work was funded by the French Agence Nationale de la Recherche under grants ANR-11-BS56-011 (AVSGeomag), ANR-13-BS06-0010 (TuDy) and ANR-14-CE33-0012 (MagLune). ISTerre is part of Labex OSUG@2020 (ANR10 LABX56). All figures were produced using matplotlib (http://matplotlib.org/) except Fig. 7 rendered with paraview (http://www.paraview.org/) and Fig. 24, drawn with inkscape (https://inkscape.org/). The contribution of AF is IPGP contribution 3859.
For our series of simulations, the computing resources needed to span a magnetic diffusion time increase by a factor of about 30 to 100 when the Ekman number is divided by 10.
Note also that supplementary meridional slices of the instantaneous velocity field can be seen at https://doi.org/10.6084/m9.figshare.5002199
REFERENCES
APPENDIX A: CODENSITY PROFILES
Time-averaged, total codensity profiles (C + C0) obtained in our simulations are represented in Fig. A2. The codensity difference across the shell is smallest in the equatorial plane, implying efficient convection takes place there.