Contents lists available at ScienceDirect

Earth and Planetary Science Letters

www.elsevier.com/locate/epsl

The transition to Earth-like torsional oscillations in magnetoconvection simulations

Robert J. Teed *, Chris A. Jones, Steven M. Tobias

Department of Applied Mathematics, University of Leeds, Leeds, LS2 9JT, UK

A R T I C L E I N F 0 A B S T R A C T

Evidence for torsional oscillations (TOs) operating within the Earth's fluid outer core has been found in the secular variation of the geomagnetic field. These waves arise via disturbances to the predominant (magnetostrophic) force balance believed to exist in the core. The coupling of the core and mantle allow TOs to affect the length-of-day of the Earth via angular momentum conservation.

Encouraged by previous work, where we were able to observe TOs in geodynamo simulations, we perform 3-D magnetoconvection simulations in a spherical shell in order to reach more Earth-like parameter regimes that proved hitherto elusive.

At large Ekman numbers we find that TOs can be present but are typically only a small fraction of the overall dynamics and are often driven by Reynolds forcing at various locations throughout the domain. However, as the Ekman number is reduced to more Earth-like values, TOs become more apparent and can make up the dominant portion of the short timescale flow. This coincides with a transition to regimes where excitation is found only at the tangent cylinder, is delivered by the Lorentz force and gives rise to a periodic Earth-like wave pattern, approximately operating on a 4 to 5 year timescale. The core travel times of our waves also become independent of rotation at low Ekman number with many converging to Earth-like values of around 4 years.

© 2015 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license

(http://creativecommons.org/licenses/by/4.0/).

CrossMark

Article history:

Received 2 September 2014

Received in revised form 23 February 2015

Accepted 28 February 2015

Available online 19 March 2015

Editor: C. Sotin

Keywords:

torsional oscillation Taylor state geodynamo

1. Introduction

The geodynamo, which operates in the Earth's iron-rich outer core and continuously replenishes its magnetic field, is believed to be one of a number of planetary dynamos that exists under a quasi-magnetostrophic regime (see, for example, Jones et al., 2011). Specifically, this suggests that, owing to the rapidly rotating nature of the Earth, the predominant force balance within the core is between Lorentz, Archimedean and Coriolis forces; this is commonly known as the MAC balance. One consequence of the magnetostrophic regime is that the averaged azimuthal Lorentz force is compelled to vanish on cylinders aligned with the rotation axis, leaving the system in a Taylor state (Taylor, 1963). Violations of Taylor's constraint manifest themselves as the acceleration of concentric cylinders with a restoring Lorentz force acting like a torsional spring that attempts to reimpose the Taylor state (Braginsky, 1970). This process ultimately leads to the excitation of torsional

* Corresponding author at: Department of Applied Mathematics and Theoretical Physics, University of Cambridge, Cambridge, CB3 0WA, UK.

E-mail addresses: R.J.Teed@damtp.cam.ac.uk (R.J. Teed), C.A.Jones@maths.leeds.ac.uk (C.A. Jones).

oscillations (TOs), which are the only Alfven waves (Alfven, 1942) in the core that act at large lengthscales, propagating in the cylindrical radial direction.

Theory has been complemented, in recent years, by improved observational data although there have been conflicting estimates of the timescale on which TOs operate. Braginsky (1984) first suggested a 60-year timescale, however, this has now been superseded by recent core flow models (Gillet et al., 2010) as well as signals in length-of-day data (Holme and de Viron, 2013; Chao et al., 2014), which indicate a 6 year operational timescale. Indeed, TOs have long been presented as an explanation for certain changes in the Earth's length-of-day via the coupling of solid and fluid regions at the core-mantle boundary (CMB) by angular momentum conservation (Jault et al., 1988; Jackson, 1997; Roberts and Aurnou, 2012), and may be connected with other features such as geomagnetic jerks (Bloxham et al., 2002; Brown et al., 2013).

Direct observation of TOs is inhibited by the inability to probe the Earth's magnetic interior beneath the CMB (Gubbins and Blox-ham, 1985), though it is believed that the TO velocity field can be seen in the recent secular variation data (Gillet et al., 2010). However, further theoretical progress can be made; in particular, improved computing resources in recent years has allowed for the

http://dx.doi.org/10.1016/j.epsl.2015.02.045

0012-821X/© 2015 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).

identification of TOs in numerical simulations (Dumberry and Blox-ham, 2003; Busse and Simitev, 2005; Wicht and Christensen, 2010; Teed et al., 2014). Wicht and Christensen (2010) provided evidence of TOs operating in the region outside the tangent cylinder (OTC) and in a recent study we (from this point on referred to as TJT, Teed et al., 2014) extended this to include the region inside the tangent cylinder (ITC) where we observed waves crossing the tangent cylinder (TC). The exact excitation mechanism of TOs, both in simulations and the Earth itself, remains unclear but progress has been made in identifying the forcing terms involved (Wicht and Christensen, 2010), which are directly calculable in numerical models. Indeed the Lorentz and Reynolds forces are shown (TJT) to be important in simulations and it is likely that they can be in Earth's core also, although other excitation mechanisms are possible (Mound and Buffett, 2006; Legaut, 2005).

As is the case with all numerical models of the Earth's core, the search for TOs is currently hindered by our inability to reach the actual parameter values of the core. The ratio of viscous to Corio-lis forces is tiny in the core and it is not possible to get close to the quasi-magnetostrophic regime for dynamo calculations. Having reached the numerical limits of our code in TJT, a new approach has been necessary to continue our investigation of TOs in simulations and this is the topic of the work presented here. Very long runs are needed to establish the growth and equilibration of a dynamo driven magnetic field. Since we are now confident that numerical dynamos giving a field similar to the geomagnetic field exist in the appropriate region of parameter space, we can make significant computational savings by imposing the dipole component of the magnetic field at the CMB. The transient field state is then much shorter, and this allows us to reach more Earthlike parameters in this work. Also, magnetoconvection simulations provide the additional benefit of allowing us to vary the imposed magnetic field strength as an input parameter to allow for a systematic exploration of the dependence of TOs on field strength.

We identify a condition for the observation of TOs in the presence of overlying convection in simulations, based on an output parameter, the short time-scale geostrophy parameter, introduced in TJT. We explore the driving forces, based on the novel approach used in TJT of separating the Lorentz force into its restoring and driving components, as well as the driving locations of the waves as we move into new parameter regimes inaccessible to us hitherto. Core travel times based on the Alfven speed are also calculated.

remaining components of the field, both at r = ro and at r = rj, we use the standard form of insulating magnetic boundary conditions (for further details see the Supplementary Material).

We previously (TJT) implemented fixed flux thermal boundary conditions; this was to promote Earth-like dynamo regimes, which, under magnetoconvection, is no longer a requirement. Therefore, in a change from TJT, we return to the simpler choice of fixed temperature conditions so that convection is instigated by differential heating. We also impose rigid (no-slip) kinematic boundary conditions at both boundaries.

The system of coupled equations for velocity, u, magnetic field, B, temperature, T, and pressure, p are:

d u Pm r -,

--+ (u • V)u =--[Vp + 2z x u - (VxB) x Bl

Pm2Ra 2

+--T r + PmV 2u,

d T Pm 2

— + (u • V)T = —V2 T, dt ( ) Pr

--V x (u x B) = V2B,

d t ( ) '

V • u = 0, V-B = 0,

(1) (2)

where we have nondimensionalised using length scale, D = ro — ri, magnetic timescale, D2/n, temperature scale, AT, and magnetic scale, VP^n.

The nondimensional parameters appearing in our equations are the Rayleigh number, Ra, Ekman number, E, Prandtl number, Pr, and magnetic Prandtl number, Pm, defined by:

gaATD 3

We take the value of the aspect ratio of the spherical shell, given by ri/ro, to be 0.35 throughout, which is apt for a model of the Earth's core.

2.2. TO theory

The observation of TOs requires the definition of various spatial and temporal averages in a cylindrical coordinate system: (s, 0, z). We continue the notation from TJT so, for any scalar field, A(t, s, 0, z), we define

2. Mathematical formulation and methods

2.1. Mathematical model

The model we employ is adapted for magnetoconvection from

that used in TJT where dynamo runs were performed. Hence, using

a spherical coordinate system (r, 0, 0), we consider a fluid filled spherical shell that rotates about the vertical with rotation rate Q (fi = Qz). Gravity acts radially inward (g = —gr) and the fluid is modelled using the Boussinesq approximation with constant values of p, v, k and n, the outer core density, kinematic viscosity, thermal diffusivity and magnetic diffusivity respectively.

The geomagnetic field is frequently decomposed into spherical harmonics, and in our numerical code the magnetic field is separated into poloidal and toroidal parts with each part being split into spherical harmonics (see, for example, Roberts, 2007). At an insulating boundary there is a matching condition between each poloidal component and its normal derivative so it is therefore convenient to implement magnetoconvection by setting the amplitude of the axial dipole component, Y^, so that this component gives B0 = (2B0 cos0, B0 sin0, 0) on the CMB at r = ro. On all the

<A)(t, s,$) = 1 h

—- *+ fAdZ + i<

A(t, s, z) = — [ AdS and 2n J 0

i —z— z+ \

Adz , (7)

Vz+ z— /

for spatial averages over 0 and z, where z+ = y7ri — s2, z— =

^ yV2 — s2 , h = 2(z+ — z—). A complication arose at this stage

in TJT, where multiple definitions of the z-average were required to account for the separation of hemispheres by the inner core and the inhomogeneous nature of the fields in the dynamo runs, when s < ri. Although the current geometry maintains this separation, we find that in practice the imposition of a dipolar magnetic field on the outer boundary creates a restriction so that the dynamics in the two hemispheres are very similar. The requirement for different definitions of the z-average OTC and ITC is therefore lifted and our z-averages in this work are calculated across the entire domain for all values of s. In TJT we found it helpful to introduce

a time t , longer than the TO travel time, over which we can average and distinguish the oscillating parts from the steady parts of the flow, so we define:

A(s,0,z) = -/ Adt and A'(t, s,0,z) = A — A,

to indicate the (temporal) mean and fluctuating parts of A respectively, and

A 11 =

A dv dt

to denote the volume average over the whole outer core, V c being the outer core volume, averaged over the whole time interval of selected run-time.

TOs oscillate in the azimuthal direction and act on cylinders and thus we consider the 0 and z-average of the 0-component of (1), noting that the Coriolis and buoyancy forces vanish during integration, to leave an equation describing the acceleration of geostrophic cylinders:

d (U0) dt

= — (0 ■ (V- uu)) + PmE—1 (0 ■ ((V X B) X B))

+ Pm(0 ■ V2u). = Fr + Fl + Fv .

As previously discussed by Wicht and Christensen (2010), the three forces on the right-hand-side of (10) can be identified as the Reynolds force, FR, the Lorentz force, FL, and the viscous force, Fv . A system operating under magnetostrophic balance (FL ^ FR, FV), subject to (10) then results in Taylor's constraint (Taylor, 1963), Fl = 0, when averaged over long timescales. Deviations to, and the reinforcement of, the Taylor state produce azimuthal oscillations of concentric cylinders as a torsional wave propagating radially; Roberts and Wu (2014) have recently investigated what they refer to as the 'modified Taylor state' where inertia, which leads to TOs, is retained in the derivation.

When analysing TOs, the general procedure is to split the velocity and magnetic field into mean and fluctuating parts in order to separate the timescale that TOs operate on from that of the convection. Some earlier work (Wicht and Christensen, 2010; Roberts and Aurnou, 2012) makes the assumption that the mean quantities are the principal parts of the Taylor state, about which TOs operate, and that the fluctuating quantities are perturbations giving rise to TOs. However, this also makes the implicit assumption that the fluctuating part of the flow is exclusively geostrophic (as noted by Taylor, 1963) since TOs operate on geostrophic cylinders. In TJT, we argued that the fluctuating part of the flow should also include an ageostrophic component since in the Earth's core, as well as in our numerical models, convection will operate, to some degree, on all timescales. This ageostrophic fluctuating velocity component could contribute to excitation and we therefore split the fluctuating part of u into separate geostrophic and ageostrophic parts by writing

u = u + u = u + sZ'(s, t)4 + ua

B = B + B',

where u is the mean part of the velocity and sZ' and ua are the geostrophic and ageostrophic parts of the fluctuating velocity, u', respectively. Additionally B and B' are the mean and fluctuating parts of the magnetic field respectively. By splitting the velocity and magnetic fields in this way we are able (see TJT) to identify the part of the Lorentz force that can act as a driving force for

TOs (separately from the restoring force) by writing FL = FLR + FLD where

flr = —- ( s3huj

1 d hs2 ds

ua = J m (B 2 K

is the time derivative of the restoring Lorentz force and U A is the Alfven speed. With the canonical torsional wave equation (Braginsky, 1970) then given by s\' = FLR we found that (10) could be written as

s'z' — Flr = Fld + Fr + Fv ; (13)

see Section 3.1 of TJT for the full details.

2.3. Methods

In this work we perform many magnetoconvection simulations using the Leeds spherical dynamo code (Willis et al., 2007; Jones et al., 2011) with the mathematical set-up described in Section 2.1. We cover a broad range of parameter space in E, Pm and B0, whilst keeping the Prandtl number and the level of supercriti-cality constant for all runs: Pr = 1 and Ra = 5Rac , with the values of Rac used being those appropriate for the onset of convection in the non-magnetic case (see, for example, Dormy et al., 2004).

Simulations are run from an initial random state, past any transient and allowed to reach statistical equilibrium before being analysed for evidence of TOs over the time, t . Further information, including a complete table of all runs performed can be found in the Supplementary Material.

2.4. Output parameters

We output several quantities from our simulations. The magnetic Reynolds number, Elsasser number and Rossby number are defined by

A =11 B2

respectively. In these definitions, the velocity and magnetic field are the nondimensional quantities and hence Rm and A can be thought of as the nondimensional velocity and square of the magnetic field because of our choice of nondimensionalisation. We additionally define further output parameters useful in our study:

\ut ■ ut

II u - u II

Roe = Ro — . n

(us )2 — (us )2 Rm2 — Rm2

E DUa (ro)

Tc = UE Tc. U A

The local Rossby number, Rof , is based on a local convective length scale (the characteristic harmonic degree of the fluid velocity, lu) rather than the global length scale and its definition here is consistent with that of previous studies such as Christensen and Aubert (2006). The geostrophy parameter, U'c, introduced in TJT, is the ratio of the root mean square fluctuating geostrophic flow to the root mean square total fluctuating flow, so it is a measure of the geostrophy of the short timescale flow and therefore indicates the dominance of TOs in the flow. The precise definition of the geostrophy parameter depends on t , but we found that provided t is long enough, as in Table A.1 in the Supplementary Material,

1 000 r

800 ■

600 ■

400 ■

200 ■

v Pm= 0.05

Pm= 0.1

A Pm= 0.5

□ Pm= i

O Pm= 5

■ o O

100 ■

50 ■

• v Pm= 0.05 No B0

• •:■ Pm= 0.1 ♦ B0= 0.1

■ 0.08 - A Pm= 0.5 0.8 •• J Ba=0.5

• □ Pm—2 0 1 " A. B0= 1

A □ □ * ♦ a □ ' O Pm= 5 B0= 10

No Ba B0=0.1 S0=0.5 B0= 1

t □ .

J ■ 5 o

■iJl a i

10"7 10"6 10"5 10"+ 10"3 1

(a) Magnetic Reynolds number, Rm

10"7 10"6 10"5 10"+ 10"3 1

(b) Elsasser number, A

Fig. 1. Parameter plots of output parameters against magnetic Ekman number encompassing all simulations. Filled points indicate runs for which TOs were observed. Symbol shapes denote Pm values and colours denote B0 values. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

the geostrophy parameter becomes insensitive to t . The quantity Uc measures the ratio of the root mean square fluctuating geostrophic flow to the total flow. This quantity is not quite so useful as U'c, but it is easier to estimate in the core, because the fluctuating geostrophic flow can be estimated from length-of-day changes. Gillet et al. (2010) give 0.4 kmyr-1 as its amplitude over the 6 year period currently expected for TOs, and the core flow estimated from secular variation studies is typically 15 kmyr-1. However, over decadal periods of time the fluctuating flow may be as large as 5 kmyr-1 (Zatman and Bloxham, 1997) and could be even larger with longer data samples. Thus when the fluctuating part is defined over a short enough timescale (for example, in this study of TOs), it is likely that only a small fraction of the core flow resides in the fluctuating part, so U'c may be significantly larger than Uc as we find in our simulations. Finally, tce is the dimensional (outer) core travel time of TOs in years and is calculated from its dimensionless counterpart, tc , which in turn is calculated using the Alfven speed, UA (s). In order to convert between the dimensionless and dimensional versions of the core travel time we have matched the Alfven speed at the outer boundary in our simulation, UA(ro), to that of the Earth at the CMB, UA & 1.784 x 10-3 ms-1, which is appropriate given the focus of this study. This conversion was used in TJT and explains the factors appearing in (17).

3. Numerical results

In this section we discuss the numerical results from our simulations by way of a series of plots against the magnetic Ekman number: En = n/QD2 = E/Pm. Each point in each plot represents a time averaged quantity from a single simulation (see Table A.1 for the numerical values) and the same symbol and colour coding is used for each plot. Different symbols indicate various values of the magnetic Prandtl number whereas points are colour-coded by B0 strength. Note that black points represent dynamo runs from TJT, which are included for reference only. In addition, filled points represent simulations for which we were able to identify TOs and the converse is true of the unfilled points.

3.1. Dynamics

The magnetic Reynolds number increases as the magnetic Ek-man number is reduced as indicated by Fig. 1a and this is to be

dq fj iJ:s...E

10"7 10"6 10"5 10"

10"J 10"

1 0"7 1 0"6 1 0"5 1 0"+ 1 (

(a) Local Rossby number, Roe (b) Geostrophy parameter, U'c

Fig. 2. Parameter plots of output parameters against magnetic Ekman number encompassing all simulations. Filled points indicate runs for which TOs were observed. Symbol shapes denote Pm values and colours denote B0 values. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

expected from the increased convective vigour at larger rotation rates arising from our choice of fixed supercriticality, Ra = 5Rac, rather than fixed Rayleigh number. Points are grouped neatly by symbol with larger Rm being permissible at larger values of Pm, and in each Pm-case a stronger imposed field lowers Rm. Runs where TOs were not observed (indicated by unfilled points) tend to be found at low values of Rm, perhaps indicating that the convec-tive driving is then too small to overcome the ohmic dissipation. There does not, however, appear to be a clear critical value of Rm below which TOs can never be found since there is evidence of the waves at quite low Rm so long as the magnetic Ekman number is also sufficiently small.

Fig. 1b shows an increase in magnetic field strength with reduced En although points with the same value of the magnetic Prandtl number are not grouped together as clearly in this case. It is not surprising to see that A also increases with B0 as more magnetic energy is imposed; this is most clearly evidenced by the grouping of blue points at large values of the Elsasser number.

Unsurprisingly, the local Rossby number, displayed in Fig. 2a, is small and is found to decrease as the magnetic Ekman number is reduced, since rotation dominates inertia at the parameters used in our suite of runs. However, the local Rossby number does provide insight in two areas. Firstly, it has been shown in previous work (Olson and Christensen, 2006) that Roi provides a good indication of the transition from multipolar to dipolar solutions with the latter being found as Roi passes below & 0.1. The values of Roi displayed in Fig. 2a suggest that had our simulations been run as dynamo, rather than magnetoconvection, simulations they would fall into the dipolar (and thus more Earth-like) regime. Secondly, although the local Rossby number is small, it is consistently (approximately) one order of magnitude larger than the Rossby number, the values of which can be found in Table A.1. Therefore, even though the role of rotation is pre-eminent both globally and locally, inertia plays a more significant role on the (local) convective length scales than the (global) system length scale. This suggests that even at small magnetic Ekman numbers inertia may remain important at small enough scales and, crucially, may influence the driving of TOs at certain locations in the domain through Reynolds forcing or ageostrophic convection acting via the Lorentz force.

Table 1

Table displaying the input and output parameter sets used for five selected simulations from our complete suite. Note that all runs have Pr = 1 and Ra = 5Rac and rigid, isothermal, magnetically insulating boundary conditions.

Run E Pm Be Ra T Rm Ro Rot a Uc U c Tc T E Tc TOs

E3P5d 10-3 5 10 1.705 x 105 0.05 80.38 0.0161 0.0339 114.80 0.113 0.545 0.010 3.44 Y

E5P2a 10-5 2 0.1 6.008 x 107 0.002 253.62 0.0009 0.0039 22.40 0.125 0.551 0.001 7.87 Y

5E6P.1d 5 x 10-6 0.1 10 1.498 x 108 0.03 13.92 0.0007 0.0084 99.75 0.155 0.848 0.018 14.84 Y

E6P.5a 10-6 0.5 0.1 6.405 x 109 0.0038 193.84 0.0005 0.0110 0.47 0.163 0.772 0.005 5.23 Y

5E5P2d 5 x 10-5 2 10 3.639 x 10' 0.005 115.33 0.0029 0.0154 104.86 0.041 0.259 0.004 3.77 N

3.2. TO identification

We identify TOs as structures in )' travelling at the Alfven speed, UA. We previously found, in TJT, that TOs could be observed in dynamo simulations when the geostrophy parameter, U'c, was greater than approximately 0.4, although U'c never exceeded & 0.65 at the Ekman numbers that could be used. We see from Fig. 2b that U'c remains a good output parameter to indicate the detectability of TOs in simulations in our current work also, where we have managed to reach lower Ekman numbers. Recall that filled symbols indicate that TOs can be identified so, with unfilled points lying almost exclusively below the line U'c & 0.4, it is helpful to think of U'c = 0.4 as a critical value that must be exceeded for a given simulation to clearly exhibit TOs. The typical value of U'c in the core is thought to be approximately 0.45 (N. Gillet, private communication) and it is therefore of a similar size to values found in many of our simulations.

Fig. 2b also shows a general trend for the short timescale flow to become more geostrophic as the magnetic Ekman number is reduced and we would therefore expect TOs to be more apparent in these runs with large U'c. This trend is noticeably more apparent if the data are plotted against En rather than against E, which is why we display these results as functions of En rather than E. This indicates that provided the magnetic diffusion is low enough, En < 10-5, TOs are always observed. For comparison, in the Earth's core En & 3 x 10-9. Unsurprisingly, at high Ekman numbers only runs with large enough values of Pm or B0 are found to have TOs with the remaining simulations falling below the critical value of U'c. As E is reduced the value of U'c increases but only when the Ekman number reaches a sufficiently small value (E < 5 x 10-7) do we consistently find TOs in runs with Pm < 0.1. However, at these small Ekman numbers, numerical constraints limit the range of Pm or B0 (and, equivalently, A). If we were able to extend our parameter space to lower Ekman numbers whilst retaining large values of Pm and B0 we might expect to find even larger values of U'c. There is no consistent dependence of U'c on magnetic field strength, and so no consistent dependence on Lehnert number X = VEA/Pm. This is slightly surprising, as it is known that convection on the short timescales considered here becomes more quasi-geostrophic (i.e. more columnar) as X is reduced (Jault, 2008; Gillet et al., 2011). However, the TOs depend on magnetic field for their existence and the magnetic field plays an important role in exciting TOs, so that although reducing X reduces the fluctuating ageostrophic components of the flow, it also reduces the TO signal too. There is therefore no simple dependence of their ratio, U'c, on X.

We have performed a wide range of simulations to display evidence of TOs that we have observed though we now focus on just five specific runs (the data of which are also shown in Table 1) with contrasting TO-qualities. We plot, in Fig. 3, the (normalised) Alfven speed for the five runs and the spatial form of UA (s) shown is typical for our suite of runs although the magnitude varies with magnetic field strength. Cylindrically averaged quantities are functions of time and cylindrical radius only and in Fig. 4 we display colour-coded density plots of (u^)' against t and s. These plots are overlaid with white curves that indicate the trajectory taken by an

Fig. 3. Plot of the normalised Alfven speed, UA (s), as a function of s, for the runs of Table 1. The true maxima for each run are: E3P5d: 225.93; E5P2a: 1688.29; 5E6P.1d: 68.10; E6P.5a: 4.65; 5E5P2d: 406.33. The dashed vertical line shows the location of the ICB.

Alfven wave given by UA (s) and thus we highlight features in (u^)' with these curves.

At large Ekman numbers (E > 5 x 10-5) TOs are often not observed but when they are the geostrophy parameter remains relatively small (U'c < 0.55) indicating that ageostrophic motions remain important. This is highlighted by Fig. 4a, for run E3P5d, where waves can be seen travelling inwards but convective motions of a similar magnitude are also occurring on the same scale. As the Ekman number is reduced, and geostrophic motions form a larger component of the short timescale flow, TOs become more evident in our runs with the plots of Figs. 4b to 4d showing three such cases. In these runs the waves are clearly excited at the TC and can propagate either inwards or outwards from this location. Fig. 4c, for run 5E6P.1d, in particular, shows a striking periodic nature of waves excited at the TC and propagating outwards to the outer boundary. This plot is reminiscent of observational data of TOs (Gillet et al., 2010) although we are only able to replicate waves OTC whereas observations appear to show waves travelling in both directions from the TC. Nevertheless we believe this is the first evidence, in simulations, of waves travelling from the TC to the outer boundary at the correct Alfven speed in a periodic Earth-like pattern. The period of the oscillation in this run is approximately 0.005 dimensionless time units and this corresponds (using (17)) to a period of around 4 to 5 years, which is of the same order as the 6 year length-of-day signal found from observations (Holme and de Viron, 2013).

The lack of flow ITC may be explained by our imposition of a dipolar field which naturally creates a strong field in this region and inhibits convection. Simulations with a large B0 were found to have very little flow ITC whereas we are sometimes able to observe convection and TOs ITC in runs with a smaller imposed field strength; for Fig. 4b, for run E5P2a, displays waves propagating ITC. In our previous work (TJT), the magnetic field strength and morphology was mostly controlled by Pm allowing for TOs and convective motions to manifest within the TC as the field was then either less strong or less dipolar. The inhibition of flow ITC at low Ekman numbers and strong field strengths also has a consequence for the direction waves are able to travel since, with the location

0.0020

(a) E3P5d

(b) E5P2a

(e) 5E5P2d

Fig. 4. Azimuthal velocity, (uj)' as a function of distance, s, from the rotation axis and time, t, in magnetic diffusion units for the runs of Table 1. White curves indicate the trajectory of a particle travelling at the Alfven speed, UA(s). For plots (a) to (d) we highlight, with these white curves, features in (uj)' that travel at the Alfven speed. In plot (e) we display a case where no features in (uj)' travelled at the Alfven speed. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

of excitation being the TC at low E, waves are only able to travel outwards, unless B0 is small. We therefore see waves travelling outwards in Figs. 4c (large B0) and 4d (small E), but inwards in Fig. 4b, with all three simulations having excitation of waves at the TC. In runs at larger Ekman numbers we find that both inwards and outwards propagating waves occur even with large imposed fields because the excitation location is not restricted to the TC and thus, in Fig. 4a, waves travel inwards within the region OTC.

Fig. 4e, for run 5E5P2d, is included as an example where TOs were not present (note that this run has U'c = 0.259 < 0.4). No structures in (uj)' were found to follow the paths taken by a wave travelling at the Alfven speed and in this plot the white curves are only included to show possible examples of wave trajectories.

The reflection of waves incident on the CMB at the equator is a complicated topic as illustrated by Schaeffer et al. (2012) who found that reflection is related to the boundary conditions imposed and also the value of Pm. Under the boundary conditions used in our simulations Schaeffer et al. (2012) predict that it would only be possible to observe reflection in runs with Pm < 0.3. Nevertheless, we do not observe any clear evidence of wave reflection in our simulations, even at Pm = 0.05.

3.3. Excitation

Having identified TOs in a majority of our runs we now consider the properties of excitation of TOs in each simulation. The

driving and dissipation of waves is controlled by the terms on the right-hand-side of (13), namely Lorentz, Reynolds and viscous forcing. In TJT we were unable to identify conclusively which of these terms was ultimately responsible for the excitation of TOs. We found that, although the Reynolds force was weak, it could be locally strong and had a better correlation with the excitation locations of TOs than the larger Lorentz force. Viscous forcing was found to be small at excitation locations but could be large within the TC and on the outer boundary and thus contribute to the dissipation of waves at those locations.

In this work we aim to clarify the driving location and driving force of TOs in various parameter regimes. To do this we survey our simulations and, based on the strength of the force and its correlation with the TO excitation location, we are able, in many runs, to identify the driving force of the TOs with a high degree of confidence. The results of this process are displayed in Fig. 5a with black and white points indicating excitation by FLD and FR, respectively. In several of our runs, however, the Lorentz and Reynolds forces are found to be of a similar order of magnitude indicating that there are two contributors to the driving; these cases are displayed as grey points in Fig. 5a.

It is clear from Fig. 5a that there is a transition in the driving force as the Ekman number is reduced, or, equivalently, as the geostrophy parameter increases. In Figs. 5a and 5b we found that it is small E rather than small Ev which gives the clearest in-

0.0 10

▼ ■

♦ a g

▼ ♦ A E

v Pm- = 0.05

Pm= = 0.1

A Pm= = 0.5

□ Pm=

O Pm= = 5

» A □

0.6 ■ ° o

• * '1

o b" 5J ■ T * - | § o • o

0.4 - s □ o

V Prn= 0.05

Pm=0.1

0.2 - A Pm= 0.5

□ Pm=2

O Pm= 5

' 10"° 10"° 10"+ 10" E

(a) Driving mechanism

10"' 10"° 10"° lO"* 10" E

(b) Driving location

Fig. 5. Parameter space plot of U'c against Ekman number for all runs showing predominant driving mechanism/location. (a): Black/white points indicate driving by Lorentz/Reynolds force. (b): Black/white points indicate driving exclusively at the TC/away from the TC. Grey points are those for which the driving mechanism/location is unclear.

dication of the transition from Lorentz driving to Reynolds stress driving. Since E depends on viscosity, this suggests that Reynolds stress driving is an artefact of the enhanced viscosity present in dynamo simulations, which is put in for numerical stability reasons rather than physical reasons. TOs driven by the Reynolds force occupy the region of parameter space where E is large and U'c small (white points), whereas as E decreases and the short timescale flow is more geostrophic there is a change to a regime where TOs are driven by the Lorentz force (black points). Between the two regimes is a region of parameter space inhabited by grey points and it is here that the two forces are of a similar magnitude, or the correlation is unclear, and we are unable to distinguish exactly which acts as the driving force.

Dynamo runs from TJT were found to be driven by the Reynolds force and, given that the Lorentz driving regime is located at smaller Ekman numbers and larger field strengths than were accessible in dynamo simulations, this also explains why we were unable to observe TOs driven by FLD previously.

Complementing the regime diagram of driving mechanisms is a similar plot for the driving location of TOs (shown in Fig. 5b), again determined by surveying our full set of simulations over their time period t . Here black points indicate runs where TOs were found to exclusively originate at the TC, white points are runs where TOs were never found to be driven at the TC, and grey points are the remaining runs (i.e. driving both at, and away from, the TC). As the Ekman number is reduced there is a transition in the driving location from one regime to the other, mirroring the property of the driving mechanism. We briefly commented in Section 3.2 that a low Ekman number favoured excitation at the TC and this is confirmed in Fig. 5b, which clearly shows this preference. At larger Ekman numbers, where U'c is also smaller, excitation can occur away from the TC and, in fact, at large enough values of E the role of the TC is greatly diminished as we were unable to find driving there at all.

In order to complement Figs. 5a and 5b we plot colour-coded density plots of FR and FLD, in Fig. 6, for the four runs of Table 1 that displayed TOs. The same white curves that were displayed in Fig. 4 are again overlaid on these plots although we now expect to see features in FR and FLD at the origin of these curves, rather than features travelling along the curves themselves.

Figs. 6a and 6b show the approximate parity in magnitude between FR and FLD in runs with large Ekman number, as well as evidence for excitation of waves throughout the region OTC. Con-

versely, the remaining plots in Fig. 6 show a clear preference for driving, exclusively at the TC, by one of the two forces. By comparing Figs. 6c and 6d it is clear from the correlation that driving is orchestrated by the Reynolds force, which is also approximately an order of magnitude larger than the Lorentz force. However, driving at the TC by FLD becomes preferable when smaller Ekman numbers are considered as evidenced by comparing Figs. 6e and 6f, as well as Figs. 6g and 6h. With the reduction in Ekman number the Reynolds force is now dwarfed by the Lorentz force throughout the domain and, in particular, at the TC where the consistent correlation between FLD and wave excitation is extremely clear in Fig. 6f.

As expected, Fig. 6 displays a change of TO-regime with reduced Ekman number that was indicated in Fig. 5. The plots of Fr and Fld help to show just how stark the transition is from a regime with TOs driven by the Reynolds (and possibly Lorentz) force throughout the core at large E to a regime where waves are driven exclusively at the TC by the Lorentz force.

3.4. core travel times

Observational evidence has suggested that the time taken for TOs to travel across the Earth's outer core is 4 years (Gillet et al., 2010) and hence we calculate the equivalent value from our simulations (tce , as defined by (17)) to enable the identification of Earth-like travel times in our work. Fig. 7 shows the core travel time, in years, against the Ekman number, separated into individual plots for each Pm to aide clarity.

There are two striking features immediately apparent from Fig. 7. First, many of our simulations permit TOs travelling across the core on Earth-like timescales, and this is true at various values of Pm and B0. Secondly, it is also clear, at least for the smaller values of Pm, that, for a given field strength and magnetic Prandtl number, the value of tce becomes independent of E at low Ekman number. This saturation of tce with E is less clear at larger values of Pm (see Figs. 7d and 7e) because these regimes are inhibited by the inability to reach small enough values of the Ekman number in our simulations. However, this intriguing result may be true at the various values of the parameters used across our set of simulations with the saturation of tce possibly becoming apparent at the relatively large Ekman numbers used for Pm = 5, see Fig. 7e. Independence from the Ekman number occurs earliest for weaker fields and hence the convergence is perhaps most clear in the B0 = 0.1 cases for Pm = 0.05, in Fig. 7a, and Pm = 0.1, in Fig. 7b, which are also aided by their ability to reach the smallest Es. Since our travel times are normalised on the Alfven speed at the CMB, this independence of E suggests that the distribution of UA as a function of s (see Fig. 3) does not change dramatically even when E is lowered to very small values.

4. Discussion

Our magnetoconvection study in this paper has allowed us to perform a large number of simulations and probe parameter space unattainable in previous dynamo studies (Wicht and Christensen, 2010; Teed et al., 2014) and therefore we have been able to explore the dynamics and excitation of TOs under new conditions. Across the set of runs we found the geostrophy parameter, U 'c, was able to give an excellent indication of the ability to observe TOs in any given simulation. Perturbations of a quasi-Taylor state leading to TOs may, on some level, always be possible, but TOs must account for ~40% of the short timescale flow, in order for the waves to be observed over convection. There was evidence of the geostrophy parameter indicating an onset value for TOs in previous work (TJT) but we have now been able to cement its role by writing U'c . & 0.4. The stability of U'c even as the Ekman num-

(a) Fr for E3P5d

0.0005 0.0010 0.0015 0.0020 t

(c) Fr for E5P2a

(e) Fr for 5E6P.ld

(g) Fr for E6P.5a

(b) Fld for E3P5d

1.8E+05

1.2E+05 o

5.9E+04 1 0

0.0E+00 « 0 8

-5.9E+04

-1.2E+05 0 n

-1.8E+05 0 0

0.0000

1.8E+05 1.2E+05 6.0E+04 0.0E+00 -6.0E + 04

.0005 0.0010 0.0015 0.0020 t

(d) FLD for E5P2a

4.1 E+04 2.8E+04 1.4E+04 0.0E+00 -1.4E+04 -2.8E+04 -4.1 E+04

(f) Fld for 5E6P.ld

0.001 0.002 t

(h) Fld for E6P.5a

Fig. 6. Reynolds force, FR (left column) and Lorentz force, FLD (right column) as a function of distance, s, from the rotation axis and time, t, in magnetic diffusion units. White curves (which are the same curves as those used in Fig. 4) indicate the trajectory of a particle travelling at the Alfven speed, UA(s). We look for features in FR and Fld that appear at the origin of these curves. (For interpretation of the references to colour in this figure, the reader is referred to the web version of this article.)

ber is reduced suggests it may continue to be a good measure of a critical value for TO onset at more Earth-like parameter regimes.

Within the suite of runs performed we have been able to identify two regimes for TOs and a transition between them. At large Ekman numbers TOs are only found when the magnetic field strength is also large; this can be achieved via large Pm, large B0, or a combination of the two. This regime is characterised by TOs generated by the Reynolds force throughout the region OTC, which are of a similar strength to the freely occurring convection. The geostrophy parameter is therefore close to the critical value for TO observation so that a small tweak of the input parameters can lead to a state without the dominance of TOs. Excitation of TOs can be controlled by the Reynolds force, regularly the largest forcing term,

because the relatively large rotation rates both contribute to larger convection velocities and inhibit large magnetic field strength leading to a relatively weak Lorentz force.

In the second regime, which exists at significantly larger rotation rates, TOs can be found in simulations at smaller values of the magnetic Prandtl number but the Ekman number needs to fall below ^ 10-6, values difficult to achieve in dynamo simulations currently. The regime is typified by a large Lorentz force at a radially thin region at the TC, which can, in several runs, lead to periodic excitation of TOs from the TC that propagate outwards through the outer core. The reduction in Ekman number leading to this regime suppresses the Reynolds force and boosts the Lorentz force so that excitation (on timescales relevant to TOs) only occurs at the TC.

15 ......... ......... ......... ......... 15 ......... 1 1

10 V 10

TT ' v o

0 0 «8

5 No B0 5 No B0

BD=0.1 B0=0.1

B0=0.5 B0=0.5

B0= 1 B0= 1

£>0=1° B0= 10

0 ........1 ........1 ........1 ........1 0 ........i .....

10~3 10~4 E

(a) Pm = 0.05

r° 10~3 10"' E

(b) Pm = 0.1

B„= 0.1

£>o= U.b

= 0.1 = 0.5 = 1 = 10

(c) Pm = 0.5

10-7 1Q-6 10-5 10-4 1C

(d) Pm = 2

„=0-5

„=10

1q-7 1Q-6 .|0-5 10-4 1Q-3

(e) Pm = 5

Fig. 7. Parameter space plots of core travel time, t^ , in years against Ekman number for various values of the magnetic Prandtl number. Filled points indicate runs for which TOs were observed. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

across the TC (Livermore and Hollerbach, 2012). From our current work it is unclear which terms are the largest in the complicated definition of FLD and therefore the exact physical mechanism, and periodic nature, of wave excitation using these terms must be examined in future work. Associated with the low Ekman number regime are core travel times that are independent of E itself and, encouragingly for many parameter sets, these travel times converge to values close to the 4 years expected for TOs in the Earth.

There also exists a region of parameter space where the two regimes overlap and in these cases the driving mechanism is often unclear because of an approximate parity between Reynolds and Lorentz forcing. However, in this region there can also be systems where driving is controlled by Reynolds forcing at the TC irrespective of the size of the Lorentz force, which, indeed, may be the larger of the two forces. This type of system, identified by a good correlation between FR and the origin of waves at the TC, was previously found in TJT and it consequently restricted the conclusions on excitation mechanisms in that work. It now appears that this type of system is at the edge of the numerically achievable parameter space of geodynamo simulations and we believe that the Lorentz forced regime, identified in our current magnetoconvec-tion simulations, would become apparent if the parameters could be pushed to more Earth-like values.

The imposition of a dipolar field on the outer boundary in our simulations may have created the effect of suppressing flow ITC due to the concentration of magnetic field near the poles. Consequently, TOs operating within, as well as crossing, the TC are rare in our current set of simulations although they were previously found in dynamo simulations. Increasing preference for generation of TOs at the TC at low E suggests that waves crossing the TC may not exist in the Earth, a hypothesis supported by observational data (Gillet et al., 2010). Nevertheless, it seems to be possible for TOs to arise at the TC and propagate inwards ITC although this was often not the case in these magnetoconvection simulations perhaps due to our choice of imposed field morphology. The simplified choice of magnetic boundary condition at the inner boundary may additionally compound the problem.

In summary, our current work has shown that as the Ekman number becomes smaller, the driving part of the Lorentz force triggers TOs at the TC with core travel times that are independent of E itself. Additionally, TOs can manifest in an Earth-like quasi-periodic fashion, with a period of around 4 to 5 years, reminiscent of observation results based on core flow models. The discovery of these key results has only been possible with the implementation of magnetoconvection since the parameter regimes necessary are currently inaccessible in dynamo models. However, it remains to be seen which aspects of our results are artefacts of our choice of a dipolar field, which is a logical harmonic to choose for a first investigation. For instance, adjusting the imposed field to one made up of different spherical harmonics, or using conducting boundary conditions could promote flow, and therefore TOs, within the TC as well as affect the convergence of core travel times at low Ekman number. These ideas should be the subject of future investigations.

Acknowledgement

Upon consideration of the form of FLD (see Eq. (13) from TJT) it appears that a large Lorentz force is generated by a combination of large magnetic fields as well as a contribution from the shearing of ageostrophic convection, possibly on various timescales. We therefore believe that excitation becomes confined to the TC because at low Ekman number it is now the only region exhibiting significant velocity shears through the discontinuity of fluid flow

This work was supported by the Natural Environment Research Council, grant NE/I012052/1.

Appendix A. Supplementary material

Supplementary material related to this article can be found online at http://dx.doi.org/10.1016/jxpsl.2015.02.045.

References

Alfven, H., 1942. Existence of electromagnetic-hydrodynamic waves. Nature 150, 405-406.

Bloxham, J., Zatman, S., Dumberry, M., 2002. The origin of geomagnetic jerks. Nature 420, 65-68.

Braginsky, S., 1970. Torsional magnetohydrodynamic vibrations in the Earth's core and variation in day length. Geomagn. Aeron. 10, 1-8.

Braginsky, S., 1984. Short-period geomagnetic secular variation. Geophys. Astrophys. Fluid Dyn. 30, 1-78.

Brown, W., Mound, J., Livermore, P., 2013. Jerks abound: an analysis of geomagnetic observatory data from 1957 to 2008. Phys. Earth Planet. Inter. 223, 62-76.

Busse, F., Simitev, R., 2005. Convection in rotating spherical fluid shells and its dynamo states. In: Soward, A., Jones, C., Hughes, D., Weiss, N. (Eds.), Mathematical Aspects of Natural Dynamos. CRC Press, pp. 359-392.

Chao, B., Chung, W., Shih, Z., Hsieh, Y., 2014. Earth's rotation variations: a wavelet analysis. Terra Nova 26 (4), 260-264.

Christensen, U., Aubert, J., 2006. Scaling properties of convection-driven dynamos in rotating spherical shells and application to planetary magnetic fields. Geophys. J. Int. 166 (1), 97-114.

Dormy, E., Soward, A., Jones, C., Jault, D., Cardin, P., 2004. The onset of thermal convection in rotating spherical shells. J. Fluid Mech. 501, 43-70.

Dumberry, M., Bloxham, J., 2003. Torque balance, Taylor's constrait and torsional oscillations in a numerical model of the geodynamo. Phys. Earth Planet. Inter. 140, 29-51.

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.

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 (3), 380-390.

Gubbins, D., Bloxham, J., 1985. Geomagnetic field analysis. Part III. Magnetic fields on the core-mantle boundary. Geophys. J. R. Astron. Soc. 80, 695-713.

Holme, R., de Viron, O., 2013. Characterization and implications of intradecadal variations in length of day. Nature 499 (7457), 202-204.

Jackson, A., 1997. Time dependence of geostrophic core-surface motions. Phys. Earth Planet. Inter. 103, 293-311.

Jault, D., 2008. Axial invariance of rapidly varying diffusionless motions in the Earth's core interior. Phys. Earth Planet. Inter. 166 (1), 67-76.

Jault, D., Gire, C., LeMouel, J.-L., 1988. Westward drift, core motion and exchanges of angular momentum between core and mantle. Nature 333, 353-356.

Jones, C., Boronski, P., Brun, A., Glatzmaier, G., Gastine, T., Miesch, M., Wicht, J., 2011. Anelastic convection-driven dynamo benchmarks. Icarus 216, 120-135.

Legaut, G., 2005. Ondes de torsion dans le noyau terrestre. Ph.D. thesis. Univ. Joseph Fourier, Grenoble, France.

Livermore, P.W., Hollerbach, R., 2012. Successive elimination of shear layers by a hierarchy of constraints in inviscid spherical-shell flows. J. Math. Phys. 53, 073104.

Mound, J., Buffett, B., 2006. Detection of a gravitational oscillation in length-of-day. Earth Planet. Sci. Lett. 243, 383-389.

Olson, P., Christensen, U., 2006. Dipole moment scaling for convection-driven planetary dynamos. Earth Planet. Sci. Lett. 250 (3), 561-571.

Roberts, P., 2007. Theory of the geodynamo. In: Olson, P. (Ed.), Treatise of Geophysics. Elsevier, pp. 67-105.

Roberts, P., Aurnou, J., 2012. On the theory of core-mantle coupling. Geophys. As-trophys. Fluid Dyn. 106 (2), 157-230.

Roberts, P., Wu, C.-C., 2014. On the modified Taylor constraint. Geophys. Astrophys. Fluid Dyn. 108 (6), 696-715.

Schaeffer, N., Jault, D., Cardin, P., Drouard, M., 2012. On the reflection of Alfven waves and its implication for Earth's core modelling. Geophys. J. Int. 191, 508-516.

Taylor, J., 1963. The magneto-hydrodynamics of a rotating fluid and the Earth's dynamo problem. Proc. R. Soc. A, Math. Phys. Eng. Sci. 274, 274-283.

Teed, R., Jones, C., Tobias, S., 2014. The dynamics and excitation of torsional waves in geodynamo simulations. Geophys. J. Int. 196 ( 2), 724-735.

Wicht, J., Christensen, U., 2010. Torsional oscillations in dynamo simulations. Geo-phys. J. Int. 181, 1367-1380.

Willis, A.P., Sreenivasan, B., Gubbins, D., 2007. Thermal core-mantle interaction: exploring regimes for 'locked' dynamo action. Phys. Earth Planet. Inter. 165 (1), 83-92.

Zatman, S., Bloxham, J., 1997. Torsional oscillations and the magnetic field with the Earth's core. Nature 388, 760-761.