# Convectively driven decadal zonal accelerations in Earth’s fluid core

Convectively driven decadal zonal accelerations in Earth’s fluid core Summary Azimuthal accelerations of cylindrical surfaces co-axial with the rotation axis have been inferred to exist in Earth’s fluid core on the basis of magnetic field observations and changes in the length-of-day. These accelerations have a typical timescale of decades. However, the physical mechanism causing the accelerations is not well understood. Scaling arguments suggest that the leading order torque averaged over cylindrical surfaces should arise from the Lorentz force. Decadal fluctuations in the magnetic field inside the core, driven by convective flows, could then force decadal changes in the Lorentz torque and generate zonal accelerations. We test this hypothesis by constructing a quasi-geostrophic model of magnetoconvection, with thermally driven flows perturbing a steady, imposed background magnetic field. We show that when the Alfvén number in our model is similar to that in Earth’s fluid core, temporal fluctuations in the torque balance are dominated by the Lorentz torque, with the latter generating mean zonal accelerations. Our model reproduces both fast, free Alfvén waves and slow, forced accelerations, with ratios of relative strength and relative timescale similar to those inferred for the Earth’s core. The temporal changes in the magnetic field which drive the time-varying Lorentz torque are produced by the underlying convective flows, shearing and advecting the magnetic field on a timescale associated with convective eddies. Our results support the hypothesis that temporal changes in the magnetic field deep inside Earth’s fluid core drive the observed decadal zonal accelerations of cylindrical surfaces through the Lorentz torque. Core, Numerical modelling, Planetary interiors 1 INTRODUCTION Decadal variations in Earth’s length-of-day (LOD) are believed to be driven by core flow variations (Gross 2015). For instance, an increase in the bulk angular velocity of the core entrains an increase in its axial angular momentum. Angular momentum conservation between the core and the mantle then implies that the latter must slow its rotation rate, thereby leading to a increase in LOD. On timescales of centuries or less, the dynamics of the outer core are expected to be dominated by a geostrophic balance between pressure gradients and the Coriolis force (e.g. Finlay et al.2010). Flows that obey such a balance have the property of being invariant, or ‘rigid’, parallel to the axis of rotation (Hough 1897; Proudman 1916; Taylor 1917). This is the Taylor–Proudman theorem. Temporal variations in the outer core’s axial angular momentum must then be carried by azimuthal accelerations in the form of rigid, co-axial cylindrical surfaces, or ‘geostrophic cylinders’, as shown in Fig.  1. Figure 1. View largeDownload slide The geometry of the zonal flows in Earth’s outer core which carry angular momentum. The depicted time varying zonal flows consist in either free Alfvén waves or forced accelerations. Figure 1. View largeDownload slide The geometry of the zonal flows in Earth’s outer core which carry angular momentum. The depicted time varying zonal flows consist in either free Alfvén waves or forced accelerations. Flows at the core–mantle boundary (CMB) may be reconstructed from the magnetic field’s secular variation observed at Earth’s surface (Holme 2015). Time-dependent accelerations of geostrophic cylinders may then be extracted from the reconstructed core flows, allowing a prediction of LOD changes to be built. A number of studies (e.g. Jault et al.1988; Jackson et al.1993) have shown that such predictions match well with the observed decadal LOD variations. Furthermore, Gillet et al. (2010) have shown that zonal accelerations of geostrophic cylinders carry changes in angular momentum which can also explain an observed 6-yr periodic LOD signal (Holme & de Viron 2013; Chao et al.2014). These studies confirm not only that the decadal and 6-yr LOD variations are due to core-mantle angular momentum exchanges, they also confirm the presence of rigid zonal accelerations in the core. However, the dynamics responsible for such zonal accelerations are not fully understood. This is the topic of our study. The forces responsible for controlling the dynamics of zonal flows can be examined by integrating the azimuthal component of the momentum equation over the surface of geostrophic cylinders. Upon doing so, the pressure, Coriolis, and buoyancy terms vanish. If one neglects inertia and viscosity, this implies that the azimuthal component of the Lorentz force must cancel out when integrated over geostrophic cylinders. In other words, the axial Lorentz torque on cylinders must vanish. This is known as Taylor’s condition (Taylor 1963), with systems obeying it said to be in a Taylor state. Reinstating inertia, perturbations in the magnetic field – and, therefore, in the Lorentz torque – can be accommodated by rigid zonal accelerations of geostrophic cylinders. When a magnetic field in a conducting fluid is distorted by the fluid’s motion, such as the differential rotation between coaxial cylinders, a current is induced. This current interacts with the original magnetic field to create a restoring force which opposes the motion that caused the initial magnetic distortion. In the language of Taylor’s condition, this restoring force nudges geostrophic cylinders back towards a Taylor state, subject to magnetic diffusion. This mechanism allows rigid zonal flows to oscillate about a Taylor state and, in doing so, support Alfvén waves. Braginsky (1970) suggested that the free modes of Alfvén waves could be responsible for decadal zonal accelerations in the outer core. In order to produce free Alfvén waves with a fundamental mode period of approximately 60 yr, corresponding to the characteristic timescale of the LOD signal (Roberts et al.2007; Gross 2015), the magnetic field strength throughout the core must be approximately 0.3 mT, similar to its observed strength at the CMB. However, modern geodynamo simulations suggest that the internal radial magnetic field is approximately ten times larger than that at the CMB (e.g. Christensen & Aubert 2006). Such a field strength should yield a fundamental Alfvén period of approximately 6 yr. It is therefore now believed that the 6-yr LOD signal is due to the propagation of free Alfvén waves (Gillet et al.2010), leaving the dynamics responsible for the decadal zonal accelerations unexplained. One possible explanation is that convective eddies in the fluid core continuously distort the magnetic field, causing spatial and temporal variations in the Lorentz force. Integrated over geostrophic cylinders, the induced Lorentz torques must be balanced by rigid zonal accelerations of the cylinders. Continual distortion of the internal magnetic field would then lead to continual zonal accelerations of the cylinders. The observed decadal zonal accelerations could therefore be a forced fluctuation about the Taylor state driven by convective flows. If this is the case, forcing must occur on timescales longer than the propagation time of free Alfvén waves. In other words, the typical convective velocity uC must be smaller than the Alfvén wave velocity uA. The relationship between uC and uA is captured by the Alfvén number   $$\mathcal {A}= \frac{u_{C}}{u_{A}},$$ (1)where   $$u_{A} = \frac{\left| \bf {B} \right|}{\sqrt{\rho \mu _{0}}},$$ (2) B is the magnetic field, ρ is the fluid density, and μ0 is the magnetic permeability of free space. Typical large-scale flow velocities in Earth’s core are of the order of uC = 10 km yr−1, or 3 × 10−4 m s−1 (Holme 2015). With a typical radial magnetic field of 3 mT (e.g. Gillet et al.2010; Buffett 2010), an outer core density of 104 kg m−3, and μ0 = 4π × 10−7 N A−2, the Alfvén wave velocity is approximately uA = 3 × 10−2 m s−1. This yields $$\mathcal {A}\approx 0.01$$. Such a small value of $$\mathcal {A}$$ supports the idea that the observed decadal zonal flows in Earth’s core may be driven by convection via time-dependent Lorentz torques. The goal of this study is to demonstrate that such a dynamical scenario is possible. One option to investigate the zonal flow dynamics in the Earth’s core is to use a 3-D numerical model of a self-generated dynamo (e.g. Christensen & Wicht 2015). However, typical Alfvén numbers in such models are $$\mathcal {A}\approx 1$$, meaning the dynamics of free Alfvén waves and convectively driven zonal accelerations are not well separated. Some recent 3-D models have been able to achieve lower values of $$\mathcal {A}\approx 0.1$$ (e.g. Aubert et al.2017; Schaeffer et al.2017), though these models are numerically expensive to run. Here, we follow a different strategy and exploit the fact that, as discussed previously, the large-scale fluid motions in the core which vary on decadal timescales are expected to be almost invariant along the rotation axis (Jault 2008). Flows of this type are often termed Quasi-Geostrophic (QG), with their existence in Earth’s core supported by the observed geomagnetic secular variation (Pais & Jault 2008; Gillet et al.2011). They also emerge from numerical models of the geodynamo (e.g. Schaeffer et al.2017). We present in this study a numerical model of the decadal timescale dynamics of Earth’s fluid core constructed within a QG framework (e.g. Aubert et al.2003). We add an induction equation to a QG model of thermal convection to follow the evolution of the magnetic field as it is sheared and advected by the flow. We likewise take into account the feedback of the Lorentz force on the flow, though in a limited way—see Section 2.7. Since we focus on the short timescale dynamics, we substitute the self-generated magnetic field resulting from dynamo action with a steady, imposed background magnetic field. We thus perform a magnetoconvection experiment, in which the perturbations of the magnetic field are tracked with respect to this imposed field. This strategy allows us to readily achieve $$\mathcal {A}< 1$$ at very modest numerical cost. 2 MODEL 2.1 Background Previous studies (e.g. Aubert et al.2003; Gillet & Jones 2006) have shown that thermally driven, Boussinesq, QG models can reproduce flow patterns similar in scale and behaviour to the ones we expect in planetary cores. The QG approximation takes advantage of the geometric constraint imposed by strong rotation: fluid motions must be dominantly rigid. Studying the dynamics of the full system is then equivalent to studying the dynamics of a slice through the system’s equatorial plane. A QG model therefore collapses 3-D convection experiments onto a 2-D domain. In our case, this domain corresponds to the shaded annulus of Fig. 2. Using the usual cylindrical coordinates (s, ϕ, z), with the ez direction aligned with the rotation axis, this annulus extends from the ‘tangent cylinder’ surrounding the equivalent of the inner core at s = s1, to the ‘equator’ at s = s2. We do not model the region inside the tangent cylinder, equivalent to the polar regions above and below the inner core. Figure 2. View largeDownload slide Geometry of our QG model. The domain of integration is the shaded 2-D annulus between s1 and s2. Figure 2. View largeDownload slide Geometry of our QG model. The domain of integration is the shaded 2-D annulus between s1 and s2. Our model is based on the QG model of thermal convection, first presented by Busse & Or (1986) and Cardin & Olson (1994) and expanded in many subsequent studies (e.g. Aubert et al.2003; Gillet & Jones 2006). In this model, 2-D equations for the momentum and thermal evolution are coupled together. To these governing equations, we add a third, two-dimension magnetic induction equation, while augmenting the momentum equation with a Lorentz force term. Although the Taylor–Proudman theorem predicts rigid motions, there are no such theoretical constraints on the temperature nor the magnetic field. However, by averaging the governing equations in the axial (ez) direction, we capture, if imperfectly, the net effect of both the temperature and magnetic field on the flow dynamics. Strategies have been devised to couple QG flows to a 3-D magnetic field (Schaeffer & Cardin 2006; Schaeffer et al.2016) and temperature (Guervilly & Cardin 2016), but here we restrict our attention to a purely 2-D model. We scale length by the radius of the outer sphere r2, time by the inverse of the angular rotational velocity Ω, temperature by the superadiabatic temperature difference ΔT between the inner and outer spheres, and the magnetic field by $$r_{2} \Omega \sqrt{\rho _{0} \mu _{0}}$$, where ρ0 is the reference density. The QG approximation remains valid so long as the Coriolis force remains dominant. In our case, this implies that both the Rossby number Ro (the ratio of inertial to Coriolis forces) and the Lehnert number λ (the ratio of magnetic to Coriolis forces) must remain ≪1. Our choice of scalings causes Ro in our model to be equivalent to |u| (the typical amplitude of the non-dimensional velocity) and λ to be equivalent to |b| (the typical amplitude of the non-dimensional magnetic field perturbation). We must therefore ensure that |u| ≪ 1 and |b| ≪ 1. 2.2 Flow and magnetic fields Because both the flow and magnetic fields are solenoidal, they may be represented in terms of vector potentials. We define the velocity u and magnetic field perturbation b as   \begin{eqnarray} \mathbf {u} = \overline{u_{\phi }}\, \mathbf {e}_{\phi }+ \frac{1}{L} \boldsymbol{\nabla }\times \left(L \psi \, \mathbf {e}_{z}\right)\, , \end{eqnarray} (3a)  \begin{eqnarray} \mathbf {b} = \overline{b_{\phi }}\, \mathbf {e}_{\phi }+ \frac{1}{L} \boldsymbol{\nabla }\times \left(L a \, \mathbf {e}_{z}\right)\, , \end{eqnarray} (3b)where ψ and a are toroidal scalars and $$L = \sqrt{1 - s^{2}}$$ is the half-column height. With an overbar denoting an azimuthal average, $$\overline{u_{\phi }}$$ and $$\overline{b_{\phi }}$$ capture the axisymmetric azimuthal (zonal) flow and zonal magnetic fields, respectively. The horizontal components of the velocity field $${\bf u}_{H} = u_{s}\mathbf {e}_{s}+ u_{\phi }\mathbf {e}_{\phi }$$ and the axial vorticity ωz are then defined as   \begin{eqnarray} u_{s} = \frac{1}{s} \frac{\partial \psi }{\partial \phi } \, , \end{eqnarray} (4a)  \begin{eqnarray} u_{\phi } = \overline{u_{\phi }}- \left(\frac{\partial }{\partial s}+ \beta \right)\psi \, , \end{eqnarray} (4b)  \begin{eqnarray} \omega _{z} = \left(2 \frac{\overline{u_{\phi }}}{s}+ \frac{\partial }{\partial s}\overline{u_{\phi }}\right)- \nabla _{H}^{2} \psi - \frac{1}{s} \frac{\partial }{\partial s}\left(s \beta \psi \right)\, . \end{eqnarray} (4c) Similarly, the horizontal components of the magnetic perturbation field $${\bf b}_{H} = b_{s}\mathbf {e}_{s}+ b_{\phi }\mathbf {e}_{\phi }$$ and the axial current jz are defined as   \begin{eqnarray} b_{s}= \frac{1}{s} \frac{\partial a}{\partial \phi } \, , \end{eqnarray} (5a)  \begin{eqnarray} b_{\phi } = \overline{b_{\phi }}- \left(\frac{\partial }{\partial s}+ \beta \right)a \, , \end{eqnarray} (5b)  \begin{eqnarray} j_{z}= \left(2 \frac{\overline{b_{\phi }}}{s} + \frac{\partial }{\partial s}\overline{b_{\phi }}\right)- \nabla _{H}^{2} a - \frac{1}{s} \frac{\partial }{\partial s}\left(s \beta a \right)\, . \end{eqnarray} (5c) Here, $$\nabla _{H}^{2}$$ specifies the cylindrical (s, ϕ) components of the Laplacian operator. The factor β, which enters the equations via the no-penetration condition on the r = r2 surface, is a measure of how L changes with s:   $$\beta = \frac{1}{L} \frac{\partial L}{\partial s} = - \frac{s}{L^{2}} \, .$$ (6) The definitions of us, uϕ and ωz according to eqs (4) follow the traditional approach in QG models (e.g. Schaeffer & Cardin 2005). Because the fluid is assumed incompressible, there is no axisymmetric radial velocity. Using an equivalent representation for bs, bϕ and jz in eqs (5) follows the strategy employed by Labbé et al. (2015). These rigid variables represent the axial averages of the truly 3-D magnetic field. Since their form is equivalent to the flow field, the implied assumption is that the magnetic field obeys the equivalent of a no-penetration condition on the spherical top and bottom boundary of a fluid column. While this cannot be rigorously justified, to first order this is approximately correct since the magnetic field near the core-mantle boundary is likely much smaller than deeper in the core. Furthermore, as we show in Appendix, representing the magnetic field as in eqs (5) is necessary to ensure conservation of angular momentum. 2.3 Momentum equation Under the QG approximation, the axially averaged, thermally driven Navier–Stokes equation reduces to an axial vorticity equation (e.g. Aubert et al.2003),   \begin{eqnarray} \frac{\partial }{\partial t}\omega _{z}&=& - Ra^{*}\frac{\partial \Theta }{\partial \phi } + E \nabla _{H}^{2} \omega _{z}+ \left(2 + \omega _{z}\right)\beta u_{s} \nonumber\\ &&- \left(u_{s}\frac{\partial }{\partial s}+ \frac{u_{\phi }}{s} \frac{\partial }{\partial \phi }\right)\omega _{z}+ F_{L} \, . \end{eqnarray} (7) Here, $$Ra^{*}= E^{2} RaP_{r}^{-1}$$ is the modified Rayleigh number, $$E = \nu (\Omega r_{2}^{2} )^{-1}$$ is the Ekman number, $$Ra= \alpha g_{0} \Delta T r_{2}^{3} \left(\nu \kappa \right)^{-1}$$ is the Rayleigh number, Pr = νκ−1 is the Prandtl number, and ν, α, g0, and κ are respectively the kinematic viscosity, thermal expansion coefficient, gravitational acceleration at r = r2, and thermal diffusivity. Θ is the local temperature perturbation from the conducting profile (see Section 2.5), and FL represents the axial component of the curl of the Lorentz force. Since we wish to study the dynamics of zonal flows, we write a separate equation for the evolution of the axisymmetric angular velocity $$\frac{\overline{u_{\phi }}}{s}$$:   $$\frac{\partial }{\partial t}\left(\frac{\overline{u_{\phi }}}{s}\right)= \Gamma _{\textrm{L}}+ \Gamma _{\textrm{R}}+ \Gamma _{\textrm{V}}\, ,$$ (8)where ΓL denotes the axial torque from Lorentz forces, ΓR the axial torque from Reynolds stresses, and ΓV the axial viscous torque. ΓR and ΓV are given by   \begin{eqnarray} \Gamma _{\textrm{R}}= - \frac{1}{s} \left(\overline{\frac{u_{s}}{s} \frac{\partial }{\partial s}s u_{\phi }} \right), \end{eqnarray} (9)  \begin{eqnarray} \Gamma _{\textrm{V}}=\phantom{-} \frac{E}{s^{3} L} \frac{\partial }{\partial s}\left(s^{3} L \frac{\partial }{\partial s}\left(\frac{\overline{u_{\phi }}}{s}\right)\right), \end{eqnarray} (10)while ΓL is given in the next section by eq. (14). Although eq. (8) is contained in eq. (7), in practice we use eq. (7) to evolve the non-axisymmetric part of ωz, and eq. (8) to evolve $$\overline{u_{\phi }}$$. We note that we have neglected in both eqs (7) and (8) viscous friction at the top and bottom spherical boundaries. These can be incorporated in the QG framework through Ekman pumping terms (e.g. Schaeffer & Cardin 2005). However, for the Ekman numbers that we can reach numerically, these friction terms would play an unrealistically dominant role in the force balance, unlike in the Earth’s core. 2.4 Induction equation and Lorentz force Since a 2-D model cannot maintain a self-sustaining dynamo (e.g. Cowling 1957; Roberts 2015), we instead impose a steady background magnetic field $${\bf B_{0}}$$ on it. For simplicity we choose a uniform field, $${\bf B_{0}}(s,\phi ,z) = B_{0s}\, \mathbf {e}_{s}$$. This uniform field can be interpreted to represent the local cumulative effect from all length scales of the background field. Two-dimensional perturbations (bs, bϕ) from that background state are then tracked via the induction equation. We form an equation for the magnetic potential a by axially averaging the s-component of the induction equation, which yields   $$\frac{\partial }{\partial t}a = - u_{\phi }B_{0s}+ \left(u_{s}b_{\phi }- u_{\phi }b_{s}\right)+ \frac{E}{P_{m}} \left(\nabla _{H}^{2} a + \frac{2 \beta a}{s} \right)\, ,$$ (11)where Pm = νη−1 is the magnetic Prandtl number and η is the magnetic diffusivity. The β factor in the last term of eq. (11) originates from the contribution of bϕ to the vector Laplacian in the s-direction. Physically, it represents an added contribution to dissipation introduced by the spherical geometry. The equation for $$\overline{b_\phi }$$ is derived from the axially averaged and azimuthally averaged ϕ-component of the induction equation,   \begin{eqnarray} \frac{\partial }{\partial t}\left(\frac{\overline{b_{\phi }}}{s}\right)&=& B_{0s}\frac{\partial }{\partial s}\left(\frac{\overline{u_{\phi }}}{s}\right)+ \frac{1}{s}\frac{1}{L} \frac{\partial }{\partial s}\Big ( L \left(\overline{u_{\phi }b_{s}} - \overline{u_{s}b_{\phi }} \right)\Big ) \nonumber\\ &&+\,\frac{1}{s^{3}} \frac{E}{P_{m}} \frac{\partial }{\partial s}\left(s^{3} \frac{\partial }{\partial s}\right)\left(\frac{\overline{b_{\phi }}}{s}\right)\, . \end{eqnarray} (12)The axially averaged Lorentz force which enters eq. (7) is   $$F_{L} = \left(\left(B_{0s}+ b_{s}\right)\frac{\partial }{\partial s}+ \frac{b_{\phi }}{s} \frac{\partial }{\partial \phi }\right)j_{z}\, .$$ (13) Although eq. (13) is the correct expression for the Lorentz force acting on flow eddies in our QG model, we note that in practice we do not compute this term, that is, we set FL = 0. The justification for this is given in Section 2.7. Meanwhile, the axial Lorentz torque in eq. (8) is   $$\Gamma _{\textrm{L}}= \Gamma _{\textrm{L}_{1}}+ \Gamma _{\textrm{L}_{2}}\, ,$$ (14)where   \begin{eqnarray} \Gamma _{\textrm{L}_{1}} = \frac{1}{s^{3}} \frac{1}{L} \frac{\partial }{\partial s}\left(s^{3} L B_{0s}\frac{\overline{b_{\phi }}}{s}\right)\, , \end{eqnarray} (15)  \begin{eqnarray} \Gamma _{\textrm{L}_{2}} = \frac{1}{s} \left(\overline{\frac{b_{s}}{s} \frac{\partial }{\partial s}s b_{\phi }} \right)\, . \end{eqnarray} (16) 2.5 Thermal equations Thermal convection is driven by the steady, superadiabatic temperature difference ΔT = T1 − T2, where T1 and T2 are the superadiabatic temperatures on the spheres r = r1 and r = r2, respectively. Taking the axial average of the resulting 3-D temperature profile T0 gives the background conducting, rigid temperature profile T as a function of cylindrical radius:   $$T = \left(\frac{r_{2} T_{2} - r_{1} T_{1}}{r_{2} - r_{1}} \right)+ \left(\frac{r_{1} r_{2}}{r_{2} - r_{1}} \right)\frac{\Delta T}{L} \ln \left(1 + L \right)\, .$$ (17)Since it is implicitly contained in the control parameter Ra*, we are free to set ΔT to an arbitrary non-zero value. We therefore choose T1 = 1, T2 = 0. These choices, along with eq. (17), are consistent with those used by Aubert et al. (2003) and Gillet & Jones (2006). Time-dependent perturbations Θ from T are tracked with the axially averaged heat equation   $$\frac{\partial \Theta }{\partial t} = -\left(\mathbf {u}_{H} \cdot \boldsymbol{\nabla }_{H} \right)\left(T + \Theta \right)+ \frac{E}{P_{r}} \nabla ^{2}_{H} \Theta \, ,$$ (18)where $$\boldsymbol{\nabla }_{H}$$ represents the (s, ϕ) components of the gradient operator. 2.6 Boundary conditions We ensure that the radial magnetic perturbation field bs drops to zero at the boundaries by imposing a = 0, and we also impose $$\frac{\overline{b_{\phi }}}{s}= 0$$. Similarly, we ensure that the temperature perturbation drops to zero on the boundary by imposing Θ = 0. To respect the no-penetration boundary condition, ψ must be constant on the inner and outer boundaries. For convenience, we use ψ = 0. We apply an additional no-slip boundary condition on the non-axisymmetric part of uϕ such that $$\frac{\partial \psi }{\partial s} = 0$$. Forcing us and uϕ to be both zero at the boundary is self-consistent with the a = 0 condition (see eq. 11) so it is the most natural choice. However, we impose a free-slip boundary condition on $$\frac{\overline{u_{\phi }}}{s}$$, such that $$\frac{\partial }{\partial s}(\frac{\overline{u_{\phi }}}{s})=0$$. The latter choice is made to minimize viscous friction effects on the zonal flows—the target of our study—which should be small in an Earth-like regime. 2.7 Parameter regime As argued in the introduction, in order to be in an Earth-like regime ($$\mathcal {A}\ll 1$$), typical convective flow speeds (uC) must be much smaller than the typical Alfvén wave speed (uA). For a given Ekman number E, convective speeds are controlled by Ra*. With our choice of non-dimensionalization, the Alfvén wave velocity in eq. (2) is   $$u_{A}= \left| {\bf B} \right| = \left| {\bf B}_{0} + \mathbf {b} \right| \, ,$$ (19)so a lower bound is set by the strength of the background magnetic field B0s. Therefore, the combination of Ra* and B0s in our model must be such that $$\mathcal {A}\ll 1$$. Furthermore, we must be in a regime where accelerations of the zonal flow are dominantly controlled by the Lorentz, rather than the Reynolds, torque. An inspection of eqs (9) and 16) shows that the Reynolds and Lorentz torques are proportional to the non-axisymmetric parts of u2 and b2, respectively, with typical scales given by their root-mean-square amplitudes $$u_{\rm {\small RMS}}= \left| \mathbf {u}- \overline{u_{\phi }} \mathbf {e}_{\phi } \right| \approx u_{C}$$ and $$b_{\rm {\small RMS}}= \left| \mathbf {b}- \overline{b_{\phi }} \mathbf {e}_{\phi } \right|$$. Hence, we must be in a regime where $$b_{\rm {\small RMS}}\gg u_{\rm {\small RMS}}$$. The parameters controlling the ratio between $$b_{\rm {\small RMS}}$$ and $$u_{\rm {\small RMS}}$$ can be established by analysing the induction equation. When our model has achieved statistical equilibrium, the source and diffusion terms must be in balance:   $${\left| \boldsymbol{\nabla }\times \left(\left({\bf B}_{0} + \mathbf {b}\right)\times \mathbf {u}\right) \right| = \left| \frac{E}{P_{m}} \boldsymbol{\nabla }_{H}^{2} \mathbf {b} \right|.}$$ (20) Using local flow length scale ℓ1 and magnetic field length scale ℓ2, eq. (20) scales as   $$\frac{\left| {\bf B} \right| \left| \mathbf {u} \right|}{\ell _{1}} = \frac{E}{P_{m}} \frac{\left| \mathbf {b} \right|}{\ell _{2}^{2}} \quad \Rightarrow \quad \frac{\left| \mathbf {b} \right|}{\left| \mathbf {u} \right|} = \frac{\left| {\bf B} \right|}{\ell _{1}} \frac{P_{m}\ell _{2}^{2}}{E} \, .$$ (21) The Alfvén timescale can be defined as τA = ℓ1/uA, and the magnetic diffusion timescale as $$\tau _{D} = P_{m}\ell _{2}^{2} / E$$. The ratio Lu = τD/τA is then the Lundquist number. Thus, by using $$\left| \mathbf {u} \right| \approx u_{\rm {\small RMS}}$$, $$\left| \mathbf {b} \right| \approx b_{\rm {\small RMS}}$$, and $$u_{A}= \left| {\bf B} \right| \approx B_{0s}$$, eq. (21) becomes   $$\frac{b_{\rm {\small RMS}}}{u_{\rm {\small RMS}}} = \frac{B_{0s}}{\ell _{1}} \frac{P_{m}\ell _{2}^{2}}{E} = \frac{\tau _{D}}{\tau _{A}} = \mathrm{Lu} \, .$$ (22) Hence, the requirement of $$b_{\rm {\small RMS}}\gg u_{\rm {\small RMS}}$$ implies that we must be in a regime where Lu ≫ 1, which in turn places constraints on the combinations of E, Pm and B0s which we may choose. For an Earth-like setup, we would ideally have Pm ≪ 1, which requires that we pick a sufficiently large ratio of B0s/E. In practice, however, numerical constraints limit both the maximum possible value of B0s and the minimum possible value of E. In addition, it is desirable to keep uA < 1, so that the Alfvén wave speed is smaller than the inertial wave speed and our model remains consistent with our QG assumption on rigid flows. Thus, for a numerically achievable ratio of B0s/E, while we may choose Pm < 1, Pm must remain sufficiently large such that solutions are in a regime characterized by Lu ≫ 1. The final aspect of the convective regime to be addressed concerns the influence of the Lorentz force FL on the vorticity equation of eq. (7). Left to evolve dynamically, flow and magnetic field lines tend to align themselves so as to limit induction by shear (e.g. Schaeffer et al.2017). However, because we impose a steady background magnetic field in our model, such a fully dynamic reorganization is not possible. As a result, for the parameters that we use (see next section), convective eddies are unduly constrained by the Lorentz force and become thin, elongated columns in the s-direction, sharing little resemblance with the large-scale flows in Earth’s core. However, when we set FL = 0 in eq. (7) and ignore the influence of the Lorentz force on the non-axisymmetric flows, we retrieve Earth-like large-scale eddies. Because our main goal is to illustrate the mechanism by which Earth-like, large-scale flows can interact with the background field to generate slow zonal accelerations, our model is a better analogue to Earth’s core with FL turned off. In other words, although convective eddies in our model are still allowed to distort the magnetic field, we do not take into account the feedback of the Lorentz force on the flow eddies. We retain the Lorentz torque ΓL in the axially symmetric torque balance of eq. (8), as it is our primary objective to show that it can be the dominant driver of rigid accelerations. 3 RESULTS 3.1 Solution scheme A semi-spectral method with 768 Fourier modes in azimuth and radial derivatives approximated by second-order finite differences between 901 points arranged on a Chebyshev grid in radius was used to discretize eqs (7), (8), (11), (12) and (18). The resulting discrete equations were then evolved in timesteps of 5 × 10−4 using a combination of a Crank–Nicolson method for the linear terms and a second-order Adams–Bashforth scheme for the nonlinear terms (e.g. He & Sun 2007). Use of a Chebyshev grid ensured fine enough spacing to resolve boundary layers (proportional in thickness to E1/2 ≈ 0.002), while maintaining reasonable computation times via coarser spacing away from the boundaries. In our setup, minimum spacing near the boundaries is 2 × 10−6, while maximum spacing in the centre of the model domain is 1 × 10−3. We set s1 = 0.35, mimicking the thickness ratio of Earth’s core, and s2 = 0.98. Limiting our domain to s2 = 0.98 instead of s2 = 1 is convenient, as it allows us to use a slightly larger grid space and timesteps (since β → ∞ as s → 1). The results that we now present come from an experiment with E = 5.0 × 10−6, Ra = 5.0 × 108, B0s = 0.15, Pm = 0.1, and Pr = 1.0. With these parameters, a numerical simulation started from small perturbations typically takes about 1000 non-dimensional time units (or about 160 full rotations) to reach statistical equilibrium. 3.2 General features of the flow and induced magnetic fields Fig. 3 shows a snapshot in time of the non-axisymmetric ωz and jz after the model’s global energy budget has reached statistical equilibrium. The presence of eddies in both the flow and the perturbed magnetic field are clearly visible. The magnetic field perturbations exhibit structures with larger wavelengths than those in the flow field, and with smoother features, a result of the more rapid magnetic to viscous diffusion (Pm = 0.1). Figure 3. View largeDownload slide Snapshots in time of (a) non-axisymmetric vorticity ωz and (b) non-axisymmetric axial current jz, as seen looking downward from the north pole. Figure 3. View largeDownload slide Snapshots in time of (a) non-axisymmetric vorticity ωz and (b) non-axisymmetric axial current jz, as seen looking downward from the north pole. Magnetic field perturbations result from the action of convective eddies shearing and advecting the sum of the background and perturbed magnetic field. Time-dependency in both the flow and the magnetic field leads to fluctuations in the Reynolds (ΓR, eq. 9) and Lorentz (ΓL, eq. 14) torques, respectively. These must be accommodated by zonal accelerations. Typical $$u_{\rm {\small RMS}}$$ and $$b_{\rm {\small RMS}}$$ values after equilibration are 0.0024 and 0.14, respectively. The amplitude of magnetic field perturbations in our numerical experiment is of the same order as the imposed background field, $$\left| \mathbf {b} \right| \approx \left| {\bf B}_{0} \right|$$. The Alfvén and Lundquist numbers of our simulation are then $$\mathcal {A}\approx 0.014$$ and Lu ≈ 70, within the region of parameter space where we expect ΓL ≫ ΓR. This is the region where the decadal timescale dynamics of zonal accelerations in Earth’s core should reside. 3.3 Time-averaged axisymmetric force balance The dashed black line of Fig. 4 shows the time-averaged zonal angular velocity. Its profile is dominated by a shear flow spanning the whole of the modelled region, retrograde at the outer boundary and prograde at the inner boundary. The amplitude of this shear flow is of the same order of magnitude as the amplitude of typical convective eddies. Figure 4. View largeDownload slide Time-averaged mean axial torques (solid lines) and time-averaged zonal angular velocity (dashed line) as a function of radius. Figure 4. View largeDownload slide Time-averaged mean axial torques (solid lines) and time-averaged zonal angular velocity (dashed line) as a function of radius. The mean (time-averaged) zonal flow of Fig. 4 differs from the one typically observed in the absence of a magnetic field. In non-magnetic convection with stress-free boundaries, the direction of the mean zonal flow is reversed – it is prograde at the equator – and its characteristic velocity is much larger than the velocities associated with convective eddies. It results from time-averaged Reynolds stresses, themselves being the product of the topographic beta effect acting on convective eddies (e.g. Cardin & Olson 1994; Christensen 2002). In the presence of a radial magnetic field, distortion of this magnetic field by the shear flow induces a restoring Lorentz force (through $$\Gamma _{\rm{L}_{1}}$$ of eq. 15) which limits the growth of the mean zonal flow. In 3-D models, this mean zonal flow is not z-invariant (e.g. Aubert 2005). In our QG model, because our magnetic field perturbation is defined with a built-in topographic beta effect identical to that of the flow, a time-averaged Maxwell stress ($$\Gamma _{\rm{L}_{2}}$$, eq. 16) is maintained in the same way as the time-averaged Reynolds stress ΓR. Since $$\Gamma _{\rm{L}_{2}}$$ has the same form, but opposite sign, as ΓR, and since the magnitude of the former typically dominates that of the latter in our model, the direction of the driven mean zonal flow is reversed to that produced in non-magnetic convection. The time-averaged $$\overline{u_\phi }$$ profile, then, is the result of a balance between the time-averaged torques. Fig. 4 illustrates this balance, showing the time-averaged Reynolds (ΓR, orange), Lorentz (ΓL, blue), and viscous (ΓV, green) torques as a function of radius. In the interior of the domain, the Reynolds and Lorentz torques largely balance one another. The viscous torque becomes more important near the boundaries, especially the outer one. The enhancement of all three torques near the outer boundary is caused by the large β-effect: both the viscous and the $$\Gamma _{\rm{L}_{1}}$$ part of the Lorentz torque depends explicitly on β (through the s-derivative of L, see eqs 10 and 15), while the Reynolds (ΓR, eq. (9)) and Maxwell ($$\Gamma _{\rm{L}_{2}}$$, eq. (16)) torques implicitly involve β through their dependence on uϕ and bϕ, respectively. We further note that both $$\Gamma _{\rm{L}_{1}}$$ and $$\Gamma _{\rm{L}_{2}}$$ are individually much larger than ΓL, as shown by Fig. 5. However, they tend to cancel one another, leaving a net Lorentz torque several orders of magnitude smaller than either one individually. This self-cancellation will be re-examined in Section 3.5. Figure 5. View largeDownload slide Time-averaged components of the mean axial Lorentz torque as a function of radius. Figure 5. View largeDownload slide Time-averaged components of the mean axial Lorentz torque as a function of radius. 3.4 Zonal accelerations Fluctuations in time with respect to the time-averaged torque balance are the main focus of our study. Fig. 6 shows the axisymmetric angular accelerations (top panel) and the time-varying parts of ΓL, ΓR, and ΓV (bottom three panels) after the system has equilibrated. Fluctuations in ΓL have a typical amplitude five times larger than those of ΓR: ∼ 6.6 × 10−5 versus ∼1.2 × 10−5. ΓV only plays a small role in the time-dependent dynamics, with RMS fluctuations of the order of ∼0.3 × 10−5. Figure 6. View largeDownload slide The angular acceleration (top panel), and the time-dependent parts of the Lorentz (second panel), Reynolds (third panel) and viscous (bottom panel) torques, as functions of cylinder radius and time. The time-averaged contribution to each torque (shown in Fig. 4) has been removed. Figure 6. View largeDownload slide The angular acceleration (top panel), and the time-dependent parts of the Lorentz (second panel), Reynolds (third panel) and viscous (bottom panel) torques, as functions of cylinder radius and time. The time-averaged contribution to each torque (shown in Fig. 4) has been removed. Thus, fluctuations in our model’s zonal accelerations are mainly controlled by the Lorentz torque. Indeed, the upper two panels of Fig. 6 suggest a strong correlation at all times and radii between $$\frac{\partial }{\partial t}(\frac{\overline{u_{\phi }}}{s})$$ and ΓL. Time on Fig. 6 has been scaled to the timescale of the fundamental Alfvén wave mode, τA = 2(s2 − s1)/uA. The joint fluctuations in the zonal accelerations and Lorentz torque cover a broad range of timescales but are dominated by periods which fall in the range of 5 to 10 times longer than τA. This slower timescale, τslow ≈ 5 − 10 τA, reflects the time-fluctuations of the magnetic field, themselves the result of induction by convective flows. These occur on a longer timescale than the Alfvén wave propagation timescale, as we expect for a regime with $$\mathcal {A}< 1$$. Scaling the temporal fluctuations on Fig. 6 to Earth’s core, taking τA ≈ 6 yr, gives a timescale τslow of 30 to 60 yr for these magnetically driven zonal accelerations, similar to the zonal accelerations inferred within Earth’s core. To further demonstrate that the slow magnetic field fluctuations originate from the underlying convective dynamics, we compute a characteristic azimuthal wave number k* (e.g. Takahashi et al.2008) at each radius s from the convolution of wavenumber k and convective speed u(k):   $$k^{*}(s) = \frac{\int _{}^{} \! k u(k) \, \mathrm{d} k}{\int _{}^{} \! u(k) \, \mathrm{d} k} \, .$$ (23) The result of this calculation as a function of radius, for the snapshot corresponding to time equal 0 on Fig. 6, is shown in panel (a) of Fig. 7. The characteristic convective length scale ℓC(s) is then calculated from the characteristic wavenumber as   $$\ell _{C}(s) = \frac{2 \pi s}{k^{*}(s)} \, .$$ (24) The characteristic convective length scale for all radial positions is shown in panel (b) of Fig. 7. Dividing these length scales by the RMS velocities at each radius yields a characteristic convective timescales τC, which is representative of the time required for the flow to create significant changes in the magnetic field:   $$\tau _{C}(s) = \frac{\ell _{C}(s)}{u_{\rm {\small RMS}}} \, .$$ (25) The dependence of τC with radius is shown in panel (c) of Fig. 7. All three plots of Fig. 7 change for different time snapshots, but are representative of the general behaviour at all times. Furthermore, they only show the dominant length and time scales of the flow, which is truly characterized by a spectrum of scales. As shown in Fig. 7(c), in the bulk of the fluid τC is approximately equal to 5 − 8 τA, within the same range as τslow. These are obviously simple order of magnitude calculations that do not capture the complex nonlinear connection between large-scale flows and time changes in the Lorentz torque. Nevertheless, the combination of Figs 6 and 7 show that it is possible to drive slow zonal accelerations by fluctuating Lorentz torques, themselves driven by underlying convective flows shearing and advecting the magnetic field. Figure 7. View largeDownload slide Characteristic (a) wave numbers, (b) length scales and (c) timescales as a function of radius at a given snapshot in time. Figure 7. View largeDownload slide Characteristic (a) wave numbers, (b) length scales and (c) timescales as a function of radius at a given snapshot in time. In addition to the slow fluctuations of Fig. 6, the convective dynamics also generate free Alfvén waves. Fig. 8 shows the second half of the zonal acceleration and Lorentz torque panels of Fig. 6 after applying a highpass filter to remove fluctuations with periods longer than 1.19 τA. This reveals the presence of periodic fluctuations with a period of approximately 1 τA. The correlation between the acceleration and Lorentz torque shows that these are indeed Alfvén waves. (Applying the same highpass filter to ΓR and ΓV produces only low-amplitude noise.) Figure 8. View largeDownload slide Output of a highpass filter applied to the second half of the zonal angular accelerations (top panel) and Lorentz torques (bottom panel) of Fig. 6. In our non-dimensional time units, the eighth-order digital Butterworth filter’s −3 dB frequency is 1/10. Figure 8. View largeDownload slide Output of a highpass filter applied to the second half of the zonal angular accelerations (top panel) and Lorentz torques (bottom panel) of Fig. 6. In our non-dimensional time units, the eighth-order digital Butterworth filter’s −3 dB frequency is 1/10. Left on their own, free Alfvén oscillations should decay away because of ohmic dissipation. In our model, they are continuously re-excited by the underlying convective dynamics, though resonant amplification remains modest and their amplitude does not rise much higher than that resulting from the forced background accelerations. Their typical RMS velocities of 1 × 10−4 are approximately 7 times smaller than the RMS velocities of 7 × 10−4 found for the slower zonal flow fluctuations. This is the reason why Alfvén waves, though present, are not apparent on Fig. 6. 3.5 Taylorization For a system in a perfect Taylor state, even though the Lorentz torque at any point may be large, cancellations occur such that the Lorentz torque integrated over a geostrophic cylinder is equal to zero. The ‘Taylorization’ is a measure of the degree of this Lorentz torque cancellation. In our model, where the magnetic field (and therefore the Lorentz torque) is assumed to be axially invariant, the Taylorization is measured by the factor $$\mathcal {T}$$:   $$\mathcal {T}(s) = \frac{\left| \oint _{}^{} \! \mathcal {M}_{\phi }(s,\phi ) \, \mathrm{d} \phi \right|}{\oint _{}^{} \! \left| \mathcal {M}_{\phi }(s,\phi ) \right| \, \mathrm{d} \phi } = \frac{ \left| \Gamma _{\textrm{L}}(s) \right| }{ \oint _{}^{} \! \left| \mathcal {M}_{\phi }(s,\phi ) \right| \, \mathrm{d} \phi } \, ,$$ (26)where $$\mathcal {M}_{\phi }(s,\phi )$$ represents the azimuthal component of the local magnetic force. A system with low Taylorization has $$\mathcal {T} \lesssim 1$$, while a system with a high Taylorization has $$\mathcal {T} \ll 1$$. Fig. 9 shows the Taylorization factor of our model as a function of cylindrical radius and time. Typical RMS values of $$\mathcal {T}$$ are approximately 3.2 × 10−4, with peak values of 2.6 × 10−3. As shown in the previous section, temporal fluctuations of the magnetic field lead to temporal fluctuations of the Lorentz torque. Thus, Taylor’s constraint is continually broken, with the Lorentz torque fluctuations being accommodated by rigid zonal accelerations. This can be observed in Fig. 6, where the times and radii of the largest zonal accelerations often coincide with the largest values of $$\mathcal {T}$$ in Fig. 9. Figure 9. View largeDownload slide Taylorization factor $$\mathcal {T}$$ as a function of radius and time. With our choice of boundary conditions, the Taylorization factor equals 1 at the boundaries. Only s ∈ [0.36, 0.97] is shown. Figure 9. View largeDownload slide Taylorization factor $$\mathcal {T}$$ as a function of radius and time. With our choice of boundary conditions, the Taylorization factor equals 1 at the boundaries. Only s ∈ [0.36, 0.97] is shown. Therefore, our model is not in a perfect Taylor state at any given moment, but instead fluctuates about an equilibrium state characterized by a low Taylorization factor. A part of the cancellation of the Lorentz torque over any given geostrophic cylinder is inherent to the form of the Maxwell torque ($$\Gamma _{\rm{L}_{2}}$$), as it involves products of the magnetic field vectors which change direction at different azimuthal points. On a time-average, $$\Gamma _{\rm{L}_{2}}$$ does not vanish but it is mostly opposed by $$\Gamma _{\rm{L}_{1}}$$, as shown by Fig. 5. Yet, $$\Gamma _{\rm{L}_{1}}$$ and $$\Gamma _{\rm{L}_{2}}$$ do not cancel exactly: the time-averaged Taylorization factor is small ($$\mathcal {T} \sim 10^{-4}$$) but not zero. In our model, the departure from a time-averaged Taylor state is dominantly caused by the persistent torque from Reynolds stresses (and from viscous forces near the domain boundaries, see Fig. 4). 4 DISCUSSION 4.1 Slow zonal accelerations and decadal timescale dynamics in Earth’s core Our model shows that in an Earth-like parameter regime characterized by $$\mathcal {A}\ll 1$$ and Lu ≫ 1, convective flows can drive fluctuations in the magnetic field which then lead to temporal fluctuations in the Lorentz torque. The latter generate zonal accelerations in the form of free Alfvén waves, but also slower forced zonal accelerations on a timescale connected to the convective flows. In order to drive large-scale, core-size, zonal accelerations, convective eddies must be large-scale themselves. Because the typical velocity of these large-scale eddies is slower than the Alfvén wave velocity, the timescale associated with the magnetic field change they produce, and thus the Lorentz torque, is slower than τA, the timescale of Alfvén wave propagation across the core. In our model, the typical convective timescale associated with large-scale eddies is 5 to 8 times longer than τA. These eddies induce Lorentz torque fluctuations over a broad range of timescales, but dominated by periods that are 5 to 10 times longer than τA. Scaling our results to Earth’s core by taking τA = 6 yr, our model produces zonal accelerations with typical timescales in the range of 30 to 60 yr, comparable to observations. Furthermore, the power spectrum of flows reconstructed from geomagnetic secular variation shows peaks at spherical harmonic degrees 8–10 (e.g. fig. 5 of Gillet et al.2015). At mid-core radius, this corresponds to a typical length scale of approximately 600 km. Taking a typical velocity of 10 km yr−1, these are associated with a typical convective timescale of approximately 60 yr. If these flows rigidly extend axially through the core, according to our mechanism they should distort the magnetic field, and thus produce fluctuations in the Lorentz torque and zonal accelerations in response, at approximately the same 60 yr timescale. Our results suggest that the inferred decadal, large-scale zonal accelerations of geostrophic cylinders in the Earth’s core can be explained by forced fluctuations of the Lorentz torque, themselves driven by the large-scale convective flows. As observed in Earth’s core, our dynamical model contains large-scale flows, slow zonal accelerations, and faster Alfvén waves. In our model, their typical non-dimensional velocities are, respectively, 2.4 × 10−3, 7 × 10−4 and 1 × 10−4, for a relative ratio between them of 34:7:1. In Earth’s core, these typical velocities are 10 km yr−1, 2 km yr−1 and 0.2 km yr−1 (e.g. Gillet et al.2015), for a ratio of 50:10:1. Therefore, although their amplitude ratios do not match those of Earth’s core exactly, our model reproduces the correct ordering between these flows. This gives us further confidence that our model is a good analogue for the decadal flow dynamics in Earth’s core. Although our dynamical model is in the correct parameter regime in terms of $$\mathcal {A}\ll 1$$ and Lu ≫ 1, other parameters remain far from Earth-like, notably the Ekman number and magnetic Prandtl number. Care must then be taken when extrapolating our results to Earth’s core. To further confirm that our results appropriately capture the decadal timescale dynamics of zonal flows in the core, it would be desirable to carry many more numerical experiments, systematically varying some of the input parameters in order to develop scaling properties for our model. We plan to do this in a future study. A different approach is to try to do a similar analysis in a 3-D geodynamo model, some of which are approaching the parameter regime of $$\mathcal {A}\ll 1$$ and Lu ≫ 1 which we have highlighted (e.g. Aubert et al.2017; Schaeffer et al.2017). Short of doing this, we may speculate on how the use of more Earth-like Ekman and magnetic Prandtl numbers would affect the dynamics of zonal flows. For the relatively large Ekman number of our numerical experiment, viscous forces remain important in the establishment of the large-scale, quasi-core-size eddies that are visible on Fig. 3 (e.g. Aurnou et al.2015). However, at E ≈ 10−15, as inferred for Earth’s core, the E1/3 scaling law suggests that the typical length scale of viscously controlled eddies should be much smaller, about 10−5 times the core size, or of the order of 10-100 m in width. Such small eddies would be inefficient at generating the core-size changes in the magnetic field that are required to drive large-scale zonal accelerations. However, large-scale eddies are present in Earth’s core, as inferred from the geomagnetic secular variation. The typical length-scale of eddies in the core must then be controlled by processes other than this viscous scaling. Two main reasons can be invoked. First, for sufficiently low Ekman numbers, small-scale structures in the core likely feed their energy to larger length scales via an inverse energy cascade driven by geostrophic turbulence (e.g. Guervilly et al.2014; Stellmach et al.2014). Second, the influence of the Lorentz force should also lead to larger convective length scales (e.g. Roberts & King 2013), a regime which is beginning to be accessible in numerical simulations (e.g. Matsui et al.2014; Yadav et al.2016; Schaeffer et al.2017). A combination of these two effects may be responsible for forming the large structures seen in the core. Both are absent in our model. Therefore, although the dynamics that sustain the large-scale eddies in our numerical model are most likely not Earth-like, the large-scale eddies themselves share some Earth-like qualities. The key ingredient for driving decadal zonal accelerations according to our mechanism is the very presence of large-scale flow eddies with decadal convective timescale, not precisely how these eddies are generated. The magnetic Prandtl number in the core is much lower than the one we have chosen in our numerical experiment. Since a lower Pm corresponds to enhanced magnetic diffusion, temporal changes in the magnetic field would be dominantly controlled by the largest length scales of the underlying convective flow. Thus, while small-scale eddies certainly exist in the core, changes in the magnetic field should still preferentially occur at the largest length scale. We therefore expect that temporal changes in the Lorentz torque can drive decadal zonal accelerations at core-size wavelengths in the radial direction, just as we observe in our numerical experiment of limited spatial resolution. 4.2 Free Alfvén waves As shown in Fig. 8, our dynamical model excites free Alfvén waves. However, their spatio-temporal properties differ from those detected in Earth’s core. In our model, they are dominated by a standing wave oscillation of the fundamental mode. In Earth’s core, they take the form of outward travelling waves (Gillet et al.2010, 2015). The reason why Alfvén waves in Earth’s core travel outward, as well as their excitation mechanism, remain unclear. Outward travelling Alfvén waves resulting from a quasi-periodic triggering near the tangent cylinder are observed in some numerical models (Teed et al.2014, 2015; Schaeffer et al.2017). When approaching an Earth-like regime, Lorentz forces are responsible for this torque, but the precise physical mechanism has not been clearly identified. It would be a valuable effort to investigate where and how Alfvén waves are excited in our model. Since we do not see a preferential propagation direction, excitation appears to be distributed evenly within the integration domain of our model. The region close to the tangent cylinder does not appear to be the seat of any form of recurrent instability, although this may be because we have not modelled the dynamics inside the tangent cylinder. Given that convective flows in our model induce changes in the magnetic field on a broad range of timescales, Alfvén waves on Fig. 8 may simply represent the resonant response to fluctuations of the Lorentz torque which occur in the vicinity of their free period range. Indeed, correlations between Figs 8 and 9 suggest that this is the case: notable increases in the Taylorization factor are often, though not always, associated with an amplitude enhancement of free Alfvén waves. This argues along the same line of a recent study which has shown that applying a stochastic forcing in the volume of the core readily excites Alfvén waves (Gillet et al.2017). Moreover, the study of Gillet et al. (2017) has also shown that electromagnetic dissipation at the CMB transforms standing Alfvén waves into outward traveling waves, with similar characteristics as those detected in the Earth’s core. Our model constitutes a dynamical realization of such a stochastic forcing and supports the idea that Alfvén waves in Earth core are simply the response to sub-decadal changes in the Lorentz torque within the bulk of the core. Our model does not include dissipation at the CMB, but we believe (though this should be tested) that adding it would also transform our standing Alfvén waves into outwardly propagating waves. 4.3 Taylorization Taylor’s constraint is continuously being broken in our model, generating zonal accelerations in response. The numerical value of the Taylorization factor associated with these fluctuations is of the order of 10−3. The time-average state (or statistical equilibrium) about which these fluctuations occur is characterized by a higher degree of Taylorization, with a Taylorization factor of the order of 10−4. Though this is small, it indicates that the equilibrium state in our model remains far from a perfect Taylor state. In the bulk interior, this is dominantly because Reynolds stresses lead to a non-zero time-averaged torque. Thus, the time-averaged Lorentz torque need not be zero (and obey a Taylor state), but only be as small as the persistent torque from Reynolds stresses. The results of our model suggest that, as was pointed out by Dumberry & Bloxham (2003), the torque from Reynolds stresses in Earth’s core may also play a leading order role in the departure from a time-averaged Taylor state. The ratio of the torque from Reynolds stresses to the Lorentz torque scales as (uC/uA)2, so it is proportional to the square of the Alfvén number. Because the Alfvén number in our model is similar to that in Earth’s core, the baseline Taylorization factor of approximately 10−4 in Fig. 9 may also be representative of that expected in Earth’s core. However, because our model is 2-D, extrapolating our results to the 3-D magnetic field of Earth is clearly not straightforward. Furthermore, Ekman friction, which we have neglected, could be important in the time-averaged torque balance, especially if turbulent processes lead to an enhanced effective viscosity. Ekman friction could balance a part of the torque from Reynolds stresses and the Taylorization factor could then be smaller than 10−4. Regardless of its exact value, the high degree of Taylorization requires that large cancellations in the equilibrium Lorentz torque over a cylinder must occur in Earth’s core. It this sense, exploring dynamo solutions in the limit of a vanishing Lorentz torque remains a worthy goal (e.g. Livermore et al.2008; Wu & Roberts 2015). This being said, one must keep in mind that the correct Taylorization factor in Earth’s core may not be asymptotically close to zero but instead be closer to 10−4. 5 CONCLUSION We have constructed a 2-D reduced model of rotationally dominated magnetoconvection capable of producing distinct short- and long-timescale accelerations in the zonal flow. The short-timescale accelerations are free Alfvén waves, while the long-timescale accelerations are magnetically forced through the evolution of the Lorentz torque. The temporal changes in the magnetic field which drive the time-varying Lorentz torque are produced by the underlying convective flows, shearing and advecting the magnetic field on a timescale associated with convective eddies. Our results provide a dynamical explanation for the rigid decadal zonal accelerations that are inferred to exist in Earth’s core on the basis of the magnetic field observations and changes in LOD. Our results offer an alternative to the recent suggestion that the decadal zonal accelerations may not reflect deep seated rigid flows but are instead free Magnetic–Archemedian–Coriolis (MAC) waves in a stratified layer at the top of the core (Buffett 2014; Buffett et al.2016; Jaupart & Buffett 2017). The zonal flows of such MAC waves are characterized by a shear in the radial direction: flow at the CMB does not reflect flow deeper in the core in the axial direction. It is then more difficult to build a prediction of LOD changes based solely on the knowledge of flows at the CMB. It is nevertheless possible that, when properly taking into account the coupling of flows in the bulk of the core with these MAC waves, a prediction of core angular momentum change may be constructed so as to match the observed LOD variations (Buffett et al.2016). Yet the very fact that a very good match exists between the observed changes in the LOD and those predicted on the basis of purely rigid zonal flows suggests that deviations from rigidity are limited. Moreover, as our dynamical model shows, convective dynamics are expected to drive temporal changes in deep seated rigid zonal flows at decadal timescales. If the top of the core is stably stratified, rigid zonal flows may drive forced MAC oscillations and/or excite free MAC waves. Hence, the zonal flows at the CMB may consist of a combination of deep seated rigid zonal flows and flows that obey a MAC balance. However, because of the good match in LOD based on purely rigid flows, we speculate that the non-rigid MAC flows, if present, make up only a fraction of the total zonal flow. ACKNOWLEDGEMENTS This work was supported by both an NSERC/CRSNG Post-Graduate Scholarship and a Discovery grant from NSERC/CRSNG. Numerical simulations were performed on computing facilities provided by WestGrid and Compute/Calcul Canada. Figures were created using the Matplotlib package for Python (Hunter 2007). We thank Nathanaël Schaeffer and an anonymous reviewer for their constructive comments and suggestions. We further express our gratitude to Nathanaël Schaeffer, who kindly shared a numerical QG code which we extended over the course of this project. REFERENCES Aubert J., 2005. Steady zonal flows in spherical shell dynamos, J. Fluid Mech. , 542, 53– 67. https://doi.org/10.1017/S0022112005006129 Google Scholar CrossRef Search ADS   Aubert J., Gillet N., Cardin P., 2003. Quasigeostrophic models of convection in rotating spherical shells, Geochem. Geophys. Geosyst. , 4, 1052, doi:10.1029/2002GC000456. https://doi.org/10.1029/2002GC000456 Google Scholar CrossRef Search ADS   Aubert J., Gastine T., Fournier A., 2017. Spherical convective dynamos in the rapidly rotating asymptotic regime, J. Fluid Mech. , 813, 558– 593. https://doi.org/10.1017/jfm.2016.789 Google Scholar CrossRef Search ADS   Aurnou J., Calkins M.A., Cheng J.S., King E.M., Nieves D., Soderlund K.M., Stellmach S., 2015. Rotating convective turbulence in Earth and planetary cores, Phys. Earth planet. Inter. , 246, 52– 71. https://doi.org/10.1016/j.pepi.2015.07.001 Google Scholar CrossRef Search ADS   Braginsky S.I., 1970. Torsional magnetohydrodynamic vibrations in the Earth’s core and variations in day length, Geomagn. Aeron. , 10, 1– 10. Buffett B.A., 2010. Tidal dissipation and the strength of the earth’s internal magnetic field, Nature , 468, 952– 954. https://doi.org/10.1038/nature09643 Google Scholar CrossRef Search ADS PubMed  Buffett B.A., 2014. Geomagnetic fluctuations reveal stable stratification at the top of the Earth’s core, Nature , 507, 484– 487. https://doi.org/10.1038/nature13122 Google Scholar CrossRef Search ADS PubMed  Buffett B.A., Knezek N., Holme R., 2016. Evidence for MAC waves at the top of earth’s core and implications for variations in length of day, Geophys. J. Int. , 204, 1789– 1800. https://doi.org/10.1093/gji/ggv552 Google Scholar CrossRef Search ADS   Busse F., Or A., 1986. Convection in a rotating cylindrical annulus: thermal Rossby waves, J. Fluid Mech. , 166, 173– 187. https://doi.org/10.1017/S0022112086000095 Google Scholar CrossRef Search ADS   Cardin P., Olson P., 1994. Chaotic thermal convection in a rapidly rotating spherical shell: Consequences for flow in the outer core, Phys. Earth planet. Inter. , 82, 235– 259. https://doi.org/10.1016/0031-9201(94)90075-2 Google Scholar CrossRef Search ADS   Chao B.F., Chung W.Y., Shih Z.R., Hsieh Y.K., 2014. Earth’s rotation variations: a wavelet analysis, Terra Nova , 26, 260– 264. https://doi.org/10.1111/ter.12094 Google Scholar CrossRef Search ADS   Christensen U., 2002. Zonal flow driven by strongly supercritical convection in rotating spherical shells, J. Fluid Mech. , 470, 115– 133. https://doi.org/10.1017/S0022112002002008 Google Scholar CrossRef Search ADS   Christensen U., Wicht J., 2015. Numerical dynamo simulations, in Treatise on Geophysics , chap. 8.10, pp. 245– 277, ed. Schubert G., Elsevier. Christensen U.R., Aubert J., 2006. Scaling properties of convection-driven dynamos in rotating spherical shells and application to planetary magnetic fields, Geophys. J. Int. , 166, 97– 114. https://doi.org/10.1111/j.1365-246X.2006.03009.x Google Scholar CrossRef Search ADS   Cowling T., 1957. The dynamo maintenance of steady magnetic fields, Q. J. Mech. Appl. Math. , 10, 129– 136. https://doi.org/10.1093/qjmam/10.1.129 Google Scholar CrossRef Search ADS   Dumberry M., Bloxham J., 2003. Torque balance, Taylor’s constraint and torsional oscillations in a numerical model of the geodynamo, Phys. Earth planet. Inter. , 140, 29– 51. https://doi.org/10.1016/j.pepi.2003.07.012 Google Scholar CrossRef Search ADS   Finlay C.C., Dumberry M., Chulliat A., Pais M., 2010. Short timescale core dynamics: theory and observations, Space Sci. Rev. , 155, 177– 218. https://doi.org/10.1007/s11214-010-9691-6 Google Scholar CrossRef Search ADS   Gillet N., Jones C.A., 2006. The quasi-geostrophic model for rapidly rotating spherical convection outside the tangent cylinder, J. Fluid Mech. , 554, 343– 370. https://doi.org/10.1017/S0022112006009219 Google Scholar CrossRef Search ADS   Gillet N., Jault D., Canet E., Fournier A., 2010. Fast torsional waves and strong magnetic field within the Earth’s core, Nature , 465, 74– 77. https://doi.org/10.1038/nature09010 Google Scholar CrossRef Search ADS PubMed  Gillet N., Schaeffer N., Jault D., 2011. Rationale and geophysical evidence for quasi-geostrophic rapid dynamics within the Earth’s outer core, Phys. Earth planet. Inter. , 187, 380– 390. https://doi.org/10.1016/j.pepi.2011.01.005 Google Scholar CrossRef Search ADS   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. , 120, 3991– 4013. https://doi.org/10.1002/2014JB011786 Google Scholar CrossRef Search ADS   Gillet N., Jault D., Canet E., 2017. Excitation of traveling torsional normal modes in an Earth’s core model, Geophys. J. Int. , 210, 1503– 1516. https://doi.org/10.1093/gji/ggx237 Google Scholar CrossRef Search ADS   Gross R., 2015. Earth rotation variations - long period, in Treatise on Geophysics , chap. 3.09, pp. 215– 261, ed. Schubert G., Elsevier. Guervilly C., Cardin P., 2016. Subcritical convection of liquid metals in a rotating sphere using a quasi-geostrophic model, J. Fluid Mech. , 808, 61– 89. https://doi.org/10.1017/jfm.2016.631 Google Scholar CrossRef Search ADS   Guervilly C., Hughes D.W., Jones C.A., 2014. Large-scale vortices in rapidly rotating Rayleigh–Bénard convection, J. Fluid Mech. , 758, 407– 435. https://doi.org/10.1017/jfm.2014.542 Google Scholar CrossRef Search ADS   He Y., Sun W., 2007. Stability and convergence of the Crank–Nicolson/Adams–Bashforth scheme for the time-dependent Navier–Stokes equations, SIAM J. Numer. Anal. , 45( 2), 837– 869. https://doi.org/10.1137/050639910 Google Scholar CrossRef Search ADS   Holme R., 2015. Large-scale flow in the core, in Treatise on Geophysics , chap. 8.04, pp. 91– 113, ed. Schubert G., Elsevier. Holme R., de Viron O., 2013. Characterization and implications of intradecadal variations in length of day, Nature , 499, 202– 205. https://doi.org/10.1038/nature12282 Google Scholar CrossRef Search ADS PubMed  Hough S.S., 1897. On the application of harmonic analysis to the dynamical theory of the tides. Part I. On Laplace’s ‘Oscillations of the First Species,’ and on the dynamics of ocean currents, Phil. Trans. R. Soc. A , 189, 201– 257. https://doi.org/10.1098/rsta.1897.0009 Google Scholar CrossRef Search ADS   Hunter J., 2007. Matplotlib: A 2D graphics environment, Comput. Sci. Eng. , 9, 90– 95. https://doi.org/10.1109/MCSE.2007.55 Google Scholar CrossRef Search ADS   Jackson A., Bloxham J., Gubbins D., 1993. Time-dependent flow at the core surface and conservation of angular momentum in the coupled core-mantle system, in Dynamics of the Earth’s Deep Interior and Earth Rotation , Vol. 72, pp. 97– 107, eds Le Mouël J.-L., Smylie D.E., Herring T., AGU Geophysical Monograph. Google Scholar CrossRef Search ADS   Jault D., 2008. Axial invariance of rapidly varying diffusionless motions in the Earth’s core interior, Phys. Earth planet. Inter. , 166, 67– 76. https://doi.org/10.1016/j.pepi.2007.11.001 Google Scholar CrossRef Search ADS   Jault D., Gire C., Le Mouël J.-L., 1988. Westward drift, core motions and exchanges of angular momentum between core and mantle, Nature , 333, 353– 356. https://doi.org/10.1038/333353a0 Google Scholar CrossRef Search ADS   Jaupart E., Buffett B., 2017. Generation of MAC waves by convection in Earth’s core, Geophys. J. Int. , 209, 1326– 1336. https://doi.org/10.1093/gji/ggx088 Google Scholar CrossRef Search ADS   Labbé F., Jault D., Gillet N., 2015. On magnetostrophic inertia-less waves in quasi-geostrophic models of planetary cores, Geophys. Astrophys. Fluid Dyn. , 109, 587– 610. https://doi.org/10.1080/03091929.2015.1094569 Google Scholar CrossRef Search ADS   Livermore P.W., Ierley G., Jackson A., 2008. The structure of Taylor’s constraint in three dimensions, Proc. R. Soc. A , 464, 3149– 3174. https://doi.org/10.1098/rspa.2008.0091 Google Scholar CrossRef Search ADS   Matsui H., King E., Buffett B., 2014. Multiscale convection in a geodynamo simulation with uniform heat flux along the outer boundary, Geochem. Geophys. Geosyst. , 15, 3212– 3225. https://doi.org/10.1002/2014GC005432 Google Scholar CrossRef Search ADS   Pais M.A., Jault D., 2008. Quasi-geostrophic flows responsible for the secular variation of the Earth’s magnetic field, Geophys. J. Int. , 173, 421– 443. https://doi.org/10.1111/j.1365-246X.2008.03741.x Google Scholar CrossRef Search ADS   Proudman J., 1916. On the motion of solids in a liquid possessing vorticity, Proc. R. Soc. A , 92, 408– 424. https://doi.org/10.1098/rspa.1916.0026 Google Scholar CrossRef Search ADS   Roberts P., 2015. Theory of the geodynamo, in Treatise on Geophysics , chap. 8.03, pp. 57– 90, ed. Schubert G., Elsevier. Roberts P., King E., 2013. On the genesis of the Earth’s magnetism, Rep. Prog. Phys. , 76, 096801, doi:10.1088/0034-4885/76/9/096801. https://doi.org/10.1088/0034-4885/76/9/096801 Google Scholar CrossRef Search ADS PubMed  Roberts P.H., Yu Z., Russell C., 2007. On the 60-year signal from the core, Geophys. Astrophys. Fluid Dyn. , 101, 11– 35. https://doi.org/10.1080/03091920601083820 Google Scholar CrossRef Search ADS   Schaeffer N., Cardin P., 2005. Quasigeostrophic model of the instabilities of the Stewartson layer in flat and depth-varying containers, Phys. Fluids , 17, 104111, doi:10.1063/1.2073547. https://doi.org/10.1063/1.2073547 Google Scholar CrossRef Search ADS   Schaeffer N., Cardin P., 2006. Quasi-geostrophic kinematic dynamos at low magnetic Prandtl number, Earth planet. Sci. Lett. , 245, 595– 604. https://doi.org/10.1016/j.epsl.2006.03.024 Google Scholar CrossRef Search ADS   Schaeffer N., Lora Silva E., Pais M.A., 2016. Can core flows inferred from geomagnetic field models explain the Earth’s dynamo?, Geophys. J. Int. , 204, 568– 877. https://doi.org/10.1093/gji/ggv488 Google Scholar CrossRef Search ADS   Schaeffer N., Jault D., Nataf H.-C., Fournier A., 2017. Turbulent geodynamo simulations: a leap towards Earth’s core, Geophys. J. Int. , 211, 1– 29. https://doi.org/10.1093/gji/ggx265 Google Scholar CrossRef Search ADS   Stellmach S., Lischper M., Julien K., Vasil G., Cheng J.S., Ribeiro A., King E.M., Aurnou J.M., 2014. Approaching the asymptotic regime of rapidly rotating convection: Boundary layers versus interior dynamics, Phys. Rev. Lett. , 113, 254501, doi:10.1103/PhysRevLett.113.254501. https://doi.org/10.1103/PhysRevLett.113.254501 Google Scholar CrossRef Search ADS PubMed  Takahashi F., Matsushima M., Honkura Y., 2008. Scale variability in convection-driven MHD dynamos at low Ekman number, Phys. Earth planet. Inter. , 167, 168– 178. https://doi.org/10.1016/j.pepi.2008.03.005 Google Scholar CrossRef Search ADS   Taylor G.I., 1917. Motion of solids in fluids when the flow is not irrotational, Proc. R. Soc. A , 93, 99– 113. https://doi.org/10.1098/rspa.1917.0007 Google Scholar CrossRef Search ADS   Taylor J.B., 1963. The magneto-hydrodynamics of a rotating fluid and the earth’s dynamo problem, Proc. R. Soc. A , 274, 274– 283. https://doi.org/10.1098/rspa.1963.0130 Google Scholar CrossRef Search ADS   Teed R., Jones C.A., Tobias S., 2014. The dynamics and excitation of torsional waves in geodynamo simulations, Geophys. J. Int. , 196, 724– 735. https://doi.org/10.1093/gji/ggt432 Google Scholar CrossRef Search ADS   Teed R.J., Jones C.A., Tobias S.M., 2015. The transition to Earth-like torsional oscillations in magnetoconvection simulations, Earth planet. Sci. Lett. , 419, 22– 31. https://doi.org/10.1016/j.epsl.2015.02.045 Google Scholar CrossRef Search ADS   Wu C.-C., Roberts P.H., 2015. On magnetostrophic mean-field solutions of the geodynamo equations, Geophys. Astrophys. Fluid Dyn. , 109, 84– 110. Yadav R.K., Gastine T., Christensen U.R., Wolk S.J., Poppenhaeger K., 2016. Approaching a realistic force balance in geodynamo simulations, Proc. Natl. Acad. Sci. USA , 113, 12 065– 12 070. https://doi.org/10.1073/pnas.1608998113 Google Scholar CrossRef Search ADS   APPENDIX: CONSERVATION OF AXIAL ANGULAR MOMENTUM The total axial angular momentum $$\mathcal {L}_{z}$$ of the outer core is contained in the concentric cylinders rotating with axisymmetric angular velocity $$\frac{\overline{u_{\phi }}}{s}$$. The time variation of $$\mathcal {L}_{z}$$ is then   $$\frac{1}{4\pi } \frac{\partial }{\partial t}\mathcal {L}_{z} = \int _{s_{1}}^{s_{2}} \! s^{3} L \frac{\partial }{\partial t}\left(\frac{\overline{u_{\phi }}}{s}\right) \, \mathrm{d} s \, .$$ (A1) As shown by eq. (8), $$\frac{\partial }{\partial t}(\frac{\overline{u_{\phi }}}{s})$$ can be decomposed into four component torques. Defining the linear operator $$\mathcal {G}(f) = \int _{s_{1}}^{s_{2}} \! s^{3} L f \, \mathrm{d} s$$, eq. (A1) may be rewritten as a sum:   \begin{eqnarray} \frac{1}{4\pi } \frac{\partial }{\partial t}\mathcal {L}_{z} &= \mathcal {G} \left(\Gamma _{\textrm{V}}\right)+ \mathcal {G} \left(\Gamma _{\mathrm{L}_{1}} \right)+ \mathcal {G} \left(\Gamma _{\mathrm{L}_{2}} \right)+ \mathcal {G} \left(\Gamma _{\textrm{R}}\right)\, . \end{eqnarray} (A2)With the input torques   \begin{eqnarray} \Gamma _{\textrm{V}}= \frac{E}{s^{3} L} \frac{\partial }{\partial s}\left(s^{3} L \frac{\partial }{\partial s}\left(\frac{\overline{u_{\phi }}}{s}\right)\right)\, , \end{eqnarray} (A3)  \begin{eqnarray} \Gamma _{\textrm{L}_{1}}= \frac{1}{s^{3} L} \frac{\partial }{\partial s}\left(s^{3} L B_{0s}\frac{\overline{b_{\phi }}}{s}\right)\, , \end{eqnarray} (A4)  \begin{eqnarray} \Gamma _{\textrm{L}_{2}}= \frac{1}{s} \left(\overline{\frac{b_{s}}{s} \frac{\partial }{\partial s}\left(s b_{\phi }\right)} \right)\, , \end{eqnarray} (A5)  \begin{eqnarray} \Gamma _{\textrm{R}}= - \frac{1}{s} \left(\overline{\frac{u_{s}}{s} \frac{\partial }{\partial s}\left(s u_{\phi }\right)} \right)\, , \end{eqnarray} (A6) the first two terms on the right-hand side of eq. (A2) reduce to $$\mathcal {G}(\Gamma _{\textrm{V}}) = s^{3} L \frac{\partial }{\partial s}(\frac{\overline{u_{\phi }}}{s})|_{s_{1}}^{s_{2}}$$ and $$\mathcal {G}(\Gamma _{\rm{L}_{1}}) = s^{3} L B_{0s}\frac{\overline{b_{\phi }}}{s}|_{s_{1}}^{s_{2}}$$. Because we assume a stress-free condition on $$\overline{u_{\phi }}$$ and $$\overline{b_{\phi }} = 0$$ at s1 and s2, both terms are zero, reducing eq. (A2) to   \begin{eqnarray} \frac{1}{4\pi } \frac{\partial }{\partial t}\mathcal {L}_{z} = \mathcal {G} \left(\Gamma _{\mathrm{L}_{2}} \right)+ \mathcal {G} \left(\Gamma _{\textrm{R}}\right)\, . \end{eqnarray} (A7) Since both u and b are solenoidal and periodic in azimuth, eqs (A5) and (A6) may be rewritten as   \begin{eqnarray} \Gamma _{\textrm{L}_{2}}= \phantom{-} \frac{1}{s} \left(\frac{1}{s^{2} L} \frac{\partial }{\partial s}\left(s^{2} L \, \overline{b_{s}b_{\phi }} \right)- \overline{\left(\beta b_{s}- \frac{\partial }{\partial z}b_{z}\right)b_{\phi }} \right)\, , \end{eqnarray} (A8)  \begin{eqnarray} \Gamma _{\textrm{R}}= - \frac{1}{s} \left(\frac{1}{s^{2} L} \frac{\partial }{\partial s}\left(s^{2} L \, \overline{u_{s}u_{\phi }} \right)- \overline{\left(\beta u_{s}- \frac{\partial }{\partial z}u_{z}\right)u_{\phi }} \right)\, . \end{eqnarray} (A9) Application of the $$\mathcal {G}$$ operator to the first term on the right-hand side of each equation again produces terms which depend only on boundary values: $$s^{2} L \left.\overline{b_{s}b_{\phi }}\right|_{s_{1}}^{s_{2}}$$ and $$s^{2} L \left.\overline{u_{s}u_{\phi }}\right|_{s_{1}}^{s_{2}}$$. Because we use no-penetration (ψ = 0 →  us = 0) and a = 0 (→  bs = 0) boundary conditions, both are zero. This leaves   \begin{eqnarray} \frac{1}{4 \pi } \frac{\partial }{\partial t}\mathcal {L}_{z} \!=\! \mathcal {G} \!\left(\overline{\left(\beta u_{s}\!-\! \frac{\partial }{\partial z}u_{z}\right)u_{\phi }} \right)\!-\! \mathcal {G} \left(\overline{\left(\beta b_{s}\!-\! \frac{\partial }{\partial z}b_{z}\right)b_{\phi }} \right).\nonumber\\ \end{eqnarray} (A10) An axial profile must be assumed for both the u and b fields. The no-penetration boundary condition at the spherical top and bottom boundaries dictates that, defining the boundary’s normal vector as $${\bf \hat{n}}$$,   $$\mathbf {u}\cdot {\bf \hat{n}} = 0 \quad \Rightarrow \quad u_{z}(L) = - \frac{s}{L} u_{s}\, .$$ (A11)Assuming that the uz velocity profile varies linearly with z, and is zero in the equatorial plane,   $$\frac{\partial }{\partial z}u_{z}= \beta u_{s}\, ,$$ (A12) which is consistent with mass conservation ($$\boldsymbol{\nabla }\cdot \mathbf {u}= 0$$), and with the definition of the flow used in eqs (3) and (4). This causes the first term on the right-hand side of eq. (A10) to vanish. Our definition of the magnetic field perturbation b in eqs (3) and (5) follows the same form as that of u:   $$\frac{\partial }{\partial z}b_{z}= \beta b_{s}\, .$$ (A13)Substituting eqs (A12) and (A13) into eq. (A10) causes the remaining terms on the right-hand side of eq. (A2) to vanish. Thus, the angular momentum is conserved:   $$\frac{\partial }{\partial t}\mathcal {L}_{z} = 0 \, .$$ (A14) © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

# Convectively driven decadal zonal accelerations in Earth’s fluid core

, Volume 213 (1) – Apr 1, 2018
13 pages

Publisher
Oxford University Press
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggx548
Publisher site
See Article on Publisher Site

### Abstract

Summary Azimuthal accelerations of cylindrical surfaces co-axial with the rotation axis have been inferred to exist in Earth’s fluid core on the basis of magnetic field observations and changes in the length-of-day. These accelerations have a typical timescale of decades. However, the physical mechanism causing the accelerations is not well understood. Scaling arguments suggest that the leading order torque averaged over cylindrical surfaces should arise from the Lorentz force. Decadal fluctuations in the magnetic field inside the core, driven by convective flows, could then force decadal changes in the Lorentz torque and generate zonal accelerations. We test this hypothesis by constructing a quasi-geostrophic model of magnetoconvection, with thermally driven flows perturbing a steady, imposed background magnetic field. We show that when the Alfvén number in our model is similar to that in Earth’s fluid core, temporal fluctuations in the torque balance are dominated by the Lorentz torque, with the latter generating mean zonal accelerations. Our model reproduces both fast, free Alfvén waves and slow, forced accelerations, with ratios of relative strength and relative timescale similar to those inferred for the Earth’s core. The temporal changes in the magnetic field which drive the time-varying Lorentz torque are produced by the underlying convective flows, shearing and advecting the magnetic field on a timescale associated with convective eddies. Our results support the hypothesis that temporal changes in the magnetic field deep inside Earth’s fluid core drive the observed decadal zonal accelerations of cylindrical surfaces through the Lorentz torque. Core, Numerical modelling, Planetary interiors 1 INTRODUCTION Decadal variations in Earth’s length-of-day (LOD) are believed to be driven by core flow variations (Gross 2015). For instance, an increase in the bulk angular velocity of the core entrains an increase in its axial angular momentum. Angular momentum conservation between the core and the mantle then implies that the latter must slow its rotation rate, thereby leading to a increase in LOD. On timescales of centuries or less, the dynamics of the outer core are expected to be dominated by a geostrophic balance between pressure gradients and the Coriolis force (e.g. Finlay et al.2010). Flows that obey such a balance have the property of being invariant, or ‘rigid’, parallel to the axis of rotation (Hough 1897; Proudman 1916; Taylor 1917). This is the Taylor–Proudman theorem. Temporal variations in the outer core’s axial angular momentum must then be carried by azimuthal accelerations in the form of rigid, co-axial cylindrical surfaces, or ‘geostrophic cylinders’, as shown in Fig.  1. Figure 1. View largeDownload slide The geometry of the zonal flows in Earth’s outer core which carry angular momentum. The depicted time varying zonal flows consist in either free Alfvén waves or forced accelerations. Figure 1. View largeDownload slide The geometry of the zonal flows in Earth’s outer core which carry angular momentum. The depicted time varying zonal flows consist in either free Alfvén waves or forced accelerations. Flows at the core–mantle boundary (CMB) may be reconstructed from the magnetic field’s secular variation observed at Earth’s surface (Holme 2015). Time-dependent accelerations of geostrophic cylinders may then be extracted from the reconstructed core flows, allowing a prediction of LOD changes to be built. A number of studies (e.g. Jault et al.1988; Jackson et al.1993) have shown that such predictions match well with the observed decadal LOD variations. Furthermore, Gillet et al. (2010) have shown that zonal accelerations of geostrophic cylinders carry changes in angular momentum which can also explain an observed 6-yr periodic LOD signal (Holme & de Viron 2013; Chao et al.2014). These studies confirm not only that the decadal and 6-yr LOD variations are due to core-mantle angular momentum exchanges, they also confirm the presence of rigid zonal accelerations in the core. However, the dynamics responsible for such zonal accelerations are not fully understood. This is the topic of our study. The forces responsible for controlling the dynamics of zonal flows can be examined by integrating the azimuthal component of the momentum equation over the surface of geostrophic cylinders. Upon doing so, the pressure, Coriolis, and buoyancy terms vanish. If one neglects inertia and viscosity, this implies that the azimuthal component of the Lorentz force must cancel out when integrated over geostrophic cylinders. In other words, the axial Lorentz torque on cylinders must vanish. This is known as Taylor’s condition (Taylor 1963), with systems obeying it said to be in a Taylor state. Reinstating inertia, perturbations in the magnetic field – and, therefore, in the Lorentz torque – can be accommodated by rigid zonal accelerations of geostrophic cylinders. When a magnetic field in a conducting fluid is distorted by the fluid’s motion, such as the differential rotation between coaxial cylinders, a current is induced. This current interacts with the original magnetic field to create a restoring force which opposes the motion that caused the initial magnetic distortion. In the language of Taylor’s condition, this restoring force nudges geostrophic cylinders back towards a Taylor state, subject to magnetic diffusion. This mechanism allows rigid zonal flows to oscillate about a Taylor state and, in doing so, support Alfvén waves. Braginsky (1970) suggested that the free modes of Alfvén waves could be responsible for decadal zonal accelerations in the outer core. In order to produce free Alfvén waves with a fundamental mode period of approximately 60 yr, corresponding to the characteristic timescale of the LOD signal (Roberts et al.2007; Gross 2015), the magnetic field strength throughout the core must be approximately 0.3 mT, similar to its observed strength at the CMB. However, modern geodynamo simulations suggest that the internal radial magnetic field is approximately ten times larger than that at the CMB (e.g. Christensen & Aubert 2006). Such a field strength should yield a fundamental Alfvén period of approximately 6 yr. It is therefore now believed that the 6-yr LOD signal is due to the propagation of free Alfvén waves (Gillet et al.2010), leaving the dynamics responsible for the decadal zonal accelerations unexplained. One possible explanation is that convective eddies in the fluid core continuously distort the magnetic field, causing spatial and temporal variations in the Lorentz force. Integrated over geostrophic cylinders, the induced Lorentz torques must be balanced by rigid zonal accelerations of the cylinders. Continual distortion of the internal magnetic field would then lead to continual zonal accelerations of the cylinders. The observed decadal zonal accelerations could therefore be a forced fluctuation about the Taylor state driven by convective flows. If this is the case, forcing must occur on timescales longer than the propagation time of free Alfvén waves. In other words, the typical convective velocity uC must be smaller than the Alfvén wave velocity uA. The relationship between uC and uA is captured by the Alfvén number   $$\mathcal {A}= \frac{u_{C}}{u_{A}},$$ (1)where   $$u_{A} = \frac{\left| \bf {B} \right|}{\sqrt{\rho \mu _{0}}},$$ (2) B is the magnetic field, ρ is the fluid density, and μ0 is the magnetic permeability of free space. Typical large-scale flow velocities in Earth’s core are of the order of uC = 10 km yr−1, or 3 × 10−4 m s−1 (Holme 2015). With a typical radial magnetic field of 3 mT (e.g. Gillet et al.2010; Buffett 2010), an outer core density of 104 kg m−3, and μ0 = 4π × 10−7 N A−2, the Alfvén wave velocity is approximately uA = 3 × 10−2 m s−1. This yields $$\mathcal {A}\approx 0.01$$. Such a small value of $$\mathcal {A}$$ supports the idea that the observed decadal zonal flows in Earth’s core may be driven by convection via time-dependent Lorentz torques. The goal of this study is to demonstrate that such a dynamical scenario is possible. One option to investigate the zonal flow dynamics in the Earth’s core is to use a 3-D numerical model of a self-generated dynamo (e.g. Christensen & Wicht 2015). However, typical Alfvén numbers in such models are $$\mathcal {A}\approx 1$$, meaning the dynamics of free Alfvén waves and convectively driven zonal accelerations are not well separated. Some recent 3-D models have been able to achieve lower values of $$\mathcal {A}\approx 0.1$$ (e.g. Aubert et al.2017; Schaeffer et al.2017), though these models are numerically expensive to run. Here, we follow a different strategy and exploit the fact that, as discussed previously, the large-scale fluid motions in the core which vary on decadal timescales are expected to be almost invariant along the rotation axis (Jault 2008). Flows of this type are often termed Quasi-Geostrophic (QG), with their existence in Earth’s core supported by the observed geomagnetic secular variation (Pais & Jault 2008; Gillet et al.2011). They also emerge from numerical models of the geodynamo (e.g. Schaeffer et al.2017). We present in this study a numerical model of the decadal timescale dynamics of Earth’s fluid core constructed within a QG framework (e.g. Aubert et al.2003). We add an induction equation to a QG model of thermal convection to follow the evolution of the magnetic field as it is sheared and advected by the flow. We likewise take into account the feedback of the Lorentz force on the flow, though in a limited way—see Section 2.7. Since we focus on the short timescale dynamics, we substitute the self-generated magnetic field resulting from dynamo action with a steady, imposed background magnetic field. We thus perform a magnetoconvection experiment, in which the perturbations of the magnetic field are tracked with respect to this imposed field. This strategy allows us to readily achieve $$\mathcal {A}< 1$$ at very modest numerical cost. 2 MODEL 2.1 Background Previous studies (e.g. Aubert et al.2003; Gillet & Jones 2006) have shown that thermally driven, Boussinesq, QG models can reproduce flow patterns similar in scale and behaviour to the ones we expect in planetary cores. The QG approximation takes advantage of the geometric constraint imposed by strong rotation: fluid motions must be dominantly rigid. Studying the dynamics of the full system is then equivalent to studying the dynamics of a slice through the system’s equatorial plane. A QG model therefore collapses 3-D convection experiments onto a 2-D domain. In our case, this domain corresponds to the shaded annulus of Fig. 2. Using the usual cylindrical coordinates (s, ϕ, z), with the ez direction aligned with the rotation axis, this annulus extends from the ‘tangent cylinder’ surrounding the equivalent of the inner core at s = s1, to the ‘equator’ at s = s2. We do not model the region inside the tangent cylinder, equivalent to the polar regions above and below the inner core. Figure 2. View largeDownload slide Geometry of our QG model. The domain of integration is the shaded 2-D annulus between s1 and s2. Figure 2. View largeDownload slide Geometry of our QG model. The domain of integration is the shaded 2-D annulus between s1 and s2. Our model is based on the QG model of thermal convection, first presented by Busse & Or (1986) and Cardin & Olson (1994) and expanded in many subsequent studies (e.g. Aubert et al.2003; Gillet & Jones 2006). In this model, 2-D equations for the momentum and thermal evolution are coupled together. To these governing equations, we add a third, two-dimension magnetic induction equation, while augmenting the momentum equation with a Lorentz force term. Although the Taylor–Proudman theorem predicts rigid motions, there are no such theoretical constraints on the temperature nor the magnetic field. However, by averaging the governing equations in the axial (ez) direction, we capture, if imperfectly, the net effect of both the temperature and magnetic field on the flow dynamics. Strategies have been devised to couple QG flows to a 3-D magnetic field (Schaeffer & Cardin 2006; Schaeffer et al.2016) and temperature (Guervilly & Cardin 2016), but here we restrict our attention to a purely 2-D model. We scale length by the radius of the outer sphere r2, time by the inverse of the angular rotational velocity Ω, temperature by the superadiabatic temperature difference ΔT between the inner and outer spheres, and the magnetic field by $$r_{2} \Omega \sqrt{\rho _{0} \mu _{0}}$$, where ρ0 is the reference density. The QG approximation remains valid so long as the Coriolis force remains dominant. In our case, this implies that both the Rossby number Ro (the ratio of inertial to Coriolis forces) and the Lehnert number λ (the ratio of magnetic to Coriolis forces) must remain ≪1. Our choice of scalings causes Ro in our model to be equivalent to |u| (the typical amplitude of the non-dimensional velocity) and λ to be equivalent to |b| (the typical amplitude of the non-dimensional magnetic field perturbation). We must therefore ensure that |u| ≪ 1 and |b| ≪ 1. 2.2 Flow and magnetic fields Because both the flow and magnetic fields are solenoidal, they may be represented in terms of vector potentials. We define the velocity u and magnetic field perturbation b as   \begin{eqnarray} \mathbf {u} = \overline{u_{\phi }}\, \mathbf {e}_{\phi }+ \frac{1}{L} \boldsymbol{\nabla }\times \left(L \psi \, \mathbf {e}_{z}\right)\, , \end{eqnarray} (3a)  \begin{eqnarray} \mathbf {b} = \overline{b_{\phi }}\, \mathbf {e}_{\phi }+ \frac{1}{L} \boldsymbol{\nabla }\times \left(L a \, \mathbf {e}_{z}\right)\, , \end{eqnarray} (3b)where ψ and a are toroidal scalars and $$L = \sqrt{1 - s^{2}}$$ is the half-column height. With an overbar denoting an azimuthal average, $$\overline{u_{\phi }}$$ and $$\overline{b_{\phi }}$$ capture the axisymmetric azimuthal (zonal) flow and zonal magnetic fields, respectively. The horizontal components of the velocity field $${\bf u}_{H} = u_{s}\mathbf {e}_{s}+ u_{\phi }\mathbf {e}_{\phi }$$ and the axial vorticity ωz are then defined as   \begin{eqnarray} u_{s} = \frac{1}{s} \frac{\partial \psi }{\partial \phi } \, , \end{eqnarray} (4a)  \begin{eqnarray} u_{\phi } = \overline{u_{\phi }}- \left(\frac{\partial }{\partial s}+ \beta \right)\psi \, , \end{eqnarray} (4b)  \begin{eqnarray} \omega _{z} = \left(2 \frac{\overline{u_{\phi }}}{s}+ \frac{\partial }{\partial s}\overline{u_{\phi }}\right)- \nabla _{H}^{2} \psi - \frac{1}{s} \frac{\partial }{\partial s}\left(s \beta \psi \right)\, . \end{eqnarray} (4c) Similarly, the horizontal components of the magnetic perturbation field $${\bf b}_{H} = b_{s}\mathbf {e}_{s}+ b_{\phi }\mathbf {e}_{\phi }$$ and the axial current jz are defined as   \begin{eqnarray} b_{s}= \frac{1}{s} \frac{\partial a}{\partial \phi } \, , \end{eqnarray} (5a)  \begin{eqnarray} b_{\phi } = \overline{b_{\phi }}- \left(\frac{\partial }{\partial s}+ \beta \right)a \, , \end{eqnarray} (5b)  \begin{eqnarray} j_{z}= \left(2 \frac{\overline{b_{\phi }}}{s} + \frac{\partial }{\partial s}\overline{b_{\phi }}\right)- \nabla _{H}^{2} a - \frac{1}{s} \frac{\partial }{\partial s}\left(s \beta a \right)\, . \end{eqnarray} (5c) Here, $$\nabla _{H}^{2}$$ specifies the cylindrical (s, ϕ) components of the Laplacian operator. The factor β, which enters the equations via the no-penetration condition on the r = r2 surface, is a measure of how L changes with s:   $$\beta = \frac{1}{L} \frac{\partial L}{\partial s} = - \frac{s}{L^{2}} \, .$$ (6) The definitions of us, uϕ and ωz according to eqs (4) follow the traditional approach in QG models (e.g. Schaeffer & Cardin 2005). Because the fluid is assumed incompressible, there is no axisymmetric radial velocity. Using an equivalent representation for bs, bϕ and jz in eqs (5) follows the strategy employed by Labbé et al. (2015). These rigid variables represent the axial averages of the truly 3-D magnetic field. Since their form is equivalent to the flow field, the implied assumption is that the magnetic field obeys the equivalent of a no-penetration condition on the spherical top and bottom boundary of a fluid column. While this cannot be rigorously justified, to first order this is approximately correct since the magnetic field near the core-mantle boundary is likely much smaller than deeper in the core. Furthermore, as we show in Appendix, representing the magnetic field as in eqs (5) is necessary to ensure conservation of angular momentum. 2.3 Momentum equation Under the QG approximation, the axially averaged, thermally driven Navier–Stokes equation reduces to an axial vorticity equation (e.g. Aubert et al.2003),   \begin{eqnarray} \frac{\partial }{\partial t}\omega _{z}&=& - Ra^{*}\frac{\partial \Theta }{\partial \phi } + E \nabla _{H}^{2} \omega _{z}+ \left(2 + \omega _{z}\right)\beta u_{s} \nonumber\\ &&- \left(u_{s}\frac{\partial }{\partial s}+ \frac{u_{\phi }}{s} \frac{\partial }{\partial \phi }\right)\omega _{z}+ F_{L} \, . \end{eqnarray} (7) Here, $$Ra^{*}= E^{2} RaP_{r}^{-1}$$ is the modified Rayleigh number, $$E = \nu (\Omega r_{2}^{2} )^{-1}$$ is the Ekman number, $$Ra= \alpha g_{0} \Delta T r_{2}^{3} \left(\nu \kappa \right)^{-1}$$ is the Rayleigh number, Pr = νκ−1 is the Prandtl number, and ν, α, g0, and κ are respectively the kinematic viscosity, thermal expansion coefficient, gravitational acceleration at r = r2, and thermal diffusivity. Θ is the local temperature perturbation from the conducting profile (see Section 2.5), and FL represents the axial component of the curl of the Lorentz force. Since we wish to study the dynamics of zonal flows, we write a separate equation for the evolution of the axisymmetric angular velocity $$\frac{\overline{u_{\phi }}}{s}$$:   $$\frac{\partial }{\partial t}\left(\frac{\overline{u_{\phi }}}{s}\right)= \Gamma _{\textrm{L}}+ \Gamma _{\textrm{R}}+ \Gamma _{\textrm{V}}\, ,$$ (8)where ΓL denotes the axial torque from Lorentz forces, ΓR the axial torque from Reynolds stresses, and ΓV the axial viscous torque. ΓR and ΓV are given by   \begin{eqnarray} \Gamma _{\textrm{R}}= - \frac{1}{s} \left(\overline{\frac{u_{s}}{s} \frac{\partial }{\partial s}s u_{\phi }} \right), \end{eqnarray} (9)  \begin{eqnarray} \Gamma _{\textrm{V}}=\phantom{-} \frac{E}{s^{3} L} \frac{\partial }{\partial s}\left(s^{3} L \frac{\partial }{\partial s}\left(\frac{\overline{u_{\phi }}}{s}\right)\right), \end{eqnarray} (10)while ΓL is given in the next section by eq. (14). Although eq. (8) is contained in eq. (7), in practice we use eq. (7) to evolve the non-axisymmetric part of ωz, and eq. (8) to evolve $$\overline{u_{\phi }}$$. We note that we have neglected in both eqs (7) and (8) viscous friction at the top and bottom spherical boundaries. These can be incorporated in the QG framework through Ekman pumping terms (e.g. Schaeffer & Cardin 2005). However, for the Ekman numbers that we can reach numerically, these friction terms would play an unrealistically dominant role in the force balance, unlike in the Earth’s core. 2.4 Induction equation and Lorentz force Since a 2-D model cannot maintain a self-sustaining dynamo (e.g. Cowling 1957; Roberts 2015), we instead impose a steady background magnetic field $${\bf B_{0}}$$ on it. For simplicity we choose a uniform field, $${\bf B_{0}}(s,\phi ,z) = B_{0s}\, \mathbf {e}_{s}$$. This uniform field can be interpreted to represent the local cumulative effect from all length scales of the background field. Two-dimensional perturbations (bs, bϕ) from that background state are then tracked via the induction equation. We form an equation for the magnetic potential a by axially averaging the s-component of the induction equation, which yields   $$\frac{\partial }{\partial t}a = - u_{\phi }B_{0s}+ \left(u_{s}b_{\phi }- u_{\phi }b_{s}\right)+ \frac{E}{P_{m}} \left(\nabla _{H}^{2} a + \frac{2 \beta a}{s} \right)\, ,$$ (11)where Pm = νη−1 is the magnetic Prandtl number and η is the magnetic diffusivity. The β factor in the last term of eq. (11) originates from the contribution of bϕ to the vector Laplacian in the s-direction. Physically, it represents an added contribution to dissipation introduced by the spherical geometry. The equation for $$\overline{b_\phi }$$ is derived from the axially averaged and azimuthally averaged ϕ-component of the induction equation,   \begin{eqnarray} \frac{\partial }{\partial t}\left(\frac{\overline{b_{\phi }}}{s}\right)&=& B_{0s}\frac{\partial }{\partial s}\left(\frac{\overline{u_{\phi }}}{s}\right)+ \frac{1}{s}\frac{1}{L} \frac{\partial }{\partial s}\Big ( L \left(\overline{u_{\phi }b_{s}} - \overline{u_{s}b_{\phi }} \right)\Big ) \nonumber\\ &&+\,\frac{1}{s^{3}} \frac{E}{P_{m}} \frac{\partial }{\partial s}\left(s^{3} \frac{\partial }{\partial s}\right)\left(\frac{\overline{b_{\phi }}}{s}\right)\, . \end{eqnarray} (12)The axially averaged Lorentz force which enters eq. (7) is   $$F_{L} = \left(\left(B_{0s}+ b_{s}\right)\frac{\partial }{\partial s}+ \frac{b_{\phi }}{s} \frac{\partial }{\partial \phi }\right)j_{z}\, .$$ (13) Although eq. (13) is the correct expression for the Lorentz force acting on flow eddies in our QG model, we note that in practice we do not compute this term, that is, we set FL = 0. The justification for this is given in Section 2.7. Meanwhile, the axial Lorentz torque in eq. (8) is   $$\Gamma _{\textrm{L}}= \Gamma _{\textrm{L}_{1}}+ \Gamma _{\textrm{L}_{2}}\, ,$$ (14)where   \begin{eqnarray} \Gamma _{\textrm{L}_{1}} = \frac{1}{s^{3}} \frac{1}{L} \frac{\partial }{\partial s}\left(s^{3} L B_{0s}\frac{\overline{b_{\phi }}}{s}\right)\, , \end{eqnarray} (15)  \begin{eqnarray} \Gamma _{\textrm{L}_{2}} = \frac{1}{s} \left(\overline{\frac{b_{s}}{s} \frac{\partial }{\partial s}s b_{\phi }} \right)\, . \end{eqnarray} (16) 2.5 Thermal equations Thermal convection is driven by the steady, superadiabatic temperature difference ΔT = T1 − T2, where T1 and T2 are the superadiabatic temperatures on the spheres r = r1 and r = r2, respectively. Taking the axial average of the resulting 3-D temperature profile T0 gives the background conducting, rigid temperature profile T as a function of cylindrical radius:   $$T = \left(\frac{r_{2} T_{2} - r_{1} T_{1}}{r_{2} - r_{1}} \right)+ \left(\frac{r_{1} r_{2}}{r_{2} - r_{1}} \right)\frac{\Delta T}{L} \ln \left(1 + L \right)\, .$$ (17)Since it is implicitly contained in the control parameter Ra*, we are free to set ΔT to an arbitrary non-zero value. We therefore choose T1 = 1, T2 = 0. These choices, along with eq. (17), are consistent with those used by Aubert et al. (2003) and Gillet & Jones (2006). Time-dependent perturbations Θ from T are tracked with the axially averaged heat equation   $$\frac{\partial \Theta }{\partial t} = -\left(\mathbf {u}_{H} \cdot \boldsymbol{\nabla }_{H} \right)\left(T + \Theta \right)+ \frac{E}{P_{r}} \nabla ^{2}_{H} \Theta \, ,$$ (18)where $$\boldsymbol{\nabla }_{H}$$ represents the (s, ϕ) components of the gradient operator. 2.6 Boundary conditions We ensure that the radial magnetic perturbation field bs drops to zero at the boundaries by imposing a = 0, and we also impose $$\frac{\overline{b_{\phi }}}{s}= 0$$. Similarly, we ensure that the temperature perturbation drops to zero on the boundary by imposing Θ = 0. To respect the no-penetration boundary condition, ψ must be constant on the inner and outer boundaries. For convenience, we use ψ = 0. We apply an additional no-slip boundary condition on the non-axisymmetric part of uϕ such that $$\frac{\partial \psi }{\partial s} = 0$$. Forcing us and uϕ to be both zero at the boundary is self-consistent with the a = 0 condition (see eq. 11) so it is the most natural choice. However, we impose a free-slip boundary condition on $$\frac{\overline{u_{\phi }}}{s}$$, such that $$\frac{\partial }{\partial s}(\frac{\overline{u_{\phi }}}{s})=0$$. The latter choice is made to minimize viscous friction effects on the zonal flows—the target of our study—which should be small in an Earth-like regime. 2.7 Parameter regime As argued in the introduction, in order to be in an Earth-like regime ($$\mathcal {A}\ll 1$$), typical convective flow speeds (uC) must be much smaller than the typical Alfvén wave speed (uA). For a given Ekman number E, convective speeds are controlled by Ra*. With our choice of non-dimensionalization, the Alfvén wave velocity in eq. (2) is   $$u_{A}= \left| {\bf B} \right| = \left| {\bf B}_{0} + \mathbf {b} \right| \, ,$$ (19)so a lower bound is set by the strength of the background magnetic field B0s. Therefore, the combination of Ra* and B0s in our model must be such that $$\mathcal {A}\ll 1$$. Furthermore, we must be in a regime where accelerations of the zonal flow are dominantly controlled by the Lorentz, rather than the Reynolds, torque. An inspection of eqs (9) and 16) shows that the Reynolds and Lorentz torques are proportional to the non-axisymmetric parts of u2 and b2, respectively, with typical scales given by their root-mean-square amplitudes $$u_{\rm {\small RMS}}= \left| \mathbf {u}- \overline{u_{\phi }} \mathbf {e}_{\phi } \right| \approx u_{C}$$ and $$b_{\rm {\small RMS}}= \left| \mathbf {b}- \overline{b_{\phi }} \mathbf {e}_{\phi } \right|$$. Hence, we must be in a regime where $$b_{\rm {\small RMS}}\gg u_{\rm {\small RMS}}$$. The parameters controlling the ratio between $$b_{\rm {\small RMS}}$$ and $$u_{\rm {\small RMS}}$$ can be established by analysing the induction equation. When our model has achieved statistical equilibrium, the source and diffusion terms must be in balance:   $${\left| \boldsymbol{\nabla }\times \left(\left({\bf B}_{0} + \mathbf {b}\right)\times \mathbf {u}\right) \right| = \left| \frac{E}{P_{m}} \boldsymbol{\nabla }_{H}^{2} \mathbf {b} \right|.}$$ (20) Using local flow length scale ℓ1 and magnetic field length scale ℓ2, eq. (20) scales as   $$\frac{\left| {\bf B} \right| \left| \mathbf {u} \right|}{\ell _{1}} = \frac{E}{P_{m}} \frac{\left| \mathbf {b} \right|}{\ell _{2}^{2}} \quad \Rightarrow \quad \frac{\left| \mathbf {b} \right|}{\left| \mathbf {u} \right|} = \frac{\left| {\bf B} \right|}{\ell _{1}} \frac{P_{m}\ell _{2}^{2}}{E} \, .$$ (21) The Alfvén timescale can be defined as τA = ℓ1/uA, and the magnetic diffusion timescale as $$\tau _{D} = P_{m}\ell _{2}^{2} / E$$. The ratio Lu = τD/τA is then the Lundquist number. Thus, by using $$\left| \mathbf {u} \right| \approx u_{\rm {\small RMS}}$$, $$\left| \mathbf {b} \right| \approx b_{\rm {\small RMS}}$$, and $$u_{A}= \left| {\bf B} \right| \approx B_{0s}$$, eq. (21) becomes   $$\frac{b_{\rm {\small RMS}}}{u_{\rm {\small RMS}}} = \frac{B_{0s}}{\ell _{1}} \frac{P_{m}\ell _{2}^{2}}{E} = \frac{\tau _{D}}{\tau _{A}} = \mathrm{Lu} \, .$$ (22) Hence, the requirement of $$b_{\rm {\small RMS}}\gg u_{\rm {\small RMS}}$$ implies that we must be in a regime where Lu ≫ 1, which in turn places constraints on the combinations of E, Pm and B0s which we may choose. For an Earth-like setup, we would ideally have Pm ≪ 1, which requires that we pick a sufficiently large ratio of B0s/E. In practice, however, numerical constraints limit both the maximum possible value of B0s and the minimum possible value of E. In addition, it is desirable to keep uA < 1, so that the Alfvén wave speed is smaller than the inertial wave speed and our model remains consistent with our QG assumption on rigid flows. Thus, for a numerically achievable ratio of B0s/E, while we may choose Pm < 1, Pm must remain sufficiently large such that solutions are in a regime characterized by Lu ≫ 1. The final aspect of the convective regime to be addressed concerns the influence of the Lorentz force FL on the vorticity equation of eq. (7). Left to evolve dynamically, flow and magnetic field lines tend to align themselves so as to limit induction by shear (e.g. Schaeffer et al.2017). However, because we impose a steady background magnetic field in our model, such a fully dynamic reorganization is not possible. As a result, for the parameters that we use (see next section), convective eddies are unduly constrained by the Lorentz force and become thin, elongated columns in the s-direction, sharing little resemblance with the large-scale flows in Earth’s core. However, when we set FL = 0 in eq. (7) and ignore the influence of the Lorentz force on the non-axisymmetric flows, we retrieve Earth-like large-scale eddies. Because our main goal is to illustrate the mechanism by which Earth-like, large-scale flows can interact with the background field to generate slow zonal accelerations, our model is a better analogue to Earth’s core with FL turned off. In other words, although convective eddies in our model are still allowed to distort the magnetic field, we do not take into account the feedback of the Lorentz force on the flow eddies. We retain the Lorentz torque ΓL in the axially symmetric torque balance of eq. (8), as it is our primary objective to show that it can be the dominant driver of rigid accelerations. 3 RESULTS 3.1 Solution scheme A semi-spectral method with 768 Fourier modes in azimuth and radial derivatives approximated by second-order finite differences between 901 points arranged on a Chebyshev grid in radius was used to discretize eqs (7), (8), (11), (12) and (18). The resulting discrete equations were then evolved in timesteps of 5 × 10−4 using a combination of a Crank–Nicolson method for the linear terms and a second-order Adams–Bashforth scheme for the nonlinear terms (e.g. He & Sun 2007). Use of a Chebyshev grid ensured fine enough spacing to resolve boundary layers (proportional in thickness to E1/2 ≈ 0.002), while maintaining reasonable computation times via coarser spacing away from the boundaries. In our setup, minimum spacing near the boundaries is 2 × 10−6, while maximum spacing in the centre of the model domain is 1 × 10−3. We set s1 = 0.35, mimicking the thickness ratio of Earth’s core, and s2 = 0.98. Limiting our domain to s2 = 0.98 instead of s2 = 1 is convenient, as it allows us to use a slightly larger grid space and timesteps (since β → ∞ as s → 1). The results that we now present come from an experiment with E = 5.0 × 10−6, Ra = 5.0 × 108, B0s = 0.15, Pm = 0.1, and Pr = 1.0. With these parameters, a numerical simulation started from small perturbations typically takes about 1000 non-dimensional time units (or about 160 full rotations) to reach statistical equilibrium. 3.2 General features of the flow and induced magnetic fields Fig. 3 shows a snapshot in time of the non-axisymmetric ωz and jz after the model’s global energy budget has reached statistical equilibrium. The presence of eddies in both the flow and the perturbed magnetic field are clearly visible. The magnetic field perturbations exhibit structures with larger wavelengths than those in the flow field, and with smoother features, a result of the more rapid magnetic to viscous diffusion (Pm = 0.1). Figure 3. View largeDownload slide Snapshots in time of (a) non-axisymmetric vorticity ωz and (b) non-axisymmetric axial current jz, as seen looking downward from the north pole. Figure 3. View largeDownload slide Snapshots in time of (a) non-axisymmetric vorticity ωz and (b) non-axisymmetric axial current jz, as seen looking downward from the north pole. Magnetic field perturbations result from the action of convective eddies shearing and advecting the sum of the background and perturbed magnetic field. Time-dependency in both the flow and the magnetic field leads to fluctuations in the Reynolds (ΓR, eq. 9) and Lorentz (ΓL, eq. 14) torques, respectively. These must be accommodated by zonal accelerations. Typical $$u_{\rm {\small RMS}}$$ and $$b_{\rm {\small RMS}}$$ values after equilibration are 0.0024 and 0.14, respectively. The amplitude of magnetic field perturbations in our numerical experiment is of the same order as the imposed background field, $$\left| \mathbf {b} \right| \approx \left| {\bf B}_{0} \right|$$. The Alfvén and Lundquist numbers of our simulation are then $$\mathcal {A}\approx 0.014$$ and Lu ≈ 70, within the region of parameter space where we expect ΓL ≫ ΓR. This is the region where the decadal timescale dynamics of zonal accelerations in Earth’s core should reside. 3.3 Time-averaged axisymmetric force balance The dashed black line of Fig. 4 shows the time-averaged zonal angular velocity. Its profile is dominated by a shear flow spanning the whole of the modelled region, retrograde at the outer boundary and prograde at the inner boundary. The amplitude of this shear flow is of the same order of magnitude as the amplitude of typical convective eddies. Figure 4. View largeDownload slide Time-averaged mean axial torques (solid lines) and time-averaged zonal angular velocity (dashed line) as a function of radius. Figure 4. View largeDownload slide Time-averaged mean axial torques (solid lines) and time-averaged zonal angular velocity (dashed line) as a function of radius. The mean (time-averaged) zonal flow of Fig. 4 differs from the one typically observed in the absence of a magnetic field. In non-magnetic convection with stress-free boundaries, the direction of the mean zonal flow is reversed – it is prograde at the equator – and its characteristic velocity is much larger than the velocities associated with convective eddies. It results from time-averaged Reynolds stresses, themselves being the product of the topographic beta effect acting on convective eddies (e.g. Cardin & Olson 1994; Christensen 2002). In the presence of a radial magnetic field, distortion of this magnetic field by the shear flow induces a restoring Lorentz force (through $$\Gamma _{\rm{L}_{1}}$$ of eq. 15) which limits the growth of the mean zonal flow. In 3-D models, this mean zonal flow is not z-invariant (e.g. Aubert 2005). In our QG model, because our magnetic field perturbation is defined with a built-in topographic beta effect identical to that of the flow, a time-averaged Maxwell stress ($$\Gamma _{\rm{L}_{2}}$$, eq. 16) is maintained in the same way as the time-averaged Reynolds stress ΓR. Since $$\Gamma _{\rm{L}_{2}}$$ has the same form, but opposite sign, as ΓR, and since the magnitude of the former typically dominates that of the latter in our model, the direction of the driven mean zonal flow is reversed to that produced in non-magnetic convection. The time-averaged $$\overline{u_\phi }$$ profile, then, is the result of a balance between the time-averaged torques. Fig. 4 illustrates this balance, showing the time-averaged Reynolds (ΓR, orange), Lorentz (ΓL, blue), and viscous (ΓV, green) torques as a function of radius. In the interior of the domain, the Reynolds and Lorentz torques largely balance one another. The viscous torque becomes more important near the boundaries, especially the outer one. The enhancement of all three torques near the outer boundary is caused by the large β-effect: both the viscous and the $$\Gamma _{\rm{L}_{1}}$$ part of the Lorentz torque depends explicitly on β (through the s-derivative of L, see eqs 10 and 15), while the Reynolds (ΓR, eq. (9)) and Maxwell ($$\Gamma _{\rm{L}_{2}}$$, eq. (16)) torques implicitly involve β through their dependence on uϕ and bϕ, respectively. We further note that both $$\Gamma _{\rm{L}_{1}}$$ and $$\Gamma _{\rm{L}_{2}}$$ are individually much larger than ΓL, as shown by Fig. 5. However, they tend to cancel one another, leaving a net Lorentz torque several orders of magnitude smaller than either one individually. This self-cancellation will be re-examined in Section 3.5. Figure 5. View largeDownload slide Time-averaged components of the mean axial Lorentz torque as a function of radius. Figure 5. View largeDownload slide Time-averaged components of the mean axial Lorentz torque as a function of radius. 3.4 Zonal accelerations Fluctuations in time with respect to the time-averaged torque balance are the main focus of our study. Fig. 6 shows the axisymmetric angular accelerations (top panel) and the time-varying parts of ΓL, ΓR, and ΓV (bottom three panels) after the system has equilibrated. Fluctuations in ΓL have a typical amplitude five times larger than those of ΓR: ∼ 6.6 × 10−5 versus ∼1.2 × 10−5. ΓV only plays a small role in the time-dependent dynamics, with RMS fluctuations of the order of ∼0.3 × 10−5. Figure 6. View largeDownload slide The angular acceleration (top panel), and the time-dependent parts of the Lorentz (second panel), Reynolds (third panel) and viscous (bottom panel) torques, as functions of cylinder radius and time. The time-averaged contribution to each torque (shown in Fig. 4) has been removed. Figure 6. View largeDownload slide The angular acceleration (top panel), and the time-dependent parts of the Lorentz (second panel), Reynolds (third panel) and viscous (bottom panel) torques, as functions of cylinder radius and time. The time-averaged contribution to each torque (shown in Fig. 4) has been removed. Thus, fluctuations in our model’s zonal accelerations are mainly controlled by the Lorentz torque. Indeed, the upper two panels of Fig. 6 suggest a strong correlation at all times and radii between $$\frac{\partial }{\partial t}(\frac{\overline{u_{\phi }}}{s})$$ and ΓL. Time on Fig. 6 has been scaled to the timescale of the fundamental Alfvén wave mode, τA = 2(s2 − s1)/uA. The joint fluctuations in the zonal accelerations and Lorentz torque cover a broad range of timescales but are dominated by periods which fall in the range of 5 to 10 times longer than τA. This slower timescale, τslow ≈ 5 − 10 τA, reflects the time-fluctuations of the magnetic field, themselves the result of induction by convective flows. These occur on a longer timescale than the Alfvén wave propagation timescale, as we expect for a regime with $$\mathcal {A}< 1$$. Scaling the temporal fluctuations on Fig. 6 to Earth’s core, taking τA ≈ 6 yr, gives a timescale τslow of 30 to 60 yr for these magnetically driven zonal accelerations, similar to the zonal accelerations inferred within Earth’s core. To further demonstrate that the slow magnetic field fluctuations originate from the underlying convective dynamics, we compute a characteristic azimuthal wave number k* (e.g. Takahashi et al.2008) at each radius s from the convolution of wavenumber k and convective speed u(k):   $$k^{*}(s) = \frac{\int _{}^{} \! k u(k) \, \mathrm{d} k}{\int _{}^{} \! u(k) \, \mathrm{d} k} \, .$$ (23) The result of this calculation as a function of radius, for the snapshot corresponding to time equal 0 on Fig. 6, is shown in panel (a) of Fig. 7. The characteristic convective length scale ℓC(s) is then calculated from the characteristic wavenumber as   $$\ell _{C}(s) = \frac{2 \pi s}{k^{*}(s)} \, .$$ (24) The characteristic convective length scale for all radial positions is shown in panel (b) of Fig. 7. Dividing these length scales by the RMS velocities at each radius yields a characteristic convective timescales τC, which is representative of the time required for the flow to create significant changes in the magnetic field:   $$\tau _{C}(s) = \frac{\ell _{C}(s)}{u_{\rm {\small RMS}}} \, .$$ (25) The dependence of τC with radius is shown in panel (c) of Fig. 7. All three plots of Fig. 7 change for different time snapshots, but are representative of the general behaviour at all times. Furthermore, they only show the dominant length and time scales of the flow, which is truly characterized by a spectrum of scales. As shown in Fig. 7(c), in the bulk of the fluid τC is approximately equal to 5 − 8 τA, within the same range as τslow. These are obviously simple order of magnitude calculations that do not capture the complex nonlinear connection between large-scale flows and time changes in the Lorentz torque. Nevertheless, the combination of Figs 6 and 7 show that it is possible to drive slow zonal accelerations by fluctuating Lorentz torques, themselves driven by underlying convective flows shearing and advecting the magnetic field. Figure 7. View largeDownload slide Characteristic (a) wave numbers, (b) length scales and (c) timescales as a function of radius at a given snapshot in time. Figure 7. View largeDownload slide Characteristic (a) wave numbers, (b) length scales and (c) timescales as a function of radius at a given snapshot in time. In addition to the slow fluctuations of Fig. 6, the convective dynamics also generate free Alfvén waves. Fig. 8 shows the second half of the zonal acceleration and Lorentz torque panels of Fig. 6 after applying a highpass filter to remove fluctuations with periods longer than 1.19 τA. This reveals the presence of periodic fluctuations with a period of approximately 1 τA. The correlation between the acceleration and Lorentz torque shows that these are indeed Alfvén waves. (Applying the same highpass filter to ΓR and ΓV produces only low-amplitude noise.) Figure 8. View largeDownload slide Output of a highpass filter applied to the second half of the zonal angular accelerations (top panel) and Lorentz torques (bottom panel) of Fig. 6. In our non-dimensional time units, the eighth-order digital Butterworth filter’s −3 dB frequency is 1/10. Figure 8. View largeDownload slide Output of a highpass filter applied to the second half of the zonal angular accelerations (top panel) and Lorentz torques (bottom panel) of Fig. 6. In our non-dimensional time units, the eighth-order digital Butterworth filter’s −3 dB frequency is 1/10. Left on their own, free Alfvén oscillations should decay away because of ohmic dissipation. In our model, they are continuously re-excited by the underlying convective dynamics, though resonant amplification remains modest and their amplitude does not rise much higher than that resulting from the forced background accelerations. Their typical RMS velocities of 1 × 10−4 are approximately 7 times smaller than the RMS velocities of 7 × 10−4 found for the slower zonal flow fluctuations. This is the reason why Alfvén waves, though present, are not apparent on Fig. 6. 3.5 Taylorization For a system in a perfect Taylor state, even though the Lorentz torque at any point may be large, cancellations occur such that the Lorentz torque integrated over a geostrophic cylinder is equal to zero. The ‘Taylorization’ is a measure of the degree of this Lorentz torque cancellation. In our model, where the magnetic field (and therefore the Lorentz torque) is assumed to be axially invariant, the Taylorization is measured by the factor $$\mathcal {T}$$:   $$\mathcal {T}(s) = \frac{\left| \oint _{}^{} \! \mathcal {M}_{\phi }(s,\phi ) \, \mathrm{d} \phi \right|}{\oint _{}^{} \! \left| \mathcal {M}_{\phi }(s,\phi ) \right| \, \mathrm{d} \phi } = \frac{ \left| \Gamma _{\textrm{L}}(s) \right| }{ \oint _{}^{} \! \left| \mathcal {M}_{\phi }(s,\phi ) \right| \, \mathrm{d} \phi } \, ,$$ (26)where $$\mathcal {M}_{\phi }(s,\phi )$$ represents the azimuthal component of the local magnetic force. A system with low Taylorization has $$\mathcal {T} \lesssim 1$$, while a system with a high Taylorization has $$\mathcal {T} \ll 1$$. Fig. 9 shows the Taylorization factor of our model as a function of cylindrical radius and time. Typical RMS values of $$\mathcal {T}$$ are approximately 3.2 × 10−4, with peak values of 2.6 × 10−3. As shown in the previous section, temporal fluctuations of the magnetic field lead to temporal fluctuations of the Lorentz torque. Thus, Taylor’s constraint is continually broken, with the Lorentz torque fluctuations being accommodated by rigid zonal accelerations. This can be observed in Fig. 6, where the times and radii of the largest zonal accelerations often coincide with the largest values of $$\mathcal {T}$$ in Fig. 9. Figure 9. View largeDownload slide Taylorization factor $$\mathcal {T}$$ as a function of radius and time. With our choice of boundary conditions, the Taylorization factor equals 1 at the boundaries. Only s ∈ [0.36, 0.97] is shown. Figure 9. View largeDownload slide Taylorization factor $$\mathcal {T}$$ as a function of radius and time. With our choice of boundary conditions, the Taylorization factor equals 1 at the boundaries. Only s ∈ [0.36, 0.97] is shown. Therefore, our model is not in a perfect Taylor state at any given moment, but instead fluctuates about an equilibrium state characterized by a low Taylorization factor. A part of the cancellation of the Lorentz torque over any given geostrophic cylinder is inherent to the form of the Maxwell torque ($$\Gamma _{\rm{L}_{2}}$$), as it involves products of the magnetic field vectors which change direction at different azimuthal points. On a time-average, $$\Gamma _{\rm{L}_{2}}$$ does not vanish but it is mostly opposed by $$\Gamma _{\rm{L}_{1}}$$, as shown by Fig. 5. Yet, $$\Gamma _{\rm{L}_{1}}$$ and $$\Gamma _{\rm{L}_{2}}$$ do not cancel exactly: the time-averaged Taylorization factor is small ($$\mathcal {T} \sim 10^{-4}$$) but not zero. In our model, the departure from a time-averaged Taylor state is dominantly caused by the persistent torque from Reynolds stresses (and from viscous forces near the domain boundaries, see Fig. 4). 4 DISCUSSION 4.1 Slow zonal accelerations and decadal timescale dynamics in Earth’s core Our model shows that in an Earth-like parameter regime characterized by $$\mathcal {A}\ll 1$$ and Lu ≫ 1, convective flows can drive fluctuations in the magnetic field which then lead to temporal fluctuations in the Lorentz torque. The latter generate zonal accelerations in the form of free Alfvén waves, but also slower forced zonal accelerations on a timescale connected to the convective flows. In order to drive large-scale, core-size, zonal accelerations, convective eddies must be large-scale themselves. Because the typical velocity of these large-scale eddies is slower than the Alfvén wave velocity, the timescale associated with the magnetic field change they produce, and thus the Lorentz torque, is slower than τA, the timescale of Alfvén wave propagation across the core. In our model, the typical convective timescale associated with large-scale eddies is 5 to 8 times longer than τA. These eddies induce Lorentz torque fluctuations over a broad range of timescales, but dominated by periods that are 5 to 10 times longer than τA. Scaling our results to Earth’s core by taking τA = 6 yr, our model produces zonal accelerations with typical timescales in the range of 30 to 60 yr, comparable to observations. Furthermore, the power spectrum of flows reconstructed from geomagnetic secular variation shows peaks at spherical harmonic degrees 8–10 (e.g. fig. 5 of Gillet et al.2015). At mid-core radius, this corresponds to a typical length scale of approximately 600 km. Taking a typical velocity of 10 km yr−1, these are associated with a typical convective timescale of approximately 60 yr. If these flows rigidly extend axially through the core, according to our mechanism they should distort the magnetic field, and thus produce fluctuations in the Lorentz torque and zonal accelerations in response, at approximately the same 60 yr timescale. Our results suggest that the inferred decadal, large-scale zonal accelerations of geostrophic cylinders in the Earth’s core can be explained by forced fluctuations of the Lorentz torque, themselves driven by the large-scale convective flows. As observed in Earth’s core, our dynamical model contains large-scale flows, slow zonal accelerations, and faster Alfvén waves. In our model, their typical non-dimensional velocities are, respectively, 2.4 × 10−3, 7 × 10−4 and 1 × 10−4, for a relative ratio between them of 34:7:1. In Earth’s core, these typical velocities are 10 km yr−1, 2 km yr−1 and 0.2 km yr−1 (e.g. Gillet et al.2015), for a ratio of 50:10:1. Therefore, although their amplitude ratios do not match those of Earth’s core exactly, our model reproduces the correct ordering between these flows. This gives us further confidence that our model is a good analogue for the decadal flow dynamics in Earth’s core. Although our dynamical model is in the correct parameter regime in terms of $$\mathcal {A}\ll 1$$ and Lu ≫ 1, other parameters remain far from Earth-like, notably the Ekman number and magnetic Prandtl number. Care must then be taken when extrapolating our results to Earth’s core. To further confirm that our results appropriately capture the decadal timescale dynamics of zonal flows in the core, it would be desirable to carry many more numerical experiments, systematically varying some of the input parameters in order to develop scaling properties for our model. We plan to do this in a future study. A different approach is to try to do a similar analysis in a 3-D geodynamo model, some of which are approaching the parameter regime of $$\mathcal {A}\ll 1$$ and Lu ≫ 1 which we have highlighted (e.g. Aubert et al.2017; Schaeffer et al.2017). Short of doing this, we may speculate on how the use of more Earth-like Ekman and magnetic Prandtl numbers would affect the dynamics of zonal flows. For the relatively large Ekman number of our numerical experiment, viscous forces remain important in the establishment of the large-scale, quasi-core-size eddies that are visible on Fig. 3 (e.g. Aurnou et al.2015). However, at E ≈ 10−15, as inferred for Earth’s core, the E1/3 scaling law suggests that the typical length scale of viscously controlled eddies should be much smaller, about 10−5 times the core size, or of the order of 10-100 m in width. Such small eddies would be inefficient at generating the core-size changes in the magnetic field that are required to drive large-scale zonal accelerations. However, large-scale eddies are present in Earth’s core, as inferred from the geomagnetic secular variation. The typical length-scale of eddies in the core must then be controlled by processes other than this viscous scaling. Two main reasons can be invoked. First, for sufficiently low Ekman numbers, small-scale structures in the core likely feed their energy to larger length scales via an inverse energy cascade driven by geostrophic turbulence (e.g. Guervilly et al.2014; Stellmach et al.2014). Second, the influence of the Lorentz force should also lead to larger convective length scales (e.g. Roberts & King 2013), a regime which is beginning to be accessible in numerical simulations (e.g. Matsui et al.2014; Yadav et al.2016; Schaeffer et al.2017). A combination of these two effects may be responsible for forming the large structures seen in the core. Both are absent in our model. Therefore, although the dynamics that sustain the large-scale eddies in our numerical model are most likely not Earth-like, the large-scale eddies themselves share some Earth-like qualities. The key ingredient for driving decadal zonal accelerations according to our mechanism is the very presence of large-scale flow eddies with decadal convective timescale, not precisely how these eddies are generated. The magnetic Prandtl number in the core is much lower than the one we have chosen in our numerical experiment. Since a lower Pm corresponds to enhanced magnetic diffusion, temporal changes in the magnetic field would be dominantly controlled by the largest length scales of the underlying convective flow. Thus, while small-scale eddies certainly exist in the core, changes in the magnetic field should still preferentially occur at the largest length scale. We therefore expect that temporal changes in the Lorentz torque can drive decadal zonal accelerations at core-size wavelengths in the radial direction, just as we observe in our numerical experiment of limited spatial resolution. 4.2 Free Alfvén waves As shown in Fig. 8, our dynamical model excites free Alfvén waves. However, their spatio-temporal properties differ from those detected in Earth’s core. In our model, they are dominated by a standing wave oscillation of the fundamental mode. In Earth’s core, they take the form of outward travelling waves (Gillet et al.2010, 2015). The reason why Alfvén waves in Earth’s core travel outward, as well as their excitation mechanism, remain unclear. Outward travelling Alfvén waves resulting from a quasi-periodic triggering near the tangent cylinder are observed in some numerical models (Teed et al.2014, 2015; Schaeffer et al.2017). When approaching an Earth-like regime, Lorentz forces are responsible for this torque, but the precise physical mechanism has not been clearly identified. It would be a valuable effort to investigate where and how Alfvén waves are excited in our model. Since we do not see a preferential propagation direction, excitation appears to be distributed evenly within the integration domain of our model. The region close to the tangent cylinder does not appear to be the seat of any form of recurrent instability, although this may be because we have not modelled the dynamics inside the tangent cylinder. Given that convective flows in our model induce changes in the magnetic field on a broad range of timescales, Alfvén waves on Fig. 8 may simply represent the resonant response to fluctuations of the Lorentz torque which occur in the vicinity of their free period range. Indeed, correlations between Figs 8 and 9 suggest that this is the case: notable increases in the Taylorization factor are often, though not always, associated with an amplitude enhancement of free Alfvén waves. This argues along the same line of a recent study which has shown that applying a stochastic forcing in the volume of the core readily excites Alfvén waves (Gillet et al.2017). Moreover, the study of Gillet et al. (2017) has also shown that electromagnetic dissipation at the CMB transforms standing Alfvén waves into outward traveling waves, with similar characteristics as those detected in the Earth’s core. Our model constitutes a dynamical realization of such a stochastic forcing and supports the idea that Alfvén waves in Earth core are simply the response to sub-decadal changes in the Lorentz torque within the bulk of the core. Our model does not include dissipation at the CMB, but we believe (though this should be tested) that adding it would also transform our standing Alfvén waves into outwardly propagating waves. 4.3 Taylorization Taylor’s constraint is continuously being broken in our model, generating zonal accelerations in response. The numerical value of the Taylorization factor associated with these fluctuations is of the order of 10−3. The time-average state (or statistical equilibrium) about which these fluctuations occur is characterized by a higher degree of Taylorization, with a Taylorization factor of the order of 10−4. Though this is small, it indicates that the equilibrium state in our model remains far from a perfect Taylor state. In the bulk interior, this is dominantly because Reynolds stresses lead to a non-zero time-averaged torque. Thus, the time-averaged Lorentz torque need not be zero (and obey a Taylor state), but only be as small as the persistent torque from Reynolds stresses. The results of our model suggest that, as was pointed out by Dumberry & Bloxham (2003), the torque from Reynolds stresses in Earth’s core may also play a leading order role in the departure from a time-averaged Taylor state. The ratio of the torque from Reynolds stresses to the Lorentz torque scales as (uC/uA)2, so it is proportional to the square of the Alfvén number. Because the Alfvén number in our model is similar to that in Earth’s core, the baseline Taylorization factor of approximately 10−4 in Fig. 9 may also be representative of that expected in Earth’s core. However, because our model is 2-D, extrapolating our results to the 3-D magnetic field of Earth is clearly not straightforward. Furthermore, Ekman friction, which we have neglected, could be important in the time-averaged torque balance, especially if turbulent processes lead to an enhanced effective viscosity. Ekman friction could balance a part of the torque from Reynolds stresses and the Taylorization factor could then be smaller than 10−4. Regardless of its exact value, the high degree of Taylorization requires that large cancellations in the equilibrium Lorentz torque over a cylinder must occur in Earth’s core. It this sense, exploring dynamo solutions in the limit of a vanishing Lorentz torque remains a worthy goal (e.g. Livermore et al.2008; Wu & Roberts 2015). This being said, one must keep in mind that the correct Taylorization factor in Earth’s core may not be asymptotically close to zero but instead be closer to 10−4. 5 CONCLUSION We have constructed a 2-D reduced model of rotationally dominated magnetoconvection capable of producing distinct short- and long-timescale accelerations in the zonal flow. The short-timescale accelerations are free Alfvén waves, while the long-timescale accelerations are magnetically forced through the evolution of the Lorentz torque. The temporal changes in the magnetic field which drive the time-varying Lorentz torque are produced by the underlying convective flows, shearing and advecting the magnetic field on a timescale associated with convective eddies. Our results provide a dynamical explanation for the rigid decadal zonal accelerations that are inferred to exist in Earth’s core on the basis of the magnetic field observations and changes in LOD. Our results offer an alternative to the recent suggestion that the decadal zonal accelerations may not reflect deep seated rigid flows but are instead free Magnetic–Archemedian–Coriolis (MAC) waves in a stratified layer at the top of the core (Buffett 2014; Buffett et al.2016; Jaupart & Buffett 2017). The zonal flows of such MAC waves are characterized by a shear in the radial direction: flow at the CMB does not reflect flow deeper in the core in the axial direction. It is then more difficult to build a prediction of LOD changes based solely on the knowledge of flows at the CMB. It is nevertheless possible that, when properly taking into account the coupling of flows in the bulk of the core with these MAC waves, a prediction of core angular momentum change may be constructed so as to match the observed LOD variations (Buffett et al.2016). Yet the very fact that a very good match exists between the observed changes in the LOD and those predicted on the basis of purely rigid zonal flows suggests that deviations from rigidity are limited. Moreover, as our dynamical model shows, convective dynamics are expected to drive temporal changes in deep seated rigid zonal flows at decadal timescales. If the top of the core is stably stratified, rigid zonal flows may drive forced MAC oscillations and/or excite free MAC waves. Hence, the zonal flows at the CMB may consist of a combination of deep seated rigid zonal flows and flows that obey a MAC balance. However, because of the good match in LOD based on purely rigid flows, we speculate that the non-rigid MAC flows, if present, make up only a fraction of the total zonal flow. ACKNOWLEDGEMENTS This work was supported by both an NSERC/CRSNG Post-Graduate Scholarship and a Discovery grant from NSERC/CRSNG. Numerical simulations were performed on computing facilities provided by WestGrid and Compute/Calcul Canada. Figures were created using the Matplotlib package for Python (Hunter 2007). We thank Nathanaël Schaeffer and an anonymous reviewer for their constructive comments and suggestions. We further express our gratitude to Nathanaël Schaeffer, who kindly shared a numerical QG code which we extended over the course of this project. REFERENCES Aubert J., 2005. Steady zonal flows in spherical shell dynamos, J. Fluid Mech. , 542, 53– 67. https://doi.org/10.1017/S0022112005006129 Google Scholar CrossRef Search ADS   Aubert J., Gillet N., Cardin P., 2003. Quasigeostrophic models of convection in rotating spherical shells, Geochem. Geophys. Geosyst. , 4, 1052, doi:10.1029/2002GC000456. https://doi.org/10.1029/2002GC000456 Google Scholar CrossRef Search ADS   Aubert J., Gastine T., Fournier A., 2017. Spherical convective dynamos in the rapidly rotating asymptotic regime, J. Fluid Mech. , 813, 558– 593. https://doi.org/10.1017/jfm.2016.789 Google Scholar CrossRef Search ADS   Aurnou J., Calkins M.A., Cheng J.S., King E.M., Nieves D., Soderlund K.M., Stellmach S., 2015. Rotating convective turbulence in Earth and planetary cores, Phys. Earth planet. Inter. , 246, 52– 71. https://doi.org/10.1016/j.pepi.2015.07.001 Google Scholar CrossRef Search ADS   Braginsky S.I., 1970. Torsional magnetohydrodynamic vibrations in the Earth’s core and variations in day length, Geomagn. Aeron. , 10, 1– 10. Buffett B.A., 2010. Tidal dissipation and the strength of the earth’s internal magnetic field, Nature , 468, 952– 954. https://doi.org/10.1038/nature09643 Google Scholar CrossRef Search ADS PubMed  Buffett B.A., 2014. Geomagnetic fluctuations reveal stable stratification at the top of the Earth’s core, Nature , 507, 484– 487. https://doi.org/10.1038/nature13122 Google Scholar CrossRef Search ADS PubMed  Buffett B.A., Knezek N., Holme R., 2016. Evidence for MAC waves at the top of earth’s core and implications for variations in length of day, Geophys. J. Int. , 204, 1789– 1800. https://doi.org/10.1093/gji/ggv552 Google Scholar CrossRef Search ADS   Busse F., Or A., 1986. Convection in a rotating cylindrical annulus: thermal Rossby waves, J. Fluid Mech. , 166, 173– 187. https://doi.org/10.1017/S0022112086000095 Google Scholar CrossRef Search ADS   Cardin P., Olson P., 1994. Chaotic thermal convection in a rapidly rotating spherical shell: Consequences for flow in the outer core, Phys. Earth planet. Inter. , 82, 235– 259. https://doi.org/10.1016/0031-9201(94)90075-2 Google Scholar CrossRef Search ADS   Chao B.F., Chung W.Y., Shih Z.R., Hsieh Y.K., 2014. Earth’s rotation variations: a wavelet analysis, Terra Nova , 26, 260– 264. https://doi.org/10.1111/ter.12094 Google Scholar CrossRef Search ADS   Christensen U., 2002. Zonal flow driven by strongly supercritical convection in rotating spherical shells, J. Fluid Mech. , 470, 115– 133. https://doi.org/10.1017/S0022112002002008 Google Scholar CrossRef Search ADS   Christensen U., Wicht J., 2015. Numerical dynamo simulations, in Treatise on Geophysics , chap. 8.10, pp. 245– 277, ed. Schubert G., Elsevier. Christensen U.R., Aubert J., 2006. Scaling properties of convection-driven dynamos in rotating spherical shells and application to planetary magnetic fields, Geophys. J. Int. , 166, 97– 114. https://doi.org/10.1111/j.1365-246X.2006.03009.x Google Scholar CrossRef Search ADS   Cowling T., 1957. The dynamo maintenance of steady magnetic fields, Q. J. Mech. Appl. Math. , 10, 129– 136. https://doi.org/10.1093/qjmam/10.1.129 Google Scholar CrossRef Search ADS   Dumberry M., Bloxham J., 2003. Torque balance, Taylor’s constraint and torsional oscillations in a numerical model of the geodynamo, Phys. Earth planet. Inter. , 140, 29– 51. https://doi.org/10.1016/j.pepi.2003.07.012 Google Scholar CrossRef Search ADS   Finlay C.C., Dumberry M., Chulliat A., Pais M., 2010. Short timescale core dynamics: theory and observations, Space Sci. Rev. , 155, 177– 218. https://doi.org/10.1007/s11214-010-9691-6 Google Scholar CrossRef Search ADS   Gillet N., Jones C.A., 2006. The quasi-geostrophic model for rapidly rotating spherical convection outside the tangent cylinder, J. Fluid Mech. , 554, 343– 370. https://doi.org/10.1017/S0022112006009219 Google Scholar CrossRef Search ADS   Gillet N., Jault D., Canet E., Fournier A., 2010. Fast torsional waves and strong magnetic field within the Earth’s core, Nature , 465, 74– 77. https://doi.org/10.1038/nature09010 Google Scholar CrossRef Search ADS PubMed  Gillet N., Schaeffer N., Jault D., 2011. Rationale and geophysical evidence for quasi-geostrophic rapid dynamics within the Earth’s outer core, Phys. Earth planet. Inter. , 187, 380– 390. https://doi.org/10.1016/j.pepi.2011.01.005 Google Scholar CrossRef Search ADS   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. , 120, 3991– 4013. https://doi.org/10.1002/2014JB011786 Google Scholar CrossRef Search ADS   Gillet N., Jault D., Canet E., 2017. Excitation of traveling torsional normal modes in an Earth’s core model, Geophys. J. Int. , 210, 1503– 1516. https://doi.org/10.1093/gji/ggx237 Google Scholar CrossRef Search ADS   Gross R., 2015. Earth rotation variations - long period, in Treatise on Geophysics , chap. 3.09, pp. 215– 261, ed. Schubert G., Elsevier. Guervilly C., Cardin P., 2016. Subcritical convection of liquid metals in a rotating sphere using a quasi-geostrophic model, J. Fluid Mech. , 808, 61– 89. https://doi.org/10.1017/jfm.2016.631 Google Scholar CrossRef Search ADS   Guervilly C., Hughes D.W., Jones C.A., 2014. Large-scale vortices in rapidly rotating Rayleigh–Bénard convection, J. Fluid Mech. , 758, 407– 435. https://doi.org/10.1017/jfm.2014.542 Google Scholar CrossRef Search ADS   He Y., Sun W., 2007. Stability and convergence of the Crank–Nicolson/Adams–Bashforth scheme for the time-dependent Navier–Stokes equations, SIAM J. Numer. Anal. , 45( 2), 837– 869. https://doi.org/10.1137/050639910 Google Scholar CrossRef Search ADS   Holme R., 2015. Large-scale flow in the core, in Treatise on Geophysics , chap. 8.04, pp. 91– 113, ed. Schubert G., Elsevier. Holme R., de Viron O., 2013. Characterization and implications of intradecadal variations in length of day, Nature , 499, 202– 205. https://doi.org/10.1038/nature12282 Google Scholar CrossRef Search ADS PubMed  Hough S.S., 1897. On the application of harmonic analysis to the dynamical theory of the tides. Part I. On Laplace’s ‘Oscillations of the First Species,’ and on the dynamics of ocean currents, Phil. Trans. R. Soc. A , 189, 201– 257. https://doi.org/10.1098/rsta.1897.0009 Google Scholar CrossRef Search ADS   Hunter J., 2007. Matplotlib: A 2D graphics environment, Comput. Sci. Eng. , 9, 90– 95. https://doi.org/10.1109/MCSE.2007.55 Google Scholar CrossRef Search ADS   Jackson A., Bloxham J., Gubbins D., 1993. Time-dependent flow at the core surface and conservation of angular momentum in the coupled core-mantle system, in Dynamics of the Earth’s Deep Interior and Earth Rotation , Vol. 72, pp. 97– 107, eds Le Mouël J.-L., Smylie D.E., Herring T., AGU Geophysical Monograph. Google Scholar CrossRef Search ADS   Jault D., 2008. Axial invariance of rapidly varying diffusionless motions in the Earth’s core interior, Phys. Earth planet. Inter. , 166, 67– 76. https://doi.org/10.1016/j.pepi.2007.11.001 Google Scholar CrossRef Search ADS   Jault D., Gire C., Le Mouël J.-L., 1988. Westward drift, core motions and exchanges of angular momentum between core and mantle, Nature , 333, 353– 356. https://doi.org/10.1038/333353a0 Google Scholar CrossRef Search ADS   Jaupart E., Buffett B., 2017. Generation of MAC waves by convection in Earth’s core, Geophys. J. Int. , 209, 1326– 1336. https://doi.org/10.1093/gji/ggx088 Google Scholar CrossRef Search ADS   Labbé F., Jault D., Gillet N., 2015. On magnetostrophic inertia-less waves in quasi-geostrophic models of planetary cores, Geophys. Astrophys. Fluid Dyn. , 109, 587– 610. https://doi.org/10.1080/03091929.2015.1094569 Google Scholar CrossRef Search ADS   Livermore P.W., Ierley G., Jackson A., 2008. The structure of Taylor’s constraint in three dimensions, Proc. R. Soc. A , 464, 3149– 3174. https://doi.org/10.1098/rspa.2008.0091 Google Scholar CrossRef Search ADS   Matsui H., King E., Buffett B., 2014. Multiscale convection in a geodynamo simulation with uniform heat flux along the outer boundary, Geochem. Geophys. Geosyst. , 15, 3212– 3225. https://doi.org/10.1002/2014GC005432 Google Scholar CrossRef Search ADS   Pais M.A., Jault D., 2008. Quasi-geostrophic flows responsible for the secular variation of the Earth’s magnetic field, Geophys. J. Int. , 173, 421– 443. https://doi.org/10.1111/j.1365-246X.2008.03741.x Google Scholar CrossRef Search ADS   Proudman J., 1916. On the motion of solids in a liquid possessing vorticity, Proc. R. Soc. A , 92, 408– 424. https://doi.org/10.1098/rspa.1916.0026 Google Scholar CrossRef Search ADS   Roberts P., 2015. Theory of the geodynamo, in Treatise on Geophysics , chap. 8.03, pp. 57– 90, ed. Schubert G., Elsevier. Roberts P., King E., 2013. On the genesis of the Earth’s magnetism, Rep. Prog. Phys. , 76, 096801, doi:10.1088/0034-4885/76/9/096801. https://doi.org/10.1088/0034-4885/76/9/096801 Google Scholar CrossRef Search ADS PubMed  Roberts P.H., Yu Z., Russell C., 2007. On the 60-year signal from the core, Geophys. Astrophys. Fluid Dyn. , 101, 11– 35. https://doi.org/10.1080/03091920601083820 Google Scholar CrossRef Search ADS   Schaeffer N., Cardin P., 2005. Quasigeostrophic model of the instabilities of the Stewartson layer in flat and depth-varying containers, Phys. Fluids , 17, 104111, doi:10.1063/1.2073547. https://doi.org/10.1063/1.2073547 Google Scholar CrossRef Search ADS   Schaeffer N., Cardin P., 2006. Quasi-geostrophic kinematic dynamos at low magnetic Prandtl number, Earth planet. Sci. Lett. , 245, 595– 604. https://doi.org/10.1016/j.epsl.2006.03.024 Google Scholar CrossRef Search ADS   Schaeffer N., Lora Silva E., Pais M.A., 2016. Can core flows inferred from geomagnetic field models explain the Earth’s dynamo?, Geophys. J. Int. , 204, 568– 877. https://doi.org/10.1093/gji/ggv488 Google Scholar CrossRef Search ADS   Schaeffer N., Jault D., Nataf H.-C., Fournier A., 2017. Turbulent geodynamo simulations: a leap towards Earth’s core, Geophys. J. Int. , 211, 1– 29. https://doi.org/10.1093/gji/ggx265 Google Scholar CrossRef Search ADS   Stellmach S., Lischper M., Julien K., Vasil G., Cheng J.S., Ribeiro A., King E.M., Aurnou J.M., 2014. Approaching the asymptotic regime of rapidly rotating convection: Boundary layers versus interior dynamics, Phys. Rev. Lett. , 113, 254501, doi:10.1103/PhysRevLett.113.254501. https://doi.org/10.1103/PhysRevLett.113.254501 Google Scholar CrossRef Search ADS PubMed  Takahashi F., Matsushima M., Honkura Y., 2008. Scale variability in convection-driven MHD dynamos at low Ekman number, Phys. Earth planet. Inter. , 167, 168– 178. https://doi.org/10.1016/j.pepi.2008.03.005 Google Scholar CrossRef Search ADS   Taylor G.I., 1917. Motion of solids in fluids when the flow is not irrotational, Proc. R. Soc. A , 93, 99– 113. https://doi.org/10.1098/rspa.1917.0007 Google Scholar CrossRef Search ADS   Taylor J.B., 1963. The magneto-hydrodynamics of a rotating fluid and the earth’s dynamo problem, Proc. R. Soc. A , 274, 274– 283. https://doi.org/10.1098/rspa.1963.0130 Google Scholar CrossRef Search ADS   Teed R., Jones C.A., Tobias S., 2014. The dynamics and excitation of torsional waves in geodynamo simulations, Geophys. J. Int. , 196, 724– 735. https://doi.org/10.1093/gji/ggt432 Google Scholar CrossRef Search ADS   Teed R.J., Jones C.A., Tobias S.M., 2015. The transition to Earth-like torsional oscillations in magnetoconvection simulations, Earth planet. Sci. Lett. , 419, 22– 31. https://doi.org/10.1016/j.epsl.2015.02.045 Google Scholar CrossRef Search ADS   Wu C.-C., Roberts P.H., 2015. On magnetostrophic mean-field solutions of the geodynamo equations, Geophys. Astrophys. Fluid Dyn. , 109, 84– 110. Yadav R.K., Gastine T., Christensen U.R., Wolk S.J., Poppenhaeger K., 2016. Approaching a realistic force balance in geodynamo simulations, Proc. Natl. Acad. Sci. USA , 113, 12 065– 12 070. https://doi.org/10.1073/pnas.1608998113 Google Scholar CrossRef Search ADS   APPENDIX: CONSERVATION OF AXIAL ANGULAR MOMENTUM The total axial angular momentum $$\mathcal {L}_{z}$$ of the outer core is contained in the concentric cylinders rotating with axisymmetric angular velocity $$\frac{\overline{u_{\phi }}}{s}$$. The time variation of $$\mathcal {L}_{z}$$ is then   $$\frac{1}{4\pi } \frac{\partial }{\partial t}\mathcal {L}_{z} = \int _{s_{1}}^{s_{2}} \! s^{3} L \frac{\partial }{\partial t}\left(\frac{\overline{u_{\phi }}}{s}\right) \, \mathrm{d} s \, .$$ (A1) As shown by eq. (8), $$\frac{\partial }{\partial t}(\frac{\overline{u_{\phi }}}{s})$$ can be decomposed into four component torques. Defining the linear operator $$\mathcal {G}(f) = \int _{s_{1}}^{s_{2}} \! s^{3} L f \, \mathrm{d} s$$, eq. (A1) may be rewritten as a sum:   \begin{eqnarray} \frac{1}{4\pi } \frac{\partial }{\partial t}\mathcal {L}_{z} &= \mathcal {G} \left(\Gamma _{\textrm{V}}\right)+ \mathcal {G} \left(\Gamma _{\mathrm{L}_{1}} \right)+ \mathcal {G} \left(\Gamma _{\mathrm{L}_{2}} \right)+ \mathcal {G} \left(\Gamma _{\textrm{R}}\right)\, . \end{eqnarray} (A2)With the input torques   \begin{eqnarray} \Gamma _{\textrm{V}}= \frac{E}{s^{3} L} \frac{\partial }{\partial s}\left(s^{3} L \frac{\partial }{\partial s}\left(\frac{\overline{u_{\phi }}}{s}\right)\right)\, , \end{eqnarray} (A3)  \begin{eqnarray} \Gamma _{\textrm{L}_{1}}= \frac{1}{s^{3} L} \frac{\partial }{\partial s}\left(s^{3} L B_{0s}\frac{\overline{b_{\phi }}}{s}\right)\, , \end{eqnarray} (A4)  \begin{eqnarray} \Gamma _{\textrm{L}_{2}}= \frac{1}{s} \left(\overline{\frac{b_{s}}{s} \frac{\partial }{\partial s}\left(s b_{\phi }\right)} \right)\, , \end{eqnarray} (A5)  \begin{eqnarray} \Gamma _{\textrm{R}}= - \frac{1}{s} \left(\overline{\frac{u_{s}}{s} \frac{\partial }{\partial s}\left(s u_{\phi }\right)} \right)\, , \end{eqnarray} (A6) the first two terms on the right-hand side of eq. (A2) reduce to $$\mathcal {G}(\Gamma _{\textrm{V}}) = s^{3} L \frac{\partial }{\partial s}(\frac{\overline{u_{\phi }}}{s})|_{s_{1}}^{s_{2}}$$ and $$\mathcal {G}(\Gamma _{\rm{L}_{1}}) = s^{3} L B_{0s}\frac{\overline{b_{\phi }}}{s}|_{s_{1}}^{s_{2}}$$. Because we assume a stress-free condition on $$\overline{u_{\phi }}$$ and $$\overline{b_{\phi }} = 0$$ at s1 and s2, both terms are zero, reducing eq. (A2) to   \begin{eqnarray} \frac{1}{4\pi } \frac{\partial }{\partial t}\mathcal {L}_{z} = \mathcal {G} \left(\Gamma _{\mathrm{L}_{2}} \right)+ \mathcal {G} \left(\Gamma _{\textrm{R}}\right)\, . \end{eqnarray} (A7) Since both u and b are solenoidal and periodic in azimuth, eqs (A5) and (A6) may be rewritten as   \begin{eqnarray} \Gamma _{\textrm{L}_{2}}= \phantom{-} \frac{1}{s} \left(\frac{1}{s^{2} L} \frac{\partial }{\partial s}\left(s^{2} L \, \overline{b_{s}b_{\phi }} \right)- \overline{\left(\beta b_{s}- \frac{\partial }{\partial z}b_{z}\right)b_{\phi }} \right)\, , \end{eqnarray} (A8)  \begin{eqnarray} \Gamma _{\textrm{R}}= - \frac{1}{s} \left(\frac{1}{s^{2} L} \frac{\partial }{\partial s}\left(s^{2} L \, \overline{u_{s}u_{\phi }} \right)- \overline{\left(\beta u_{s}- \frac{\partial }{\partial z}u_{z}\right)u_{\phi }} \right)\, . \end{eqnarray} (A9) Application of the $$\mathcal {G}$$ operator to the first term on the right-hand side of each equation again produces terms which depend only on boundary values: $$s^{2} L \left.\overline{b_{s}b_{\phi }}\right|_{s_{1}}^{s_{2}}$$ and $$s^{2} L \left.\overline{u_{s}u_{\phi }}\right|_{s_{1}}^{s_{2}}$$. Because we use no-penetration (ψ = 0 →  us = 0) and a = 0 (→  bs = 0) boundary conditions, both are zero. This leaves   \begin{eqnarray} \frac{1}{4 \pi } \frac{\partial }{\partial t}\mathcal {L}_{z} \!=\! \mathcal {G} \!\left(\overline{\left(\beta u_{s}\!-\! \frac{\partial }{\partial z}u_{z}\right)u_{\phi }} \right)\!-\! \mathcal {G} \left(\overline{\left(\beta b_{s}\!-\! \frac{\partial }{\partial z}b_{z}\right)b_{\phi }} \right).\nonumber\\ \end{eqnarray} (A10) An axial profile must be assumed for both the u and b fields. The no-penetration boundary condition at the spherical top and bottom boundaries dictates that, defining the boundary’s normal vector as $${\bf \hat{n}}$$,   $$\mathbf {u}\cdot {\bf \hat{n}} = 0 \quad \Rightarrow \quad u_{z}(L) = - \frac{s}{L} u_{s}\, .$$ (A11)Assuming that the uz velocity profile varies linearly with z, and is zero in the equatorial plane,   $$\frac{\partial }{\partial z}u_{z}= \beta u_{s}\, ,$$ (A12) which is consistent with mass conservation ($$\boldsymbol{\nabla }\cdot \mathbf {u}= 0$$), and with the definition of the flow used in eqs (3) and (4). This causes the first term on the right-hand side of eq. (A10) to vanish. Our definition of the magnetic field perturbation b in eqs (3) and (5) follows the same form as that of u:   $$\frac{\partial }{\partial z}b_{z}= \beta b_{s}\, .$$ (A13)Substituting eqs (A12) and (A13) into eq. (A10) causes the remaining terms on the right-hand side of eq. (A2) to vanish. Thus, the angular momentum is conserved:   $$\frac{\partial }{\partial t}\mathcal {L}_{z} = 0 \, .$$ (A14) © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.

### Journal

Geophysical Journal InternationalOxford University Press

Published: Apr 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations