Kinematic validation of a quasi-geostrophic model for the fast dynamics in the Earth’s outer core

Maffei, S.; Jackson, A.

doi:10.1093/gji/ggx263

Summary

We derive a quasi-geostrophic (QG) system of equations suitable for the description of the Earth’s core dynamics on interannual to decadal timescales. Over these timescales, rotation is assumed to be the dominant force and fluid motions are strongly invariant along the direction parallel to the rotation axis. The diffusion-free, QG system derived here is similar to the one derived in Canet et al. but the projection of the governing equations on the equatorial disc is handled via vertical integration and mass conservation is applied to the velocity field. Here we carefully analyse the properties of the resulting equations and we validate them neglecting the action of the Lorentz force in the momentum equation. We derive a novel analytical solution describing the evolution of the magnetic field under these assumptions in the presence of a purely azimuthal flow and an alternative formulation that allows us to numerically solve the evolution equations with a finite element method. The excellent agreement we found with the analytical solution proves that numerical integration of the QG system is possible and that it preserves important physical properties of the magnetic field. Implementation of magnetic diffusion is also briefly considered.

Core, Rapid time variations, Numerical modelling

1 INTRODUCTION

The Earth’s magnetic field shows oscillations and variations on a wide range of temporal scales (Hulot et al.2015). This variability is the observable consequence of the rich dynamics taking place in the Earth’s outer core, where the geomagnetic field is generated and continuously altered by the complex interplay between fluid flows and the geomagnetic field itself (Jackson & Finlay 2015; Roberts 2015).

Numerical geodynamo models (Christensen & Wicht 2015) can capture qualitative features of the geomagnetic secular variation and of the spatial structure of the geomagnetic field at the core–mantle boundary (CMB) but the whole temporal spectrum is impossible to cover due to the prohibitive computational costs required. For the same reason, current numerical simulations operate in a parameter regime that is far away from what is thought to be the one in which the geodynamo operates. For instance, the ratio of the rotation to viscous timescales, as measured by the Ekman number Ek, is supposed to be extremely small in the Earth’s core for which Ek = O(10⁻¹⁵). Recent numerical simulations (Schaeffer et al.2017) are capable of reaching Ek = O(10⁻⁷) suggesting that they might not operate in the same regime as the true geodynamo. Conversely, according to Aubert et al. (2017), these simulations might be in the correct asymptotic regime. In the same paper, the authors present the results from simulations where the value Ek = 10⁻⁸ has been reached with hyperdiffusivity applied to the smallest length scales. The ratio of magnetic to viscous dissipation timescales, as measured by the magnetic Prandtl number Pm is also orders of magnitude higher in numerical simulations than in the core: Pm = O(10⁻⁵–10⁻⁶) in the core while Pm > 3 × 10⁻² in recent simulations (Aubert et al.2017; Schaeffer et al.2017). These numbers indicate that the enormous spatial and temporal scale separation that is thought to be characteristic of the dynamics in the core is not present in numerical models. Unfortunately, by lowering Pm it becomes more difficult to obtain self-sustained dynamo simulation, as clarified by the relationship Rm = PmRe between the Reynolds number Re (a measure of the vigour of the flow velocities) and the magnetic Reynolds number Rm. Since the latter has to overcome a certain critical value for dynamo action to take place, one can see how lowering Pm to Earth-like values tends to result in the generation of weaker magnetic fields. Therefore stronger fields require more vigorous flows and consequently additional computational burden. As a consequence, geodynamo simulations typically have a lower ratio of magnetic to kinetic energy (less than 10) than what it is believed to be relevant for the Earth (where this value is about 100). See Schaeffer et al. (2017) for a summary of recent geodynamo simulations and for their comparison with the Earth’s dynamo.

A possible route towards the development of numerical models capable of operating in more realistic parameter regimes is given by a class of models that attempts to represent only the physics that is thought to be relevant for the Earth’s core. In particular, models based on the quasi-geostrophic (QG) approximation take advantage of the smallness of the Ekman and Rossby (Ro) numbers, the latter being a measure of the importance of inertia over rotational forces. These non-dimensional numbers are an indication of how rotation (via the Coriolis force) is dominant in the force budget of the core. As such the Earth’s core is said to be in rapid rotation. The fluid flows characteristic of such a system are essentially 2-D, with vertical component strongly inhibited by the presence of rotation, as predicted by the Proudman–Taylor theorem (Greenspan 1968; Jacobs 1987). Even when magnetic and thermal effects are included (by including in the momentum equation the Lorentz force and buoyancy, respectively) flows retain their columnar structure (vertically invariant), as long as the main balance remains geostrophic, that is, a balance between pressure gradient and the Coriolis force. The remaining forces, most notably inertia, Lorentz force and buoyancy, enter the balance at subdominant order, together with the Coriolis and pressure forces due to ageostrophic flows (Gillet et al.2011; Calkins et al.2015). An account of the experimental and numerical evidence for bidimensionality in rapidly rotating flows is given in Cardin & Olson (2007) and Williams et al. (2010).

These theoretical arguments, together with experimental and numerical evidence suggest that it is acceptable to impose a flow structure that is almost invariant along the vertical direction. This approach has been originally developed for studies of thermal convection (Cardin & Olson 1994; Aubert et al.2003; Gillet & Jones 2006; Guervilly & Cardin 2016) and has subsequently been used for studies of kinematic dynamos (Schaeffer & Cardin 2006; Schaeffer et al.2016), the study of interannual and decadal dynamics in a data assimilation framework (Canet et al.2009) and the study of magneto-hydrodynamic oscillations in the Earth’s core (Canet et al.2014; Labbé et al.2015). The columnar flow approximation has also been successfully applied to the inversion of geomagnetic data for the retrieval of the flows at the surface of the core (Pais & Jault 2008; Gillet et al.2009, 2015).

The use of the columnar flow approximation for the study of interannual and decadal dynamics has been justified in Jault (2008) and Gillet et al. (2011) where an impulsive perturbation of the inner core boundary (ICB) propagated into the exterior shell of a 3-D numerical model of a rotating fluid permeated by a magnetic field. These transient motions are shown to be highly columnar provided that magnetic forces do not become dominant with respect to rotation and that the disturbances propagate on timescales much shorter than the magnetic diffusion timescale and much longer than the rotation period of the shell. In the Earth’s core the propagation timescales of magnetic disturbances (the Alfvén velocity) across the core is thought to be of the order of 6–8 yr (Gillet et al.2010, 2015) and the magnetic diffusion timescale can be estimated to be at least 10⁴ times this value (see Canet et al. (2014) and references therein). It is safe to assume that, in the Earth’s outer core, interannual and decadal flows and fluctuations are columnar. Examples of such motions include Alfvén torsional waves (Braginsky 1970) that propagate on the Alfvén timescales mentioned above, and slow magneto-hydrodynamic oscillations (Hide 1966; Malkus 1967; Finlay 2011) that, depending on the geometry and intensity of the magnetic field in the core, could move on decadal to centennial timescales and therefore be a possible mechanism for the westward drift observed in magnetic field models (Hide 1966; Finlay & Jackson 2003; Jackson 2003). The detection of these motions in 3-D geodynamo numerical simulation has proven challenging, mainly due to the high values of Ek and Pm (Teed et al.2014; Hori et al.2015). Additionally, care has to be taken in relating the observation of interannual and decadal magnetic driven oscillations to their counterpart detected in geodynamo simulations, as the timescales are likely to be not in the correct ratio with rotational or advective timescales. As both input parameters (like Ek and Pm) and diagnostic quantities (such as the ratio of magnetic to kinetic energy) in current simulations are not the same as for the Earth’s core, so are the relative timescales (Glatzmaier 2002) of relevant phenomena. In particular the ratio of the Alfvén timescale to the convective turnover or to the rotational timescale tend to be smaller in simulations than in the Earth so that equating one timescale to Earth’s values, makes other ones wrong. These issues may be also tackled by QG models currently under development.

In the present study we consider a variation of the QG model of Canet et al. (2009). The key aspect of that study is that the magnetic field is considered through vertical averages of quadratic combinations of its equatorial components. In this way both the momentum and the induction equations are projected on the equatorial plane and a fully 2-D model is obtained. The properties of the model of Canet et al. (2009) have not yet been fully characterised through a numerical simulation of the system, but data assimilation experiments have been performed in simplified settings. Subsequent QG models of magnetohydrodynamic oscillations (Canet et al.2014; Labbé et al.2015) have assumed that the magnetic field too could be approximated as columnar, with the vertical component of the field being either negligible or linearized along the vertical direction. This approximation is more restrictive than the general description of Canet et al. (2009), in which no assumption about the vertical structure of the magnetic field is made, but it allows the consistent inclusion of magnetic diffusion and the natural imposition of insulating boundary conditions at the CMB. Both these issues are addressed in the present study.

We perform kinematic forward numerical simulations of a system of equation similar to the one of Canet et al. (2009), the most notable differences being the use of a velocity formulation that enforces mass conservation (Schaeffer & Cardin 2005) and the use of vertically integrated magnetic quantities instead of vertically averaged ones to describe the magnetic field. In addition, the presence of the inner core is neglected from now on. We made this choice to simplify our study and to focus attention on the significant aspects of the QG model that we now derive. In Canet et al. (2009) and Canet et al. (2014) the magnetohydrodynamic equations are solved only outside the tangent cylinder, the imaginary cylinder aligned with the rotation axis of the Earth and tangent to the inner core and separating the regions above and below the inner core from the rest of the outer core. Whether the geostrophic balance or the columnar flow approximation hold in the fluid regions inside the tangent cylinder is an open question. In Canet et al. (2009, 2014), these regions are completely excluded from the QG dynamics on the basis that, if the flows are indeed vertically invariant outside the tangent cylinder little feedback is expected from the fluid regions above and below the inner core. A full sphere represents the simplest canonical system to study and the correct geometry for an early Earth, in which the inner core is not yet present, and some experimental configurations. As prescribed velocity field we chose a time-independent zonal flow with a cylindrical radial dependence. As we show later in the paper, this choice allows significant simplifications in the governing equations and the derivation of a relatively simple analytical solution. We do not expect the fundamental conclusion of our study to vary for more complex flows (such as non-axysimmetric flows). In Section 2, the system of equations under the QG approximation is derived and its peculiar characteristics discussed. In Section 3, we illustrate the ideal kinematic setup and the criteria chosen to validate the QG equations. By imposing a time-invariant velocity field we can study specific properties of the model in a highly simplified setup. The results of the validation are given in Section 4. In Sec-tion 5, we propose an alternative to magnetic diffusion, not possible to implement consistently in the model we derive here. Comments and conclusions are given in Sections 6 and 7.

2 GOVERNING EQUATIONS

We describe the Earth’s outer core as a rotating fluid sphere of radius r_c, whose axis of rotation is |$\boldsymbol{\Omega }$| (assuming that it passes through the centre of the sphere), with homogeneous density ρ, kinematic viscosity ν and magnetic diffusivity η. The velocity u, the non-hydrostatic pressure p and the magnetic field B are described in a cylindrical coordinate system with origin at the centre of the sphere, coordinates (s, ϕ, z) and unit vectors (e_s, e_ϕ, e_z). The vertical direction e_z is parallel to the rotation vector |$\boldsymbol{\Omega }$|⁠, z is the distance from the equatorial plane, s is the distance from the rotation axis and ϕ is the longitude. The radial distance from the centre of the sphere is therefore

\begin{eqnarray} r=\sqrt{s^2+z^2}. \end{eqnarray}

(1)

The boundary of the sphere (the CMB) is r = r_c and, in cylindrical coordinates, it is located at a height ±H above the equatorial plane, where:

\begin{eqnarray} H=\sqrt{r_c^2-s^2}. \end{eqnarray}

(2)

The geometry of the system is depicted in Fig. 1. Note that H = 0 at s = r_c.

Figure 1.

Geometry of a quasi-geostrophic column in a full sphere. H is the half height of the column, x = s cos ϕ and y = s sin ϕ, s is the distance from the rotation axis and ϕ is the azimuth.

Open in new tab Download slide

The equations describing the combined evolution of u and B are, under the hypothesis of incompressibility:

\begin{eqnarray} \frac{\partial {\bf u}}{\partial t} + ({\bf u}\cdot \nabla ){\bf u} + 2\mathbf {\boldsymbol\Omega } \times {\bf u} = -\frac{\nabla p}{\rho } + \frac{({\bf B}\cdot \nabla ){\bf B}}{\mu _0\rho } + \nu \nabla ^2{\bf u} \end{eqnarray}

(3)

\begin{eqnarray} \frac{\partial {\bf B}}{\partial t} = \nabla \times ({\bf u}\times {\bf B}) + \eta \nabla ^2 {\bf B} \end{eqnarray}

(4)

\begin{eqnarray} \nabla \cdot {\bf u}=0 \end{eqnarray}

(5)

\begin{eqnarray} \nabla \cdot {\bf B} =0 \end{eqnarray}

(6)

where the magnetic pressure term |B|²(2μ₀)⁻¹ has for simplicity being incorporated in the pressure p and μ₀ is the vacuum magnetic permeability constant. These equations are, respectively, the momentum equation per unit mass, the induction equation and the divergence free condition for the velocity field and the magnetic field and are commonplace in the study of the Earth’s core (Jacobs 1987). We consider the non-dimensional version of these equation by introducing |$\mathcal {B}$|⁠, |$\mathcal {U}$|⁠, r_c, |$r_c/\mathcal {U}$|⁠, |$\mathcal {P}=\rho \mathcal {U}^2$| as the typical scales for magnetic field, velocity, length, time and pressure respectively. Then the non-dimensional MHD equations are obtained by making the following substitutions:

\begin{eqnarray} {\bf u} &\rightarrow& \frac{{\bf u}}{\mathcal {U}}\nonumber\\ p&\rightarrow& \frac{p}{\mathcal {P}}\nonumber\\ {\bf B} & \rightarrow& \frac{{\bf B}}{\mathcal {B}}\nonumber\\ \frac{\partial }{\partial t} & \rightarrow& \frac{r_c}{\mathcal {U}}\frac{\partial }{\partial t}\nonumber\\ \nabla & \rightarrow& r_c \nabla . \end{eqnarray}

(7)

The governing equations written in non-dimensional form are

\begin{eqnarray} && \frac{\partial {\bf u}}{\partial t} + ({\bf u}\cdot \nabla ){\bf u} + \frac{2 \Omega r_c}{\mathcal {U}} {\bf e}_z\times {\bf u}\nonumber\\ &&\quad = - \nabla p + \frac{\mathcal {B}^2}{\rho \mu _0 \mathcal {U}^2}({\bf B}\cdot \nabla ){\bf B} + \frac{\nu }{ \mathcal {U} r_c} \nabla ^2 {\bf u} \end{eqnarray}

(8)

\begin{eqnarray} \frac{\partial {\bf B}}{\partial t} = \nabla \times ({\bf u}\times {\bf B}) + \frac{\eta r_c}{\mathcal {U}} \nabla ^2{\bf B} \end{eqnarray}

(9)

\begin{eqnarray} \nabla \cdot {\bf u}=\nabla \cdot {\bf B} =0. \end{eqnarray}

(10)

Here |$\Omega =|\boldsymbol{\Omega }|$| is the rate of rotation. We introduce the Alfvén velocity:

\begin{eqnarray} V_a=\frac{\mathcal {B}}{\sqrt{\rho \mu _0}} \end{eqnarray}

(11)

which is the velocity of propagation of magnetic disturbances in an electrically conducting fluid (Alfvén 1942). We then define the Lehnert, Ekman and Lundquist numbers as in Jault (2008):

\begin{eqnarray} Le=\frac{V_a}{r_c\Omega }=\frac{\mathcal {B}}{r_c\Omega \sqrt{\rho \mu _0}}, \end{eqnarray}

(12)

\begin{eqnarray} Ek=\frac{\nu }{\Omega r_c^2} \end{eqnarray}

(13)

and

\begin{eqnarray} Lu =\frac{r_c \mathcal {B}}{\eta (\mu \rho )^{1/2}}, \end{eqnarray}

(14)

respectively. More specifically, Le measures the ratio of the rotation period over the Alfvén timescale, which is the propagation time of a magnetic disturbance of velocity V_a over the distance r_c. Assuming the same values adopted in Canet et al. (2014), its value in the Earth’s core is about 10⁻⁴. This value indicates that motions propagating on the timescale of magnetic disturbances are highly influenced by rotation. The importance of the diffusivities is measured by the Lundquist and the Ekman numbers. They are, respectively, the ratio of the magnetic diffusion timescale over the Alfvén timescale and the ratio of the rotation period over the viscous diffusion timescale. Since Lu ≃ 10⁵ and Ek ≃ 10⁻¹⁵, it is safe to assume that, on the short time scales we are interested in, we could entirely neglect both diffusivities (especially viscosity).

In what follows we revisit the derivation of the governing equations for a QG model that follows closely that of Canet et al. (2009). The quasi-z invariance of the flow suggests that the dynamics can be projected on the equatorial plane and a 2-D system can be solved in place of the original 3-D described by eqs (8)–(10). The treatment of the magnetic field requires special care as in general its vertical structure is unconstrained.

2.1 Columnar flow assumption and vertically integrated equations

Given the values of Le, Ek and Lu the conclusions of Jault (2008) apply and we consider the flow u to have a columnar structure. Following Schaeffer & Cardin (2005), Schaeffer & Cardin (2006) and Canet et al. (2014) we impose the velocity field to be of the form:

\begin{eqnarray} {\bf u}(s,\phi ,z) = \frac{1}{H}\nabla \times (\Psi {\bf e}_z) - \frac{z}{H^3} \frac{\partial \Psi }{\partial \phi }{\bf e}_z. \end{eqnarray}

(15)

This formulation describes a columnar flow that satisfies the constraint (5) and that is entirely defined by the 2-D scalar potential Ψ(s, ϕ). Furthermore, together with the condition

\begin{eqnarray} \Psi (s=1)=0, \end{eqnarray}

(16)

the flow (15) satisfies the non-penetration boundary condition

\begin{eqnarray} {\bf e}_r\cdot {\bf u}(r=1) = 0 \end{eqnarray}

(17)

everywhere on the CMB.

Given the flow structure (15) it is natural to project the momentum eq. (8) on the equatorial plane. Following Canet et al. (2014) and references therein we take the curl of eq. (8) and consider its vertical component (the axial vorticity equation):

\begin{eqnarray} &&\frac{\partial \omega _0}{\partial t} + ({\bf u}_e\cdot \nabla _e)\omega _0 +\frac{2}{Le}\frac{s}{H^2}u_s\nonumber\\ &&\quad= {\bf e}_z\cdot \nabla \times [({\bf B}\cdot \nabla ){\bf B}]_e + \frac{E}{Le}\nabla ^2\omega _0. \end{eqnarray}

(18)

Here []_e denotes the equatorial part of the vectorial quantity enclosed in square brackets and β = H⁻¹∂_sH = −sH⁻². By inserting the columnar flow definition (15), the axial vorticity ω₀ is

\begin{eqnarray} \omega _0 = {\bf e}_z\cdot \nabla \times {\bf u}_e = \frac{1}{H}\left[ -\nabla ^2 + \beta \frac{\partial }{\partial s} \right]\Psi . \end{eqnarray}

(19)

Note that the first term on the right-hand side of eq. (18) is not independent of z. In Canet et al. (2014) and Labbé et al. (2015), the magnetic field B is assumed to have a structure equivalent to eq. (15) and eq. (18) is easily projected on the equatorial plane. In Canet et al. (2009), no assumption is made regarding the vertical structure of B and eq. (18) is projected by considering its vertical average. Here we consider vertical integrals. We introduce the following notation

\begin{eqnarray} \lbrace \cdot \rbrace =\int _{-H}^H \cdot \ {\rm d}z \end{eqnarray}

(20)

to indicate a vertical integration and apply it to the axial vorticity eq. (18). The vertical average operation is then (2H)⁻¹{ · }. Both integration and average transform a 3-D quantity in a 2-D one by integrating along the vertical direction. The difference is that at s = 1 the value of the averaged quantities is equal to the value of the integrand there, while the integrated ones, as long as the original integrand is not singular at the boundary, vanishes.

As the terms related to the vorticity and velocity fields are not a function of z we obtain:

\begin{eqnarray} &&{ 2H\left[\frac{\partial \omega _0}{\partial t} + ({\bf u}_e\cdot \nabla _e)\omega _0 +\frac{2}{Le}\frac{s}{H^2}u_s - \frac{E}{Le}\nabla ^2\omega _0\right]}\nonumber\\ &&\quad= \lbrace {\bf e}_z\cdot \nabla \times [({\bf B}\cdot \nabla ){\bf B}]_e \rbrace . \end{eqnarray}

(21)

We now manipulate the right-hand side as described in Canet et al. (2009). By making use of the solenoidal nature of B and of Leibniz’ integration rule we obtain:

\begin{eqnarray} \&&{ \lbrace {\bf e}_z\cdot \nabla \times [({\bf B}\cdot \nabla ){\bf B}]_e\rbrace}\nonumber\\ &&{\quad = \left(\frac{\partial ^2}{\partial s^2} + \frac{3}{s}\frac{\partial }{\partial s} - \frac{1}{s^2}\frac{\partial ^2}{\partial \phi ^2}\right) \lbrace B_sB_\phi \rbrace} \nonumber\\ &&\quad-\,\frac{1}{s}\left(\frac{1}{s}\frac{\partial }{\partial \phi } + \frac{\partial ^2}{\partial s\partial \phi }\right)\left(\left\lbrace B_s^2\right\rbrace -\left\lbrace B_\phi ^2\right\rbrace \right) \end{eqnarray}

(22)

where terms evaluated on the surface r = 1 have been neglected by assuming that the magnetic field at the CMB is much weaker than the magnetic field in the interior of the domain. This assumption is supported by observational studies of torsional waves (Gillet et al.2010, 2015) and by geodynamo simulations (Aubert et al.2009). For convenience, we introduce the quadratic magnetic quantities:

\begin{eqnarray} a &=& \left\lbrace B_s^2\right\rbrace\nonumber \\ b &=& \left\lbrace B_\phi ^2\right\rbrace\nonumber \\ c &=& \left\lbrace B_s B_\phi \right\rbrace . \end{eqnarray}

(23)

Neglecting the nonlinear advection term and the viscosity in eq. (21) we obtain:

\begin{eqnarray} 2H\left(\frac{\partial \omega _0}{\partial t} +\frac{2}{Le}\frac{s}{H^2}u_s\right)&=& \left(\frac{\partial ^2}{\partial s^2} + \frac{3}{s}\frac{\partial }{\partial s} - \frac{1}{s^2}\frac{\partial ^2}{\partial \phi ^2}\right) c\nonumber\\ && -\,\frac{1}{s}\left(\frac{1}{s}\frac{\partial }{\partial \phi } + \frac{\partial ^2}{\partial s\partial \phi }\right)(a-b). \end{eqnarray}

(24)

By neglecting the effect of magnetic diffusion as well, evolution equations for a, b and c can be obtained from the induction eq. (9)¹:

\begin{eqnarray} \frac{\partial a}{\partial t}=-H \left( {\bf u}\cdot \nabla _e \right)\left( \frac{a}{H} \right) + 2a\frac{\partial u_s}{\partial s} +\frac{2c}{s}\frac{\partial u_s}{\partial \phi } \end{eqnarray}

(25)

\begin{eqnarray} \frac{\partial b}{\partial t} =-\frac{1}{H}\left({\bf u}\cdot \nabla _e\right)\left( Hb\right) +2sc\frac{\partial }{\partial s}\left(\frac{u_\phi }{s}\right) - 2b\frac{\partial u_s}{\partial s} \end{eqnarray}

(26)

\begin{eqnarray} \frac{\partial c}{\partial t} =-\left({\bf u}\cdot \nabla _e\right) c+sa\frac{\partial }{\partial s}\left( \frac{u_\phi }{s} \right) +\frac{b}{s}\frac{\partial u_s}{\partial \phi }. \end{eqnarray}

(27)

Note that in the derivation of these equations, no assumption on the intensity of the magnetic field is necessary. Surface terms have to be neglected in the derivation of eq. (22) for the system (24)–(27) to be closed in the quantities a, b and c. For the same reason, magnetic diffusion needs to be neglected from the induction equation. It has to be pointed out that the vertical component of B has not been neglected in any way (apart from the assumption of the field at the CMB being negligible). Thanks to the solenoidal nature of the magnetic field the component B_z is implicitly considered in the quantities a, b, c and in their derivatives. Because of the vertical integration, it is impossible to retrieve the original vector B from the quadratic magnetic variables. The differences between eqs (25)–(27) and eqs (17)–(19) of Canet et al. (2009) are due to the mass conservation here considered in the definition (15).

The system (24)–(27) has to be solved on the equatorial plane z = 0, with s ≤ 1. On this domain, the boundary conditions on a, b and c naturally arise given their nature of integral quantities. Since H vanishes on s = 1 and given that the quantities |$B_s^2$|⁠, |$B_\phi ^2$| and B_sB_ϕ are expected to remain finite everywhere:

\begin{eqnarray} a(s=1)=b(s=1)=c(s=1)=0. \end{eqnarray}

(28)

Physically, if we consider the mantle to be an insulator, we should enforce continuity of B with the exterior field B_e, valid for r > 1 which satisfies:

\begin{eqnarray} \nabla \times {\bf B}_e=0. \end{eqnarray}

(29)

However, we were not able to write this condition for a, b and c. In the more restrictive formulation employed in Canet et al. (2014) and Labbé et al. (2015) the insulating boundary condition can be naturally imposed at the CMB. The boundary conditions for the stream function Ψ is given in eq. (16). More specifically we have to make sure that

\begin{eqnarray} \Psi (s\rightarrow 1) = \mathcal {O}(H^3) \end{eqnarray}

(30)

in order for u_s to be zero and u_ϕ to be finite at s = 1.

3 VALIDATION UNDER THE KINEMATIC APPROXIMATION

The QG equations derived in the previous section can in principle be evolved forward in time from a given initial condition. Apart from simple experiments performed in Canet et al. (2009), there has been no characterization of the property of the system (24)–(27) through forward numerical simulations. In the remainder of the paper we present results of forward numerical models under the kinematic approximation. Although the QG model discussed here is targeted to magnetic fields and flows with interannual and decadal temporal variability, we impose a steady velocity field so that we can ignore the momentum equation and focus on solving the induction eqs (25)–(27). As we are neglecting the momentum equation, we follow Jones (2007) and consider typical values of velocity in the outer core to be U = 5 × 10⁻⁴ m s⁻¹ ≃ 15 km yr⁻¹ which result in a non-dimensional timescale of |$r_c/\mathcal {U}=220$| yr. These timescales lie at the upper limit of validity of the QG model developed here.

We note the following important properties of the quantities a, b and c. First, while a and b are bound to be positive, c can assume positive or negative sign. Second, Canet et al. (2009) noted that the Cauchy–Schwarz inequality applies:

\begin{eqnarray} ab \ge c^2. \end{eqnarray}

(31)

In particular we can introduce the quantity

\begin{eqnarray} q^2 = ab-c^2 \ge 0 \end{eqnarray}

(32)

that satisfies the following evolution equation:

\begin{eqnarray} \frac{\partial q^2}{\partial t}=-({\bf u}_e\cdot \nabla _e)q^2. \end{eqnarray}

(33)

3.1 Initial field and prescribed flow

The initial magnetic field used in this study is a poloidal magnetic field defined in spherical coordinates (r, θ, ϕ) by the following potential:

\begin{eqnarray} P_1^1=\frac{r^2}{2}\left( (r^2+k)\sin (\theta )\sin (\phi ) \right) \end{eqnarray}

(34)

where k = −5/3 is a constant. The magnetic field in the same system is

\begin{eqnarray} {\bf B}_{(r,\theta ,\phi )} = \nabla \times \nabla \times \left(P_1^1{\bf e}_r\right). \end{eqnarray}

(35)

This corresponds to a dipole whose axis is aligned along the y direction and satisfies insulating boundary conditions at the CMB (see Fig. 2). The corresponding quantities a, b and c are shown in Fig. 3. In this figure, the initial fields are interpolated on the finest mesh and second largest Lagrange polynomial degree we used for the numerical experiments (see the next section). The mathematical expressions of the initial conditions are, in cylindrical coordinates:

\begin{eqnarray} a&=&\frac{2}{45} H (61-62s^2+21s^4)\sin ^2\phi \nonumber\\ b&=&\frac{2}{45} H (61-152s^2+96s^4)\cos ^2\phi \nonumber\\ c&=&\frac{1}{45} H (61-107s^2+36s^4)\sin (2\phi ). \end{eqnarray}

(36)

Fig. 3(d) shows the quantity:

\begin{eqnarray} q^2 = ab-c^2 \end{eqnarray}

(37)

which, according to the Cauchy–Schwarz inequality (31) has to be positive. The analytical form of q² for the initial condition (36) is:

\begin{eqnarray} q^2= \frac{16}{45}s^4(1-s^2)^3\sin ^2(2\phi ). \end{eqnarray}

(38)

Figure 2.

Initial magnetic field for the kinematic analysis of the QG model. The left and right panels show a projection on the x, y plane (the equatorial plane) and the right panel is a meridional cut in the z, y plane, respectively.

Open in new tab Download slide

Initial conditions for (in order from panels a to d) a, b, c and q2 calculated from the initial magnetic field of Fig. 2. Maximum and minimum values are shown at the top and at the bottom of the colour scale for each quantity.

Figure 3.

Initial conditions for (in order from panels a to d) a, b, c and q² calculated from the initial magnetic field of Fig. 2. Maximum and minimum values are shown at the top and at the bottom of the colour scale for each quantity.

Open in new tab Download slide

The velocity field we impose is a simple azimuthal flow:

\begin{eqnarray} {\bf u}=s(s^2-1){\bf e}_\phi . \end{eqnarray}

(39)

Its radial functional form is shown in Fig. 4. As mentioned earlier more involved flows can be prescribed, as long as they can be described by the columnar formulation (15). With this velocity field the evolution eqs (25)–(27) are greatly simplified:

\begin{eqnarray} \frac{\partial a}{\partial t}=-\frac{u_\phi }{s}\frac{\partial a}{\partial \phi } \end{eqnarray}

(40)

\begin{eqnarray} \frac{\partial b}{\partial t} =-\frac{u_\phi }{s}\frac{\partial b}{\partial \phi } +2sc\frac{\partial }{\partial s}\left(\frac{u_\phi }{s}\right) \end{eqnarray}

(41)

\begin{eqnarray} \frac{\partial c}{\partial t} =-\frac{u_\phi }{s}\frac{\partial c}{\partial \phi }+sa\frac{\partial }{\partial s}\left( \frac{u_\phi }{s} \right) \end{eqnarray}

(42)

\begin{eqnarray} \frac{\partial q^2}{\partial t}=-\frac{u_\phi }{s}\frac{\partial q^2}{\partial \phi } . \end{eqnarray}

(43)

Here we also included the evolution equation for the control variable q². Note that a and q² evolve according to a simple advection equation. Therefore the analytical solution for these two quantities is

\begin{eqnarray} a_T(s,\phi ,t)=a_0\left(s,\phi -\frac{u_\phi }{s}t\right) \end{eqnarray}

(44)

where a₀(s, ϕ) = a(t = 0) is the initial condition and a_T indicates the true solution. Substitution of this form in the evolution equation for a will prove that this is indeed the solution. A similar theoretical solution holds for q²:

\begin{eqnarray} q^2_T(s,\phi ,t)=q^2_0\left(s,\phi -\frac{u_\phi }{s}t\right). \end{eqnarray}

(45)

Figure 4.

Radial structure of the zonal flow (39).

Open in new tab Download slide

Given the solution for a we can then solve (42):

\begin{eqnarray} c_T (s,\phi ,t) &=&\frac{2}{45} H \sin \left(\phi -\frac{u_\phi}{s} t\right) \bigg[2 (21 s^4-62 s^2+61) s^2 t \nonumber\\ &&\times\,\sin \left(\phi -\frac{u_\phi }{s}t\right) + (36 s^4-107 s^2+61)\nonumber\\ &&\times\cos \left(\phi -\frac{u_\phi}{s}t\right)\bigg] \end{eqnarray}

(46)

and finally, combining (44)–(46) with (37) we obtain a solution for b:

\begin{eqnarray} b_T(s,\phi ,t) &=& \frac{1}{45} H \bigg\lbrace 4 (36 s^4-107 s^2+61) s^2 t \sin \left[2 \left(\phi -\frac{u_\phi }{s} t\right)\right] \nonumber\\ && +\, (-84 s^8 t^2+248 s^6 t^2+s^4 (96-244 t^2)\nonumber\\ &&-\,152 s^2+61) \cos \left[2 \left(\phi -\frac{u_\phi }{s} t\right)\right] +\, 84 s^8 t^2\nonumber\\ &&-\,248 s^6 t^2+244 s^4 t^2+96 s^4-152 s^2+61\bigg\rbrace \end{eqnarray}

(47)

and a fully analytical description of the solution is at hand. Note that while a_T and q_T are given by pure advection of the initial condition, the fields c_T and b_T grow linearly and quadratically with time, respectively.

3.2 Numerical experiments setup

To assess the quality of the forward modelling, we evolve eqs (40)–(42) forward in time with the finite element method, commercially available software COMSOL Multiphysics.² The reason for this choice is twofold. First of all, we are interested in validating the QG equations for which there is almost no phenomenological knowledge. Therefore, since at this stage high performance is not required, we decided not to develop a customized code but to use an available, flexible software preferably with the possibility of testing various types of spatial discretization and time-stepping techniques. The second reason is that the finite element method is local in nature. We are mainly concerned with the positivity of a and b and, assuming that eqs (40)–(42) actually can preserve this property with time, there is always a risk that numerical noise generates spurious negative values. This is particularly true considering that we neglected diffusion. A local method should ensure that any numerical artefacts remain localized in space, while in a global method, such as a spectral one, there is no such guarantee.

In particular we are interested in preserving positivity of a, b and q². The quantities we compute in order to establish the performance of the model are:

The equatorial (non-dimensional) magnetic energy
\begin{eqnarray} 2E_m=\int _{{\rm core}}\left(B_s^2+B_\phi ^2\right) {\rm d}V= \int _{{\rm disc}} (a+b)\,{\rm d}S. \end{eqnarray}
(48)
The integral on the right-hand side is calculated over the surface of the equatorial disc. In general, we do not have information about the vertical component of the magnetic field, so we can only calculate the equatorial component of the magnetic energy. However, in the simple case we are studying here the vertical component of the magnetic field does not contribute to the evolution in time of the 3-D magnetic energy evaluated over the whole core. So E_m represents the time varying part of the total magnetic energy of the system and its evolution equation is
\begin{eqnarray} \frac{\partial E_m}{\partial t}=\int _{\rm disc} s c \frac{\partial }{\partial s}\left(\frac{u_\phi }{s}\right) {\rm d}S. \end{eqnarray}
(49)
Knowing the theoretical solutions for q² we can compute the domain integrals of the misfit
\begin{eqnarray} {\rm mean}\left(q^2-q^2_T\right)=\int _{\rm disc} |q^2-q^2_T| {\rm d}S \end{eqnarray}
(50)
where the theoretical solutions are of the form (45). A similar calculation can be performed for a, b and c. Here, however, we mainly focus on the deviations of the numerical solution of q² from its analytical form since, as we show below, this is the quantity that is most prone to numerical noise.
The minimum values of q², a and b. Since these must be always positive quantities their minimum value should always be greater than zero.

In COMSOL, we must choose several parameters to tune the numerical method, the most significant of which are the number of mesh elements used to discretize the spatial domain and the type and degree of the polynomials we use to discretize the functions representing the variables solved on the mesh. The mesh refinements range from extremely coarse to extremely fine (see Table 1), each of which is assigned an index: the lower the index, the finer the mesh. As for the discretization inside each element, we use Lagrange polynomials with variable degree from 1 to 3. The models are identified with the nomenclature LlNn where L is the degree of the polynomials used and N is the mesh identifier. So the model run with L = 2 and N = 1 will be called L2N1. Since by choosing a coarser mesh and a high degree L the finite element method tends towards a spectral method, we decide to keep L low and use fine meshes.

3.3 Cartesian formulation

As COMSOL allows great flexibility in the choice of the domain, its algorithms are formulated in Cartesian coordinates. As such any discretized physical field in the appropriate coordinate system is always smooth and differentiable everywhere on the domain. However, considering the geometry of the system, eqs (25)–(27) have been formulated in polar coordinates and an artificial singularity has been introduced at s = 0 (Lewis & Bellan 1990). This is evident in Figs 3 where the values of a, b and c are non-unique in ϕ as s vanishes. In Cartesian coordinates, sharp gradients develop near the origin and evolving eqs (40)–(42) in COMSOL results in significant numerical noise around s = 0 (not shown here). As a result, spikes of positive and negative values of similar magnitude appear to be generated in the field q², creating a violation of the Cauchy–Schwarz inequality that increases with time. At t = π/2 the quantity q² locally reaches the minimum value of −1.17. Considering that q² should always be positive and given its maximum expected value (see Fig. 3d), this violation of positivity is totally unacceptable.

These observations motivated us to develop an alternative formulation suited for numerical integration of the evolution equations with a local method that discretizes the variables in Cartesian coordinates, in the same way the finite element method is implemented by COMSOL. Let us introduce the quantities:

\begin{eqnarray} A&=&\left\lbrace B_x^2\right\rbrace =\int _{-H}^H B_x^2 {\rm d}z\nonumber \\ B&=&\left\lbrace B_y^2\right\rbrace =\int _{-H}^H B_y^2 {\rm d}z\nonumber \\ C&=&\left\lbrace B_xB_y\right\rbrace =\int _{-H}^H B_xB_y {\rm d}z. \end{eqnarray}

(51)

These are the analogues to a, b and c but in Cartesian coordinates. Here the cylindrical and Cartesian coordinates are related by the following transformations:

\begin{eqnarray} x &=& s \cos \phi\nonumber \\ y &=& s \sin \phi \end{eqnarray}

(52)

and z is not modified in the coordinate transformation. The relationship between the cylindrical quantities a, b and c and the Cartesian variations A, B and C is described by the formulae

\begin{eqnarray} A&=&a\cos ^2\phi + b\sin ^2\phi - 2c\sin \phi \cos \phi\nonumber \\ B &=& a\sin ^2\phi + b\cos ^2\phi + 2c\sin \phi \cos \phi \nonumber\\ C&=&(a-b)\sin \phi \cos \phi + c(\cos ^2\phi - \sin ^2\phi ) \end{eqnarray}

(53)

and

\begin{eqnarray} a&=&A\cos ^2\phi + B\sin ^2\phi + 2C\sin \phi \cos \phi \nonumber\\ b &=& A\sin ^2\phi + B\cos ^2\phi - 2C\sin \phi \cos \phi \nonumber\\ c&=&(B-A)\sin \phi \cos \phi + C(\cos ^2\phi - \sin ^2\phi ). \end{eqnarray}

(54)

The Cauchy–Schwarz inequality holds for the Cartesian quantities too:

\begin{eqnarray} AB\ge C^2. \end{eqnarray}

(55)

It can be easily proven that ab − c² = AB − C² and therefore the quantity q² is an invariant of the transformation. The evolution equations for A, B and C can be derived from (40)–(42) and the relationships (53):

\begin{eqnarray} \frac{\partial A}{\partial t}=-({\bf u}_e\cdot \nabla _e)A +A\left(2\frac{\partial u_x}{\partial x} + \frac{\partial u_z}{\partial z} \right) +2C \frac{\partial u_x}{\partial y} \end{eqnarray}

(56)

\begin{eqnarray} \frac{\partial B}{\partial t}=-({\bf u}_e\cdot \nabla _e)B +B\left(2\frac{\partial u_y}{\partial y} + \frac{\partial u_z}{\partial z}\right) +2C \frac{\partial u_y}{\partial x} \end{eqnarray}

(57)

\begin{eqnarray} \frac{\partial C}{\partial t}=-({\bf u}_e\cdot \nabla _e)C +A\frac{\partial u_y}{\partial x} +B\frac{\partial u_x}{\partial y} \end{eqnarray}

(58)

\begin{eqnarray} \frac{\partial q^2}{\partial t}=-({\bf u}_e\cdot \nabla _e)q^2. \end{eqnarray}

(59)

The boundary conditions on A, B, C are, as for a, b and c:

\begin{eqnarray} A(x^2+y^2=1)=B(x^2+y^2=1)=C(x^2+y^2=1)=0. \end{eqnarray}

(60)

The velocity field (39) too is expressed in the Cartesian coordinates:

\begin{eqnarray} u_x &=& -y (x^2+y^2-1)\nonumber\\ u_y &=& x(x^2 + y^2 -1) \end{eqnarray}

(61)

where u_ϕ is defined in eq. (39). The evolution equation for q² is unchanged. The initial conditions (36) can be re-written in the Cartesian formulation:

\begin{eqnarray} A&=& 2x^2y^2\sqrt{1-x^2-y^2}\nonumber\\ B&=& \frac{2}{45}\left( 61+96x^4-62y^2+21y^4+8x^2(-19+9y^2) \right)\nonumber\\ &&\times\sqrt{1-x^2-y^2}\nonumber\\ C&=&-\frac{2}{3}xy(-3+4x^2+y^2)\sqrt{1-x^2-y^2}\nonumber\\ q^2&=&A B -C^2= \frac{64}{45}x^2y^2(1-x^2-y^2)^3. \end{eqnarray}

(62)

These fields are shown in Fig. 5. Clearly the Cartesian formulation has removed the issue of regularity at the origin. We can now perform numerical integration of the eqs (56)–(58), subject to (60) and with initial conditions (62). The quadratic quantities a, b and c can be calculated via (54). As described in the previous sections, we can compare the numerical results for a and q² with their analytical solutions to assess the quality of the simulation. The numerical setup (namely the spatial discretization and time stepping) is as described in the Section 3.2 and handled modifying the relevant parameters in COMSOL Multiphysics.

Figure 5.

Initial conditions for (in order from panels a to d) A, B, C and q² calculated from the initial magnetic field of Fig. 2. Maximum and minimum values are shown at the top and at the bottom of the colour scale for each quantity.

Open in new tab Download slide

4 RESULTS

Different simulations have been run over the window 0 ≤ t ≤ π. The simulations have been performed changing the mesh refinement (mesh 1 or mesh 0) and the degree of the Lagrange polynomials. The comparison between some of these simulations is shown in Fig. 6. The magnetic energy shows that the numerical simulations have reached convergence with very small variations detectable only toward the end of the computation time interval. In Figs 7 and 8, the solution for the model L2N0 is plotted at t = 3π/4. The quadratic fields A, B and C are growing in time, each reaching values of the order 10 at the end of the simulation. Remember that q² is not a variable solved in the simulation, it is a derived value: q² = AB − C². Namely, it is a small quantity (of the order 10⁻⁴) derived from a difference from bigger quantities. In this sense it is a ’weak’ point of the simulation, difficult to keep smooth, as the snapshots of Fig. 8(d) shows. The relative amplitude of the numerical fluctuations is anyway small if we compare it to the magnitude of the quadratic fields themselves, of the order of 0.001/10 = 10⁻⁴ in the L2N0 simulation, and even smaller in the L3N0 ones. The same order of magnitude holds for the negative value in q² developing at the end of the simulation. This is the amount by which the Cauchy–Schwarz inequality is violated. Looking at A and B one can see that this is also the ratio between the most negative and positive values, expressing somehow the relative strength of positivity violation. In Fig. 9, we plot the logarithm of the absolute differences between the numerically calculated a, b, c and q² and their analytical solutions given in eqs (44)–(47). As already pointed out the departures from the analytical solutions are very small compared to the magnitudes of the quadratic fields calculated.

The magnitudes of the misfits and violations of the positivity and Cauchy–Schwarz inequality depend on the spatial discretization adopted, as is clear from Fig. 6. Provided that the numerical model has enough degrees of freedom we argue that we managed to successfully solve the kinematic problem, with acceptable departures from the expected solution. Further refinement of the spatial grid and more accurate choices of the time stepping technique might result in more refined models capable of better performance than the ones illustrated here. We were not interested in obtaining a truly advanced numerical model but only in studying the possibility of integrating the kinematic equations forward in time. Our attention was focused on the open questions highlighted in Section 3 regarding the conservation of positivity and of the Cauchy–Schwarz inequality. Interestingly we found that the eqs (25)–(27) might not be adequate to solve the kinematic problem using a local numerical strategy due to numerical noise developing at the point s = 0. The reformulation of the problem in Cartesian coordinates, leading to eqs (56)–(58) removes this issue, providing numerically acceptable solutions. In the spherical shell geometry considered by Canet et al. (2009), where the equations are not solved near the origin, this strategy might not be necessary.

Figure 6.

Performance comparison between different models in the ideal case. The nomenclature is LlNn where l is the degree of the Lagrange polynomial used as basis function and n is the number associated with the mesh. (a) Magnetic energy (48). (b) Zoom-in of magnetic energy towards the end of the computational interval. (c) Averaged difference between the numerically calculated q² and the analytical solution |$q_T^2$|⁠, according to eq. (50). (d) Minimum value of q².

Open in new tab Download slide

Figure 7.

Numerical solution for the kinematic study for the diffusion free case. Variables (in order from panels a to c) A, B and C, respectively, in the diffusion free case at the instant t = 3π/4. The flow is given in eq. (39) and the initial condition in Fig. 5. The spatial mesh is the N = 0 mesh of Table 1 and the Lagrange polynomials are of degree 2. Maximum and minimum values are shown at the top and at the bottom of the colour scale for each quantity.

Open in new tab Download slide

Numerical solution for the kinematic study for the diffusion free case. Variables (in order from panels a to d) a, b, c and q2 respectively, in the diffusion free case at the instant t = 3π/4. The fields have been calculated from the numerical solutions for A, B and C according to the formulae (54). The flow is given in eq. (39) and the initial condition in Fig. 3. The spatial mesh is the N = 0 mesh of Table 1 and the Lagrange polynomials are of degree 2. Maximum and minimum values are shown at the top and at the bottom of the colour scale for each quantity.

Figure 8.

Numerical solution for the kinematic study for the diffusion free case. Variables (in order from panels a to d) a, b, c and q² respectively, in the diffusion free case at the instant t = 3π/4. The fields have been calculated from the numerical solutions for A, B and C according to the formulae (54). The flow is given in eq. (39) and the initial condition in Fig. 3. The spatial mesh is the N = 0 mesh of Table 1 and the Lagrange polynomials are of degree 2. Maximum and minimum values are shown at the top and at the bottom of the colour scale for each quantity.

Open in new tab Download slide

Departure of a, b, c and q2 from the true analytical solutions aT, bT, cT and $q_T^2$, respectively, for the diffusion free case. Shown is the logarithm of the absolute difference as a function of space for the instant t = 3π/4. The fields a and q2 are calculated from the fields A, B and C according to the formulae (54). The spatial mesh is the N = 0 mesh of Table 1 and the Lagrange polynomials are of degree 2. Maximum and minimum values are shown at the top and at the bottom of the colour scale for each quantity.

Figure 9.

Departure of a, b, c and q² from the true analytical solutions a_T, b_T, c_T and |$q_T^2$|⁠, respectively, for the diffusion free case. Shown is the logarithm of the absolute difference as a function of space for the instant t = 3π/4. The fields a and q² are calculated from the fields A, B and C according to the formulae (54). The spatial mesh is the N = 0 mesh of Table 1 and the Lagrange polynomials are of degree 2. Maximum and minimum values are shown at the top and at the bottom of the colour scale for each quantity.

Open in new tab Download slide

5 IMPLEMENTING DIFFUSION

As discussed in Section 2.1, it is not possible to implement diffusion in eqs (25)–(27). Diffusion cannot be expressed solely in terms of a, b and c. The same limitation holds for the Cartesian formulation (56)–(58). However, mostly for numerical stability reasons, it might be desirable to retain some form of diffusion or damping in eqs (25)–(27). Here we propose to introduce a linear damping term in the induction eq. (4) and re-write it as

\begin{eqnarray} \frac{\partial {\bf B}}{\partial t}=\nabla \times ({\bf u}\times {\bf B}) - \hat{\eta }{\bf B} \end{eqnarray}

(63)

where |$\hat{\eta }$| is a positive damping constant. To write this equation in non-dimensional form in the case of the kinematic problem we introduce the following definition for the magnetic Reynolds number:

\begin{eqnarray} Rm^{-1}=\frac{r_c\hat{\eta }}{U} \end{eqnarray}

(64)

As before U = 5 × 10⁻⁴ ms⁻¹ ≃ 15 km yr⁻¹ is the characteristic value for the velocity field and r_c is the radius of the outer core. If we choose |$\hat{\eta }$| so that it gives the same decay rate as the large scale solution to the free decay problem (Jacobs 1987) we have |$\hat{\eta }^{-1} \simeq 30$| kyr and:

\begin{eqnarray} Rm \simeq 100. \end{eqnarray}

(65)

The non-dimensional induction equation becomes

\begin{eqnarray} \frac{\partial {\bf B}}{\partial t}=\nabla \times ({\bf u}\times {\bf B}) - Rm^{-1}{\bf B}. \end{eqnarray}

(66)

Eqs (25)–(27) are modified accordingly:

\begin{eqnarray} \frac{\partial a}{\partial t}=-H \left( {\bf u}\cdot \nabla _e \right)\left( \frac{a}{H} \right) + 2a\frac{\partial u_s}{\partial s} +\frac{2c}{s}\frac{\partial u_s}{\partial \phi }-2Rm^{-1}a \end{eqnarray}

(67)

\begin{eqnarray} \frac{\partial b}{\partial t} =-\frac{1}{H}\left({\bf u}\cdot \nabla _e\right)\left( Hb\right) +2sc\frac{\partial }{\partial s}\left(\frac{u_\phi }{s}\right) - 2b\frac{\partial u_s}{\partial s}-2Rm^{-1}b\nonumber\\ \end{eqnarray}

(68)

\begin{eqnarray} \frac{\partial c}{\partial t} =-\left({\bf u}\cdot \nabla _e\right) c+sa\frac{\partial }{\partial s}\left( \frac{u_\phi }{s} \right) +\frac{b}{s}\frac{\partial u_s}{\partial \phi }-2Rm^{-1}c. \end{eqnarray}

(69)

And the quantity q² satisfies the equation:

\begin{eqnarray} \frac{\partial q^2}{\partial t}=-({\bf u}_e\cdot \nabla _e)q^2 -4Rm^{-1}q^2. \end{eqnarray}

(70)

In this case too a, b and q² are positive quantities. The system (56)–(58) is also modified according to (66):

\begin{eqnarray} \frac{\partial A}{\partial t}=-({\bf u}_e\cdot \nabla _e)A +A\left(2\frac{\partial u_x}{\partial x} + \frac{\partial u_z}{\partial z} - \frac{2}{Rm}\right) +2C \frac{\partial u_x}{\partial y} \end{eqnarray}

(71)

\begin{eqnarray} \frac{\partial B}{\partial t}=-({\bf u}_e\cdot \nabla _e)B +B\left(2\frac{\partial u_y}{\partial y} + \frac{\partial u_z}{\partial z} - \frac{2}{Rm}\right) +2C \frac{\partial u_y}{\partial x} \end{eqnarray}

(72)

\begin{eqnarray} \frac{\partial C}{\partial t}=-({\bf u}_e\cdot \nabla _e)C +A\frac{\partial u_y}{\partial x} +B\frac{\partial u_x}{\partial y} - \frac{2}{Rm} C. \end{eqnarray}

(73)

Under the action of the same zonal flow (39) the analytical solutions for a, b and c are a generalization of the ideal case considered above where solutions (44), (46) and (47) are multiplied by a decay factor:

\begin{eqnarray} a_T(s,\phi ,t)=a_0\left(s,\phi -\frac{u_\phi }{s}t\right)\textrm{e}^{-2Rm^{-1}t}. \end{eqnarray}

(74)

In a similar way analytical solutions for b and c in presence of damping are obtained. The generalization of the solution for q² is similar, but the damping factor is doubled:

\begin{eqnarray} q^2_T(s,\phi ,t)=q^2_0\left(s,\phi -\frac{u_\phi }{s}t\right)\textrm{e}^{-4Rm^{-1}t}. \end{eqnarray}

(75)

We performed simulations with Rm = 10¹, 10² and 10³. As for the diffusion free case, we integrate in time eqs (71)–(73) and we calculate the quantities a and q² from A, B and C. The results for models L1N0 and L2N0 are shown in Figs 10 and 11. As expected the magnetic energy growth with time is slower with lower values of Rm, corresponding to higher values of the damping coefficient. However, the difference is modest for Rm up to 100. The solution is visually similar to the diffusion free case, with the damping introducing the expected effect of lowering the magnitude of both the kinematic solution and of the numerical artefacts. The misfit between q² and |$q^2_T$| shown in Figs 10(c) and 11(c) decreases with increasing damping because all fields are damped, and so is the error. The relative errors and the relative magnitude of the negative features with respect to the amplitude of the solution is similar to the diffusion free case and so is the location of the artefacts.

Figure 10.

Comparison between different kinematic numerical solutions for the L1N0 model in presence of damping. Different values of Rm are considered, Rm = ∞ corresponding to case with no damping. (a) Magnetic energy (48). (b) Zoom in of magnetic energy toward the end of the computational interval. (c) Averaged difference between the numerically calculated q² and the analytical solution |$q_T^2$|⁠, according to (50). (d) Minimum value of q².

Open in new tab Download slide

Figure 11.

Comparison between different kinematic numerical solutions for the L2N0 model in the presence of damping. Different values of Rm are considered, Rm = ∞ corresponding to case with no damping. (a) Magnetic energy (48). (b) Zoom in of magnetic energy towards the end of the computational interval. (c) Averaged difference between the numerically calculated q² and the analytical solution |$q_T^2$|⁠, according to eq. (50). (d) Minimum value of q².

Open in new tab Download slide

6 DISCUSSION

The preservation of positivity for the quantities a, b and q² is of fundamental importance for the QG model (24)–(27). It has been pointed out how on this ground, a spectral expansion has been ruled out to discretize the quadratic variables a, b and c on the equatorial plane. It is however possible to further manipulate eqs (25)–(27) so that positivity is preserved independently of the details of the numerical algorithm. For example under the following representation:

\begin{eqnarray} a=\left\lbrace B_s^2\right\rbrace =\textrm{e}^\gamma \sqrt{c^2+q^2} \end{eqnarray}

(76)

\begin{eqnarray} b=\left\lbrace B_\phi ^2\right\rbrace = \textrm{e}^{-\gamma }\sqrt{c^2+q^2} \end{eqnarray}

(77)

the Cauchy–Schwarz inequality (31) is identically satisfied. Evolution equations for γ, q and c can be derived introducing (76), (77) in (25)–(27). However the resulting system of equation is more involved than (25)–(27) and we chose not to investigate its properties further.

We now discuss the relationship between a, b, c and the physical magnetic field in the core. From the quantity a we can calculate the cylindrical radial propagation speed of torsional Alfvén waves in the core (Braginsky 1970), which, in non-dimensional units is

\begin{eqnarray} V_s = \sqrt{\frac{1}{4 \pi H}\int _{-H}^H \int _0^{2\pi } B_s^2 {\rm d}z \ {\rm d}\phi }= \sqrt{\frac{1}{4 \pi H}\int _0^{2\pi } a\, {\rm d}\phi }. \end{eqnarray}

(78)

This is the average of |$B_s^2$| on cylinders aligned with the rotation axis (the geostrophic cylinders). Different papers have been devoted with the calculation of V_s based on torsional waves observation (Hide et al.2000; Gillet et al.2010, 2015) as its quantification offers a precious constraint on the intensity of the magnetic field inside the core. The latest results suggest a value of 3 mT for the average intensity of the magnetic field in the interior of the core, based on estimates of V_s. In a similar way the quantity b can be quantified from the observation of slow hydrodynamic oscillations propagating in the azimuthal direction around the core (Hide 1966; Malkus 1967; Canet et al.2014; Hori et al.2015). However the detection of such oscillations is more challenging than for torsional oscillations and an observational estimate of b is still missing. The quantity c, to the best of our knowledge, is not directly related to any detectable oscillation and in the context of the QG model derived here, it is a purely mathematical quantity that we need to close the system (24)–(27). It is needed to complete the coupling of the momentum equation with the induction equation and is related to a and b via the Cauchy-Schwarz inequality (31).

7 CONCLUSIONS

In this study we presented a kinematic study of a QG model that is tailored to the study of the interannual and decadal dynamics of the Earth’s outer core. In Section 2, we discussed how the QG formalism arises from the 3-D magneto-hydrodynamic governing equations. Under the hypothesis of fast rotation, the columnar flow assumption is the key to developing realistic numerical models that have the potential to reach parameter regimes much closer to the real Earth than 3-D numerical models. The key operation is the projection of the governing equations on the equatorial disc. This operation is far from trivial when applied to the momentum equation, as it requires assumptions and delicate mathematic manipulations on the terms describing the body forces (Lorentz force and buoyancy in geodynamo simulations). The formulation presented in Sec-tion 2.1 has the advantage, compared to the QG formulation of Canet et al. (2014) and Labbé et al. (2015), of handling the magnetic field under very general assumptions. This is achieved by representing the magnetic field via the quadratic quantities |$B_s^2$|⁠, |$B_\phi ^2$| and B_sB_ϕ projected on the equatorial disc. We call this the quadratic formulation. The resulting equations are however much more complicated and magnetic diffusion is impossible to accommodate. Compared to the vertically averaged formulation of Canet et al. (2009), the equations presented here, based on a vertically integrated formulation, have the advantage of providing a natural boundary condition for the quantities a, b and c. The equations for the vertically averaged formalism can be derived from (24)–(27) dividing them by 2H.

The QG model developed here, although based on the formulation proposed in Canet et al. (2009), has never been integrated in time. Our kinematic formulation highlighted some deficiencies in the equations in polar coordinates that we removed when developing a novel Cartesian formulation.

We investigated the effect of damping as an alternative to diffusion. The simple damping term introduced in the induction eq. (66) leads to a system of equation that is now closed in a, b and c. We find that this simple solution is efficient in lowering the amplitude by which the positivity of a, b and q² is violated.

In future studies it is necessary to characterize the full system (24)–(27) where the Lorentz force provides a link between the evolutions of the magnetic field and the fluid flows. A first step could be the calculation of magnetohydrodynamic normal modes in a similar study to those of Canet et al. (2014) and Labbé et al. (2015). In such studies care has to be taken in testing the effect of the approximations that led to the derivation of (24), namely the assumption that surface terms can be neglected compared to the vertically integrated ones. The discarded surface terms could be of importance at the equator, where the height H of the domain vanishes and singularities may compromise the derivation if not treated carefully.

Finally we point out how the use of the axial vorticity eq. (18) could not be the best option to describe the evolution of the stream function Ψ under the columnar flow approximation (15). In Labbé et al. (2015), an alternative formulation is proposed that improves the predictions of the momentum equation in the equatorial regions s ≃ 1. This technique allowed Labbé and co-authors to calculate hydrodynamic solutions that are in better agreement with 3-D calculations than previously known QG calculations (Canet et al.2014). In future work, the axial vorticity eq. (24) should therefore be replaced with an equation that is derived following the methodology of Labbé et al. (2015).

The final goal of the QG studies such the one presented here is the derivation of a numerical model capable of simulating the Earth’s core dynamics in a realistic parameter regimes. Such a model will allow the implementation of data assimilation systems that, making use of modern geomagnetic observations, could open a window on the core that has never been achieved with 3-D models. A first step toward this goal is to perform closed loop data assimilation experiment the kinematic system (25)–(27), similarly to the ones presented in Canet et al. (2009). After the momentum eq. (24) has been confidently integrated forward in time, data assimilation experiments similar to the ones performed in Fournier et al. (2013) can be attempted.

Table 1.

Open in new tab

Types of meshes used in COMSOL Multiphysics for the kinematic study. The domain to be discretized is the equatorial disc. ‘tri’ indicates the number of elements used to discretize the bulk of the domain. ‘edg’ indicates the number of edge elements, or curved elements used to represent the boundary of the domain.

Mesh identifier (N)	Mesh name	tri	edg
0	Custom 1	192 732	900
1	Extremely fine	24 934	316
2	Extra fine	6534	160

Table 1.

Open in new tab

Types of meshes used in COMSOL Multiphysics for the kinematic study. The domain to be discretized is the equatorial disc. ‘tri’ indicates the number of elements used to discretize the bulk of the domain. ‘edg’ indicates the number of edge elements, or curved elements used to represent the boundary of the domain.

Mesh identifier (N)	Mesh name	tri	edg
0	Custom 1	192 732	900
1	Extremely fine	24 934	316
2	Extra fine	6534	160

Acknowledgments

This study has been supported by the ETH grant number ETH-1412-1 and by the SNF grant number 200020_143596. We wish to thank Dr Dominique Jault (Université Grenoble Alpes) and Dr Elisabeth Canet (former ETH) for fruitful discussions and useful suggestions. We also wish to express our gratitude to an anonymous reviewer whose comments and suggestions greatly improved the quality of the original manuscript.

1

These equations correct errors in the equations given in Jault & Finlay (2015). The vertically averaged formulation employed there can be retrieved by dividing eqs (25)–(27) by 2H. Note that eqs (66) of Jault & Finlay (2015) would be formally correct for the vertically integrated formulation of the present study. The interested reader is referred to section 4.2.4 of Maffei (2016) for additional details.

2

https://www.comsol.com

REFERENCES

Alfvén

H.

,

1942

.

Existence of electromagnetic-hydrodynamic waves

,

Nature

,

150

(

3805

),

405

–

406

.

Aubert

J.

,

Gillet

N.

,

Cardin

P.

,

2003

.

Quasigeostrophic models of convection in rotating spherical shells

,

Geochem. Geophys. Geosyst.

,

4

(

7

),

1052

,

doi:10.1029/2002GC000456

.

Aubert

J.

,

Labrosse

S.

,

Poitou

C.

,

2009

.

Modelling the palaeo-evolution of the geodynamo

,

Geophys. J. Int.

,

179

(

3

),

1414

–

1428

.

Aubert

J.

,

Gastine

T.

,

Fournier

A.

,

2017

.

Spherical convective dynamos in the rapidly rotating asymptotic regime

,

J. Fluid Mech.

,

813

,

558

–

593

.

Braginsky

S.

,

1970

.

Torsional magnetohydrodynamic vibrations in the Earth’s core and variations in day length

,

Geomagn. Aeron.

,

10

,

1

–

8

.

Google Scholar

OpenURL Placeholder Text

WorldCat

Calkins

M.A.

,

Julien

K.

,

Tobias

S.M.

,

Aurnou

J.M.

,

2015

.

A multiscale dynamo model driven by quasi-geostrophic convection

,

J. Fluid Mech.

,

780

,

143

–

166

.

Canet

E.

,

Fournier

A.

,

Jault

D.

,

2009

.

Forward and adjoint quasi-geostrophic models of the geomagnetic secular variation

,

J. geophys. Res.

,

114

,

B11101

,

doi:10.1029/2008JB006189

.

Canet

E.

,

Finlay

C.C.

,

Fournier

A.

,

2014

.

Hydromagnetic quasi-geostrophic modes in rapidly rotating planetary cores

,

Phys. Earth planet. Inter.

,

229

,

1

–

15

.

Cardin

P.

,

Olson

P.

,

1994

.

Chaotic thermal convection in a rapidly rotating spherical shell: consequences for flow in the outer core

,

Phys. Earth planet. Inter.

,

82

(

3–4

),

235

–

259

.

Cardin

P.

,

Olson

P.

,

2007

.

8.11—Experiments on Core Dynamics

, in

Treatise on Geophysics

, pp.

319

–

343

, ed.

Schubert

G.

,

Elsevier

.

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Christensen

U.R.

,

Wicht

J.

,

2015

.

8.10—Numerical Dynamo Simulations

, in

Treatise on Geophysics

, 2nd edn, pp.

245

–

277

, ed.

Schubert

G.

,

Elsevier

.

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Finlay

C.C.

,

2011

.

Waves in the presence of magnetic fields, rotation and convection

, in

Dynamos: Lecture Notes of the Les Houches Summer School 2007, vol. 88 of Les Houches Summer School Proceedings

, pp.

403

–

450

, eds

Cardin

P.

,

Cugliandolo

L.F.

,

Elsevier

.

Finlay

C.C.

,

Jackson

A.

,

2003

.

Equatorially Dominated Magnetic Field Change at the Surface of Earth’s Core

,

Science

,

300

(

5628

),

2084

–

2086

.

Fournier

A.

,

Nerger

L.

,

Aubert

J.

,

2013

.

An ensemble Kalman filter for the time-dependent analysis of the geomagnetic field

,

Geochem. Geophys. Geosyst.

,

14

,

4035

–

4043

.

Gillet

N.

,

Jones

C.A.

,

2006

.

The quasi-geostrophic model for rapidly rotating spherical convection outside the tangent cylinder

,

J. Fluid Mech.

,

554

,

343

–

369

.

Gillet

N.

,

Pais

M.A.

,

Jault

D.

,

2009

.

Ensemble inversion of time-dependent core flow models

,

Geochem. Geophys. Geosyst.

,

10

,

Q06004

,

doi:10.1029/2008GC002290

.

Gillet

N.

,

Jault

D.

,

Canet

E.

,

Fournier

A.

,

2010

.

Fast torsional waves and strong magnetic field within the Earth’s core

,

Nature

,

465

(

7294

),

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.

,

202–203

,

78

–

88

.

Google Scholar

OpenURL Placeholder Text

WorldCat

Gillet

N.

,

Jault

D.

,

Finlay

C.C.

,

2015

.

Planetary gyre, time-dependent eddies, torsional waves, and equatorial jets at the Earth’s core surface

,

J. geophys. Res.

,

120

(

6

),

3991

–

4013

.

Glatzmaier

G.A.

,

2002

.

Geodynamo Simulations—How Realistic Are They?

,

Annu. Rev. Earth Planet. Sci.

,

30

(

1

),

237

–

257

.

Greenspan

H.P.

,

1968

.

The Theory of Rotating Fluids

,

CUP Archive

.

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Guervilly

C.

,

Cardin

P.

,

2016

.

Subcritical convection of liquid metals in a rotating sphere using a quasi-geostrophic model

,

J. Fluid Mech.

,

808

,

61

–

89

.

Hide

R.

,

1966

.

Free hydromagnetic oscillations of the Earth’s core and the theory of the geomagnetic secular variation

,

Phil. Trans. R. Soc. A

,

259

(

1107

),

615

–

647

.

Hide

R.

,

Boggs

D.H.

,

Dickey

J.O.

,

2000

.

Angular momentum fluctuations within the Earth’s liquid core and torsional oscillations of the core–mantle system

,

Geophys. J. Int.

,

143

(

3

),

777

–

786

.

Hori

K.

,

Jones

C.A.

,

Teed

R.J.

,

2015

.

Slow magnetic Rossby waves in the Earth’s core

,

Geophys. Res. Lett.

,

42

(

16

),

6622

–

6629

.

Hulot

G.

,

Sabaka

T.J.

,

Olsen

N.

,

Fournier

A.

,

2015

.

5.02 - The Present and Future Geomagnetic Field

, in

Treatise on Geophysics

, 2nd edn, pp.

33

–

78

, ed.

Schubert

G.

,

Elsevier

.

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Jackson

A.

,

2003

.

Intense equatorial flux spots on the surface of the Earth’s core

,

Nature

,

424

(

6950

),

760

–

763

.

Jackson

A.

,

Finlay

C.

,

2015

.

5.05—Geomagnetic secular variation and its applications to the core

, in

Treatise on Geophysics

, 2nd edn, pp.

137

–

184

, ed.

Schubert

G.

,

Elsevier

.

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Jacobs

J.A.

,

1987

.

Geomagnetism

,

vol. 2

,

Academic Press

.

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Jault

D.

,

2008

.

Axial invariance of rapidly varying diffusionless motions in the Earth’s core interior

,

Phys. Earth planet. Inter.

,

166

(

1–2

),

67

–

76

.

Jault

D.

,

Finlay

C.C.

,

2015

.

8.09—Waves in the core and mechanical core–mantle interactions

, in

Treatise on Geophysics

, 2nd edn, pp.

225

–

244

, ed.

Schubert

G.

,

Elsevier

.

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Jones

C.

,

2007

.

Thermal and compositional convection in the outer core

, in

Treatise on Geophysics

, pp.

131

–

185

, ed.

Schubert

G.

,

Elsevier

.

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Labbé

F.

,

Jault

D.

,

Gillet

N.

,

2015

.

On magnetostrophic inertia-less waves in quasi-geostrophic models of planetary cores

,

Geophys. Astrophys. Fluid Dyn.

,

109

(

6

),

587

–

610

.

Lewis

H.R.

,

Bellan

P.M.

,

1990

.

Physical constraints on the coefficients of Fourier expansions in cylindrical coordinates

,

J. Math. Phys.

,

31

(

11

),

2592

–

2596

.

Maffei

S.

,

2016

.

Quasi-geostrophic models for fast dynamics in the earth’s outer core

,

PhD thesis

,

ETH Zurich

.

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Malkus

W.V.R.

,

1967

.

Hydromagnetic planetary waves

,

J. Fluid Mech.

,

28

(

04

),

793

–

802

.

Pais

M.A.

,

Jault

D.

,

2008

.

Quasi-geostrophic flows responsible for the secular variation of the Earth’s magnetic field

,

Geophys. J. Int.

,

173

(

2

),

421

–

443

.

Roberts

P.H.

,

2015

.

8.03—Theory of the geodynamo

, in

Treatise on Geophysics

, 2nd edn, pp.

57

–

90

, ed.

Schubert

G.

,

Elsevier

.

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Schaeffer

B.

,

Cardin

P.

,

2005

.

Quasi-geostrophic model of the instabilities of the stewartson layer

,

Phys. Fluids

,

17

,

104

–

111

.

Schaeffer

N.

,

Cardin

P.

,

2006

.

Quasi-geostrophic kinematic dynamos at low magnetic Prandtl number

,

Earth planet. Sci. Lett.

,

245

(

3–4

),

595

–

604

.

Schaeffer

N.

,

Silva

E.L.

,

Pais

M.A.

,

2016

.

Can core flows inferred from geomagnetic field models explain the Earth’s dynamo?

,

Geophys. J. Int.

,

204

(

2

),

868

–

877

.

Schaeffer

N.

,

Jault

D.

,

Nataf

H.-C.

,

Fournier

A.

,

2017

.

Geodynamo simulations with vigorous convection and low viscosity

,

arXiv:1701.01299

.

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Teed

R.J.

,

Jones

C.A.

,

Tobias

S.M.

,

2014

.

The dynamics and excitation of torsional waves in geodynamo simulations

,

Geophys. J. Int.

,

196

,

724

–

735

.

Williams

P.D.

,

Read

P.L.

,

Haine

T.W.N.

,

2010

.

Testing the limits of quasi-geostrophic theory: application to observed laboratory flows outside the quasi-geostrophic regime

,

J. Fluid Mech.

,

649

,

187

–

203

.

Download all slides

Month:	Total Views:
June 2017	63
July 2017	15
August 2017	20
September 2017	18
October 2017	18
November 2017	9
December 2017	16
January 2018	8
February 2018	7
March 2018	10
April 2018	11
May 2018	2
June 2018	7
August 2018	6
September 2018	7
October 2018	8
November 2018	3
December 2018	5
January 2019	4
February 2019	1
March 2019	2
April 2019	3
May 2019	1
June 2019	4
August 2019	2
September 2019	2
October 2019	3
November 2019	7
December 2019	5
February 2020	7
March 2020	4
April 2020	6
May 2020	2
June 2020	15
July 2020	2
September 2020	6
October 2020	8
November 2020	9
December 2020	2
January 2021	8
February 2021	7
March 2021	9
April 2021	13
May 2021	3
June 2021	8
July 2021	12
August 2021	6
September 2021	2
October 2021	8
November 2021	11
December 2021	4
January 2022	5
February 2022	5
March 2022	1
April 2022	9
May 2022	5
June 2022	1
July 2022	13
August 2022	5
September 2022	6
October 2022	4
November 2022	5
December 2022	4
January 2023	2
February 2023	4
March 2023	5
May 2023	4
June 2023	5
July 2023	3
August 2023	4
September 2023	3
October 2023	2
November 2023	9
December 2023	8
January 2024	10
February 2024	17
March 2024	23
April 2024	8
May 2024	12
June 2024	5
July 2024	14
August 2024	6
September 2024	4
October 2024	3
November 2024	3
December 2024	8
January 2025	11
February 2025	3
March 2025	8
April 2025	7
May 2025	3

Article Contents

Kinematic validation of a quasi-geostrophic model for the fast dynamics in the Earth’s outer core

Summary

1 INTRODUCTION

2 GOVERNING EQUATIONS

2.1 Columnar flow assumption and vertically integrated equations

3 VALIDATION UNDER THE KINEMATIC APPROXIMATION

3.1 Initial field and prescribed flow

3.2 Numerical experiments setup

3.3 Cartesian formulation

4 RESULTS

5 IMPLEMENTING DIFFUSION

6 DISCUSSION

7 CONCLUSIONS

Acknowledgments

REFERENCES

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

Article Contents

Kinematic validation of a quasi-geostrophic model for the fast dynamics in the Earth’s outer core Free

Summary

1 INTRODUCTION

2 GOVERNING EQUATIONS

2.1 Columnar flow assumption and vertically integrated equations

3 VALIDATION UNDER THE KINEMATIC APPROXIMATION

3.1 Initial field and prescribed flow

3.2 Numerical experiments setup

3.3 Cartesian formulation

4 RESULTS

5 IMPLEMENTING DIFFUSION

6 DISCUSSION

7 CONCLUSIONS

Acknowledgments

REFERENCES

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

This Feature Is Available To Subscribers Only

Kinematic validation of a quasi-geostrophic model for the fast dynamics in the Earth’s outer core