The interplay of fast waves and slow convection in geodynamo simulations nearing Earth's core conditions
TThe interplay of fast waves and slow convection ingeodynamo simulations nearing Earth’s coreconditions
Julien Aubert and Nicolas Gillet Universit´e de Paris, Institut de physique du globe de Paris, CNRS, F-75005 Paris,France. CNRS, ISTerre, University of Grenoble Alpes, Grenoble, France
Summary
Ground observatory and satellite-based determinations of temporal variations in the geo-magnetic field probe a decadal to annual time scale range where Earth’s core slow, inertia-less convective motions and rapidly propagating, inertia-bearing hydromagnetic waves arein interplay. Here we numerically model and jointly investigate these two important featureswith the help of a geodynamo simulation that (to date) is the closest to the dynamical regimeof Earth’s core. This model also considerably enlarges the scope of a previous asymptoticscaling analysis, which in turn strengthens the relevance of the approach to describe Earth’score dynamics. Three classes of hydrodynamic and hydromagnetic waves are identified inthe model output, all with propagation velocity largely exceeding that of convective ad-vection: axisymmetric, geostrophic Alfv´en torsional waves, and non-axisymmetric, quasi-geostrophic Alfv´en and Rossby waves. The contribution of these waves to the geomagneticacceleration amounts to an enrichment and flattening of its energy density spectral profileat decadal time scales, thereby providing a constraint on the extent of the f − range ob-served in the geomagnetic frequency power spectrum. As the model approaches Earth’score conditions, this spectral broadening arises because the decreasing inertia allows forwaves at increasing frequencies. Through non-linear energy transfers with convection un-derlain by Lorentz stresses, these waves also extract an increasing amount of energy fromthe underlying convection as their key time scale decreases towards a realistic value. Theflow and magnetic acceleration energies carried by waves both linearly increase with theratio of the magnetic di ff usion time scale to the Alfv´en time scale, highlighting the domi-nance of Alfv´en waves in the signal and the stabilising control of magnetic dissipation atnon-axisymmetric scales. Extrapolation of the results to Earth’s core conditions supportsthe detectability of Alfv´en waves in geomagnetic observations, either as axisymmetric tor-sional oscillations or through the geomagnetic jerks caused by non-axisymmetric waves.In contrast, Rossby waves appear to be too fast and carry too little magnetic energy to bedetectable in geomagnetic acceleration signals of limited spatio-temporal resolution. Key words:
Dynamo: theories and simulations; satellite magnetics; Rapid time variations.
Preprint submitted to Geophysical Journal International 15 February 2021 a r X i v : . [ phy s i c s . g e o - ph ] F e b Introduction
The geomagnetic signal emanating from Earth’s convecting fluid outer core con-tains energy over a wide range of time scales, from the longest paleomagnetic fea-tures on hundreds of million years down to the shortest variations observed overyears and less at observatories and in satellite surveys. Throughout this range, timeseries of the palaeomagnetic dipole (Constable & Johnson, 2005; Panovska et al.,2013) as well as series of the field components at ground observatories (De Santiset al., 2003; Lesur et al., 2018) support the existence of several power-law ranges inthe geomagnetic frequency power spectrum. At centennial periods and longer, syn-thetic time series produced by numerical simulations of the geodynamo have beenused to relate the properties of this spectrum to key features of Earth’s core geody-namics and magnetohydrodynamics (Olson et al., 2012) and to calibrate stochasticevolution models describing these processes (Bu ff ett & Matsui, 2015; Meduri &Wicht, 2016). At shorter time scales, within the range of geomagnetic secular vari-ation that spans periods from centuries down to years, it has however been chal-lenging to follow the same kind of approach because of the traditional limitationsof geodynamo simulations.The situation has improved recently with a new generation of models (see a re-cent review in Wicht & Sanchez, 2019) working in a regime closer to the physicalconditions of the core. A key physical time scale in the secular variation range isthe core overturn time τ U = D / U ≈
130 years, estimated using the core thick-ness D = U =
17 km / yr.Comparison with numerical simulations (Bouligand et al., 2016; Aubert, 2018) hasrevealed that this time scale relates to a corner frequency in the geomagnetic fieldpower spectrum, and marks the start of a high-frequency range with a spectrum fol-lowing the power law P ∼ f − (with f the frequency) similar to that inferred fromobservatory time series (De Santis et al., 2003). This range may extend up to fre-quencies of about 1 yr − (Lesur et al., 2018), though the separation of the internallygenerated geomagnetic field from external sources is already di ffi cult at this point.At such high frequencies, fast hydromagnetic waves are expected to be present inthe geomagnetic signal. Recent advances in geomagnetic field modelling have forinstance enabled the retrieval of torsional Alfv´en waves in the core (Gillet et al.,2010, 2015), at a fundamental period close to 6 years. Other types of Alfv´en wavesalso appear to play a key role in the description of geomagnetic jerks (Aubert &Finlay, 2019). The main goal of this study is to explore the role of these waves instructuring the geomagnetic power spectrum at high frequencies. More generally,we wish to understand the interplay between fast hydromagnetic waves and theslower convection of Earth’s core. This is also essential for linking the rapid geo-magnetic signals now routinely observed with satellites to the physical propertiesof Earth’s core, and for the prospect of developing better geomagnetic predictionsat horizons of human interest. 2his problem may be viewed from a di ff erent angle by associating the time scalesrelevant to the geomagnetic variations to a hierarchy of force balances in the core. Ithas long remained di ffi cult to grasp the place taken by the Lorentz force in this hier-archy, because of the non-linear nature of this force and the self-sustained characterof the magnetic field. The situation has however been largely clarified recently, fol-lowing theoretical advances (Davidson, 2013; Calkins et al., 2015; Aurnou & King,2017; Calkins, 2018) and advancing numerical explorations of the geodynamo sim-ulation parameter space (Aubert et al., 2017; Schae ff er et al., 2017; Aubert, 2019;Schwaiger et al., 2019, 2021). Throughout this parameter space indeed, the rapidplanetary rotation implies that at the system scale, the leading-order equilibriumin amplitude occurs between the Coriolis and pressure forces. In most models, thisquasi-geostrophic (QG) equilibrium is perturbed at the next order with a triple MACbalance between the buoyancy, Lorentz and residual ageostrophic component of theCoriolis force. In conditions approaching those of Earth’s core, the MAC balanceis exactly satisfied at a rather large scale d ⊥ ∼ D /
10, where magnetic energy isfed into the system from the buoyancy force (Aubert, 2019). Writing the QG-MACbalance yields two other independent estimates of the core overturn time τ U (e.g.Davidson, 2013) τ U = DU ≈ ρ Ω d ⊥ g o C ≈ ρµ Ω d ⊥ B . (1)Using typical values for the rotation rate Ω = .
29 10 − s − , density ρ = / m ,magnetic permeability µ = π − H / m, core surface gravity g o =
10 m / s , con-vective density anomaly C = − kg / m (e.g. Jones, 2015; Aubert, 2020) andmagnetic field strength B = τ U = O (10 ) yr, consistent with the first direct estimate mentioned above. Throughthe QG-MAC balance, the core overturn time scale is therefore that of the slowconvective evolution of the velocity field, adjusting in an inertia-less manner to thecreation and advection of density anomaly and magnetic field structures. A rangeof length scales and amplitudes is obviously expected for these structures, suchthat a broader range of convective time scales is expected around τ U , which ap-pears open-ended on the long side but bound to decadal periods on the short side(Aubert, 2018).While recent advanced dynamo models largely preserve the dynamics described byearlier models on time scales of centuries and longer (Aubert et al., 2017; Aubert,2018), their strength and purpose is to render the interannual and decadal range oftime scales much shorter than τ U , where the dynamics can become inertial againand take the form of rapid responses to transient disruptions in the QG-MAC bal-ance. At conditions approaching those of Earth’s core, inertia indeed comes severalorders of magnitude below the MAC forces (Aubert, 2019), such that the flows andmagnetic signals associated with these disruptions are expected to remain smallrespectively to the convective flow and main magnetic field. This should naturallylead to quasi-linear perturbations in the governing equations and, in the presence ofrestoring forces, wave propagation. The rapid, inertia-bearing (I) waves that have3een theoretically investigated (e.g. Finlay, 2008) also mainly involve the magnetic(M) and Coriolis (C) forces (with buoyancy sustaining rapid waves only in contextsof stable stratification). At the axisymmetric level, the Coriolis force identicallyvanishes on axial cylindrical surfaces, leading to purely MI torsional waves thatrepresent a purely geostrophic case of Alfv´en waves. These have been identifiedin Earth’s core by Gillet et al. (2010, 2015) and have become increasingly clear innumerical simulations as the model parameters have become more realistic (Wicht& Christensen, 2010; Teed et al., 2014; Schae ff er et al., 2017; Aubert, 2018). Therelevant time scale for torsional waves is the Alfv´en time τ A = D √ ρµ B , (2)which, using the estimates quoted above, indeed satisfies τ A ≈ (cid:28) τ U in Earth’score. At the non-axisymmetric level, the Coriolis force in principle takes over theinertial force (Lehnert, 1954; Hide, 1966; Finlay et al., 2010; Labb´e et al., 2015),leading to slow MC modes that do not contribute to rapid dynamics. In interestingrecent developments, it has however been shown that the influence of the Corio-lis force can be mitigated in a variety of ways, leading to fast, quasi-geostrophic(axially columnar), non-axisymmetric dynamics. Starting their survey from low-frequency MC modes, Gerick et al. (2020) have identified modes approaching theAlfv´en frequency 1 /τ A from below as their radial complexity is increased. Theirkinetic to magnetic energy ratio also increases, meaning that the influence of iner-tia is gradually restored. In advanced numerical dynamo simulations that may beseen as an extreme case of spatial complexity, Aubert (2018) have highlighted aspatial segregation of the force balance. In regions where the magnetic field is rel-atively homogeneous, localised and nearly MI Alfv´en wave dynamics is observedwith frequencies also closely approaching the Alfv´en frequency from below. TheMC part of the balance remains confined close to regions of strong magnetic fieldheterogeneity. Despite the constraints set by the Coriolis force, MI Alfv´en wavestherefore appear to be possible at any spatial scale, but detecting these in the ge-omagnetic signal is a challenge as they evolve close to, or at kinetic to magneticenergy equipartition. Their magnetic signature is therefore nominally small respec-tively to the main magnetic field produced by convection, the energy of which ex-ceeds the total kinetic energy in Earth’s core by a factor A − = ( τ U /τ A ) ≈ Ω (on the order of an Earth day), CI inertial waves populate a discrete and densefrequency spectrum (Rieutord & Valdettaro, 1997). The issue of detectability maybe seen as even worse here as the kinetic to magnetic energy ratio of inertial wavesfurther increases with increasing frequency (Gerick et al., 2020), such that theyweakly interact with the magnetic field. There remains however the possibility of4low, quasi-geostrophic ’Rossby’ modes (Zhang et al., 2001; Busse et al., 2005;Canet et al., 2014) with a pulsation ω (cid:28) Ω approaching the Alfv´en frequency fromabove as the spatial complexity is increased, which enter the realm of interannualgeomagnetic signals and also carry higher amounts of magnetic energy. The secondgoal of this study is to use our models to better quantify the detectability of Alfv´enand slow Rossby waves in Earth’s core.We report here on a new numerical geodynamo model that enables a joint explo-ration of all these waves and convective features at an unprecedentedly realisticlevel of separation between the key time scales τ U , τ A and τ Ω . Being to date theclosest to Earth’s core conditions, this model also considerably enlarges the param-eter range over which an asymptotic scaling analysis relevant to Earth’s core can beperformed, which forms the third goal of this study. The manuscript is organised asfollows: section 2 presents the existing path theory and the new numerical model.Results are presented in section 3 and discussed in section 4. We use a standard numerical model for Boussinesq convection, thermochemicaldensity anomaly transport, and magnetic induction in the magnetohydrodynamicapproximation. Full details on the equation set, boundary conditions and numericalmethod can be found in Aubert et al. (2013, 2017) and Aubert (2018), where thespecific configuration is denoted as ’CE’ (Coupled Earth). The unknown fields arethe velocity field u , magnetic field B and density anomaly field C . The fluid domainis an electrically conducting and rotating spherical shell of thickness D = r o − r i and aspect ratio r i / r o = .
35 representing the Earth’s outer core. The shell takesplace between an inner core of radius r i , and a solid mantle between radii r o and r E = . r o (the surface of the Earth), both of which are conducting and electro-magnetically coupled to the fluid shell. The electrical conductivity of the inner coreis set to the same value σ as that of the outer core, while the mantle features a thinconducting layer at its base of thickness ∆ and conductivity σ m such that the ratioof conductances is ∆ σ m / D σ = − . The inner core is furthermore gravitationallycoupled to the mantle, and the three layers can present axial di ff erential rotationswith respect to each other, with the constant axial rotation of the ensemble definingthe planetary rotation vector Ω = Ω e z . Moments of inertia for these three layers re-spect their Earth-like proportions. A homogeneous mass anomaly flux F imposedat radius r = r i drives convection from below, and the mass anomaly flux vanishesat r = r o (neutral buoyancy beneath the outer surface). A volumetric buoyancy sinkis present in the volume to ensure mass conservation. On top of this homogeneousbuoyancy distribution, lateral heterogeneities in the mass anomaly flux are super-5mposed at r = r i and r = r o , respecting the ’CE’ setup. Stress-free mechanicalboundary conditions are also imposed on the fluid shell. This removes the need toresolve extremely thin viscous boundary layer that have a negligible influence onthe solution (Aubert et al., 2017). As in our previous work, all models presentedhere produced a self-sustained, dipole dominated magnetic field of Earth-like ge-ometry (Christensen et al., 2010) that did not reverse polarity during the integrationtime.We consider a spherical coordinate system ( r , θ, ϕ ) with unit vectors e r , e θ , e ϕ as-sociated to the cylindrical system ( s , z , ϕ ) of unit vectors e s , e z , e ϕ . The numer-ical implementation involves a decomposition of u , B and C in spherical har-monics up to degree and order (cid:96) max and a discretisation in the radial directionon a second-order finite-di ff erencing scheme with NR grid points. We use thespherical harmonics transform library SHTns (Schae ff er, 2013) freely availableat https://bitbucket.org /nschaeff/shtns . We use a second-order, semi-implicit time stepping scheme. The solution is approximated using smoothly ramp-ing hyperdi ff usivity applied on the velocity and density anomaly fields, but not onthe magnetic field that remains fully resolved. The functional form of hyperdi ff u-sivity involves a cut-o ff (cid:96) h =
30 below which hyperdi ff usivity is not applied, andan exponential ramping parameter q h such that ( ν e ff , κ e ff ) = ( ν, κ ) q (cid:96) − (cid:96) h h , describingthe increase of hyperdi ff usivity with spherical harmonic degree (cid:96) for (cid:96) ≥ (cid:96) h . Whencompared against fully resolved references (Aubert, 2019), the approximated solu-tions adequately preserve the QG-MAC balance and the large-scale morphology ofthe solution. They enable the computation of long temporal sequences that wouldotherwise not be feasible and are therefore tailored towards the present analysis ofdynamics in the time domain. Values of NR , (cid:96) max and q h used in our models arereported in Table 1. We recall here the four main parameters of the model, the flux-based Rayleigh,Ekman, Prandtl and magnetic Prandtl numbers: Ra F = g o F πρ Ω D , E = ν Ω D , Pr = νκ , Pm = νη . (3)Aside of the already introduced planetary rotation rate Ω , core thickness D , coresurface gravity field g o , core density ρ , bottom-driven mass anomaly flux F , theseexpressions involve the viscous, thermo-chemical and magnetic di ff usivities ν , κ and η (with η = /µσ ). Details on the correspondance between flux-based andcanonical Rayleigh numbers may be found in Christensen & Aubert (2006). Ourmain model case in this study uses the path theory (Aubert et al., 2017) that bridgesthe 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 single6 able 1Input parameters of numerical models (see text for definitions). All models have Pr = (cid:15) Pathposition(percent) runlength( τ U units) E Ra F Pm NR (cid:96) max q h −
29 332 3 10 − . − .
33 10 −
50 196 10 − − −
71 117 3 10 − . − . − + o ff -path 109 10 − . − + Pm- o ff -path 68 10 − . − + o ff -path 37 10 − − variable (cid:15) : Ra F ( (cid:15) ) = (cid:15) Ra F (1) , E ( (cid:15) ) = (cid:15) E (1) , Pr ( (cid:15) ) = , Pm ( (cid:15) ) = √ (cid:15) Pm (1) . (4)Here Ra F (1) = . − , E (1) = − and Pm (1) = . Pr = ff usivity of equalstrength is applied to the velocity and density anomaly fields. Previously (Aubertet al., 2017; Aubert, 2018), we have explored this path down to (cid:15) = .
33 10 − ,halfway in a logarithmic progression scale between the start ( (cid:15) =
1) and the Earth’score conditions located at (cid:15) = − . Here we mainly report on a model at (cid:15) = − ,at 71 percent of this path on the same logarithmic scale (from hereafter the 71pmodel). Table 1 reports on the input parameters of this main model case, and of theother cases used in this study. To date, the 71p model is the closest to Earth’s coreconditions that has been reached in a numerical geodynamo simulation, in partic-ular regarding the nearly inviscid behaviour achieved at large scales, as witnessedby the Ekman number E = − . This is of course made possible by the hyper-di ff usive approximation and the use of stress-free boundaries. For comparison, thefully resolved numerical simulation currently closest to core conditions is that ofSchae ff er et al. (2017), operating at E = − .In our previous analyses that stopped at 50 percent of the path (the so-called Mid-path model, Aubert et al., 2017; Aubert, 2018; Aubert & Finlay, 2019), it was foundthat many of the important diagnostics related to waves evolve rather subtly, andscale with weak powers of (cid:15) . To reach our main goals, it is therefore necessaryto cover a wide range of this parameter. This is the main reason that motivated adirect leap from 50 to 71 percent of the path here, rather than the progression in7alf-decades of (cid:15) that we favoured previously. The 71p model has been initialisedby taking a state from the Midpath model. Initial transients could almost be elimi-nated by applying a new approach. First we applied the path scaling laws (Aubertet al., 2017) to cast the amplitudes of u , B , and C towards predictions close to theiractual values at 71 percent of the path. The numerical mesh was then refined toachieve NR = (cid:96) max = ff usivity param-eters (cid:96) h =
30 and q h = .
14. After 9.5 overturn times τ U , the lateral resolution hasbeen increased to (cid:96) max = ff usivitystrength down to q h = .
09 for the rest of the computation (about 107.5 overturns).Three secondary models have also been derived from the Midpath model by chang-ing its Rayleigh number (label Ra + ) or its magnetic Prandtl number (labels Pm +/ -).These models complement the series computed along the path with a limited ex-ploration across the path around the Midpath position.The root-mean-squared velocity and magnetic field amplitudes U and B measuredin the shell over the course of the simulation give access to the overturn and Alfv´entime scales τ U = D / U and τ A = D √ ρµ/ B . These are presented in Table 2a asdimensionless ratios taken respectively to the two other important time scales in thesystem, the rotational time scale τ Ω = / Ω and the magnetic di ff usion time τ η = D /η . As we advance along the path towards Earth’s core conditions (or decrease (cid:15) ), the numerical models preserve an approximately constant value of the magneticReynolds number Rm = τ η /τ U ≈ τ η , τ U with the Alfv´en and rotational time scales. Time scale separationcan be measured using the Alfv´en number A = τ A /τ U , Lundquist number S = τ η /τ A , Rossby number Ro = τ Ω /τ U and Lehnert number λ = τ Ω /τ A . Along thelines introduced by Jault (2008), Gillet et al. (2011) and Aubert et al. (2017), theasymptotic regime relevant to the geodynamo may be formalised through the twoconditions λ (cid:28) A (cid:28)
1, which in turn imply Ro (cid:28) Rm (cid:29) S (cid:29)
1. In Aubert (2018) the entry into this regime has been located at 29 percentof the path, from the consideration of flow and magnetic acceleration signatures.The 71p model provides the largest and most realistic separation between the fourimportant time scales τ η , τ U , τ A , τ Ω that has been obtained to date in a geodynamosimulation. The unidimensional nature of the path however implies that all timescales are covarying, which may lead to di ffi culties in disentangling the exact e ff ectof each of these scales. For this reason, the three other secondary models takenacross the path close to Midpath conditions will be useful in the analysis presentedin section 3. The computation of the 71p model also provides an opportunity to update the scal-ing analysis initially presented in Aubert et al. (2017), which stopped at 50 percent8 able 2Output parameters, presented in a. as dimensionless ratios of time scales, and in b. asdimensional time scales obtained by setting a dimensional value to the magnetic di ff usiontime τ η (see text). Estimates for Earth’s core are also reported in panel b (taken from Gilletet al., 2010; Lhuillier et al., 2011; Christensen et al., 2012; Aubert et al., 2017; Aubert,2018; Aubert & Finlay, 2019). a.Case Ro = τ Ω τ U λ = τ Ω τ A Rm = τ η τ U A = τ A τ U S = τ η τ A τ τ η τ τ η .
26 10 − .
15 10 − .
12 10 − .
11 10 − Midpath* 2 .
40 10 − .
10 10 − .
06 10 − .
72 10 −
71p 4 .
31 10 − .
82 10 − .
89 10 − .
88 10 − + .
45 10 − .
16 10 − .
65 10 − .
78 10 − + Pm- 4 .
78 10 − .
10 10 − . − .
09 10 − + .
22 10 − .
37 10 − .
32 10 − .
62 10 − b.Case τ η (yr) τ U (yr) τ A (yr) 2 πτ Ω (days) τ τ + + Pm- 135000 123 19 134 419 12.350pPm + − ≈ ≈ ≈ ≈ Table 3Additional dimensionless outputs of the 71p model.
Ohmic fraction f ohm ff usion length scale d min / D .
82 10 − Taylorisation level on axial cylinders T . −
9f the path. For this purpose, here we list and report on the few additional time-averaged dimensionless outputs that are needed (see Aubert et al., 2017, for fulldefinitions and Table 3 for values): the fraction f ohm of convective power that isdissipated in Ohmic losses, the magnetic di ff usion length scale d min / D , or squareroot of the ratio between time-averaged magnetic energy and Ohmic losses, and thetime-averaged level of enforcement T of the Taylor constraint on axial cylindricalsurfaces in the shell. As done at previous positions along the path (Aubert et al.,2017; Aubert & Finlay, 2019), we also determine a scale-dependent representationof the force balance in the 71p model, where the root-mean-squared amplitude ofeach force is represented as a function of the spherical harmonic degree (cid:96) . Thisscale-dependent force balance tool has previously revealed its strength in decipher-ing the complex hierarchy of forces present in the system, including the non-trivialroles of hydrodynamic and magnetic pressure, and in associating the successivebalances coming along this hierarchy to physically grounded length scales (Aubertet al., 2017; Aubert, 2019; Schwaiger et al., 2019, 2021). For the purpose of geophysical applications, it has become common practice to castthe dimensionless outputs of the numerical model back into the dimensional world,in a way that is rationalised by the physical equilibria preserved along the path toEarth’s core (see e.g. Aubert, 2020, for a detailed discussion). Length scales are di-mensioned by setting D = ff usion time to τ η = ff usivity η = . / s at the mid-point of current es-timates (Aubert et al., 2017). The resulting dimensional time scales values are sum-marised in Table 2b. Together with the constant value Rm ≈ τ η ensures a geophysically realistic value τ U ≈
130 yr forthe convective overturn time. This also leads to geophysically realistic values fortwo other important time scales pertaining to convectively driven geomagnetic sig-nals, the secular variation and acceleration time scales τ ≈
400 yr and τ ≈
10 yr(see Lhuillier et al., 2011; Christensen et al., 2012; Aubert, 2018, for definitions anddiscussion). This implies that all models taken along the path are adjusted for anadequate simulation of slow inertialess convection and convection-driven geomag-netic variations in the decadal to secular range. Turning now to the time scalespertaining to fast hydromagnetic waves, the progress made by the 71p model overprevious e ff orts becomes tangible when we consider the dimensional value of theAlfv´en time scale τ A = . πτ Ω = . τ U , τ and τ . This is straightforwardly achieved by adjusting τ η in proportionof the magnetic Reynolds numbers achieved by these models. Model 50pRa + Pm-uses a combination of Rayleigh and magnetic Prandtl number such that Rm keeps avalue close to that of models along the path, so the same time basis τ η = + and 50pPm + is adjusted accordingly tothe ratio between their magnetic Reynolds number and that of model 50pRa + Pm-,leading to the respective time bases τ η = τ η = τ η and Rm , particularly concerning the possibility of high thermal andelectrical conductivity in the core (e.g. Pozzo et al., 2012).As in our previous studies (Aubert, 2018; Aubert & Finlay, 2019; Aubert, 2020),other outputs of the model are dimensioned by using the invariants of the parame-ter space path, particularly the stability of the magnetic Reynolds number and theQG-MAC force balance. The velocity field u and Alfv´en velocity B / √ ρµ are firstexpressed in a dimensionless manner in units of η/ D , to respectively give local val-ues of the magnetic Reynolds number Rm and Lundquist number S . The results arethen multiplied by the dimensional value of η/ D obtained from the choices madeabove. As we noted previously (Aubert, 2018), from B / √ ρµ a dimensional valueof B can be obtained, but this value is not realistic unless we are at the end of thepath. At 71 percent of the path, this procedure yields B = . B = B = C we take advantage ofthe preservation of the balance between buoyancy and Coriolis force along the path(Aubert, 2020). The dimensionless density anomaly field is therefore expressed inunits of ρ Ω η/ g o D and the result is then multiplied with Earth’s core dimensionalestimate for ρ Ω η/ g o D obtained with our choices for η and the other values providedin section 1. The physical time span integrated for the 71p model is 117 core overturns τ U (Ta-ble 2a), which corresponds to 14000 years of physical time, using the rescalingprocedures introduced in section 2.4. It should be emphasized that achieving thisrun length while simulating strongly scale-separated dynamical features comes ata sizeable numerical cost. The numerical time step of the computation has indeedbeen capped to about a hundredth of the rotational time τ Ω . This is needed becauseof the explicit treatment of the Coriolis force in the time stepping scheme. In the11ater part of the computation, the time step for instance corresponds to 0.3 hours ofphysical time. The 71p run required about 246 million numerical time steps and 15million single-core CPU hours spent over the course of 2 years. These figures high-light that the 71p model is an extreme computational endeavour, despite the use ofapproximations to reduce its spatial complexity. To tackle this challenge, signifi-cant code optimisations have been performed during the computation, to reach awall time per iteration of 0.06 seconds on 2560 cores of the AMD-Rome partitionof GENCI-TGCC in France.A number of output time series are extracted from the model in addition to thetime-averaged diagnostics introduced above. As time scales shorter than the nu-merical day 2 πτ Ω (11.8 days in the 71p model) contain a negligible amount ofenergy (Aubert, 2018), the temporal spacing between samples has been set to thisvalue. Time series of the core surface magnetic field B ( r o , θ, ϕ, t ) and velocity field u ( r o , θ, ϕ, t ) have been systematically recorded up to spherical harmonic degree 30.As magnetic time series in our earlier models (Aubert, 2018; Aubert & Finlay,2019) were only recorded up to degree 13, here we further truncate B ( r o , θ, ϕ, t )after degree 13 for consistency. Time derivatives ˙ u = ∂ u /∂ t (the flow accelera-tion), ˙ B (the magnetic variation) and ¨ B (the magnetic acceleration) are computedby finite-di ff erencing. We define time series for the energies E SV , E SA of the mag-netic variation and acceleration as E SV ( t ) = S (cid:90) S ˙ B ( r o , θ, ϕ, t ) d S , (5) E SA ( t ) = S (cid:90) S ¨ B ( r o , θ, ϕ, t ) d S . (6)Here S is the spherical surface at r = r o . These quantities can also be evaluatedat Earth’s surface S E located at r E = . r o , where B can be upward-continued byassuming an insulating mantle. As in our previous study (Aubert & Finlay, 2019),an Earth surface jerk energy time series is defined in a way that factors in the limitedresolution of geomagnetic observations, as a sliding energy di ff erence of magneticacceleration averaged over consecutive 3-years windows: E J ( t ) = S E (cid:90) S E (cid:18)(cid:104) ¨ B (cid:105) t + t − (cid:104) ¨ B (cid:105) tt − (cid:19) ( r E , θ, ϕ, t ) d S , (7)where the square brackets stand for time averaging. We define the energy associatedto the core surface flow acceleration ˙ u ( r o , t ) up to spherical harmonic degree 30 as K ( t ) = S (cid:90) S ˙ u ( r o , θ, ϕ, t ) d S . (8)We label as ˙ u S ( r o , θ, ϕ, t ) the equator-symmetric component of ˙ u , with the associ-ated energy K S . In the asymptotic regime reached beyond 29 percent of the path,12he large-scale flow acceleration is dominantly equator-symmetric (Aubert, 2018),as confirmed here by ratios K S / K increasing from 0.7 to 0.84 between 29 and 71percent of the path. We further denote as ˙ u ZS ϕ ( r o , θ, ϕ, t ) the axisymmetric (zonal)average of ˙ u S · e ϕ , and the corresponding energy as K ZS . By using a standard Thom-son multitaper method with concentration half-bandwidth W = / ∆ t (where ∆ t is the run length), we also decompose the Earth-surface magnetic accceleration¨ B ( r E , θ, ϕ ), the core surface flow accelerations ˙ u ( r o , θ, ϕ ) and ˙ u ZS ϕ ( r o , θ, ϕ ) in the fre-quency domain to calculate the respective energy density spectra ˆ E SA ( f ), ˆ K ( f ) andˆ K ZS ( f ). Given the duration of our runs (see Table 1) this ensures that the overturnfrequency 1 /τ U is well resolved in all model cases.For reasons of storage space, full three-dimensional outputs of the simulation stateat native resolution are available only in selected portions of the numerical com-putation. In these portions and as in Aubert (2018), we introduce the azimuthalacceleration on geostrophic cylinders (whose axis is aligned with the rotation axis e z ): ˙ u g ( s , t ) = π ( z + − z − ) (cid:90) π (cid:90) z + z − ˙ u ( s , ϕ, z , t ) · e ϕ d z d ϕ, (9)as well as the Alfv´en speed given by the part of poloidal magnetic field permeatingthese cylinders: c A ( s , t ) = π ( z + − z − ) ρµ (cid:90) π (cid:90) z + z − ( B ( s , ϕ, z , t ) · e s ) d z d ϕ. (10)Here z − and z + are the ordinates along e z of the intersections between geostrophiccylinders of radius s and the spherical boundaries of the shell. Our first task is to update our previous analysis of the path in parameter space(Aubert et al., 2017; Aubert, 2018; Aubert & Finlay, 2019) in the light of the newdata point acquired at 71 percent of this path. Figure 1 presents a scale-dependentforce balance diagram obtained with a snapshot of the 71p model. Similarly toour previous results at earlier positions along the path (Aubert et al., 2017; Aubert,2018; Aubert & Finlay, 2019), a leading-order QG balance is still observed betweenthe pressure and Coriolis forces. At the next order follows a MAC balance betweenthe buoyancy, Lorentz and ageostrophic Coriolis forces. The cross-over betweenthe Lorentz and buoyancy forces defines the harmonic degree (cid:96) ⊥ = π D / d ⊥ wherethe triple balance is exactly respected, with a value (cid:96) ⊥ =
13 that has only slightlyevolved since Midpath conditions where (cid:96) ⊥ =
12 (Aubert, 2019). This scale, at13 -1 -2 -3 -4 r e l a t i v e r . m . s . a m p li t ude -5 Pressure gradientCoriolis force Viscous force ageo s t r oph i c C o r i o li s Buoyancy forceInertiaLorentz force ⊥ Fig. 1. Root-mean-squared amplitude of the forces in a snapshot of the 71p model, pre-sented as functions of the spherical harmonic degree (cid:96) (see section 2.3). All forces arenormalised respectively to the maximum reached by the pressure gradient. The snapshotwas taken from the later part of the run where (cid:96) max =
170 and q h = .
09 (see Fig. 3). which convective energy is injected into the magnetic field, is therefore confirmedto remain stable along the path. At the second order in amplitude, inertial forces arenow located two orders of magnitude below the MAC forces. As previously shownin Schwaiger et al. (2019), the widening of the gap between MAC and inertial forcesalong the path scales with the inverse squared Alfv´en number A − , in line with theevolution observed between Midpath (see Fig. 1 in Aubert, 2018) and 71 percent ofpath conditions. In the 71p model, viscosity comes five order of magnitude belowthe Coriolis force at large scales. Our choice for spatial resolution and strength ofhyperdi ff usivity in the later part of the run (from where the snapshot was taken)result in a negligible influence of viscosity up to spherical harmonic (cid:96) ≈
80, afterwhich its strength ramps up and exceeds that of inertia, as was previously the casein the Midpath model.The stability of the force balance implies that previously determined scaling lawsfor the main model outputs are essentially unchanged when adding the 71p modeldata. We have previously shown (Aubert et al., 2017) that along the path and for thehyperdi ff usive, large-eddy (LES) simulations, the QG-MAC balance, together withthe constancy of the magnetic Reynolds number Rm imply the following scalinglaws for the Rossby number and ohmic fraction-corrected Lehnert number: Ro ∼ (cid:15) / , (11) λ/ f / ∼ (cid:15) / . (12)The way the path models adhere to these scalings is essentially unchanged whenadding the 71p model data points (Fig. 2a,b). From this, the scaling A = Ro /λ ∼ (cid:15) / o DNSLESDNS fit, ε LES fit, ε λ / f oh m / A ε f oh m d m i n -4 -3 -2 -3 -2 DNS fit, ε LES fit, ε -1 DNS fit, ε LES fit, ε -5 -4 -3 -2 -1 -5 -4 -3 -2 -1 DNS fit, ε LES fit, ε -5 -4 -3 -2 -1 -2 -1 DNS fit, ε LES fit, ε a b cfed ε ε Fig. 2. Evolution of the Rossby number Ro (a), Lehnert number λ (b), Alfv´en number A (c),Ohmic fraction f ohm (d), magnetic dissipation length scale d min (e) and Taylor constraintenforcement level T (f) along the parameter space path, as functions of the path parame-ter (cid:15) (Earth’s core conditions are towards the left of the graphs). Shown are results fromhyperdi ff usive solutions along the path (LES or large-eddy-simulations) up to 71 percentof the path ( (cid:15) = − ), and the direct numerical simulations (DNS) up to 21 percent ofthe path ( (cid:15) = .
33 10 − ). See Tables 2a and 3 for data between 29 and 71 percent of thepath, and Aubert et al. (2017); Aubert (2019) for other numerical values. Green and orangelines respectively show the results of least-squares power-law fits to the LES and DNS data,which closely adhere to the theories of Aubert et al. (2017); Davidson (2013), respectively.Shaded grey areas represent the ± is also respected (Fig. 2c). Relative to the Midpath model ( (cid:15) = .
33 10 − ), theOhmic fraction has increased at 71 percent of the path to reach f ohm = .
89 (Fig.2d), meaning that despite the use of hyperdi ff usivity, dissipation is essentially ofOhmic nature in this model. The magnetic dissipation length scale d min remainsessentially invariant along the path (Fig. 2e, note the linear ordinate axis), a con-sequence of the constant magnetic Reynolds number and constrained spatial struc-ture of the solution. Finally, the level to which the solution approaches a Taylorstate (Taylor, 1963) can be measured by the Taylorisation level T (Fig. 2f), withthe 71p model exhibiting the highest levels of taylorisation observed to date alongthe path (lowest value of T ). Somewhat surprisingly, this model appears to respectthe Taylor constraint even more strongly than the expectation that could be drawnfrom the previous models. For reference, the scaling results from fully resolved di-rect numerical simulations (DNS) are also reported on Fig. 2. These simulationshave so far been computed only down to 21 percent of the path ( (cid:15) = .
33 10 − ,Aubert, 2019). In this previous study, we have shown that these adhere to the powerlaws predicted by the QG-MAC theory (Davidson, 2013; Aubert et al., 2017) i.e.15 ig. 3. Magnetic variation and acceleration properties of the 71p model. a,c: Hammerprojections of the time-averaged, root-mean squared radial magnetic variation ˙ B · e r (a)and acceleration ¨ B · e r (c) at the core surface, from magnetic field output up to sphericalharmonic degree 13. b,d: Time series of the magnetic variation and acceleration energies E SV (b) and E SA (d) at the core surface. Red lines in (b,d) mark the time-averaged values[ E SV ] and [ E SA ]. The blue and green bands between (b) and (d) locate the notable eventsrelated to spatial and temporal resolution changes in the model computation. Ro ∼ (cid:15) / , λ/ f / ∼ (cid:15) / , d min ∼ (cid:15) / . The di ff erence between path and DNS the-ories scalings essentially stems from the fact that the length scales are constrainedby hyperdi ff usivity in the former situation, while they are weakly evolving in thelatter case.In Fig. 3 we next turn to key aspects of magnetic variation and acceleration inthe 71p model, which we analyse in a manner similar to Fig. 3 of Aubert (2018).Consistently with the defining principles of the path, and with the approximatelyconstant values of τ U , τ observed along this path (Table 2), the properties of mag-netic variation at 71 percent of the path (Fig. 3a,b) remain essentially similar tothose observed since the start of the path, with an Earth-like geographic localisa-tion that is essentially dictated by the heterogeneous mass anomaly flux imposedat the inner boundary (Aubert et al., 2013). The pattern of magnetic acceleration isalso similar to that previously obtained from 29 percent of the path onwards, withthis acceleration being mainly observed within an equatorial band in the Easternhemisphere (0 o E − o E). The 71p model confirms that the relative dominanceof equatorial over polar magnetic acceleration continues to increase as we advancealong the path. As was the case in previous models, the 71p model features inter-mittent pulses in the magnetic acceleration energy, with a clear increase in strengthand frequency relative to the Midpath conditions. Relatively to the midpath Model,the time-averaged acceleration energy [ E SA ] of the 71p model (red line in Fig. 3d)consequently increases by about 25 percent, while [ E SV ] is almost unchanged (Fig.16 ig. 4. a: Time series of the Earth surface jerk energy E J (see section 2.5). b: Distributionof jerk energy as a function of the jerk recurrence time, for models previously computedin the range 29 percent to 50 percent of the parameter space path (Aubert & Finlay, 2019)and the 71p model. Shown in dashed line is the expected distribution at 71 percent of thepath from an extrapolation E J ∼ (cid:15) − . previously obtained from the behaviours between29 and 50 percent of the path (Aubert & Finlay, 2019). c: Distribution of jerk energy inthe 71p model (sampled each year) as a function of the latitude of the amplitude maximaobserved in the radial magnetic acceleration at the core surface. Jerk statistics in b,c arecomputed in the time window [4200 yr , τ ≈ √ [ E SV ] / [ E SA ] = . . τ A . This leads to geomagnetic signatures that reproducethe characteristics of recent geomagnetic jerks (Aubert & Finlay, 2019). To char-acterise the evolving properties of magnetic acceleration pulses and jerks along the17arameter space path, we previously found useful to analyse time series of the jerkenergy E J (Fig. 4a and equation 7). This quantity mirrors the occurence of acceler-ation pulses in a way that is more straightforwardly comparable to satellite-baseddeterminations of the geomagnetic acceleration. New to the 71p model is the occur-rence of extremely fast magnetic acceleration pulses (for instance at times 11070and 11740 in Fig. 3d) that are smoothed by the 3-year averaging process involvedin the definition of E J (see corresponding times in Fig. 4a). Indeed, while thisobservation-based averaging window was clearly shorter than the typical Alfv´entime over which the jerks develop in the models at 29 to 50 percent of the path,in the 71p model it is now only slightly shorter than τ A = . E J and representing the resulting estimates for jerk recurrencetime as a function of E J . This representation confirms that the 71p model producessignificantly more energetic jerks than the Midpath model at any recurrence time.However, the energy levels that are obtained are slightly recessed relative to theextrapolation that could be made in Aubert & Finlay (2019) from models span-ning 29 to 50 percent of the path. This obviously relates to the smoothing e ff ectdescribed above. It is finally useful to analyse the distribution of E J with the lati-tude of maxima in the associated magnetic acceleration pulses at the core surface(Fig. 4c). Consistently with the root-mean-squared pattern seen in Fig. 3c, magneticacceleration pulses occur preferably at high latitudes (near the projection on thecore-mantle boundary of the geostrophic cylinder tangent to the inner core, here-after called the tangent cylinder) and the equator, with equatorial dominance for thestrongest events. The residual equatorial asymmetry observed in Fig. 4c for strongevents links with the rarity of such events and the limited duration of the sequence.In the 71p model, about 80 percent of the events of magnitude E J ≥
100 nT / yr have core surface foci at latitudes less than 20 degrees away from the equator. The71p model therefore pursues the trend previously reported in Aubert (2018) towardsequatorial localisation of pulses and jerks as we advance along the path. In Fig. 5 we present the energy density spectra ˆ E SA ( f ) and ˆ K ( f ) of the Earth surfacemagnetic acceleration energy and core surface flow acceleration energy for mod-els located between 29 and 71 percent of the path as well as across the path at 50percent. At frequencies up to that of the overturn, f U = /τ U = U / D , Fig. 5a-d con-firms that all models feature similar spectral energy densities, an indication of theinvariant dynamics caused by slow secular convection. At frequencies beyond f U ,all models also feature a range where the energy density profile is flattened. Consid-ering that each successive time derivative multiplies the energy density by f , thisflat range corresponds to a f − range in the energy spectral density of the observedgeomagnetic variation, and to a f − range in that of the magnetic energy, both of18 ig. 5. Frequency domain spectral density ˆ E SA ( f ) of the magnetic acceleration energy atEarth’s surface (a,c,e), and ˆ K ( f ) of the flow acceleration energy at the core surface (b,d,f),for models taken along (a,b) and across (c,d) the parameter space path. a,b: spectra of the 29percent, 50 percent and 71 percent of path models, plotted together in a. with the additionalcontribution of core surface magnetic di ff usion in the 71p model. b,d: spectra of the 50percent of path, 50pRa + Pm-, 50pRa + and 50pPm + models. In a-d the dashed vertical linemarks the overturn frequency f U = /τ U . Flat and slanted grey lines respectively markthe power-laws f and f − , the crossing of which defines the cut-o ff frequency f c . Therelative positions and numerical values of the time scales τ η , τ U and τ A are also recalled.e,f: respective compounds of a,c and b,d, with frequency axis rescaled for each model by( S / S ) − / , where S = τ η /τ A and S is the Lundquist number of the Midpath model. which have received reasonable observational support (De Santis et al., 2003; Lesuret al., 2018). The flat range for flow acceleration corresponds to a f − range for flowvelocity, which has previously been assumed as part of stochastic core flow mod-elling strategies (Gillet et al., 2015, 2019). In the vicinity of the overturn frequency f U , it is logical to associate this range of approximately flat spectra energy densityto the signature of convective motions. At higher frequencies towards the decadalrange, decay spectral tails following a f − power law are observed both in flow andmagnetic acceleration energies. A cut-o ff frequency f c that marks the start of thedecay range is defined from the crossing between the flat ( f ) and slanted ( f − )grey lines in Fig. 5a-d. The decay spectral tails can no longer be associated with19onvection because they are not invariant as we progress along the path (Fig. 5a,b),as now better illustrated with the new data brought by the 71p model and the asso-ciated increase of the cut-o ff frequency f c . The combination of models along (Fig.5a,b) and across the path (Fig. 5c,d) enables us to pinpoint the main control on f c as well as on the energy density within the decay range. The influence of the rota-tional time scale τ Ω can be immediately ruled out because the spectra of cross-pathmodels with constant Ro or λ (see table 2a) do not superimpose. The key time scaletherefore appears to be the Alfv´en time τ A , but the question remains whether itshould be considered in proportion of the overturn time τ U (i.e by using the Alfv´ennumber A as a control) or magnetic di ff usion time τ η (i.e by using S ).A least-squares fit performed on the f c values extracted from Figs. 5a-d yields f c ∼ S . ± . , indicating that all spectra can be collapsed by rescaling the frequencyaxis in proportion of S − / , as further verified in Fig. 5e,f. The alternative scaling f c ∼ A − / can be readily discarded as the Alfv´en numbers for the two cases 50pRa + and 50pPm + with magnetic Reynolds numbers Rm = τ η /τ U ≥ Rm ≈ f U is constant prior to re-scaling the frequency axis, thescaling for the cut-o ff frequency hence writes f c ≈ f U ( S / S ) / , (13)with S ≈
900 from our results. The cut-o ff frequency f c and f − decay tails aresimilar for both the flow and magnetic acceleration energy spectra. In the decayrange f ≥ f c , the functional form of these energy densities is thereforeˆ E SA ( f ) ≈ E maxSA (cid:32) ff c (cid:33) − ≈ E maxSA f U S S f − , (14)ˆ K ( f ) ≈ K max (cid:32) ff c (cid:33) − ≈ K max f U S S f − . (15)Here E maxSA , K max are the values of the plateaus reached in the flat energy densityrange. Within the decay tails, the flow and magnetic energy densities therefore lin-early increase with the number S of Alfv´en wave periods contained within a mag-netic di ff usion time. This strongly suggests that these decay tails f ≥ f c are domi-nated by the contribution from Alfv´en waves. The range f ≤ f U is the convectiverange, and the interval [ f U , f c ] where energy density presents a plateau is the rangeof interplay between waves and convection. As the Lundquist number S increases,the plateau broadens because of the elevated contribution from Alfv´en waves at itsrightmost edge, with S in equation (13) representing the minimal Lundquist num-ber needed for this e ff ect to occur, or more generally for waves to be significant inthe simulation. Writing the time-derivative of the magnetic induction equation¨ B = ∇ × ( ˙ u × B ) + ∇ × (cid:16) u × ˙ B (cid:17) + η ∆ ˙ B , (16)20he similar spectral shapes (14,15) also indicate that in the wave range at f ≥ f c ,the first term in the right-hand-side dominates the production of magnetic acceler-ation, as previously anticipated for rapid dynamics (Lesur et al., 2010; Christensenet al., 2012; Aubert, 2018). A direct calculation of the contribution by core sur-face magnetic di ff usion to ˆ E SA ( f ) (Fig. 5a) confirms that the term η ∆ ˙ B in (16) issubdominant. The di ff usive e ff ects observed through the control of S on ˆ E SA ( f )and ˆ K ( f ) therefore originate in the bulk of the core. Incidentally, the functionalforms (14,15) may also be rewritten by involving the ratio [( U / f ) /δ ] , where U / f is a length scale excited by convection at frequency f and δ = √ ητ A is the bulkmagnetic di ff usion length scale at the Alfv´en time. We can pursue our analysis by examining either the flow or magnetic accelera-tion signal carried by the waves. The former is more convenient, because rapidazimuthal flow acceleration presents a highly columnar (and therefore equator-symmetric) structure at advanced positions along the path. Torsional waves are astraightforward and interesting case study of the Alfv´en waves present in the sys-tem, because they can be easily isolated in the axially-columnar, axisymmetric,azimuthal part ˙ u g of ˙ u defined by equation (9). Fig. 6a,b shows time-cylindricalradius representation of ˙ u g during a 100-yr sequence of the 71p model. As in Fig.10 of Aubert (2018), the presence and dominance of torsional waves in this signalis attested using the similarity of the patterns with ray-tracing tracks that representpropagation at the theoretical speed c A ( s , t ) (equation 10) commensurate with theone-dimensional Alfv´en velocity D / √ τ A ≈
225 km / yr. Compared to Fig. 10 ofAubert (2018), the signal outside the tangent cylinder (Fig. 6b) has logically in-creased in propagation speed (because of the decrease of τ A from Midpath to 71pconditions), but less trivially also in amplitude. The waves in the northern part ofthe tangent cylinder (Fig. 6a) are even faster (as they sample a stronger magneticfield), and also more intense. As the model features a jump in c A of about 15 per-cent across the tangent cylinder at s = ffi cientrelative to Midpath conditions (Fig. 10 of Aubert, 2018), where reflection was elu-sive. Given that our simulations operate with stress-free boundary conditions, thisresults appears to be at variance with the plane layer theory laid out in Schae ff eret al. (2012); Schae ff er & Jault (2016); Gillet et al. (2017), which predicts that thereflection coe ffi cient should not depend on the imput parameters varied along thepath (most notably Pm ), and should even decrease due to the presence of a con-ducting outer layer and the increasing Lundquist number S . We anticipate that the21 ig. 6. Time-cylindrical radius diagrams of the geostrophic flow acceleration ˙ u g (equation9) sampling the interior of the fluid in the Northern part of the tangent cylinder (a), andthe region outside the tangent cylinder (b), at native resolution in the 71p model. c: time–cylindrical radius diagram obtained at the core surface and outside the tangent cylinder byrepresenting the equator-symmetric, axisymmetric azimuthal flow acceleration ˙ u ZS ϕ up tospherical harmonic degree 30. Note the larger extremes of the color map in a. relativelyto b,c. Green curves in a.,b. correspond to outwards and inwards ray-tracing theoreticalpropagation tracks at the column-averaged Alfv´en speed c A ( s , t ) (equation 10). discrepancy ties with the treatment of the singularity present for torsional waves atthe equator of the core-mantle boundary, which is probably a ff ected by hyperdi ff u-sivity as the wave pattern shrinks down to small latitudinal length scales. The 71pmodel is less a ff ected by this problem than the Midpath model as it operates withhigher native resolution and reduced hyperdi ff usivity (see table 1). With stress-freeboundaries and the conducting layer, the plane layer theory predicts a reflectioncoe ffi cient R = (1 − Q ) / (1 + Q ) (Schae ff er & Jault, 2016), where the quality fac-tor writes Q = S ( ∆ σ m / D σ )( B r ( r o ) / B ). The conductance ratio is ∆ σ m / D σ = − for all models along the path, and the ratio B r ( r o ) / B ≈ / ig. 7. Frequency domain spectral density ˆ K ZS ( f ) of the core surface azimuthal, equa-tor-symmetric and axisymmetric acceleration energy for models taken (a) along and (b)across the parameter space path. a: spectra of the start, 29 percent, 50 percent and 71 per-cent of path models. b: spectra of the 50 percent of path, 50pRa + Pm-, 50pRa + and 50pPm + models. The numerical values of the time scales τ U and τ A are also recalled in a,b. In a., thearrows locate the frequency of the fundamental Alfv´en wave f = ( √ r o τ A ) − for modelsalong the path, and the slanted grey line marks the power law f − . c: compound of a.,b.,with frequency axis rescaled for each model by A / A , where A = τ A /τ U and A is theAlfv´en number of the Midpath model. therefore predicts Q = . R = .
5, qualitatively in line with Fig. 6b.At the conditions of low inertia and viscosity where our models operate, torsionalwaves are excited by perturbations of the Lorentz force averaged over axial cylin-ders, with these perturbations being linked to the underlying convection (Teed et al.,2014, 2019). This process naturally leads to excitation of large-scale nature. Thisexplains why almost all of the signal present in ˙ u g and evaluated from bulk dataat native resolution (Fig. 6b) can be retrieved in a time-cylindrical radius diagramof ˙ u ZS ϕ (Fig. 6c), the equator-symmetric zonal flow acceleration at the core surfaceup to spherical harmonic degree 30. In Fig. 7 we further analyse the associated23nergy density spectra ˆ K ZS ( f ) in the frequency domain. Beyond 29 percent of thepath (Fig. 7a), most of the signal corresponds to the contribution from torsionalwaves, with the invariant contribution from axisymmetric thermal winds (estimatedby representing the start of path model in Fig. 7a) being subdominant. In these ad-vanced models, ˆ K ZS ( f ) features a plateau centered around the fundamental wavefrequency f = ( √ r o τ A ) − (the factor √ f − decay associated with the limited length scale content of the excitation source,and much steeper than that observed for the full core surface flow in Fig. 5b. Theevolution of this spectral shape along the path reflects a wave content that remainsat large and invariant spatial scales while moving towards higher frequencies as theAlfv´en velocity is increased. Using cross-path models (Fig. 7b) again together withalong-path cases, the cut-o ff frequency is indeed confirmed to extend linearly withthe inverse Alfv´en number 1 / A = τ U /τ A (Fig. 7c), and a possible scaling with S canbe discarded. The influence of magnetic di ff usion is marginal here, because of thelarge length scales at which the waves are excited, and because radial core surfacedi ff usion is subdominant (recall Fig. 5a). E ff ects of di ff usion are only seen in thedecay tails in Fig. 7c, with the less di ff usive models 50pPm + and 50pRa + showinga slightly less steep decay than the other models. The simple form of the spectraimplies that the time averaged energy [ K ZS ] (the integral of the profiles along thefrequency axis) scales as [ K ZS ] ∼ A − . (17)This means that the waves are able to draw an increasing amount of energy from theunderlying convection as the separation between τ A and τ U increases. The leadingcontribution to this increase of [ K ZS ] comes from the additional frequencies that aremade available to wave motion as the fundamental Alfv´en frequency f is increased.Resonant forcing by convection (e.g. Teed et al., 2019) does not appear to be thedominant excitation mechanism, because the gap widens along the path betweenthe low and constant frequencies of thermal winds (again represented by the start ofpath model in Fig. 7a) and the increasing frequency f . Significant energy levels arehowever preserved in the intermediate frequency range separating the two, wheremotion is neither the consequence of convection nor of waves, but rather resultsfrom a non-linear interaction between the two, with the dominant non-linearity forenergy transfers being the Lorentz force in this case (Teed et al., 2014; Aubert et al.,2017). The observation of a broadening flat energy density profile in this interactionrange is therefore supportive of a non-linear transfer of energy from convection towaves underlain by the Lorentz stresses. Though the stabilising control of magneticdissipation is not present at the axisymmetric scale, as we have shown above, itis expected to indirectly come from the non-axisymmetric contributions to theseLorentz stresses. 24 ig. 8. State of the 71p model at time 4749.8 yr of the simulation. Core surface Hammerprojections (a,c) and partial equatorial plane views (b,d) of the azimuthal flow accelera-tion. In a,b, the total azimuthal acceleration ˙ u ϕ is presented at native spatial and temporalresolution. In c,d, the fast and large-scale component of the equator-symmetric azimuthalacceleration ˙ u S ϕ is presented. The fast component is obtained after removal of a 5-year run-ning average. The large-scale component is obtained in c. by truncation after sphericalharmonic degree 30, and in d. by lateral truncation after azimuthal wavenumber 30 andapplication of a running average of thickness 40 km in the radial direction. For clarity also,accelerations larger than 4 km . yr − are not represented in d. The partial equatorial views inb,d, are taken between the grey meridians delineated in a,c. In these equatorial views also,the density anomaly (green-brown, maximal values ± . − kg / m ) and radial magneticfield (orange-purple, maximal values ± . Searching for non-axisymmetric waves is intrinsically more di ffi cult as the signa-ture of convection is stronger than for axisymmetric motion, and can attain shorttime scales at small spatial scales (see equation 1). The wave signature can there-fore not be straightforwardly disentangled from that of convection, for instance innative resolution snapshots of the azimuthal flow acceleration ˙ u ϕ at the core sur-face (Fig. 8a) and within the equatorial plane (Fig. 8b). Focusing on the large-scalecomponent (up to degree 30), where the acceleration is overwhelmingly equator-symmetric i.e. ˙ u ϕ ≈ ˙ u S ϕ (see section 2.5) clarifies the picture by partly removing the25onvective signal originating at small spatial scales. For the purpose of illustration,in Fig. 8c,d we further consider the rapid part of the large-scale signal, obtainedafter removal of a 5-year running average. This processing highlights two otherclasses of rapidly propagating features (located by labels in Fig. 8c,d) in addition tothe already described axisymmetric torsional waves: transverse quasi-geostrophicAlfv´en (QGA) waves and longitudinal Rossby waves. Non-axisymmetric, axially columnar QGA waves propagate along magnetic fieldlines of arbitrary orientation perpendicular to the rotation axis. The typical distri-bution of field lines in the model (see e.g. Fig. 8 in Aubert, 2019) mainly promotescylindrical-radially propagating waves in the bulk of the fluid that are carried bythe azimuthal flow u ϕ (similarly to torsional waves). When tracked in the equato-rial plane (Fig. 8d), these are laterally limited by concentrations of radial magneticfield, as previously described in Aubert (2018). Near concentrations of azimuthalmagnetic field found at low latitudes, the propagation direction can also acquire anazimuthal component (see e.g. Fig. 8d near America). While the QGA waves werepreviously mostly exhibited in the vicinity of magnetic acceleration pulse events,they are now significantly stronger in the 71p model (see section 3.5) and we canconfirm that they are ubiquitous and constantly sent out by the deep convection.To better characterise the propagation of QGA wave features at the core surfacein the 71p model, in Fig. 9 we now specifically focus on the non-axisymmetric,azimuthal and equator-symmetric flow acceleration ˙ u S ϕ − ˙ u ZS ϕ and examine time-cylindrical radius diagrams at a specific longitude delineated by an arrow in Fig.8c,d. Representing this quantity at native spatial and temporal resolution (Fig. 9a)mainly reveals the outwards advection of flow by convection at speeds commen-surate with the one-dimensional convective velocity D / √ τ U =
11 km / yr. Onthis representation, however, faint signals propagating outwards at higher speedscan already be seen. In the large-scale component (Fig. 9b) that removes the sig-nature of fast convective signatures, we now observe together the slow convec-tive and faster wave signals. Singling out these latter signals by again remov-ing a 5-year running average in time (Fig. 9c) reveals QGA waves with outwardpropagation speeds close to, but lower than the one-dimensional Alfv´en velocity D / √ τ A =
225 km / yr. Comparing Figs. 9c and 6c also shows that QGA wavescan reach significantly smaller time scales and hence shorter radial wavelengthsthan torsional waves.While the Coriolis force identically vanishes from the magneto-inertial balance thatdrives axisymmetric torsional waves, it is necessarily present at the non-axisymmetriclevel and fundamentally modifies the wave equation (e.g. Finlay, 2008; Finlay et al.,2010). This leads in particular to magneto-Coriolis waves with phase velocitiesmuch smaller than the Alfv´en velocity. These are furthermore dispersive, with the26 ig. 9. Time-cylindrical radius diagrams of the non-axisymmetric, equator-symmetric coresurface azimuthal flow acceleration ˙ u S ϕ − ˙ u ZS ϕ at longitude 112 . o E (located by arrows inFig. 8c,d) in the 71p model, presented at (a) native spatial and temporal resolution, (b)large-scale (spherical harmonic degree up to 30), and (c) large-scale for the fast componentobtained after removal of a 5-year running average in time. The sets of slanted green linesdenote advection at the one-dimensional convective velocity D / √ τ U =
11 km / yr, andpropagation at the one-dimensional Alfv´en wave velocity D / √ τ A =
225 km / yr. wave velocity increasing as the wavelength decreases. The e ff ects of the Corio-lis force can however be mitigated for axially-invariant waves of su ffi ciently smallwavelength, that propagate in a direction perpendicular to the rotation vector. In thiscase, the wave velocity can approach the Alfv´en velocity from below (Gerick et al.,2020), as the influence of inertia is restored within the magneto-Coriolis balance.The axially columnar structure, propagation direction perpendicular to the rotationvector, velocities lower than D / √ τ A and short radial wavelengths of the QGAwaves observed here can also be understood within this framework. Our configu-ration however di ff ers from that of Gerick et al. (2020) in the sense that the back-ground magnetic field is highly complex, with relatively homogeneous regions sur-rounded by strongly heterogeneous field line concentrations corresponding to sharpgradients (e.g. Fig. 8b,d, see also Fig. 8 in Aubert, 2019). This situation promotesa spatial segregation of the force balance (Aubert, 2018). The magneto-Coriolispart of this balance is mostly observed at the edges of the QGA wavefronts whichcorrespond to the locations of strongest magnetic field heterogeneity. In the more27agnetically homogeneous regions, however, the Coriolis force is mitigated to thepoint where wave motion becomes locally magneto-inertial again. The residual in-fluence of the Coriolis force nevertheless tends to promote smaller radial lengthscales for QGA waves relative to torsional waves, which also rationalise the dif-ferent behaviours of energy density spectra ˆ K ( f ) and ˆ K ZS ( f ) (compare Figs. 5b,d,fand 7) as regards the steepness of their high-frequency decay and the influence ofbulk magnetic di ff usion (or lack thereof) on this decay. We also expect this residualinfluence to increase as the waves approach the core surface at equatorial position,because of the increasing slope of the spherical boundary. This leads to additionalwave slowdown compared to torsional waves and hence (the wave period beingpreserved) further reduction of the radial length scale (compare again Figs. 9c and6c.) The Rossby waves observed in Fig. 8c propagate eastwards in the vicinity of theequatorial plane. These are slow inertial waves in the sense that their pulsation ω ismuch smaller than the rotation rate Ω , thereby ensuring an evolution under a quasi-geostrophic force balance and an axially columnar structure (Zhang et al., 2001;Busse et al., 2005; Canet et al., 2014). From a geomagnetic standpoint, they arehowever fast in the sense that ω is larger than the typical pulsation 2 π/τ A of Alfv´enwaves. As a result, they are only slightly modified by the presence of the magneticfield (e.g. Finlay et al., 2010) and typically feature magnetic to kinetic energy ratiossmaller than unity (e.g. Gerick et al., 2020). They hence do not bear a significantsignature in the magnetic acceleration signals of Fig. 5, where a possible controlfrom τ Ω is elusive.In a full sphere (Zhang et al., 2001) (a configuration that reasonably applies to oursituation where the Rossby waves are confined at low to mid-latitudes, Fig. 8c), andwhen ω/ Ω <<
1, the normalised pulsation ω/ Ω of hydrodynamic Rossby wavescan be approximated by the expression ω Ω ≈ m + (cid:115) + m ( m + N (2 N + m + − , (18)which involves the azimuthal mode number m as well as a number N describingtheir level of complexity in the radial direction. Fig. 10a presents frequency-domainenergy density profiles of the equatorial acceleration seen in Fig. 8a, broken downinto contributions from distinct azimuthal mode numbers m . The peaks that can beobserved at various mode numbers are nearly synchronised in the vicinity of severaldistinct normalised pulsations. We focus on the three normalised pulsations ω/ Ω ≈ . , . , . ff erent values of m .The eastward propagation velocity of the corresponding waves can be determinedby performing a Radon transform on time-longitude diagrams of the fast equatorial28 a s t w a r d w a v e v e l o c i t y a t C M B ( k m . y r - ) m -1 ω Ω ≈ speeds identified on Radon spectraRossby modes theoretical speeds no r m a li s ed ene r g y den s i t y normalised pulsation ω/Ω m=1m=5m=6m=10 ω Ω ≈ ω Ω ≈
13 16 18 19 20 22 23 24N=6 8 10 11 12 13 1415 16 16N=7 7 8 8 89 99 m=2
N=11
Fig. 10. a: Frequency-domain energy density profiles of ˙ u ϕ at the equator in a 130-yearsequence of the 71p model. The acceleration signal is broken down into contributions fromdistinct azimuthal mode numbers m and the energy is normalised relatively to the peakat ω/ Ω ≈ . m at each of the three normalised pulsations ω/ Ω ≈ . , . , .
025 corresponding to the peaks in panel a (see arrows of samecolor). Also reported as functions of m are the closest matching phase velocities c = ω r o / m (squares) obtained by adjusting the radial complexity number N in equation (18). The greyline represents a m − law. acceleration signal, again broken down by mode number m . At each value of m , Fig.10b reports the propagation velocities identified at distinct Radon energy peaks.Next to observed velocities, Fig. 10b also reports theoretical phase velocities c = ω/ k , where k = m / r o and ω is determined from equation (18) by finding the radialcomplexity level N that best matches the observed velocity. The agreement betweenobserved and theoretical velocities is excellent throughout the range investigatedfor ω/ Ω and m , and provides an unambiguous characterisation of the Rossby wavesseen in Fig. 8a (and more clearly in Fig. 8c after filtering out the low-frequencycontent). The synchronisation of ω/ Ω across distinct mode numbers m is possiblebecause despite their discrete character, the possible values given by equation (18)form a ensemble that is dense enough to o ff er ( m , N ) couples yielding a pulsationclose to ω/ Ω ≈ . , . , . ω/ Ω ,the phase speed c = ω r o / m of each wave category in Fig. 10b naturally scales like m − . Unlike Alfv´en waves, Rossby waves are indeed highly dispersive in nature,with small-scale waves being significantly slower than large-scale waves. Whileall the waves observed here have periods within the interannual time scale range(see arrows in Fig. 10a), we note that at the largest scales, their velocities up to29 cc e l e r a t i on ene r g y ( k m . y r - ) -5 -4 -3 -2 -1 -2 -1 path parameter ε
71% of path 50% of path 29% of path start of path ε -0.25 [K][K ZS ]1 [K]-[K ] ε -0.13 Fig. 11. a: Evolution with the path parameter (cid:15) of the time-averaged, total core surfaceflow acceleration energy [ K ] up to spherical harmonic degree 30, plotted together with twodiagnostics of the evolution of wave energy along the path (Earth’s core conditions are to-wards the left of the graph). The energy of axisymmetric torsional waves is represented bythe axisymmetric, azimuthal and equator-symmetric part [ K ZS ], and the energy of non-ax-isymmetric waves is represented by the di ff erence [ K ] − [ K ] of [ K ] to its value at 29percent of the path. The grey lines mark the (cid:15) . and (cid:15) . power-law behaviours. km / yr are the fastest observed in the 71p model. In Fig. 11 we present the evolution of flow acceleration energy along the parameterspace path for the various waves that we have isolated. As measured by [ K ZS ], theenergy of torsional waves remains subdominant relatively to the total accelerationenergy [ K ], but gains two orders of magnitude along the path to reach [ K ZS ] / [ K ] ≈ . K ZS ] ∼ A − ∼ (cid:15) − / from equation (17) is well approached by our results. in contrast, [ K ] asymptoti-cally evolves less steeply (with [ K ] ∼ (cid:15) − . from our results) because it containsan invariant contribution from convection. This contribution is approximately re-moved by subtracting to [ K ] the acceleration energy [ K ] at the entry into theasymptotic regime, and [ K ] − [ K ] may be seen as the total wave contributionthat includes non-axisymmetric motion. The spectral form (15) suggests the scal-ing [ K ] − [ K ] ∼ S ∼ A − Rm ∼ (cid:15) − / , again in fair agreement (albeit being clearlyvalidated by two points only) with our numerical data at advanced path positions.At 71 percent of the path, the total wave contribution [ K ] − [ K ] amounts to abouthalf the total flow acceleration energy [ K ], and is five times stronger than the energy[ K ZS ] of axisymmetric torsional waves. If we denote as ˙ U W the typical accelerationof waves in the total acceleration ˙ U , in our advanced models this means that˙ U W / ˙ U ≈ (cid:32) [ K ] − [ K ][ K ] (cid:33) / = O (1) . (19)30e terminate our analysis by estimating the contributions of waves to the core flowand magnetic field (and no longer to the time derivatives of these quantities, as wehave done so far). Because these are relatively small, as we shall see, and becauseof the di ffi culties in disentangling the waves from convection that we have alreadyalluded to, this cannot be straightforwardly achieved from the numerical modeloutput and we need to take an alternative route. Assuming that the waves mainlyevolve with time scale τ A , while the convection evolves with time scale τ U (with A = τ A /τ U ), and denoting again as U W the contribution of waves to the typicalvelocity field amplitude U , we estimate that U W / U = A ˙ U W / ˙ U = O ( A ) . (20)We can further assume that the typical magnetic field amplitude B W carried by theAlfv´en waves obeys energy equipartition i.e. B W / U W = √ ρµ . Using the definition A = √ ρµ U / B , we then also estimate that B W / B = B W U W U W U UB = O ( A ) (21)This estimate is in principle valid at any altitude above the core-mantle boundary,if we assume that waves and convection share a similar length scale content. Equa-tions (20,21) show that the wave signatures in the total flow and magnetic fields arevery subdominant since A (cid:28)
1, and confirm that di ff erentiating the velocity fieldonce in time and the magnetic field twice in time is essential to highlighting thesecontributions. In numerical models sampling the parameter space path introduced in Aubert et al.(2017), dynamics of convective and wave origin occupy di ff erent, but overlappingfrequency ranges, as seen for instance in spectral energy density profiles of the mag-netic and flow acceleration (Fig. 5). Of particular interest is the secular to decadalrange located between the (approximately constant) overturn frequency f U = /τ U and the cut-o ff frequency f c , where energy density presents a plateau. This is whereslow convection and rapid waves are in potential interplay. We have seen that slowconvection accounts for the leftmost part of this plateau, the invariance of whichindeed mirrors the kinematic invariance observed along the path. It is possiblethat this plateau relates with vorticity equivalence in the spatial domain (David-son, 2013; Aubert, 2019), which states that vorticity tends to be evenly distributedby magnetic turbulence in the range [ d min , d ⊥ ] between the scale of magnetic dis-sipation and the dominant length scale at which convection powers the magnetic31eld. Assuming a simple linear correspondance between frequencies and lengthscales, this would indeed imply u ( f ) ∼ f − and ˙ u ( f ) ∼ f . This reasoning is ap-pealing for explaining the leftmost end of this plateau but falls short of interpretingthe rightmost end, the flattening of which along the path is a consequence of an ele-vating contribution coming primarily from hydromagnetic Alfv´en waves. This wasconfirmed in particular by the control on the cut-o ff frequency f c of the Lundquistnumber S (equation 13), which measures the number of Alfv´en wave periods oc-curring over the course of a magnetic di ff usion time. The case of torsional waves(Fig. 7) helped to illustrate the spectral broadening that can be expected as the timescale separation between the waves and convection increases, because of the higherwave frequencies that are made available and the non-linear wave-convection inter-actions underlain by Lorentz stresses. We hence conclude that the flat magneticacceleration energy density profile in the range [ f U , f c ] (corresponding to the f − range observed in energy density spectra of the geomagnetic field, De Santis et al.,2003; Lesur et al., 2018) can be ascribed to the combined e ff ect of non-linearlyinteracting slow convection and rapid Alfv´en waves.At higher frequencies f ≥ f c , most of the wave range is characterised by a f − decay power law. In this range, the wave acceleration energy is controlled by theLundquist number S = τ η /τ A (equations 14,15), again highlighting the importanceof time scale separation. This also underlines the stabilising control of magneticdi ff usion felt by waves in the bulk of the fluid (and not at the core surface, Fig. 5).Di ff usion is essentially active at non-axisymmetric scales, where waves typicallyreach small radial length scales mitigating the influence of the Coriolis acceleration.In terms of length scales also, the dependence f c ∼ √ S (equation 13) highlightsthe importance of the bulk magnetic di ff usion length scale at the Alfv´en time scale δ = √ ητ A = D / √ S . Waves at the largest scales such as axisymmetric torsionalwaves are immune to di ff usion (Fig. 7) but are nevertheless limited in amplitudebecause they are powered by non-linear couplings that involve non-axisymmetricscales where di ff usion-limited Alfv´en waves are found. As we advance along theparameter space path, and despite the tightening constraints set by the QG-MACbalance, the availability of higher wave frequencies and the non-linear energy trans-fers cause an increase in the total wave acceleration energy (Fig. 11). At 71 percentof the path we evaluate the total contribution from waves to the flow accelerationenergy to be equivalent to that of convection when considering the signal up tospherical harmonic degree 30. At a fifth of the wave acceleration energy, the con-tribution from torsional waves is subdominant.Using values τ η = τ U =
130 yr and τ A = S ≈ and f c ≈ .
07 yr − from the scaling (13).The f range of interplay between waves and convection should therefore ex-tend over a decade only between f U = − yr − and f c . The Alfv´en frequency1 /τ A = . − is predicted to lie well into the f − range, at an acceleration energydensity level of about ( f c τ A ) ≈ /
50 below the plateau (from equations 14,15).Assuming Rm ≈ τ η ) rather than the value Rm ≈ ff frequency to f c ≈ . − and the rel-ative attenuation at the Alfv´en frequency to 1 /
25. These last results fall short ofaccounting for the observed extent of the f acceleration ( f − magnetic) range upto frequency f ≈ − (Lesur et al., 2018), such that more fundamental modifi-cations of the modelling set-up appear to be needed. The path models are possiblyunderpowered at high frequencies because their design involves a neutrally buoy-ant core-mantle boundary. Exploring the e ff ects of convective instability near thecore surface on the high-frequency content of magnetic acceleration may thereforeprovide important geophysical constraints. It is unlikely that any degree of upperouter core stratification would help to bring the numerical results close to the ge-omagnetic observation, as this would further hinder the short-timescale dynamics(Aubert & Finlay, 2019). The stratification would also degrade the morphologicalresemblance of the model core surface magnetic field to the present-day field (Gas-tine et al., 2020). From the observational standpoint, this discussion also highlightsthe importance of e ff orts aiming at cleaning the geomagnetic acceleration signalfrom external contributions close to annual frequencies (see e.g. Finlay et al., 2017). Equation (21) suggests that the relative contribution B W from waves in the mag-netic field amplitude B is B W / B = O ( A ) i.e. B W ≈ A = . − of the core. Though this may appear extremelyweak, it remains within the typical resolution of satellite-based observations (seee.g. Finlay et al., 2016), and we have seen (Figs. 5,11, equation 19) that di ff eren-tiating the magnetic field twice in time, or the velocity field once in time, leadsto sizeable wave contributions to the total acceleration energy. This underlines thecrucial importance of extracting the geomagnetic acceleration from ground obser-vatory and satellite data at a good level of spatial and temporal resolution. Thisalso motivates research towards elaborate approaches aiming at extracting a reli-able flow acceleration from the magnetic acceleration signal (see a recent review inGillet, 2019). Furthermore, the increase of wave energy along the path (Fig. 11), theubiquitous character of waves (Figs. 8,9) as well as the large-scale content of tor-sional waves (Fig. 6) appear to promote some optimism as regards the detectabilityof waves in the geomagnetic signal emanating from Earth’s core.Among the signal caused by waves in the magnetic acceleration, hydromagneticAlfv´en waves are dominant because (unlike high-frequency Rossby waves) theyachieve equipartition between the kinetic and magnetic energies that they carry(Fig. 5). These should hence be most straightforwardly detectable in the geomag-netic acceleration signal. Concerning the typical flow amplitude U ZS of torsionalwaves, following the estimation strategy leading to equation (20) we infer U ZS / U ≈ A (cid:112) [ K ZS ] / [ K ], with an extrapolation to the end of the path yielding [ K ZS ] / [ K ] ≈ .
17. Using U = D /τ U ≈
17 km / yr and A = . − , this finally leads to33 ZS ≈ . / yr, somewhat smaller than (but of the same order of magnitude as)the value ≈ . / yr obtained from modelling of the rapidly evolving core flow(see Fig. 13 in Gillet et al., 2015). Likewise, at the end of the path the total waveamplitude U W including axisymmetric and non-axisymmetric contributions shouldbe higher. From equations (19,20) we infer U W / U ≈ A √ ([ K ] − [ K ]) / [ K ] ≈ A and therefore U W ≈ .
25 km / yr, again comparable to the value ≈ . / yr ob-tained by Gillet et al. (2015) for the non-zonal flow. We finally note that the value[ K ZS ] / ([ K ] − [ K ZS ]) ≈ . ff ects (Aubert & Finlay, 2019). We have found stronger andmore frequent jerks in the 71p model than at earlier path positions, though somejerks events may become too fast for being noticeable in determinations of the ge-omagnetic acceleration with limited temporal resolution (Figs. 3,4). Previously, wehad ascribed the increase of jerk energy E J along the path to an increased levelof wave radial shoaling, but the present results rather incite us to simply associateit with the linear increase of wave acceleration energy with S that we have docu-mented in this study.The case for detectability of hydrodynamic Rossby waves in the geomagnetic sig-nal appears considerably less obvious, as we did not find a conclusive influence ofthe rotational time scale τ Ω in the magnetic acceleration energy spectra of our mod-els (Fig. 5). Only the slowest Rossby modes with periods approaching the Alfv´entime scale τ A can in principle be detected from geomagnetic observations, becausefaster modes carry less magnetic than kinetic energy (e.g. Gerick et al., 2020).Although these modes also feature the slowest eastward propagation speeds (Fig.10), these speeds are still up to several thousands of kilometers per year in the71p model and should even be ten times faster in the core as they scale linearlywith 1 /τ Ω (equation 18). Seen from the standpoint of typical geomagnetic accel-eration timescales, these would therefore amount to almost instantaneous signalsand would hence be very hard to isolate as propagating features, unless the signalis considered at the smallest spatial scales where the velocity of these dispersivewaves is considerably slower. Eastward-propagating, equatorial geomagnetic sig-nals at speeds in the range 500 − / yr and wavenumbers m = − N (equation 18) was allowed to reach high values N >
60, which is beyond thelatitudinal resolution available for geomagnetic acceleration.34 .3 The path theory in the light of the 71p model and future prospects towardsreaching Earth’s core conditions
Following the path theory introduced in Aubert et al. (2017), this study has intro-duced a model located at 71 percent of this path. This has considerably enlargedthe asymptotic portion (beyond 29 percent) where similar models are available, andthe 71p model could again provide a complete validation to this theory. The leadingorder QG and first-order MAC force balances are stable, while the amplitude of in-ertia and viscosity continually decrease along the path (compare Fig. 1 to Fig. 1 ofAubert, 2018). Power-driven, di ff usivity-free scaling laws proposed in Aubert et al.(2017) could be confirmed (Fig. 2). In the 71p model this leads to a state where anoverwhelming fraction of the injected convective power is Ohmically dissipated,and where the Taylor constraint is enforced at a high level (Table 3). The dynamois also in a strong-field state, with the ratio of magnetic to kinetic energy beingdirectly given by 1 / A = τ U are still invariant and carried over from the start ofpath, while the short-timescale dynamics is gradually enriched in a way that hasbeen documented here in detail (Fig. 3,5). Here we have demonstrated that most,and therefore probably all of the parameter space path is devoid of abrupt physicaltransitions. This further rationalises the relevance of earlier dynamo models lo-cated close to the start of the path to describe the geodynamo, and this also furtherstrengthens the likeliness of core dynamics being in a similar dynamical regime asthat observed at 71 percent of the path.With a computation carried out over the course of two years, involving 246 milliontime steps and 15 million single-core CPU hours, the model at 71 percent of the pathmay be seen as an extreme endeavour that apparently obscures the prospect of beingable to terminate the exploration of this path. A positive note is that the computationof this model provided the opportunity to perform several low-level optimisationsin the numerical code towards a faster execution at given number of cores andalso towards the possibility to use more cores while maintaining a good strongparallel scaling. Together with an improved generation of supercomputers providedby GENCI in France, we could achieve at least a fivefold increase in computationspeed for the path models. This was used to increase the spatial resolution anddecrease the hyperdi ff usivity at which the 71p model runs, and also opened the wayto an exploration of cross-path models at 50 percent of the path in a way that couldnot have been feasible only two years ago, at the time of our earlier study (Aubert,2018). We have also seen that it is no longer necessary to advance in half-decadesof (cid:15) along the path, and that much larger leaps can be achieved through a simplere-scaling procedure of checkpoint files. The next stop along the road is thereforemost probably the final one, i.e. being able to simulate the geodynamo exactly inEarth’s core conditions. This will presumably imply another round of optimisations35nd another generation of Tier-1 supercomputers, but we can in principle foreseethe completion of this challenge in the coming decade. Acknowledgements
The authors thank two anonymous referees for comments and Dominique Jault,Nathana¨el Schae ff er and Thomas Gastine for discussions and help in code optimi-sation. JA acknowledges support from the Fondation Simone et Cino Del Duca ofInstitut de France (2017 Research Grant). This project has also been funded by ESAin the framework of EO Science for Society, through contract 4000127193 / / NL / IA(SWARM +
4D Deep Earth: Core). NG was partially supported by the French Cen-tre National d’Etudes Spatiales (CNES) for the study of Earth’s core dynamics inthe context of the Swarm mission of ESA. Numerical computations were performedat S-CAPAD, IPGP and using HPC resources from GENCI-TGCC and GENCI-CINES (Grant numbers A0060402122 and A0080402122).
Data availability
The numerical code and simulation data analysed in this study are available fromthe corresponding author upon reasonable request. The core surface data from the71 percent of path model is also available online at the URLhttps: // / data. References
Alfv´en, H. & F¨althammar, C., 1963.
Cosmical Electrodynamics , chap. 3.4.5, Ox-ford Univ. Press, 2nd edn.Aubert, J., 2018. Geomagnetic acceleration and rapid hydromagnetic wave dy-namics in advanced numerical simulations of the geodynamo,
Geophys. J. Int. , (1), 531–547.Aubert, J., 2019. Approaching Earth’s core conditions in high-resolution geody-namo simulations, Geophys. J. Int. , (S1), S137–S151.Aubert, J., 2020. Recent geomagnetic variations and the force balance in Earth’score, Geophys. J. Int. , (1), 378–393.Aubert, J. & Finlay, C. C., 2019. Geomagnetic jerks and rapid hydromagneticwaves focusing at Earth’s core surface, Nature Geosci. , (5), 393–398.Aubert, J., Finlay, C. C., & Fournier, A., 2013. Bottom-up control of geomagneticsecular variation by the Earth’s inner core, Nature , , 219–223.36ubert, J., Gastine, T., & Fournier, A., 2017. Spherical convective dynamos in therapidly rotating asymptotic regime, J. Fluid. Mech. , , 558–593.Aurnou, J. M. & King, E. M., 2017. The cross-over to magnetostrophic convectionin planetary dynamo systems, Proc. Roy. Soc. A , (2199), 20160731.Bouligand, C., Gillet, N., Jault, D., Schae ff er, N., Fournier, A., & Aubert, J., 2016.Frequency spectrum of the geomagnetic field harmonic coe ffi cients from dynamosimulations, Geophys. J. Int. , (2), 1142–1157.Bu ff ett, B. & Matsui, H., 2015. A power spectrum for the geomagnetic dipolemoment, Earth Plan. Sci. Lett. , , 20 – 26.Busse, F. H., Zhang, K., & Liao, X., 2005. On slow inertial waves in the solarconvection zone, Ap. J. , (2), L171–L174.Calkins, M. A., 2018. Quasi-geostrophic dynamo theory, Phys. Earth. Planet. Int. , , 182 – 189.Calkins, M. A., Julien, K., Tobias, S. M., & Aurnou, J. M., 2015. A multiscaledynamo model driven by quasi-geostrophic convection, J. Fluid. Mech. , ,143–166.Canet, E., Finlay, C., & Fournier, A., 2014. Hydromagnetic quasi-geostrophicmodes in rapidly rotating planetary cores, Phys. Earth Planet. Int. , , 1 – 15.Chi-Dur´an, R., Avery, M. S., Knezek, N., & Bu ff ett, B. A., 2020. Decomposi-tion of geomagnetic secular acceleration into traveling waves using complexempirical orthogonal functions, Geophys. Res. Lett. , (17), e2020GL087940,e2020GL087940 10.1029 / Geophys. J. Int. , , 97–114.Christensen, U. R., Aubert, J., & Hulot, G., 2010. Conditions for Earth-like geody-namo models, Earth. Plan. Sci. Let. , (3-4), 487–496.Christensen, U. R., Wardinski, I., & Lesur, V., 2012. Timescales of geomagneticsecular acceleration in satellite field models and geodynamo models, Geophys.J. Int. , (1), 243–254.Chulliat, A., Alken, P., & Maus, S., 2015. Fast equatorial waves propagating at thetop of the Earth’s core, Geophys. Res. Lett. , (9), 3321–3329.Constable, C. & Johnson, C., 2005. A paleomagnetic power spectrum, Phys. EarthPlanet. Int. , (1), 61 – 73, Studies of the Earth’s Deep Interior.Davidson, P. A., 2013. Scaling laws for planetary dynamos, Geophys. J. Int. , (1), 67–74.De Santis, A., Barraclough, D., & Tozzi, R., 2003. Spatial and temporal spectraof the geomagnetic field and their scaling properties, Phys. Earth Planet. Int. , (2), 125 – 134, Magnetic Field Modelling.Finlay, C. C., 2005. Hydromagnetic waves in Earth’s core and their influence ongeomagnetic secular variation , Ph.D. thesis, School of Earth and Environment,University of Leeds.Finlay, C. C., 2008. Waves in the presence of magnetic fields, rotation and con-vection, in
Lecture notes on Les Houches Summer School: Dynamos , vol. 88,chap. 8, pp. 403–450, eds Cardin, P. & Cugliandolo, L. F., Elsevier.37inlay, C. C., Olsen, N., Kotsiaros, S., Gillet, N., & Tø ff ner-Clausen, L., 2016.Recent geomagnetic secular variation from Swarm and ground observatories asestimated in the CHAOS-6 geomagnetic field model, Earth, Planets and Space , (1), 112.Finlay, C. C., Lesur, V., Th´ebault, E., Vervelidou, F., Morschhauser, A., & Shore,R., 2017. Challenges handling magnetospheric and ionospheric signals in inter-nal geomagnetic field modelling, Space. Sci. Rev. , (1), 157–189.Finlay, C. C., Dumberry, M., Chulliat, A., & Pais, M. A., 2010. Short TimescaleCore Dynamics: Theory and Observations, Space. Sci. Rev. , (1-4), 177–218.Gastine, T., Aubert, J., & Fournier, A., 2020. Dynamo-based limit to the extent ofa stable layer atop Earth’s core, Geophys. J. Int. , (2), 1433–1448.Gerick, F., Jault, D., & Noir, J., 2020. Fast Quasi-Geostrophic Magneto-CoriolisModes in the Earth’s core, Geophys. Res. Lett. , , 2020GL090803, doi:10.1029 / Geomagnetism, Aeronomy andSpace Weather : a Journey from the Earth’s Core to the Sun , chap. 9, eds Man-dea, M., Korte, M., Petrovsky, E., & Yau, A., International Association of Geo-magnetism and Aeronomy.Gillet, N., Jault, D., & Finlay, C. C., 2015. Planetary gyre, time-dependent eddies,torsional waves and equatorial jets at the Earth’s core surface,
J. Geophys. Res. , , 3991–4013.Gillet, N., Jault, D., & Canet, E., 2017. Excitation of travelling torsional normalmodes in an earth’s core model, Geophys. J. Int. , (3), 1503–1516.Gillet, N., Huder, L., & Aubert, J., 2019. A reduced stochastic model of core surfacedynamics based on geodynamo simulations, Geophys. J. Int. , (1), 522–539.Gillet, N., Jault, D., Canet, E., & Fournier, A., 2010. Fast torsional waves andstrong magnetic field within the Earth’s core, Nature , (7294), 74–77.Gillet, N., Schae ff er, N., & Jault, D., 2011. Rationale and geophysical evidencefor quasi-geostrophic rapid dynamics within the Earth’s outer core, Phys. EarthPlanet. Int. , (3-4, SI), 380–390.Hide, R., 1966. Free hydromagnetic oscillations of the earth’s core and the theoryof the geomagnetic secular variation, Phil. Trans. Roy. Soc. A , (1107), 615–647.Jault, D., 2008. Axial invariance of rapidly varying di ff usionless motions in theEarth’s core interior, Phys. Earth Planet. Int. , (1-2), 67–76.Jones, C., 2015. 8.05 - thermal and compositional convection in the outer core,in Treatise on Geophysics (Second Edition) , pp. 115 – 159, ed. Schubert, G.,Elsevier, Oxford, second edition edn.Labb´e, F., Jault, D., & Gillet, N., 2015. On magnetostrophic inertia-less waves inquasi-geostrophic models of planetary cores,
Geophys. Astrophys. Fluid Dyn. , (6), 587–610.Lehnert, B., 1954. Magnetohydrodynamic Waves Under the Action of the CoriolisForce., Ap. J. , , 647.Lesur, V., Wardinski, I., Asari, S., Minchev, B., & Mandea, M., 2010. Modelling38he earth’s core magnetic field under flow constraints, Earth, Plan. Space , (6),503–516.Lesur, V., Wardinski, I., Baerenzung, J., & Holschneider, M., 2018. On the fre-quency spectra of the core magnetic field gauss coe ffi cients, Phys. Earth Planet.Int. , , 145–158.Lhuillier, F., Fournier, A., Hulot, G., & Aubert, J., 2011. The geomagnetic secular-variation timescale in observations and numerical dynamo models, Geophys.Res. Lett. , , L09306.Meduri, D. G. & Wicht, J., 2016. A simple stochastic model for dipole momentfluctuations in numerical dynamo simulations, Front. Earth. Sci. , , 38.Olson, P., Christensen, U., & Driscoll, P., 2012. From superchrons to secular vari-ation: A broadband dynamo frequency spectrum for the geomagnetic dipole, Earth. Plan. Sci. Let. , , 75–82.Panovska, S., Finlay, C., & Hirt, A., 2013. Observed periodicities and the spectrumof field variations in holocene magnetic records, Earth Plan. Sci. Lett. , , 88 –94.Pozzo, M., Davies, C. J., Gubbins, D., & Alf`e, D., 2012. Thermal and electricalconductivity of iron at Earth’s core conditions, Nature , (7398), 355–358.Rieutord, M. & Valdettaro, L., 1997. Inertial waves in a rotating spherical shell, Journal of Fluid Mechanics , , 77–99.Schae ff er, N., 2013. E ffi cient spherical harmonic transforms aimed at pseudospec-tral numerical simulations, Geophys. Geochem. Geosystems. , (3), 751–758.Schae ff er, N. & Jault, D., 2016. Electrical conductivity of the lowermost mantleexplains absorption of core torsional waves at the equator, Geophys. Res. Lett. , (10), 4922–4928.Schae ff er, N., Jault, D., Cardin, P., & Drouard, M., 2012. On the reflection of alfv´enwaves and its implication for earth’s core modelling, Geophys. J. Int. , (2),508–516.Schae ff er, N., Jault, D., Nataf, H.-C., & Fournier, A., 2017. Turbulent geodynamosimulations: a leap towards Earth’s core, Geophys. J. Int. , (1), 1–29.Schwaiger, T., Gastine, T., & Aubert, J., 2019. Force balance in numerical geody-namo simulations: a systematic study, Geophys. J. Int. , (S1), S101–S114.Schwaiger, T., Gastine, T., & Aubert, J., 2021. Relating force balances and flowlength scales in geodynamo simulations, Geophys. J. Int. , (3), 1890–1904.Taylor, J., 1963. Magneto-hydrodynamics of a rotating fluid and Earths dynamoproblem, Proc. Roy. Soc. A , , 274–283.Teed, R. J., Jones, C. A., & Tobias, S. M., 2014. The dynamics and excitation oftorsional waves in geodynamo simulations, Geophys. J. Int. , (2), 724–735.Teed, R. J., Jones, C. A., & Tobias, S. M., 2019. Torsional waves driven by con-vection and jets in Earth’s liquid core, Geophys. J. Int. , (1), 123–129.Wicht, J. & Christensen, U. R., 2010. Torsional oscillations in dynamo simulations, Geophys. J. Int. , (3), 1367–1380.Wicht, J. & Sanchez, S., 2019. Advances in geodynamo modelling, Geophys. As-trophys. Fluid Dyn. , (1-2), 2–50.Zhang, K., Earnshaw, P., Liao, X., & Busse, F. H., F. H., 2001. On inertial waves39n a rotating fluid sphere, Journal of Fluid Mechanics ,437