TY - JOUR AU - Tobias, S. M. AB - Abstract The predominant force balance in rapidly rotating planetary cores is between Coriolis, pressure, buoyancy and Lorentz forces. This magnetostrophic balance leads to a Taylor state where the spatially averaged azimuthal Lorentz force is compelled to vanish on cylinders aligned with the rotation axis. Any deviation from this state leads to a torsional oscillation, signatures of which have been observed in the Earth's secular variation and are thought to influence length of day variations via angular momentum conservation. In order to investigate the dynamics of torsional oscillations (TOs), we perform several 3-D dynamo simulations in a spherical shell. We find TOs, identified by their propagation at the correct Alfvén speed, in many of our simulations. We find that the frequency, location and direction of propagation of the waves are influenced by the choice of parameters. Torsional waves are observed within the tangent cylinder and also have the ability to pass through it. Several of our simulations display waves with core traveltimes of 4–6 yr. We calculate the driving terms for these waves and find that both the Reynolds force and ageostrophic convection acting through the Lorentz force are important in driving TOs. Earth rotation variations, Dynamo: theories and simulations, Rapid time variations 1 INTRODUCTION Rapidly rotating planetary dynamos, including the geodynamo, are believed to be operating under the magnetostrophic regime (see, for example, Jones 2011). In this regime, although the Lorentz force may be locally strong, the averaged azimuthal Lorentz force must vanish on geostrophic cylinders (Taylor 1963). A dynamo with a magnetic field organized in such a way is said to be in a Taylor state, which provides a severe constraint for dynamo generated fields. Any violation of the state can be represented as an acceleration of the cylinders and stretches radial magnetic field into azimuthal field. The resultant Lorentz force acts like a torsional spring in an attempt to restore the Taylor state (Braginsky 1970) and leads to the driving of torsional oscillations (TOs) of the cylinders. These oscillations, which are dependent only on cylindrical radius and time, are a type of Alfvén wave (Alfvén 1942). Torsional waves are believed to be continually driven in the Earth's core and are traceable in observational data. However, there has been some ambiguity as to the period for the fundamental modes of the TOs. Early observational data (Braginsky 1984) inferred a decadal timescale; however, more recent data obtained from core flow models by Gillet et al. (2010) show a much shorter period of around 4–6 yr. Previous work (Jault et al.1988; Jackson 1997; Zatman & Bloxham 1997; Bloxham et al.2002; Buffett et al.2009) has suggested that TOs may be responsible for various observed features of the Earth's dynamics; these include changes in length-of-day variations (Jault et al.1988; Jackson 1997) and geomagnetic jerks (Bloxham et al.2002). Additionally, it may be possible to infer information about the magnetic field within the core via core flow models (Zatman & Bloxham 1997; Buffett et al.2009; Gillet et al.2010). This is useful since geomagnetic data from the Earth's surface can only be reliably transferred down as far as the core-mantle boundary (CMB; Gubbins & Bloxham 1985). The specific driving mechanism of torsional waves in the Earth has not been identified, with several possible explanations of their excitation. Reynolds and Lorentz forces are often thought to play a role (Wicht & Christensen 2010), but other mechanisms also exist such as gravitational coupling between the inner core and the mantle (Mound & Buffett 2006) and induction due to variations in the external magnetospheric field (Légaut 2005). Numerical simulations are an obvious tool to analyse the dynamics of torsional waves; however, difficulties arise owing to the inability to reach appropriate Earth-like parameter values. Previous efforts (Dumberry & Bloxham 2003; Busse & Simitev 2005; Wicht & Christensen 2010) to locate torsional waves in simulations have been undertaken with Wicht & Christensen (2010) providing the most clear evidence yet of their observation in the region outside the tangent cylinder (OTC). Recent studies by Schaeffer et al. (2012) and Cox et al. (2013) have focused on the reflection of Alfvén waves at boundaries and their dispersion. Schaeffer et al. (2012) suggest that simulations run with rigid boundary conditions cannot exhibit wave reflection when the viscosity is too large. We investigate torsional wave production and dynamics in numerical simulations. We employ a systematic exploration of available parameter space and include analysis of the region inside the tangent cylinder (ITC) which has often been omitted in previous studies. This allows us to attempt to observe not only torsional waves ITC, but also the propagation of such waves across the tangent cylinder (TC). We estimate core traveltimes for the oscillations and, by bandpass filtering our data, we are able to determine whether the timescales that identified TOs operate on are correct. We also explore possible excitation mechanisms by calculating the relevant driving terms. In particular, we take a novel approach by separating the Lorentz force into its constituent parts: a restoring force and a driving force arising from ageostrophic convection. 2 MATHEMATICAL FORMULATION We adapt the model described by Jones et al. (2011) to incompressible systems (using the Boussinesq approximation). We shall extend to the compressible parameter space in future work. Our geometry is based on the Earth's core using a spherical polar coordinate system, (r, θ, ϕ). We consider a spherical shell that is radially bounded above at r = ro by an electrically insulating mantle and below at r = ri by an electrically insulating inner core. The system rotates about the vertical (z-axis) with rotation rate Ω and gravity acts radially inward so that g = −g r. The fluid is assumed to have constant values of ρ, ν, κ and η, the outer core density, kinematic viscosity, thermal diffusivity and magnetic diffusivity, respectively. Both the inner core and mantle are fixed in our model. Several recent papers (Sakuraba & Roberts 2009; Christensen et al.2010; Hori et al.2010) have argued that allowing for internal heat sources (or sinks) and imposing fixed heat flux (as opposed to fixed temperature) thermal boundary conditions in models may significantly influence the generation of solutions with Earth-like magnetic field morphologies. We, therefore, employ the use of fixed flux thermal boundary conditions for our work; specifically, we set zero flux on the CMB and the flux entering at the ICB is then balanced by a sink term, −ϵ, in the temperature equation. This mathematical setup is, in a physical sense, representative of a model for compositional convection, though we still refer to the buoyancy variable as the temperature. We shall be considering timescales shorter than any variation in the mean temperature so that the internal heating, with zero flux at the CMB, must satisfy the heat flux equation: \begin{eqnarray} -\!\frac{4\pi }{3}\epsilon \left(r_o^3 - r_i^3\right) = 4\pi \kappa r_i^2\left.\frac{\partial T}{\partial r}\right|_{r=r_i}, \end{eqnarray} (1) where T is the temperature. We non-dimensionalize the basic system of equations on the length scale, D = ro − ri, magnetic timescale, D2/η, temperature scale, ϵD2/η, and magnetic scale, |$\sqrt{\rho \mu _0\Omega \eta }$|⁠. The relevant system of coupled equations for velocity, u, magnetic field, B, temperature, T and pressure, p are: \begin{eqnarray} \frac{\partial {\bf u}}{\partial t} + ({\bf u\cdot \nabla }){\bf u} &= &-\frac{Pm}{E}\left[\nabla p + 2{\bf \hat{z}}\times {\bf u} - (\nabla \times {\bf B})\times {\bf B}\right] \nonumber \\ && +\; \frac{Pm^2Ra}{Pr} T{\bf r} + Pm\nabla ^2{\bf u}, \end{eqnarray} (2) \begin{equation} \frac{\partial {T}}{\partial t} + ({\bf u\cdot \nabla }){T} = \frac{Pm}{Pr}\nabla ^2 T - 1, \end{equation} (3) \begin{equation} \frac{\partial {\bf B}}{\partial t} - \nabla \times ({\bf u}\times {\bf B}) = \nabla ^2{\bf B}, \end{equation} (4) \begin{equation} \nabla \cdot {\bf u} = 0, \end{equation} (5) \begin{equation} \nabla \cdot {\bf B} = 0. \end{equation} (6) Eqs (2)–(4) are the incompressible Navier-Stokes, temperature and induction equations, respectively, and (5) and (6) describe the solenoidal conditions for velocity and magnetic field. The non-dimensional parameters appearing in our equations are the Rayleigh number, Ra, Ekman number, E, Prandtl number, Pr, and magnetic Prandtl number, Pm, defined by: \begin{eqnarray} Ra=\frac{g\alpha |\epsilon | D^5}{\nu \kappa \eta }, \quad E = \frac{\nu }{\Omega D^2}, \quad Pr = \frac{\nu }{\kappa }, \quad Pm = \frac{\nu }{\eta }. \end{eqnarray} (7) The radius ratio, β = ri/ro, is an additional parameter but in this work we restrict ourselves to the value appropriate to the Earth's core, namely β = 0.35. Note that under the non-dimensionalization chosen, the sink term in (3) has been scaled to unity. The magnitude of ϵ appears only in the definition of the Rayleigh number where it occupies the driving role usually taken by the temperature difference across the domain which appears in the classical definition of the Rayleigh number. 3 THEORY AND METHODS 3.1 Taylor's constraint and TOs The analysis of TOs requires consideration of the forces on geostrophic cylinders and hence the introduction of a cylindrical polar coordinate system, (s, ϕ, z), is beneficial. Averages over ϕ and z are required and hence for any scalar field A we define \begin{eqnarray} \skew4\bar{A}(t,s,z) = \frac{1}{2\pi }\int _0^{2\pi }A\mathrm{d}\phi , \quad \langle A\rangle (t,s,\phi ) = \frac{1}{h}\int _{z_-}^{z_+} A \mathrm{d}z. \end{eqnarray} (8) Here, h(s) = z+(s) − z−(s) and OTC we simply have that |$z_\pm =\pm \sqrt{r_o^2-s^2}$|⁠. Within the TC the definition of z± may remain the same if an average over the entire z domain is desired. However, ITC we may wish to average over the two hemispheres separately, which we refer to as ITCN and ITCS for north and south of the inner core, respectively. For ITCN (ITCS) we then have that |$z_+=\sqrt{r_o^2-s^2}$| and |$z_-=\sqrt{r_i^2-s^2}$||$(z_+=-\sqrt{r_i^2-s^2}$| and |$z_-=-\sqrt{r_o^2-s^2})$|⁠. For later convenience, we also define two further quantities for a scalar, or vector, field A. The first of these quantities, |$\skew4\tilde{A}$|⁠, is simply the time average of A over some time period, τ. The second quantity, A′, is the fluctuating part of A. Therefore, we define |$\skew4\tilde{A}$| and A′ by \begin{eqnarray} \skew4\tilde{A}(s,\phi ,z) = \frac{1}{\tau } \int _0^\tau A\mathrm{d}t \quad {\rm and } \quad A^\prime (t,s,\phi ,z)=A-\skew4\tilde{A}, \end{eqnarray} (9) respectively. A′ is useful because it removes from A the mean background state which only varies on a long timescale. Standard torsional oscillation theory relies on the ability to separate the timescales in this way successfully. The ϕ and z averages of the ϕ-component of (2) illustrate the forces that can accelerate geostrophic cylinders. Three such forces can be identified (Wicht & Christensen 2010); namely the Reynolds force, Lorentz force and viscous force leading to the equation \begin{eqnarray} \frac{\partial \langle \overline{u_\phi }\rangle }{\partial t} &=& -\langle \overline{\boldsymbol{\hat{\phi }}\cdot (\nabla \cdot {\bf uu})}\rangle + PmE^{-1}\langle \overline{\boldsymbol{\hat{\phi }}\cdot ((\nabla \times {\bf B})\times {\bf B})}\rangle \nonumber \\ & &+\, Pm\langle \overline{\boldsymbol{\hat{\phi }}\cdot \nabla ^2{\bf u}}\rangle \nonumber \\ & \equiv & F_{\rm R} + F_{\rm L} + F_{\rm V}. \end{eqnarray} (10) The Coriolis and buoyancy forces have vanished during the integration process since in the former there is no net flow across the cylinder and no ϕ-component in the latter. This has consequences in the core where the fluid is believed, at leading order, to be in magnetostrophic balance (between Lorentz, Coriolis and Archimedean forces). Taylor (1963) noted that in systems where the force balance is magnetostrophic the constraint \begin{eqnarray} F_{\rm L} = 0, \end{eqnarray} (11) arises. The Lorentz force can be partially integrated (see, for example, Wicht & Christensen 2010) to give \begin{eqnarray} F_{\rm L} = \frac{Pm}{E}\frac{1}{hs^2}\frac{\partial }{\partial s}s^2h\langle \overline{B_sB_\phi }\rangle + \frac{Pm}{E}\frac{1}{h}\left[\frac{s}{z}\overline{B_sB_\phi }+\overline{B_zB_\phi }\right]_{z_-}^{z_+}. \end{eqnarray} (12) We are able to neglect the magnetic coupling terms in this expression at this stage due to our use of insulating boundary conditions at both the CMB and the inner core boundary (ICB; Jones et al.2011). However, if one were to allow for a conducting inner core (or mantle), the contribution from these surface terms would be non-zero resulting in an additional forcing in the system that is not discussed further here. For discussion of how this coupling term arises see Roberts & Aurnou (2012). Upon consideration of the time derivative of the expression for FL in (12), we find that we require expressions for the time derivatives of components of the magnetic field. We substitute from the induction equation and retain all terms on the right-hand side of (4), to determine that \begin{eqnarray} \skew4\dot{F}_{\rm L}\ =\ \frac{Pm}{E}\frac{1}{hs^2}\frac{\partial }{\partial s}s^2h\langle\; \overline{\dot{B}_sB_\phi + B_s\dot{B}_\phi }\;\rangle \end{eqnarray} (13) \begin{eqnarray} \hphantom{\skew4\dot{F}_{\rm L}\;}&=& \ \frac{Pm}{E}\frac{1}{hs^2}\frac{\partial }{\partial s}s^2h\bigg\lbrace \left\langle\; \overline{sB_s({\bf B}\cdot \nabla )\frac{u_\phi }{s}}\;\right\rangle \nonumber \\ && +\; \left\langle\; \overline{\frac{B_\phi }{s}({\bf B}\cdot \nabla )(su_s)}\;\right\rangle \nonumber \\ && -\; \left\langle\; \overline{\left({\bf u}\cdot \nabla +\frac{2}{s^2}\right)(B_sB_\phi )}\;\right\rangle \nonumber \\ && +\; \left\langle\; \overline{B_s\nabla ^2B_\phi + B_\phi \nabla ^2 B_s}\;\right\rangle \bigg\rbrace . \end{eqnarray} (14) In order to make further progress, we use the definitions (9) to split the velocity and magnetic field into mean and fluctuating parts. Previous studies (Wicht & Christensen 2010; Roberts & Aurnou 2012) have essentially assumed that the mean quantities, |$\tilde{\bf u}$| and |$\tilde{\bf B}$|⁠, are the principal parts of the Taylor state and that the fluctuating quantities, u′ and B′, are perturbations associated with the TOs. However, this is not the full picture since it requires the assumption that u′ is purely geostrophic as explicitly stated by Taylor (1963). In reality, the convection will be operating, to some degree, on all timescales and this phenomenon is likely to be an important driving mechanism. Hence, rather than assuming a geostrophic form for our velocity fluctuation we develop a different technique and instead split u′ into geostrophic (sζ′) and ageostrophic parts (⁠|${\bf u}_a^{\prime }$|⁠) so that \begin{eqnarray} {\bf u} = \tilde{\bf u} + {\bf u}^{\prime } = \tilde{\bf u} + s\zeta ^{\prime }(s,t)\boldsymbol{\hat{\phi }} + {\bf u}_a^{\prime }, \qquad {\bf B} = \tilde{\bf B} + {\bf B}^{\prime }. \end{eqnarray} (15) Upon substitution of these forms into our expression for |$\skew4\dot{F}_{\rm L}$|⁠, we find that ζ′ only appears in the first term on the right-hand side of (14). Considering only the mean magnetic field parts of this term and calling it |$\skew4\dot{F}_{\rm LR}$| gives \begin{eqnarray} \skew4\dot{F}_{\rm LR} = \frac{1}{hs^2}\frac{\partial }{\partial s}\left(s^3h U_{\rm A}^2 \frac{\partial \zeta }{\partial s}\right), \qquad U_{\rm A}=\sqrt{\frac{Pm}{E}\langle \overline{\tilde{B}_s^2}\rangle }, \end{eqnarray} (16) where we have defined the Alfvén speed, UA. Eq. (14) can then be written as \begin{eqnarray} \skew4\dot{F}_{\rm L} = \skew4\dot{F}_{\rm LR} + \skew4\dot{F}_{\rm LD}, \end{eqnarray} (17) where |$\skew4\dot{F}_{\rm LD}$| is a complicated expression made up of the remaining terms on the right-hand side of (14). Thus, it involves terms containing the components of |$\tilde{\bf B}$|⁠, B′, |$\tilde{\bf u}$|⁠, |${\bf u}_a^{\prime }$|⁠, as well as ζ′. If we now take the time derivative of (10) and use the result of (17) we find that \begin{eqnarray} s\ddot{\zeta }^{\prime } = \skew4\dot{F}_{\rm LR} + \skew4\dot{F}_{\rm LD} + \skew4\dot{F}_{\rm R} + \skew4\dot{F}_{\rm V}, \end{eqnarray} (18) noting that |$\langle \overline{\boldsymbol{\hat{\phi }}\cdot {\bf u}_a^{\prime }}\rangle =0$| by definition. By writing the expression for |$\ddot{\zeta }^{\prime }$| in this way, we have been able to separate the term involved in the balance of the torsional wave equation from the remaining terms. The standard canonical wave equation as found in previous work (see, for example, Braginsky 1970) is represented by |$s\ddot{\zeta }^{\prime } = \skew4\dot{F}_{\rm LR}$|⁠. Consequently, if we time integrate (18) to acquire \begin{eqnarray} s\dot{\zeta }^{\prime } - F_{\rm LR} = F_{\rm LD} + F_{\rm R} + F_{\rm V}, \end{eqnarray} (19) we find that FLR is the restoring force whereas FLD, FR and FV are driving forces. Torsional waves in the core must be driven and dissipated by some mechanism(s) and hence the terms on the right-hand side of (19), namely FR, FV and FLD, fulfil this role. They are driving (and dissipative) forces which are able to create, destroy and alter the nature of propagating torsional waves. When performing diagnostics on our simulations, one of our interests will be analysing the terms on the right-hand side of (19). This will allow us to identify which forces are able to act as excitation mechanisms at various points in the domain. We look at this in Section 4.5. 3.2 Output parameters In addition to quantities described in Section 3.1, we also output several other parameters from our simulations. The magnetic Reynolds number, Elsasser number, Rossby number, relative dipole moment, Taylorization parameter and velocity correlation parameters are defined by \begin{equation} Rm = \frac{UD}{\eta }, \end{equation} (20) \begin{equation} \Lambda = \frac{|B|^2}{\rho \mu \eta \Omega }, \end{equation} (21) \begin{equation} Ro = \frac{U}{\Omega D}, \end{equation} (22) \begin{equation} f_{\rm dip} = \left(\frac{E_M^{(1,0)}(r_o)}{\sum _{l=1}^{12}\sum _{m=0}^{l}E_M^{(l,m)}(r_o)}\right)^{1/2}, \end{equation} (23) \begin{equation} \mathcal {T} = \frac{\langle \overline{\boldsymbol{\hat{\phi }}\cdot (\nabla \times {\bf B})\times {\bf B}}\rangle }{\langle \overline{|\boldsymbol{\hat{\phi }}\cdot (\nabla \times {\bf B})\times {\bf B}|}\rangle }, \end{equation} (24) \begin{equation} U_{\rm C} = \frac{\sqrt{|\langle \overline{u_\phi }\rangle ^2-\widetilde{\langle \overline{u_\phi }\rangle ^2}|}}{Rm}, \end{equation} (25) \begin{equation} U^{\prime }_{{\rm C}} = \sqrt{\left|\frac{\langle \overline{u_\phi }\rangle ^2-\widetilde{\langle \overline{u_\phi }\rangle ^2}}{Rm^2-\widetilde{Rm^2}}\right|}, \end{equation} (26) respectively. Here, U and B are (dimensional) rms values of velocity and magnetic field, respectively, and |$E_M^{(l,m)}(r)$| represents the magnetic energy in the (l, m) harmonic at radius r. The parameter UC is a measure of the relative amplitude of the short timescale azimuthal velocity and therefore gives an indication of the strength of TOs compared with the underlying flow. Similarly, |$U^{\prime }_{\rm C}$| is a ratio of the short timescale geostrophic flow to the total short timescale flow. Owing to our choice of non-dimensionalization, the magnetic Reynolds and Elsasser numbers can be identified with the non-dimensional velocity and square of the magnetic field, respectively. All parameters can be averaged over time (as well as space for |$\mathcal {T}$|⁠, UC and |$U^{\prime }_{\rm C}$|⁠). The parameters defined in eqs (20)–(24) give an indication of the sort of regime that the dynamo is in, a point we address in Section 4.1. 3.3 Methods We perform several simulations, using the Leeds spherical dynamo code (Jones et al.2011) which uses a pseudo-spectral numerical scheme with finite differences in the radial direction. We run the code at parameter regimes and with boundary conditions that facilitate the production of Earth-like dynamos. Guided by previous work (Hori et al.2010), we therefore employ the use of fixed flux thermal boundary conditions for all of our simulations. Rigid kinematic boundary conditions are primarily used, although one set of simulations is repeated with stress-free boundaries as way of comparison. In parameter space, we perform simulations at a range of Ekman numbers since the existence of TOs requires the dynamo to be near magnetostrophic balance, which in turn is dependent on a small Ekman number. Thus, by decreasing the Ekman number over the range 10−4 to 10−6 TOs should become more apparent. We focus on Pr = 1 and each simulation is at the same value of criticality; that is Ra/Rac ≃ 8.32 for all runs. However, we do vary the magnetic Prandtl number, Pm ∈ [1, 5], in order to allow for a range in the magnetic field strength. The values of Rac used are for the onset of non-magnetic convection (see, for example, Dormy et al.2004). Table 1 displays the input parameters for the set of runs performed as well as the kinetic boundary conditions employed. Table 1. Input and output parameter sets used for the various simulations. Note that all runs have fixed flux thermal boundary conditions with zero flux on the outer boundary and an internal heat sink. Run . E . Ra . Pr . Pm . BCs . τ . UA(ro) . τE . Rm . Λ . Ro . fdip . |$\mathcal {T}$| . UC . |$U^{\prime }_{\rm C}$| . TOs . 4R1 10−4 4.937 × 106 1 1 NS 0.02 0.07 0.05 98.12 0.896 0.010 0.890 0.586 0.03 0.25 – 4R2 10−4 4.937 × 106 1 2 NS 0.02 1.44 1.16 135.60 1.888 0.007 0.867 0.587 0.03 0.20 – 4R3 10−4 4.937 × 106 1 3 NS 0.02 15.67 12.65 152.39 5.672 0.005 0.847 0.267 0.03 0.22 – 4R4 10−4 4.937 × 106 1 4 NS 0.014 22.26 12.57 183.97 10.358 0.005 0.776 0.198 0.07 0.45 OTC 4R5 10−4 4.937 × 106 1 5 NS 0.014 29.38 16.60 217.05 15.621 0.004 0.741 0.168 0.13 0.46 OTC, ITC 5R1 10−5 1 × 108 1 1 NS 0.006 5.02 1.21 128.54 0.319 0.001 0.924 0.389 0.12 0.41 OTC 5R2 10−5 1 × 108 1 2 NS 0.006 14.28 3.46 203.35 1.740 0.001 0.904 0.373 0.11 0.55 OTC, ITC 5R3 10−5 1 × 108 1 3 NS 0.006 90.07 21.80 330.52 16.197 0.001 0.722 0.222 0.10 0.61 OTC, ITC 5R4 10−5 1 × 108 1 4 NS 0.003 90.27 10.93 355.91 17.433 0.001 0.713 0.205 0.11 0.61 OTC, ITC 5R5 10−5 1 × 108 1 5 NS 0.003 123.90 15.00 437.07 19.252 0.001 0.742 0.118 0.13 0.63 OTC 6.5R1 5 × 10−6 2.493 × 108 1 1 NS 0.004 7.77 1.26 155.43 0.325 0.001 0.917 0.381 0.09 0.40 OTC, ITC 6.5R2 5 × 10−6 2.493 × 108 1 2 NS 0.004 22.08 3.56 267.72 2.400 0.001 0.955 0.298 0.10 0.58 OTC, ITC 6.5R3 5 × 10−6 2.493 × 108 1 3 NS 0.004 29.17 4.71 383.57 3.631 0.001 0.946 0.298 0.11 0.60 OTC, ITC 6.5R4 5 × 10−6 2.493 × 108 1 4 NS 0.002 259.22 20.92 575.84 23.637 0.001 0.752 0.161 0.12 0.63 OTC 6.5R5 5 × 10−6 2.493 × 108 1 5 NS 0.002 243.47 19.65 599.00 20.080 0.001 0.752 0.168 0.13 0.65 OTC, ITC 6R1 10−6 2.132 × 109 1 1 NS 0.002 15.66 1.26 372.87 0.561 <0.001 0.918 0.193 0.13 0.47 OTC, ITC 5F1 10−5 1.265 × 108 1 1 SF 0.008 5.09 1.64 172.71 0.368 0.002 0.918 0.350 0.09 0.51 OTC 5F2 10−5 1.265 × 108 1 2 SF 0.005 16.59 3.35 226.40 2.164 0.001 0.955 0.304 0.10 0.41 OTC, ITC 5F3 10−5 1.265 × 108 1 3 SF 0.003 94.57 11.45 336.97 18.817 0.001 0.676 0.218 0.10 0.48 OTC, ITC 5F4 10−5 1.265 × 108 1 4 SF 0.003 89.94 10.89 402.81 18.578 0.001 0.738 0.175 0.10 0.62 OTC, ITC 5F5 10−5 1.265 × 108 1 5 SF 0.002 109.47 8.83 560.84 23.636 0.001 0.719 0.214 0.11 0.60 OTC, ITC Run . E . Ra . Pr . Pm . BCs . τ . UA(ro) . τE . Rm . Λ . Ro . fdip . |$\mathcal {T}$| . UC . |$U^{\prime }_{\rm C}$| . TOs . 4R1 10−4 4.937 × 106 1 1 NS 0.02 0.07 0.05 98.12 0.896 0.010 0.890 0.586 0.03 0.25 – 4R2 10−4 4.937 × 106 1 2 NS 0.02 1.44 1.16 135.60 1.888 0.007 0.867 0.587 0.03 0.20 – 4R3 10−4 4.937 × 106 1 3 NS 0.02 15.67 12.65 152.39 5.672 0.005 0.847 0.267 0.03 0.22 – 4R4 10−4 4.937 × 106 1 4 NS 0.014 22.26 12.57 183.97 10.358 0.005 0.776 0.198 0.07 0.45 OTC 4R5 10−4 4.937 × 106 1 5 NS 0.014 29.38 16.60 217.05 15.621 0.004 0.741 0.168 0.13 0.46 OTC, ITC 5R1 10−5 1 × 108 1 1 NS 0.006 5.02 1.21 128.54 0.319 0.001 0.924 0.389 0.12 0.41 OTC 5R2 10−5 1 × 108 1 2 NS 0.006 14.28 3.46 203.35 1.740 0.001 0.904 0.373 0.11 0.55 OTC, ITC 5R3 10−5 1 × 108 1 3 NS 0.006 90.07 21.80 330.52 16.197 0.001 0.722 0.222 0.10 0.61 OTC, ITC 5R4 10−5 1 × 108 1 4 NS 0.003 90.27 10.93 355.91 17.433 0.001 0.713 0.205 0.11 0.61 OTC, ITC 5R5 10−5 1 × 108 1 5 NS 0.003 123.90 15.00 437.07 19.252 0.001 0.742 0.118 0.13 0.63 OTC 6.5R1 5 × 10−6 2.493 × 108 1 1 NS 0.004 7.77 1.26 155.43 0.325 0.001 0.917 0.381 0.09 0.40 OTC, ITC 6.5R2 5 × 10−6 2.493 × 108 1 2 NS 0.004 22.08 3.56 267.72 2.400 0.001 0.955 0.298 0.10 0.58 OTC, ITC 6.5R3 5 × 10−6 2.493 × 108 1 3 NS 0.004 29.17 4.71 383.57 3.631 0.001 0.946 0.298 0.11 0.60 OTC, ITC 6.5R4 5 × 10−6 2.493 × 108 1 4 NS 0.002 259.22 20.92 575.84 23.637 0.001 0.752 0.161 0.12 0.63 OTC 6.5R5 5 × 10−6 2.493 × 108 1 5 NS 0.002 243.47 19.65 599.00 20.080 0.001 0.752 0.168 0.13 0.65 OTC, ITC 6R1 10−6 2.132 × 109 1 1 NS 0.002 15.66 1.26 372.87 0.561 <0.001 0.918 0.193 0.13 0.47 OTC, ITC 5F1 10−5 1.265 × 108 1 1 SF 0.008 5.09 1.64 172.71 0.368 0.002 0.918 0.350 0.09 0.51 OTC 5F2 10−5 1.265 × 108 1 2 SF 0.005 16.59 3.35 226.40 2.164 0.001 0.955 0.304 0.10 0.41 OTC, ITC 5F3 10−5 1.265 × 108 1 3 SF 0.003 94.57 11.45 336.97 18.817 0.001 0.676 0.218 0.10 0.48 OTC, ITC 5F4 10−5 1.265 × 108 1 4 SF 0.003 89.94 10.89 402.81 18.578 0.001 0.738 0.175 0.10 0.62 OTC, ITC 5F5 10−5 1.265 × 108 1 5 SF 0.002 109.47 8.83 560.84 23.636 0.001 0.719 0.214 0.11 0.60 OTC, ITC Open in new tab Table 1. Input and output parameter sets used for the various simulations. Note that all runs have fixed flux thermal boundary conditions with zero flux on the outer boundary and an internal heat sink. Run . E . Ra . Pr . Pm . BCs . τ . UA(ro) . τE . Rm . Λ . Ro . fdip . |$\mathcal {T}$| . UC . |$U^{\prime }_{\rm C}$| . TOs . 4R1 10−4 4.937 × 106 1 1 NS 0.02 0.07 0.05 98.12 0.896 0.010 0.890 0.586 0.03 0.25 – 4R2 10−4 4.937 × 106 1 2 NS 0.02 1.44 1.16 135.60 1.888 0.007 0.867 0.587 0.03 0.20 – 4R3 10−4 4.937 × 106 1 3 NS 0.02 15.67 12.65 152.39 5.672 0.005 0.847 0.267 0.03 0.22 – 4R4 10−4 4.937 × 106 1 4 NS 0.014 22.26 12.57 183.97 10.358 0.005 0.776 0.198 0.07 0.45 OTC 4R5 10−4 4.937 × 106 1 5 NS 0.014 29.38 16.60 217.05 15.621 0.004 0.741 0.168 0.13 0.46 OTC, ITC 5R1 10−5 1 × 108 1 1 NS 0.006 5.02 1.21 128.54 0.319 0.001 0.924 0.389 0.12 0.41 OTC 5R2 10−5 1 × 108 1 2 NS 0.006 14.28 3.46 203.35 1.740 0.001 0.904 0.373 0.11 0.55 OTC, ITC 5R3 10−5 1 × 108 1 3 NS 0.006 90.07 21.80 330.52 16.197 0.001 0.722 0.222 0.10 0.61 OTC, ITC 5R4 10−5 1 × 108 1 4 NS 0.003 90.27 10.93 355.91 17.433 0.001 0.713 0.205 0.11 0.61 OTC, ITC 5R5 10−5 1 × 108 1 5 NS 0.003 123.90 15.00 437.07 19.252 0.001 0.742 0.118 0.13 0.63 OTC 6.5R1 5 × 10−6 2.493 × 108 1 1 NS 0.004 7.77 1.26 155.43 0.325 0.001 0.917 0.381 0.09 0.40 OTC, ITC 6.5R2 5 × 10−6 2.493 × 108 1 2 NS 0.004 22.08 3.56 267.72 2.400 0.001 0.955 0.298 0.10 0.58 OTC, ITC 6.5R3 5 × 10−6 2.493 × 108 1 3 NS 0.004 29.17 4.71 383.57 3.631 0.001 0.946 0.298 0.11 0.60 OTC, ITC 6.5R4 5 × 10−6 2.493 × 108 1 4 NS 0.002 259.22 20.92 575.84 23.637 0.001 0.752 0.161 0.12 0.63 OTC 6.5R5 5 × 10−6 2.493 × 108 1 5 NS 0.002 243.47 19.65 599.00 20.080 0.001 0.752 0.168 0.13 0.65 OTC, ITC 6R1 10−6 2.132 × 109 1 1 NS 0.002 15.66 1.26 372.87 0.561 <0.001 0.918 0.193 0.13 0.47 OTC, ITC 5F1 10−5 1.265 × 108 1 1 SF 0.008 5.09 1.64 172.71 0.368 0.002 0.918 0.350 0.09 0.51 OTC 5F2 10−5 1.265 × 108 1 2 SF 0.005 16.59 3.35 226.40 2.164 0.001 0.955 0.304 0.10 0.41 OTC, ITC 5F3 10−5 1.265 × 108 1 3 SF 0.003 94.57 11.45 336.97 18.817 0.001 0.676 0.218 0.10 0.48 OTC, ITC 5F4 10−5 1.265 × 108 1 4 SF 0.003 89.94 10.89 402.81 18.578 0.001 0.738 0.175 0.10 0.62 OTC, ITC 5F5 10−5 1.265 × 108 1 5 SF 0.002 109.47 8.83 560.84 23.636 0.001 0.719 0.214 0.11 0.60 OTC, ITC Run . E . Ra . Pr . Pm . BCs . τ . UA(ro) . τE . Rm . Λ . Ro . fdip . |$\mathcal {T}$| . UC . |$U^{\prime }_{\rm C}$| . TOs . 4R1 10−4 4.937 × 106 1 1 NS 0.02 0.07 0.05 98.12 0.896 0.010 0.890 0.586 0.03 0.25 – 4R2 10−4 4.937 × 106 1 2 NS 0.02 1.44 1.16 135.60 1.888 0.007 0.867 0.587 0.03 0.20 – 4R3 10−4 4.937 × 106 1 3 NS 0.02 15.67 12.65 152.39 5.672 0.005 0.847 0.267 0.03 0.22 – 4R4 10−4 4.937 × 106 1 4 NS 0.014 22.26 12.57 183.97 10.358 0.005 0.776 0.198 0.07 0.45 OTC 4R5 10−4 4.937 × 106 1 5 NS 0.014 29.38 16.60 217.05 15.621 0.004 0.741 0.168 0.13 0.46 OTC, ITC 5R1 10−5 1 × 108 1 1 NS 0.006 5.02 1.21 128.54 0.319 0.001 0.924 0.389 0.12 0.41 OTC 5R2 10−5 1 × 108 1 2 NS 0.006 14.28 3.46 203.35 1.740 0.001 0.904 0.373 0.11 0.55 OTC, ITC 5R3 10−5 1 × 108 1 3 NS 0.006 90.07 21.80 330.52 16.197 0.001 0.722 0.222 0.10 0.61 OTC, ITC 5R4 10−5 1 × 108 1 4 NS 0.003 90.27 10.93 355.91 17.433 0.001 0.713 0.205 0.11 0.61 OTC, ITC 5R5 10−5 1 × 108 1 5 NS 0.003 123.90 15.00 437.07 19.252 0.001 0.742 0.118 0.13 0.63 OTC 6.5R1 5 × 10−6 2.493 × 108 1 1 NS 0.004 7.77 1.26 155.43 0.325 0.001 0.917 0.381 0.09 0.40 OTC, ITC 6.5R2 5 × 10−6 2.493 × 108 1 2 NS 0.004 22.08 3.56 267.72 2.400 0.001 0.955 0.298 0.10 0.58 OTC, ITC 6.5R3 5 × 10−6 2.493 × 108 1 3 NS 0.004 29.17 4.71 383.57 3.631 0.001 0.946 0.298 0.11 0.60 OTC, ITC 6.5R4 5 × 10−6 2.493 × 108 1 4 NS 0.002 259.22 20.92 575.84 23.637 0.001 0.752 0.161 0.12 0.63 OTC 6.5R5 5 × 10−6 2.493 × 108 1 5 NS 0.002 243.47 19.65 599.00 20.080 0.001 0.752 0.168 0.13 0.65 OTC, ITC 6R1 10−6 2.132 × 109 1 1 NS 0.002 15.66 1.26 372.87 0.561 <0.001 0.918 0.193 0.13 0.47 OTC, ITC 5F1 10−5 1.265 × 108 1 1 SF 0.008 5.09 1.64 172.71 0.368 0.002 0.918 0.350 0.09 0.51 OTC 5F2 10−5 1.265 × 108 1 2 SF 0.005 16.59 3.35 226.40 2.164 0.001 0.955 0.304 0.10 0.41 OTC, ITC 5F3 10−5 1.265 × 108 1 3 SF 0.003 94.57 11.45 336.97 18.817 0.001 0.676 0.218 0.10 0.48 OTC, ITC 5F4 10−5 1.265 × 108 1 4 SF 0.003 89.94 10.89 402.81 18.578 0.001 0.738 0.175 0.10 0.62 OTC, ITC 5F5 10−5 1.265 × 108 1 5 SF 0.002 109.47 8.83 560.84 23.636 0.001 0.719 0.214 0.11 0.60 OTC, ITC Open in new tab Each run is initially time integrated from a random state for at least one-tenth of a magnetic diffusion time apart from run 6R1 which is run for a shorter period due to resolution constraints. We have also checked that the system has achieved statistical equilibrium and thus we have integrated past any transient stage of the solution. In order to search for TOs we then analyse a period of time, τ, of every run. The value of τ for each run, indicated in Table 1, is run dependent and varies between 0.002 and 0.02 of a diffusion time. In order to estimate core traveltimes of TOs, we must choose a method of scaling our non-dimensional time back to a dimensional time. Consideration of the diffusion timescale reveals that it is not ideal for conversion in our study of TOs since our fields in these units are often too strong. Therefore, we choose to convert by matching the Alfvén speed at the CMB. Using 0.2 mT as the magnetic field strength in the cylindrical radial direction at the CMB (Gillet et al.2010) and ρ = 1 × 104 kg m−3 (as well as μ0 = 4π × 10−7 Hm−1) this gives an approximate physical Earth value of the Alfvén speed at the CMB: |$U_{\rm A}^E=1.784\times 10^{-3}$| m s−1. We use the non-dimensional values of UA(ro) (the Alfvén speed at the equator at the CMB) given in Table 1, as well as D = 2.27 × 106 m, to calculate the dimensional version of τ for each run, using \begin{equation} \tau ^E = \frac{DU_{\rm A}}{U_{\rm A}^E}\tau . \end{equation} (27) This value is also displayed, in units of years, in Table 1 and we shall discuss core traveltimes in Section 4.3. By including the region ITC in our analysis, we present ourselves with a complication since it is not obvious how to deal with the regions north and south of the inner core. For example, when performing averages over z do we average over the entire vertical from pole to pole or instead retain the distinction between the hemispheres? Consequently, there is also the issue of how to treat waves propagating across the TC since they may originate (or terminate) in either hemisphere. These issues were not present in the previous work on torsional wave analysis in dynamo simulations (Wicht & Christensen 2010) where the region ITC was omitted. We choose to allow for both scenarios by performing both sets of averages. Therefore, in our analysis we average over the entire region ITC, but also perform averages over each hemisphere separately (that is over ITCN and ITCS). For the region OTC, averages are always performed across all z-space. 4 NUMERICAL RESULTS 4.1 Field strength and morphology Time-averaged output parameters calculated from our numerical results are displayed in Table 1. The values of |$\mathcal {T}$|⁠, UC and |$U^{\prime }_{{\rm C}}$| have also been averaged spatially in s so that they are indications of the degree of Taylorization and relative TO strength for the entire shell. In Table 1, we also indicate, for each run, whether TOs were identified and if so, the region(s) of the shell that they were observed. Within our full set of simulations we are able to identify two major magnetohydrodynamic regimes for which the fluid in each run can organize itself. The weak field regime has Λ ∼ O(1) whereas the strong field regime has a much larger Elsasser number and, as one would expect, the latter regime is found at larger values of the magnetic Prandtl number. Velocity structures are larger in the strong field regime. However, it should be noted that even in the weak field regime the convection is not as small scale as one may expect for such a rapidly rotating system. This is due to the employment of fixed flux thermal boundary conditions, which have been found to significantly affect the size of velocity structures (Hori et al.2010). With current estimates that Rm ≈ 1000 for the Earth's outer core, Table 1 indicates that only our high Pm, low E runs begin to approach Earth-like magnetic Reynolds numbers since larger Rayleigh numbers are required for more Earth-like dynamos (Christensen et al.2010). Simulations in the strong field regime produce Elsasser numbers too large for the Earth where Λ ∼ O(1), however this may be a result of the way in which the magnetic field is scaled. Christensen (2010) points out that the Elsasser number may decrease as the Ekman number becomes small. The converse is true of the dipolarity, which decreases to near Earth-like values for our larger Pm runs. In Figs 1 and 2, we plot Br, truncated at harmonic degree 12, at the CMB for runs at two different values of Pm. Although both figures show dipolar fields, the dipolarity is visibly stronger in Fig. 1 than Fig. 2, which has patches of reversed flux. These plots are representative of the radial magnetic field for the two different regimes seen across all of our runs. As we shall discuss later, the two regimes will also have implications on where and what sort of TOs can be found. Figure 1. Open in new tabDownload slide The radial magnetic field at the CMB for the run 5R2. Figure 1. Open in new tabDownload slide The radial magnetic field at the CMB for the run 5R2. Figure 2. Open in new tabDownload slide The radial magnetic field at the CMB for the run 6.5R5. Figure 2. Open in new tabDownload slide The radial magnetic field at the CMB for the run 6.5R5. The values of |$\mathcal {T}$| displayed in Table 1 show that the degree of Taylorization is well correlated with the Elsasser number and therefore a strong field regime is beneficial for Taylorization. Our values for the Taylorization parameter broadly match those found in previous literature (Wicht & Christensen 2010) for the parameters used. The velocity correlation parameter, UC, shows that TOs typically contribute approximately 10 per cent of the total flow, which is slightly larger than that found by Wicht & Christensen (2010) but is a similar value to the contribution they find from geostrophic cylinders. The relative size of the TOs increases slowly as the system becomes more Taylorized and we were unable to detect TOs in simulations where UC ≲ 0.07. The second velocity correlation parameter, |$U^{\prime }_{\rm C}$|⁠, which measures the geostrophic fraction of the short timescale flow, again shows the tendency for TOs to become more prevalent at lower Ekman number. This measure shows an increase from ≈45 to ≈63 per cent as E is reduced from 10−4 to 5 × 10−6. However, even at their most dominant, TOs only contribute ≈65 per cent of the total flow on short timescales suggesting that ageostrophic motions remain significant and thus may contribute to excitation of TOs via the Lorentz force. 4.2 Identification of TOs In a similar vein to Wicht & Christensen (2010) we identify TOs by structures in the azimuthal fields moving radially in s with the local Alfvén speed. In order to observe features operating on short timescales we analyse the fields with the time average removed; that is we consider |$u_\phi ^\prime$| and its spatial average relevant to the problem in hand. For each run we evaluate the quantity |$\langle \overline{\tilde{B}_s^2}\rangle$| for use in the definition of UA. Figs 3 and 4 show UA as a function of s for the two runs 6.5R2 and 6.5R5, respectively. Blue and red curves indicate a z-average over the northern and southern hemisphere, respectively whereas the black curve is an average performed over all z-space. These plots are typical for all runs with the same values of Pm so we do not present further plots of UA here. The form of UA is broadly similar in the two cases: increasing rapidly from the origin (but not identically zero at s = 0), reaching a peak at the TC (clearly located at s ≈ 0.538) and generally decreasing OTC as the equatorial region at the CMB is approached. The main difference is an increase in the magnitude of the Alfvén speed as the magnetic Prandtl number is increased. This is to be expected owing to the dependence of UA on Pm shown in (16). The only major difference in the form of UA at different magnetic Prandtl numbers is that runs with lower Pm tend to retain their peak Alfvén speed for a significant region OTC. Conversely, at higher Pm the Alfvén speed, as a function of s, decreases more or less immediately and monotonically from the TC to the CMB at the equator. Figure 3. Open in new tabDownload slide Alfvén speed, as a function of s, for the run 6.5R2. Figure 3. Open in new tabDownload slide Alfvén speed, as a function of s, for the run 6.5R2. Figure 4. Open in new tabDownload slide Alfvén speed, as a function of s, for the run 6.5R5. Figure 4. Open in new tabDownload slide Alfvén speed, as a function of s, for the run 6.5R5. In Figs 5– 9, we display colour-coded density plots of |$\langle \overline{u_\phi }\rangle ^\prime$| in ts-space for several runs. For these figures we have chosen runs from both regimes described in Section 4.1. Each of the figures contains three plots which display the different possible averaging domains ITC. The top/middle plot is for ITCN/ITCS whereas the bottom plot takes the average over the entire z-domain. Each plot contains the same data OTC. Overlaying each plot are several white curves that display trajectories that features take when travelling at the Alfvén speed, UA. Note that these curves do not have a constant gradient since the Alfvén speed is a function of s. Figure 5. Open in new tabDownload slide Azimuthal velocity, |$\langle \overline{u_\phi }\rangle ^\prime$|⁠, for the run 4R5, as a function of distance, s, from the rotation axis and time, t, in magnetic diffusion units. Figure 5. Open in new tabDownload slide Azimuthal velocity, |$\langle \overline{u_\phi }\rangle ^\prime$|⁠, for the run 4R5, as a function of distance, s, from the rotation axis and time, t, in magnetic diffusion units. Figure 6. Open in new tabDownload slide Azimuthal velocity, |$\langle \overline{u_\phi }\rangle ^\prime$|⁠, for the run 5R2. Figure 6. Open in new tabDownload slide Azimuthal velocity, |$\langle \overline{u_\phi }\rangle ^\prime$|⁠, for the run 5R2. Figure 7. Open in new tabDownload slide Azimuthal velocity, |$\langle \overline{u_\phi }\rangle ^\prime$|⁠, for the run 5R5. Figure 7. Open in new tabDownload slide Azimuthal velocity, |$\langle \overline{u_\phi }\rangle ^\prime$|⁠, for the run 5R5. Figure 8. Open in new tabDownload slide Azimuthal velocity, |$\langle \overline{u_\phi }\rangle ^\prime$|⁠, for the run 6.5R2. Figure 8. Open in new tabDownload slide Azimuthal velocity, |$\langle \overline{u_\phi }\rangle ^\prime$|⁠, for the run 6.5R2. Figure 9. Open in new tabDownload slide Azimuthal velocity, |$\langle \overline{u_\phi }\rangle ^\prime$|⁠, for the run 6.5R5. Figure 9. Open in new tabDownload slide Azimuthal velocity, |$\langle \overline{u_\phi }\rangle ^\prime$|⁠, for the run 6.5R5. The first run that we display plots for is a run with Pm = 5 and E = 10−4, which is the largest value of the Ekman number considered. We were unable to detect TOs during runs in the weak field regime at this large an Ekman number which, given the large value of |$\mathcal {T}$| and small value of UC for these runs, is not a surprising result. In Fig. 5, for run 4R5, several structures in |$\langle \overline{u_\phi }\rangle ^\prime$| can be identified as torsional waves since they follow a trajectory predicted by UA. These features appear regularly and can be seen to originate at various locations of the domain indicating that the waves can, but are not obliged to appear from the TC. Within the TC a wave propagates inwards from the TC in the northern hemisphere (at t ≃ 0.011); the only feature to do so in this run. In Figs 6 (for a weak field regime at Pr = 2) and 7 (for a strong field regime at Pr = 5) the Ekman number has been reduced by an order of magnitude compared with Fig. 5. In both sets of plots several TOs are again immediately apparent. Features in |$\langle \overline{u_\phi }\rangle ^\prime$| travel slower in the lower Pm case owing to the smaller magnetic field strength generated at lower magnetic Prandtl number. However, it is certainly notable, from the timescale on the plots alone, that waves are propagating significantly faster at lower Ekman number, as expected from (16). There is evidence of an inward propagating wave passing through the TC (at s ≈ 0.538) in Fig. 6 shortly after t = 0.002. It is clear from the top and middle plots that this wave continues to propagate in the southern hemisphere ITC but does not ITCN. At t ≃ 0.005 a second structure again appears to pass through the TC, this time in both hemispheres. Run 5R2 also has an approximately similar number of inward and outward propagating waves. Conversely, run 5R5 is dominated by two structures originating at the TC and moving radially outwards towards the equator at the CMB. Neither inwards propagating TOs nor TOs within the TC were identified in this run. When the Ekman number is reduced further to E = 5 × 10−6, for runs 6.5R2 and 6.5R5, we continue to observe faster moving waves with lower Ekman number. Other than the difference in the speed of the waves, run 6.5R2 is rather similar to run 5R2 since Fig. 8 displays several oscillations propagating both inwards and outwards as well as persistence through the TC. There are TOs propagating from the TC in run 6.5R5 as well as possible evidence of waves ITC propagating in either direction. However, several of the features highlighted with white curves in Fig. 9 will become more apparent when we apply bandpass filtering and thus we retain further discussion until Section 4.4. Figs 10 and 11 show a series of snapshots of |$\overline{u_\phi }^\prime$| in a meridional section for two runs. In the first set of snapshots, for run 5R2, we see that the azimuthal velocity is very columnar both inside and outside the TC. However, it proves difficult to see evidence of propagation of these columns either inwards or outwards. Analysis of a movie shows occasional propagation of columns but for the most part the oscillations are reminiscent of standing waves. This is to be expected because we observed from Fig. 6 that this run contains both inwards and outwards moving waves in approximately equal numbers. Therefore, it is tricky to distinguish between the two directions of travel, which only become apparent when averaging over both ϕ and z. Figure 10. Open in new tabDownload slide Series of snapshots of |$\overline{u_\phi }^\prime$| for the run 5R2. Panels from left- to right-hand side are at the following times: t = 0.0004, t = 0.0008, t = 0.0012 and t = 0.0016. Figure 10. Open in new tabDownload slide Series of snapshots of |$\overline{u_\phi }^\prime$| for the run 5R2. Panels from left- to right-hand side are at the following times: t = 0.0004, t = 0.0008, t = 0.0012 and t = 0.0016. Figure 11. Open in new tabDownload slide Series of snapshots of |$\overline{u_\phi }^\prime$| for the run 4R5. Panels from left- to right-hand side are at the following times: t = 0.009, t = 0.010, t = 0.012 and t = 0.014. Figure 11. Open in new tabDownload slide Series of snapshots of |$\overline{u_\phi }^\prime$| for the run 4R5. Panels from left- to right-hand side are at the following times: t = 0.009, t = 0.010, t = 0.012 and t = 0.014. Although the columnar structure of Fig. 11, for run 4R5, is less striking, we are able to observe features moving radially outwards. Between t = 0.009 and t = 0.010, a positive (red) structure in |$\overline{u_\phi }^\prime$| propagates towards the equator and by t = 0.012 it has dissipated at the boundary. This is shortly followed by a negative (blue) structure that at t = 0.009 resides in the centre of the region OTC but by t = 0.014 has moved to the equator as a newly formed positive structure now dominates OTC. These outwards propagating positive and negative features can be directly matched with those of Fig. 5 for the section of time from t = 0.009 to t = 0.014. The plots displayed, and more generally the runs considered, in this section are representative of other runs from Table 1 that are in neighbouring regions of parameter space. The general features observed in the figures can be extrapolated to the runs for which we have not displayed plots. For example, runs with Pm = 1 are found to have an even more columnar structure with even fewer propagating waves compared with the Pm = 2 cases. Additionally, we find that repeating runs with stress-free boundary conditions do not appear to alter our findings from the rigid case since various plots of the data for the runs 5F1 to 5F5 broadly match those of runs 5R1 to 5R5. This is, perhaps, not surprising when reflecting on the similarity of the output parameters from these two sets of runs (Table 1). One feature of TOs that we have not observed is the possible reflection of waves at the equator. This is true not only for the runs for which we have displayed plots, but, more generally, is the case across all of our simulations. Our results are therefore in agreement with Schaeffer et al. (2012) who suggest that the observation of wave reflection in dynamo simulations with insulating no-slip boundary conditions is not possible due to a small reflection coefficient. However, Schaeffer et al. (2012) also highlight the crucial role of the magnetic Prandtl number and, in particular, the special case of Pm = 1, where total reflection or absorption is possible depending on the kinematic boundary conditions. We were unable to identify reflection at the CMB even in the special case of stress-free boundaries and Pm = 1 (run 5F1). 4.3 Core traveltimes We are able to estimate the traveltime for our observed waves to cross the outer core. However, such estimates must be treated with a degree of caution since the parameter regimes used to produce these simulations are not that of the Earth resulting in a difficulty in identifying the timescale to use when converting back from our non-dimensional time to physical time. Consequently, there are various ways of re-scaling the dimensionless results to Earth-like values and, as discussed in Section 3.3, we have converted to dimensional time by matching the Alfvén speed at the CMB resulting in the values of τE given in Table 1. Since TOs are operating on timescales less than τE, we find that the outer core traveltime of the TOs in our simulations ranges from ≈1 to ≈15 yr. TOs in the core are currently believed (Gillet et al.2010; Wicht & Christensen 2010) to operate on a 4–10-yr timescale and, from our set of simulations, the runs in the strong field regime fare best at operating on or near to this timescale. In particular, runs 4R5, 5R3, 5R4, 5R5, 6.5R4 and 6.5R5 have all shown TOs with core crossing traveltimes in the 4–10-yr range. TOs also make up a larger component of the overall flow in the strong field regime as evidenced by larger values of UC and |$U^{\prime }_{\rm C}$|⁠. 4.4 Bandpass filtering In order to observe TOs more clearly in our simulation data, we perform bandpass filtering on our ts-data from Section 4.2. Hence, we perform a Fourier transform on our velocity in the t-direction to acquire data in frequency-space defined on frequency points, fn, given by fn = n/τ for n ≥ 0. We then filter out certain modes, fn, using a step function and inverse transform the data back to obtain the velocity that operates on timescales corresponding to the unfiltered frequencies only. This is a similar analysis to that performed by Gillet et al. (2010) albeit on our synthetic data rather than observational data. Figs 12– 15 show the azimuthal velocity in ts-space, filtered of certain frequencies, for several of our simulations. The plots in each figure follow the same layout as previous figures so from top to bottom: data for ITCN, ITCS and the average over the entire z-average, respectively. In all of our runs, we find that filtering out higher frequencies allows us to better identify the TOs in our data. Fig. 12, for run 5R5, further highlights the two TOs that were identified in this data previously (cf. Fig. 7). These data have been filtered of frequency modes above f4 (as well as the mean) and thus corresponds to timescales between 3.749 and 7.499 yr. If we instead filter these low frequency modes out of the data, we remove the structures travelling at the correct Alfvén speed. This is notable in Fig. 13, again for run 5R5, where all but frequency modes f6 to f8 (corresponding to timescales between 1.875 and 2.500 yr) are filtered. The structures present in |$\langle \overline{u_\phi }\rangle$| no longer follow the trajectories given by the white curves and instead move outwards at a faster rate. Figure 12. Open in new tabDownload slide |$\langle \overline{u_\phi }\rangle$| bandpass filtered over modes f2 to f4, for the run 5R5. Figure 12. Open in new tabDownload slide |$\langle \overline{u_\phi }\rangle$| bandpass filtered over modes f2 to f4, for the run 5R5. Figure 13. Open in new tabDownload slide |$\langle \overline{u_\phi }\rangle$| bandpass filtered over modes f6 to f8, for the run 5R5. Figure 13. Open in new tabDownload slide |$\langle \overline{u_\phi }\rangle$| bandpass filtered over modes f6 to f8, for the run 5R5. Figure 14. Open in new tabDownload slide |$\langle \overline{u_\phi }\rangle$| bandpass filtered over modes f2 to f4, for the run 6.5R2. Figure 14. Open in new tabDownload slide |$\langle \overline{u_\phi }\rangle$| bandpass filtered over modes f2 to f4, for the run 6.5R2. Figure 15. Open in new tabDownload slide |$\langle \overline{u_\phi }\rangle$| bandpass filtered over modes f2 to f4, for the run 6.5R5. Figure 15. Open in new tabDownload slide |$\langle \overline{u_\phi }\rangle$| bandpass filtered over modes f2 to f4, for the run 6.5R5. Further bandpass filtered plots for |$\langle \overline{u_\phi }\rangle$|⁠, also over the frequency modes f2 to f4, for runs 6.5R2 and 6.5R5 are presented in Figs 14 and 15, respectively. These modes correspond to timescales between 0.891 and 1.782 yr and 4.912 and 9.823 yr for 6.5R2 and 6.5R5, respectively. We have omitted plots filtered of higher frequencies for runs 6.5R2 and 6.5R5 due to their similarity to the plots of Fig. 13. All data filtered over ranges other than approximately modes f2 to f4 only show structures moving at rates inconsistent with the TO Alfvén speed. Fig. 14 allows us to identify a complicated structure of inwards and outwards propagating waves OTC near to the TC, which was not immediately obvious in the earlier unfiltered plots (cf. Fig. 8). It is clear that some inwards moving waves propagate through the TC and often into one hemisphere only. For example, the earliest instance of an inwards propagating wave in Fig. 14 reaches the TC at t ≃ 0.0006 and passes through into the region ITCS but not in the northern hemisphere. Filtering all but low frequency structures again highlights the previously identified TOs in Fig. 15, for run 6.5R5 (cf. Fig. 9). In fact, several of the features previously identified have only become clear upon filtering. We can clearly see the structures propagating outwards from the TC, as well as inwards from the TC in the northern hemisphere. Conversely, the structures ITC in the southern hemisphere propagate outwards and through into the region OTC. This run, in particular, highlights the complicated nature of waves incident on the TC. The sensitivity in the bandpass filtering and preference for low frequency modes draws our attention to two points. First, it validates our choice of τ for each run since TOs appearing at low frequencies implies that they do indeed operate approximately on the τ timescale. Secondly, the lack of TOs appearing at higher frequencies also suggests that TOs do not operate on timescale much smaller than τ. This was not immediately obvious from our unfiltered data. 4.5 Excitation mechanisms We now explore the role various forces have in the driving of the torsional waves observed in Sections 4.2 and 4.4. In Section 3.1, we discussed how there were three possible driving forces in our system and hence we plot quantities appearing on the right-hand side of (19). Since we aim to find correlation between these forcing terms and the origins of TOs we retain, on our plots throughout this section, the white curves from the associated azimuthal velocity plots of Sections 4.2 and 4.4. However, in our ts-contour plots for FL, FV and FLD, we do not expect features to be travelling along the white curves; rather we expect to find features at the origins of the curves. From Fig. 16, displaying forcing terms for run 5R5 (for the regions OTC and ITCS only), we can make several observations. All three forces are weak for most of the region OTC except at the TC itself. The viscous dissipation and the Lorentz forcing are also strong at the equator, where the rapid changes in velocity due to the CMB boundary layer have a significant effect. Within the TC all three forces, but especially FV and FLD, are larger. However, one of the most striking features of these plots in the context of TO driving is the excellent correlation between large Reynolds force at the TC and the excitation of waves represented by the origin of the two curves. Although the Reynolds force is clearly weaker than the Lorentz forcing (by approximately a factor of three), its correlation is superior since there are regions of large Lorentz force that do not coincide with TO initiation. Conversely, whenever the Reynolds force is large at the TC, a TO is produced. Figure 16. Open in new tabDownload slide Forcing terms for ITCS and OTC for the run 5R5. From top to bottom: FR, FV and FLD. Figure 16. Open in new tabDownload slide Forcing terms for ITCS and OTC for the run 5R5. From top to bottom: FR, FV and FLD. In Fig. 17, we again plot forcing terms, this time for run 6.5R5. The plots for the three forces are broadly similar to the 5R5 case OTC. Once again the locations of the origin of identified TOs are well correlated with large regions of Reynolds force, this time ITC. A lack of correlation of large FR at the TC with the waves propagating outwards there suggests that the waves ITC do indeed traverse the TC and thus do not require an excitation mechanism at the TC in this case. Evidence for correlation between Reynolds forcing and TO excitation comes not only from Figs 16 and 17, but from a series of snapshots from our runs, too numerous to display here. Figure 17. Open in new tabDownload slide Forcing terms for ITCS and OTC for the run 6.5R5. From top to bottom: FR, FV and FLD. Figure 17. Open in new tabDownload slide Forcing terms for ITCS and OTC for the run 6.5R5. From top to bottom: FR, FV and FLD. 5 DISCUSSION Through our numerical simulations we have observed TOs at a range of Ekman numbers including the relatively large E = 10−4. These oscillations are able to propagate either inwards or outwards in the cylindrical radial direction. The torsional waves travel fastest under parameter regimes that promote the production of strong magnetic fields. Thus, large magnetic Prandtl number and rapidly rotating regimes produce the quickest oscillations. TOs are often found to propagate from the TC, both inwards and outwards. Hence, we have observed waves ITC, a region of the spherical shell not considered in previous work. Although waves are mostly found to originate at the TC, it is possible for excitation to occur at other locations in the shell. This indicates a complicated non-uniform excitation mechanism with various processes likely to excite oscillations at the different locations. Within our set of simulations, we identified two dynamo regimes for which a given system is able to organize itself. Whether the dynamo is in a weak or strong field regime has implications on the torsional waves observed. Weak field regimes found at Pm ∈ [1, 3] for a range of Ekman numbers are able to produce approximately equal numbers of inward and outward propagating waves. Conversely, strong field regimes found at Pm ∈ [3, 5] are dominated by waves of outwards propagation. Plots (and movies) of meridional sections of |$\overline{u_\phi }^\prime$| are able to show the outwards propagation of columns in strong field runs whereas the same graphics show features more reminiscent of standing waves in the weak field runs. The speed of waves is found to best match that predicted for the Earth in the strong field regime with a core traveltime of between 4 and 10 yr. This may not be surprising, however, given that we believe the strong field regime to be closer to that of the Earth. Oscillations observed ITC almost exclusively originate at the TC and thus move radially inwards. This is either via an excitation mechanism at the TC or by a wave propagating across the TC from OTC. Additionally, weak field regimes are more likely to promote TOs within the TC. If waves are being excited at the TC then the weak field regime, with its greater ability to promote inwards propagation, is naturally preferred for disturbances ITC. Conversely, the preference for outwards movement in the strong field regime leads to disturbances at the TC commonly travelling through the region OTC towards the equator. We are unable to observe a consistent periodic appearance of waves in our simulations, even over the relatively short time periods we have considered; the maximum number of consecutive waves is approximately three. However, the oscillation period of waves depends on the magnitude of the damping in the system (as well as Bs) and our highly forced, highly damped system may not permit a periodic nature of waves similar to that seen in geomagnetic data (Gillet et al.2010). One of the most intriguing results from our simulations is the apparent ability of waves to cross the TC, a phenomenon also discussed by Jault & Légaut (2005). Waves can cross in either direction; however, waves entering the region ITC often dissipate quickly, probably owing to the large viscous dissipation there. Features propagating from OTC are often absorbed into only one hemisphere ITC suggesting that conditions and flow patterns have to be desirable, in a given hemisphere, for a crossing of the TC to take place in this direction. The crossing of waves in the opposite direction is possible but rarer. The likelihood of movement of oscillations into the region OTC is increased if waves are found to be approaching the TC in each hemisphere approximately concurrently. Since the regions north and south of the inner core effectively act independently, propagation from ITC to OTC is a random and often infrequent phenomenon resulting in the scarcity of such events. One of our most studied simulations (6.5R5) was one of the few to display propagation of waves from ITC to OTC. We have been able to investigate the excitation mechanisms of torsional waves within our simulations. We split these into three categories, the damping due to viscous forces, the Reynolds forces and the Lorentz forces. We have shown that the Lorentz force can be usefully divided into that part which gives the restoring force of the torsional oscillation itself, and the part that comes from the ageostrophic convection. Although the convection is relatively small-scale, the Lorentz force it produces does not vanish when averaged over the Taylor cylinder, and may be an important excitation mechanism for TOs. Despite the Reynolds force consistently being the weakest of the three forces, correlation with TO propagation from the TC leads us to conclude that it is also an important excitation mechanism in our simulations. At lower, more Earth-like, Ekman numbers the Reynolds forcing will inevitably become small relative to the Lorentz force and may play a diminished role. However, the thin region near the TC may well become thinner at low Ekman number, so the velocity gradients driving the Reynolds force might be sufficient to have an effect even though the velocity itself is small in magnitude. Indeed, Livermore & Hollerbach (2012) recently suggested excitation at the TC may occur via a discontinuity at the TC being approached as the Ekman number tends to zero. The scaling of the relative size of the Reynolds and Lorentz contribution with Ekman number needs to be explored further, but this will require a new approach, as reducing the Ekman number is notoriously expensive in full geodynamo simulations. The Lorentz force excited by ageostrophic convection, which seems particularly strong inside the TC, is currently the preferred explanation of TO excitation in the Earth's core. Viscous forces were found to be significant near the CMB equator and inside the TC in our models, though we expect their impact to be much reduced at the very low Ekman numbers of the Earth's core. Their damping effect may be replaced by electromagnetic coupling with the mantle and the inner core, which has not yet been included in our model. Several of the observations from our results highlight a common problem in numerical geodynamo simulations: we are restricted by limited computing resources when attempting to reach a parameter regime that can quantitatively replicate many of the geodynamo's features, including TOs. A reduction of geometric complexity by considering, for example, magnetoconvection in an annulus would help to alleviate this problem by allowing one to perform simulations at more realistic Ekman numbers. Alternatively, spherical geometry could be retained and a lower Ekman number achieved by performing simulations of magnetoconvection where the magnetic field strength can be systematically adjusted and the requirement of a long period of time integration to ensure a dynamo state is found is not necessary. Additionally, our current study omits electromagnetic coupling to the inner core via conducting magnetic boundary conditions on the ICB, which may play a role in torsional wave dynamics. These topics are the subject of future work. This work was supported by the Natural Environment Research Council, grant NE/I012052/1. REFERENCES Alfvén H. . Existence of electromagnetic-hydrodynamic waves , Nature , 1942 , vol. 150 (pg. 405 - 406 ) Google Scholar Crossref Search ADS WorldCat Bloxham J. , Zatman S. , Dumberry M. . The origin of geomagnetic jerks , Nature , 2002 , vol. 420 (pg. 65 - 68 ) Google Scholar Crossref Search ADS PubMed WorldCat Braginsky S. . Torsional magnetohydrodynamic vibrations in the Earth's core and variation in day length , Geomag. Aeron. , 1970 , vol. 10 (pg. 1 - 8 ) Google Scholar OpenURL Placeholder Text WorldCat Braginsky S. . Short-period geomagnetic secular variation , Geophys. astrophys. Fluid Dyn. , 1984 , vol. 30 (pg. 1 - 78 ) Google Scholar Crossref Search ADS WorldCat Buffett B. , Mound J. , Jackson A. . Inversion of torsional oscillations for the structure and dynamics of Earth's core , Geophys. J. Int. , 2009 , vol. 177 (pg. 878 - 890 ) Google Scholar Crossref Search ADS WorldCat Busse F. , Simitev R. . Soward A. , Jones C. , Hughes D. , Weiss N. . Convection in rotating spherical fluid shells and its dynamo states , Mathematical Aspects of Natural Dynamos , 2005 CRC Press (pg. 359 - 392 ) Google Scholar OpenURL Placeholder Text WorldCat Christensen U. . Dynamo scaling laws and applications to the planets , Space Sci. Rev. , 2010 , vol. 152 (pg. 565 - 590 ) Google Scholar Crossref Search ADS WorldCat Christensen U. , Aubert J. , Hulot G. . Conditions for Earth-like geodynamo models , Earth planet. Sci. Lett. , 2010 , vol. 296 (pg. 487 - 496 ) Google Scholar Crossref Search ADS WorldCat Cox G. , Livermore P. , Mound J. . , Geophys. J. Int. , 2013 Forward models of torsional waves: dispersion and geometric effects doi: 10.1093/gji/ggt414 Dormy E. , Soward A. , Jones C. , Jault D. , Cardin P. . The onset of thermal convection in rotating spherical shells , J. Fluid Mech. , 2004 , vol. 501 (pg. 43 - 70 ) Google Scholar Crossref Search ADS WorldCat Dumberry M. , Bloxham J. . Torque balance, Taylor's constrait and torsional oscillations in a numerical model of the geodynamo , Phys. Earth planet. Inter. , 2003 , vol. 140 (pg. 29 - 51 ) Google Scholar Crossref Search ADS WorldCat Gillet N. , Jault D. , Canet E. , Fournier A. . Fast torsional waves and strong magnetic field within the Earth's core , Nature , 2010 , vol. 465 (pg. 74 - 77 ) Google Scholar Crossref Search ADS PubMed WorldCat Gubbins D. , Bloxham J. . Geomagnetic field analysis. Part III. Magnetic fields on the core-mantle boundary , Geophys. J. R. astr. Soc. , 1985 , vol. 80 (pg. 695 - 713 ) Google Scholar Crossref Search ADS WorldCat Hori K. , Wicht J. , Christensen U. . The effect of thermal boundary conditions on dynamos driven by internal heating , Phys. Earth planet. Inter. , 2010 , vol. 182 (pg. 85 - 97 ) Google Scholar Crossref Search ADS WorldCat Jackson A. . Time dependence of geostrophic core-surface motions , Phys. Earth planet. Inter. , 1997 , vol. 103 (pg. 293 - 311 ) Google Scholar Crossref Search ADS WorldCat Jault D. , Gire C. , LeMouel J.-L. . Westward drift, core motion and exchanges of angular momentum between core and mantle , Nature , 1988 , vol. 333 (pg. 353 - 356 ) Google Scholar Crossref Search ADS WorldCat Jault D. , Légaut G. . Soward A.M. , Jones C.A. , Hughes D.W. , Weiss N.O. . Alfvén waves within the Earth's core , Fluid Dynamics and Dynamos in Astrophysics and Geophysics , 2005 CRC Press (pg. 277 - 294 ) chap. 9 Google Scholar OpenURL Placeholder Text WorldCat Jones C. . Planetary magnetic fields and fluid dynamos , Ann. Rev. Fluid Mech. , 2011 , vol. 43 (pg. 583 - 614 ) Google Scholar Crossref Search ADS WorldCat Jones C. , Boronski P. , Brun A. , Glatzmaier G. , Gastine T. , Miesch M. , Wicht J. . Anelastic convection-driven dynamo benchmarks , Icarus , 2011 , vol. 216 (pg. 120 - 135 ) Google Scholar Crossref Search ADS WorldCat Légaut G. . Ondes de torsion dans le noyau terrestre , PhD thesis , 2005 Grenoble, France Univ. Joseph Fourier Google Scholar OpenURL Placeholder Text WorldCat Livermore P.W. , Hollerbach R. . Successive elimination of shear layers by a hierarchy of constraints in inviscid spherical-shell flows , J. Math. Phys. , 2012 , vol. 53 pg. 073104 Google Scholar Crossref Search ADS WorldCat Mound J. , Buffett B. . Detection of a gravitational oscillation in length-of-day , Earth planet. Sci. Lett. , 2006 , vol. 243 (pg. 383 - 389 ) Google Scholar Crossref Search ADS WorldCat Roberts P. , Aurnou J. . On the theory of core-mantle coupling , Geophys. astrophys. Fluid Dyn. , 2012 , vol. 106 2 (pg. 157 - 230 ) Google Scholar Crossref Search ADS WorldCat Sakuraba A. , Roberts P. . Generation of a strong magnetic field using uniform heat flux at the surface of the core , Nature Geosci. , 2009 , vol. 2 (pg. 802 - 805 ) Google Scholar Crossref Search ADS WorldCat Schaeffer N. , Jault D. , Cardin P. , Drouard M. . On the reflection of Alfvén waves and its implication for Earth's core modelling , Geophys. J. Int. , 2012 , vol. 191 (pg. 508 - 516 ) Google Scholar Crossref Search ADS WorldCat Taylor J. . The magneto-hydrodynamics of a rotating fluid and the Earth's dynamo problem , Proc. Roy. Soc. A , 1963 , vol. 274 (pg. 274 - 283 ) Google Scholar Crossref Search ADS WorldCat Wicht J. , Christensen U. . Torsional oscillations in dynamo simulations , Geophys. J. Int. , 2010 , vol. 181 (pg. 1367 - 1380 ) Google Scholar OpenURL Placeholder Text WorldCat Zatman S. , Bloxham J. . Torsional oscillations and the magnetic field with the Earth's core , Nature , 1997 , vol. 388 (pg. 760 - 761 ) Google Scholar Crossref Search ADS WorldCat © The Author 2013. Published by Oxford University Press on behalf of The Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. © The Author 2013. Published by Oxford University Press on behalf of The Royal Astronomical Society. TI - The dynamics and excitation of torsional waves in geodynamo simulations JF - Geophysical Journal International DO - 10.1093/gji/ggt432 DA - 2014-02-01 UR - https://www.deepdyve.com/lp/oxford-university-press/the-dynamics-and-excitation-of-torsional-waves-in-geodynamo-t0ReLI0Qgd SP - 724 EP - 735 VL - 196 IS - 2 DP - DeepDyve ER -