The magnetic non-hydrostatic shallow-water model

Abstract

We consider the dynamics of a set of reduced equations describing the evolution of a magnetised, rotating stably stratified fluid layer, atop a stagnant dense, perfectly conducting layer. We consider two closely related models. In the first, the layer has, above it, relatively light fluid where the magnetic pressure is much larger than the gas pressure, and the magnetic field is largely force-free. In the second model, the magnetic field is constrained to lie within the dynamical layer by the implementation of a model diffusion operator for the magnetic field. The model derivation proceeds by assuming that the horizontal velocity and the horizontal magnetic field are independent of the vertical coordinate, whilst the vertical components in the layer have a linear dependence on height. The full system comprises evolution equations for the magnetic field, horizontal velocity and height field together with a linear elliptic equation for the vertically integrated non-hydrostatic pressure. In the magneto-hydrostatic limit, these equations simplify to equations of shallow-water type. Numerical solutions for both models are provided for the fiducial case of a Gaussian vortex interacting with a magnetic field. The solutions are shown to differ negligibly. We investigate how the interaction of the vortex changes in response to the magnetic Reynolds number ${Rm}$, the Rossby deformation radius $L_D$, and a Coriolis buoyancy frequency ratio $f/N$ measuring the significance of non-hydrostatic effects. The magneto-hydrostatic limit corresponds to $f/N\to 0$.

1. Introduction

The modelling of stably stratified regions can play an important role in explaining the dynamics of many geophysical and astrophysical objects. On Earth, examples include the stratosphere, the oceans and the outer layer of the Earth's fluid outer core (sometimes known as the Earth's inner ocean; Braginsky 1999). Astrophysical examples include the radiative zones of certain stars and the outermost layers of certain solar system planets and exoplanets. There are many circumstances in which the stably stratified layer is electrically conducting and the fluid interacts with magnetic fields. These include the solar tachocline (see e.g. Tobias 2005; Christensen-Dalsgaard & Thompson 2007; Spiegel 2007), the region of stably stratified, sheared, magnetised plasma at the base of the solar convection zone, which is believed to be the seat of the solar dynamo, and solar activity that takes the form of active regions at the solar surface, as well as the outer layers of certain exoplanets known as hot Jupiters (see e.g. Rogers & Komacek 2014; Yadav & Thorngren 2017; Fortney, Dawson & Komacek 2021).

Modelling of these stable layers, which are often thin – in the sense that the horizontal length scales for the dynamics are much longer than typical radial (or vertical) length scales – often involves the derivation and study of reduced models. These reduced models are commonly derived from the full ‘parent’ system of equations by making use of a small parameter (sometimes the ratio of the vertical to horizontal length scales) or by considering leading-order balances (e.g. geostrophic and hydrostatic). Of these models, in hydrodynamics the most famous is the hydrostatic single-layer shallow-water model. Here, one considers a free-surface flow of uniform density in the presence of gravity and, for the cases of geophysical and astrophysical relevance, background rotation.

The shallow-water model can be derived by assuming that the horizontal flow is independent of height (columnar motion), and additionally making the hydrostatic approximation. The latter is valid for disturbances with large horizontal scales, but becomes increasingly inappropriate for small-scale disturbances in the horizontal (see e.g. Dritschel & Jalali (2019, 2020), and references therein). An approach with wider applicability is to relax the hydrostatic approximation while retaining the assumption of height independence of the horizontal velocity. This hydrodynamic extension, often termed the ‘Green–Naghdi’ model (Green & Naghdi 1976), has been considered by a number of authors (for a short review, see Jalali & Dritschel 2021) and has been shown explicitly to be substantially more accurate than the shallow-water model in describing the dynamics of shallow rotating free-surface flows (Dritschel & Jalali (2020); see also Nadiga, Margolin & Smolarkiewicz (1996) for unidirectional non-rotating flow over topography).

Just as for the hydrodynamic case, magneto-hydrodynamic (MHD) models can be reduced from the full three-dimensional system. The most severe approximation for stably stratified systems is to consider a strictly two-dimensional system with a magnetic field and flow constrained to lie in a plane. Here, it has been shown that even an extremely weak mean magnetic field can lead to a significant change in the dynamics through the action of the Lorentz force breaking the material conservation of potential vorticity (see e.g. Tobias, Diamond & Hughes 2007; Dritschel, Diamond & Tobias 2018). However, this geometry is overly restrictive, and more sophisticated reduced models are needed.

Significant progress was made by Gilman (2000), who extended the shallow-water equations to include MHD effects. He considered a perfectly conducting thin layer of magnetised fluid, which is an excellent test-bed for studying waves (see e.g. Schecter, Boyd & Gilman 2001; Zaqarashvili, Oliver & Ballester 2009; Márquez-Artavia, Jones & Tobias 2017), joint instabilities of the magnetic field and differential rotation (Dikpati, Gilman & Rempel 2003), and the nonlinear evolution of such instabilities (see e.g. Dikpati et al. 2021) in magnetised, rotating stratified layers. However, this formalism does not allow any self-consistent procedure for including the effects of magnetic diffusion and so is less useful for nonlinear turbulence problems. We note here that in both astrophysical settings discussed above, the magnetic Prandtl number ${Pm} = \nu /\eta$ (where $\nu$ is the viscous diffusion, and $\eta$ is the magnetic diffusivity) is small, so magnetic diffusion is the primary dissipative process. The magnetic field therefore dissipates on a length scale much larger than that of the velocity – motivating the construction of models where the magnetic diffusion is included explicitly but the viscosity is not (Dritschel et al. 2018). It is possible, however, to consider models where the diffusive layer is bounded above by an idealised perfectly conducting fluid, allowing surface currents to form at the interface. In this case, one can devise dissipative terms that enforce the condition that the field remains tangent to the free surface in the shallow-water limit (Gilbert, Griffiths & Hughes 2022); this is one of the models that we will consider here.

The MHD shallow-water equations were subsequently extended to include dispersive effects by Dellar (2003). These equations are the magnetised versions of the Green–Naghdi equations, though again only the case with no magnetic diffusivity (and the properties of waves therein) was considered. A generalisation of this model has been developed recently by Alonso-Orán (2020), though seemingly independent of Dellar (2003).

In the present paper, we derive reduced models of rotating, stratified magnetised fluid layers that allow for the presence of magnetic diffusion (though these models are still inviscid). In § 2, we set up the model, with a full derivation. Numerical results for the various models are described in § 3, with conclusions offered in § 4.

2. Set-up of the model and derivation

We consider a magnetised, conducting system comprising three layers in a local Cartesian domain. The bottom layer ($z<0$) is either a perfectly conducting solid or a very dense perfectly conducting fluid, with density $\rho _b$, such that the interface between it and the fluid layer above remains fixed at $z=0$. The layer of primary interest is of uniform density $\rho _l \ll \rho _b$ lying above this flat bottom at $z=0$ and extending to a free surface at $z=h(x,y,t)$. A final layer of much smaller density $\rho _u$ (and hence negligible gas pressure) lies above the free surface. The fluid is assumed to be rotating uniformly about the vertical $z$ axis at rate $\varOmega$, and gravity $g$ acts downwards in $z$.

Let $\rho$ denote a fiducial density. Then the governing three-dimensional incompressible inviscid, MHD equations that we consider in the middle layer are given by (see e.g. Chandrasekhar 1981)

(2.1)$$\begin{gather} \rho\left(\frac{\partial{{\boldsymbol{u}}}}{\partial{t}}+{\boldsymbol{u}}\boldsymbol{\cdot}\boldsymbol{\nabla}{\boldsymbol{u}}+f{\boldsymbol{k}}\times{\boldsymbol{u}}+g{\boldsymbol{k}}\right) ={-}\boldsymbol{\nabla} p + {\boldsymbol{j}}\times{\boldsymbol{B}}, \end{gather}$$

(2.2)$$\begin{gather}\frac{\partial{{\boldsymbol{B}}}}{\partial{t}} = \boldsymbol{\nabla}\times({\boldsymbol{u}}\times{\boldsymbol{B}}) + \eta\,\nabla^2{\boldsymbol{B}}, \end{gather}$$

(2.3a,b)$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{u}} = 0,\quad \boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{B}} = 0, \end{gather}$$

where ${\boldsymbol {u}}$ is the velocity field, ${\boldsymbol {B}}$ is the magnetic field, ${\boldsymbol {j}}=\boldsymbol {\nabla }\times {\boldsymbol {B}}/\mu$ is the current density (with $\mu$ the magnetic permeability, assumed constant throughout), $p$ is the pressure, ${f=2\varOmega}$ is the Coriolis frequency, ${\boldsymbol {k}}$ is the (upwards) vertical unit vector, and $\eta =1/(\sigma \mu )$ is the magnetic diffusivity (with $\sigma$ the conductivity, also assumed constant).

We next scale pressure and the magnetic variables to remove non-essential constants. If the substitutions

(2.4a–c)\begin{equation} \frac{p}{\rho_l} \to p,\quad \frac{{\boldsymbol{B}}}{\sqrt{\rho_l\mu}} \to {\boldsymbol{B}}, \quad \frac{\mu{\boldsymbol{j}}}{\sqrt{\rho_l\mu}} \to {\boldsymbol{j}} \end{equation}

are made, then the magnetic field is measured in terms of the Alfvén velocity, while ${\boldsymbol {j}}=\boldsymbol {\nabla }\times {\boldsymbol {B}}$ is analogous to the vorticity relation $\boldsymbol {\omega }=\boldsymbol {\nabla }\times {\boldsymbol {u}}$. Equations (2.2) and (2.3a,b) are unchanged, while the momentum equation (2.1) is rescaled to

(2.5)\begin{equation} \frac{\partial{{\boldsymbol{u}}}}{\partial{t}}+{\boldsymbol{u}}\boldsymbol{\cdot}\boldsymbol{\nabla}{\boldsymbol{u}}+f{\boldsymbol{k}}\times{\boldsymbol{u}}+g{\boldsymbol{k}} ={-}\boldsymbol{\nabla} p + {\boldsymbol{j}}\times{\boldsymbol{B}}. \end{equation}

Traditional shallow-water models (both hydrodynamic and magneto-hydrodynamic) are derived by making two assumptions. They suppose that the layer of fluid $0 \le z \le h(x,y,t)$ has a characteristic horizontal scale $L$ much larger than the fluid depth $h$. Further assuming that the horizontal velocity is independent of $z$, and making the hydrostatic approximation (ignoring the vertical acceleration ${\rm D}{w}/{\rm D}{t}$), leads to the well-known shallow-water model in the absence of a magnetic field first derived by Saint-Venant (1871). More recently, this model (with both assumptions on spatial scale and magneto-hydrostatic balance) was extended to a perfectly conducting magnetised fluid by Gilman (2000). There has been much development of these MHD shallow-water models, including the investigation of instabilities and wave phenomena (see e.g. Schecter et al. 2001; Zaqarashvili et al. 2009; Márquez-Artavia et al. 2017).

In the hydrodynamical context, the shallow-water models were extended to include non-hydrostatic effects by Serre (1953) but are often attributed to Green & Naghdi (1976), who explained how the equations can be derived simply by vertically averaging the parent three-dimensional system when the horizontal velocity is assumed to be independent of height (see also the discussion in § 3 of Dritschel & Jalali 2020). Notably, the vertically averaged equations possess all of the material and integral invariants of the parent hydrodynamic system, and moreover these invariants are simply the vertical averages of their three-dimensional counterparts (Miles & Salmon 1985; Dritschel & Jalali 2020). Direct comparisons with the parent three-dimensional free-surface model in Dritschel & Jalali (2020) demonstrate that the non-hydrostatic model (including rotation) is substantially more accurate than the hydrostatic one, despite the fact that neither model captures deep-water waves whose velocity field varies strongly with depth (see e.g. §§ 4 and figure 16 in Dritschel & Jalali 2019). Similar improvements were reported by Nadiga et al. (1996) in the non-rotating case for unidirectional flow (independent of one coordinate) over topography.

The magnetised version of the non-hydrostatic model with no magnetic diffusivity was first derived formally by Dellar (2003), who also demonstrated that the wave properties of this model capture more accurately those of the full three-dimensional parent system, especially at small scales. In this paper, we derive and study the dynamics of two systems with two different models of magnetic diffusion. In the first, resistive diffusion allows the magnetic field to leave and enter the domain, but in the layer above, the magnetic pressure dominates the gas pressure, which can be assumed negligible – the layer above is assumed to reach rapidly a nearly force-free configuration. In the second scenario, an approximate model of diffusion is considered that does not allow flux to leave the layer (Gilbert et al. 2022), and the fluid above is again considered to have an approximately zero pressure acting on the free surface. We believe that the first of these models is more relevant to the outer layers of planets and stars that lie below a highly compressible magnetised atmosphere dominated by the magnetic pressure (high plasma $\beta$). The second of these models might be preferred for the interior of stars (in regions such as the solar tachocline). As we will see, the difference between these models is negligible in the astrophysically relevant regime of high magnetic Reynolds number $Rm$, with the difference between them appearing only at ${O} (Rm^{-1})$.

2.1. Model derivation

Following Dellar (2003) and Dritschel & Jalali (2020), hereafter all three-dimensional vectors and vector operators are indicated by a subscript 3, while corresponding two-dimensional quantities have no subscript. Thus the velocity field is ${\boldsymbol {u}}_3=({\boldsymbol {u}},w)$, where $w$ is the vertical component. The magnetic field is ${\boldsymbol {B}}_3=({\boldsymbol {B}},B_z)$, where $B_z$ is the vertical component. The horizontal components of ${\boldsymbol {u}}$ are $(u,v)$, while those of ${\boldsymbol {B}}$ are $(B_x,B_y)$. The horizontal position vector is ${\boldsymbol {x}}$, while the vertical coordinate is $z$. In our first derivation we focus on the model with regular magnetic diffusion given by $\eta \,\nabla ^2_3{\boldsymbol {B}}_3$.

To simplify notation, we use

(2.6a,b)\begin{equation} D\equiv\frac{\partial{}}{\partial{t}}+{\boldsymbol{u}}\boldsymbol{\cdot}\boldsymbol{\nabla},\quad D_3\equiv\frac{\partial{}}{\partial{t}}+{\boldsymbol{u}}_3\boldsymbol{\cdot}\boldsymbol{\nabla}_3, \end{equation}

to denote material derivatives in two dimensions and three dimensions, respectively.

In this notation, the scaled three-dimensional equations (2.2), (2.3a,b) and (2.5) read

(2.7)$$\begin{gather} D_3{\boldsymbol{u}}_3+f{\boldsymbol{k}}\times{\boldsymbol{u}}_3+g{\boldsymbol{k}} ={-}\boldsymbol{\nabla}_3p+{\boldsymbol{j}}_3\times{\boldsymbol{B}}_3, \end{gather}$$

(2.8)$$\begin{gather}D_3{\boldsymbol{B}}_3-{\boldsymbol{B}}_3\boldsymbol{\cdot}\boldsymbol{\nabla}_3{\boldsymbol{u}}_3 = \eta\,\nabla^2_3{\boldsymbol{B}}_3, \end{gather}$$

(2.9a,b)$$\begin{gather}\boldsymbol{\nabla}_3\boldsymbol{\cdot}{\boldsymbol{u}}_3 = 0,\quad \boldsymbol{\nabla}_3\boldsymbol{\cdot}{\boldsymbol{B}}_3 = 0, \end{gather}$$

where the induction equation (2.8) has been rewritten making use of (2.9a,b). The boundary conditions are $w=B_z=0$ on $z=0$, $w=Dh$ on $z=h$, and $p=0$ on ${z=h}$. There is much literature on the construction of the external field for $z>h$ (see Appendix B). In practice, however, we do not need to do this. The model derived below is entirely independent of its continued solution outside of the middle layer ($0\leq z\leq h$), provided that the gas pressure in the uppermost layer is considered negligible (see Appendix A).

We are guided by Green & Naghdi (1976), who argued that it is necessary to assume only that the horizontal velocity ${\boldsymbol {u}}$ is independent of $z$, then average vertically (see also Miles & Salmon 1985). No perturbation expansion is needed, and in fact, Green & Naghdi (1976) warn against such expansions as they do not guarantee conservation (see § 3.1 in Dritschel & Jalali 2020). Here, we simply make the same assumption for the magnetic field, taking the horizontal part ${\boldsymbol {B}}$ to be independent of $z$.

Starting with incompressibility, $\boldsymbol {\nabla }_3\boldsymbol {\cdot }{\boldsymbol {u}}_3 = 0$, we have

(2.10)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{u}}+\frac{\partial{w}}{\partial{z}}=0 \quad\Rightarrow\quad w={-}z\delta, \end{equation}

where $\delta \equiv \boldsymbol {\nabla }\boldsymbol {\cdot }{\boldsymbol {u}}$ is the horizontal divergence (and is independent of $z$). This already satisfies $w=0$ on $z=0$. On $z=h$, this gives the usual mass continuity equation for $h$:

(2.11)\begin{equation} Dh={-}h\delta \quad\Rightarrow\quad \frac{\partial{h}}{\partial{t}}+\boldsymbol{\nabla}\boldsymbol{\cdot}(h{\boldsymbol{u}})=0\,. \end{equation}

Next, the divergence-free condition on ${\boldsymbol {B}}_3$ implies, with $\tau =\boldsymbol {\nabla }\boldsymbol {\cdot }{\boldsymbol {B}}$,

(2.12)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{B}}+\frac{\partial{B_z}}{\partial{z}}=0 \quad\Rightarrow\quad B_z={-}z\tau, \end{equation}

using the boundary condition $B_z=0$ on $z=0$. This boundary condition arises as the lower boundary at $z=0$ is assumed to be perfectly conducting (see Appendix B). It is possible to relax this assumption and allow the field to penetrate the lower boundary, but this would entail coupling the magneto-hydrodynamics of the fluid layer with a diffusing three-dimensional magnetic field in the lowest layer, making the model significantly more complicated.

Turning next to the induction equation (2.8), the horizontal part is

(2.13)\begin{equation} D{\boldsymbol{B}}-{\boldsymbol{B}}\boldsymbol{\cdot}\boldsymbol{\nabla}{\boldsymbol{u}} = \eta\,\nabla^2{\boldsymbol{B}}. \end{equation}

This is already independent of $z$, so it is already vertically averaged. The vertical part of (2.8) is

(2.14)\begin{equation} \left(D-z\delta\,\frac{\partial{}}{\partial{z}}\right)({-}z\tau) -\left({\boldsymbol{B}}\boldsymbol{\cdot}\boldsymbol{\nabla}+({-}z\tau)\,\frac{\partial{}}{\partial{z}}\right)({-}z\delta)= \eta\,\nabla^2({-}z\tau). \end{equation}

All terms are proportional to $z$, and cancelling the common factor $-z$ gives

(2.15)\begin{equation} D\tau-{\boldsymbol{B}}\boldsymbol{\cdot}\boldsymbol{\nabla}\delta=\eta\,\nabla^2\tau. \end{equation}

This equation is, however, implied by the horizontal divergence of (2.13), so it is not new.

Now consider the momentum equation (2.7). Following Dritschel & Jalali (2020), we divide the pressure $p$ into a hydrostatic part $p_h=g(h-z)$ and a remaining non-hydrostatic part $p_n=p-p_h$. Then vertically averaging the horizontal part of (2.7), we obtain

(2.16)\begin{equation} D{\boldsymbol{u}}+f{\boldsymbol{k}}\times{\boldsymbol{u}} ={-}g\,\boldsymbol{\nabla}{h}-\frac{\boldsymbol{\nabla}\bar{p}_n}{h} + \bar{\boldsymbol{F}}_b, \end{equation}

(see § 3.1 in Dritschel & Jalali 2020), where

(2.17)\begin{equation} \bar{p}_n\equiv\int_0^h p_n\,{\rm d}{z} \end{equation}

is the vertically integrated non-hydrostatic pressure, and

(2.18)\begin{equation} \bar{\boldsymbol{F}}_b\equiv\frac{1}{h}\int_0^h \left(\,j_z{\boldsymbol{k}}\times{\boldsymbol{B}}-B_z\,\boldsymbol{\nabla}{B_z}\right){\rm d}{z} \end{equation}

is the vertically averaged horizontal part of the Lorentz force ( ${\boldsymbol {j}}_3\times {\boldsymbol {B}}_3$), while ${j_z=\partial {B_y}/}\partial {x}-\partial {B_x}/\partial {y}$ is the $z$ component of the current density and ${\boldsymbol {k}}\times {\boldsymbol {B}}=(-B_y,B_x)$. In (2.18), the first term in the integrand is independent of $z$, but $B_z=-z\tau$ is proportional to $z$. Carrying out the integration over $z$, we obtain

(2.19)\begin{equation} \bar{\boldsymbol{F}}_b=j_z{\boldsymbol{k}}\times{\boldsymbol{B}}-\frac{1}{6} h^2\,\boldsymbol{\nabla}{\tau^2}. \end{equation}

The vertically integrated non-hydrostatic pressure $\bar {p}_n$ in (2.17) is determined from the vertical component of (2.7) as follows. Since ${\boldsymbol {k}}\boldsymbol {\cdot }(\,{\boldsymbol {j}}_3\times {\boldsymbol {B}}_3)={\boldsymbol {B}}\boldsymbol {\cdot }\boldsymbol {\nabla }{B_z}$ owing to the fact that ${\boldsymbol {B}}$ is independent of $z$, the vertical component of (2.7) is

(2.20)\begin{equation} \left(D-z\delta\,\frac{\partial{}}{\partial{z}}\right)({-}z\delta) ={-}\frac{\partial{p_n}}{\partial{z}} + ({\boldsymbol{B}}\boldsymbol{\cdot}\boldsymbol{\nabla})({-}z\tau). \end{equation}

Expanding and simplifying, we find

(2.21)\begin{equation} \frac{\partial{p_n}}{\partial{z}}=z(D\delta-\delta^2-{\boldsymbol{B}}\boldsymbol{\cdot}\boldsymbol{\nabla}\tau). \end{equation}

As this is a first-order equation in $z$ for $p_n$, we need only one boundary condition, namely $p_n=0$ at $z=h$, to determine $p_n$ everywhere. Note that this boundary condition can be achieved in two ways. If the field is confined to the fluid layer (via an appropriate choice of diffusion), then this is achieved through the usual hydrodynamic approximations, i.e. that the fluid in the layer above has such a small density that it exerts a small constant pressure. If the field does diffuse into the layer above, then the gas pressure can still be small compared with the magnetic tension and pressure. For a low-$\beta$ (compressible) plasma above, the magnetic pressure dominates and the magnetic field is believed to adjust rapidly through a sequence of nearly force-free equilibria (see e.g. Wiegelmann & Sakurai 2021; Guo et al. 2016; and Appendices A and B).

Integrating (2.21) in $z$, we thus find

(2.22)\begin{equation} p_n=\tfrac{1}{2}(z^2-h^2)(D\delta-\delta^2-{\boldsymbol{B}}\boldsymbol{\cdot}\boldsymbol{\nabla}\tau). \end{equation}

A further integration, now from $z=0$ to $h$, yields the vertically integrated non-hydrostatic pressure $\bar {p}_n$ (see (2.17)) appearing in the horizontal momentum equation (2.16):

(2.23)\begin{equation} \bar{p}_n={-}\tfrac{1}{3} h^3(D\delta-\delta^2-{\boldsymbol{B}}\boldsymbol{\cdot}\boldsymbol{\nabla}\tau). \end{equation}

This completes the model equations, which consist of the evolution equations (2.11) for $h$, (2.13) for ${\boldsymbol {B}}$ and (2.16) for ${\boldsymbol {u}}$, alongside the relations (2.19) and (2.23). Except for the diffusive term in (2.13) and the inclusion of rotation, the equations are equivalent to those first derived by Dellar (2003). Note, however, using (2.23) in (2.16) leads to an implicit equation for ${\boldsymbol {u}}$, since a time derivative of $\delta =\boldsymbol {\nabla }\boldsymbol {\cdot }{\boldsymbol {u}}$ occurs in the $D\delta$ term in $\bar {p}_n$. Nevertheless, this is the standard form of the so-called Green–Naghdi equations (when there is no magnetic field).

Pearce & Esler (2010) deal with this problem (when ${\boldsymbol {B}}={\boldsymbol {0}}$) by forming an evolution equation for $\delta$ and combining part of the implicit term in $\bar {p}_n$ with $\partial {\delta }/\partial {t}$, so that the left-hand side of the evolution equation has the form $(1-\frac 13 H^2\,\nabla ^2)\,\partial {\delta }/\partial {t}$, where $H$ is the mean fluid depth. Yet, part of the implicit term associated with depth variations remains on the right-hand side, so the equation is still implicit and one cannot fully benefit from the smoothing effect arising from the inversion of the $(1-\frac 13 H^2\,\nabla ^2)$ operator. Alternatively, Holm (1988) propose using a ‘momentum variable’ ${\boldsymbol {m}}=h{\boldsymbol {u}}-\frac 13\boldsymbol {\nabla }(h^3\boldsymbol {\nabla }\boldsymbol {\cdot }{\boldsymbol {u}})$ to remove the implicitness, then deduce ${\boldsymbol {u}}$ by inverting the definition of ${\boldsymbol {m}}$ (see also Dellar 2003).

Another approach, arguably simpler, is to remove this implicitness by forming an elliptic equation directly for $\bar {p}_n$ that does not involve any time derivatives. For this, we follow Dritschel & Jalali (2020). Taking the divergence of the horizontal momentum equation (2.16), a little rearrangement yields

(2.24)\begin{equation} D\delta-\delta^2=\varXi-\boldsymbol{\nabla}\boldsymbol{\cdot}(h^{{-}1}\,\boldsymbol{\nabla}\bar{p}_n), \end{equation}

where

(2.25)\begin{equation} \varXi \equiv f\zeta-g\,\nabla^2 h + 2(J(u,v)-\delta^2)+\boldsymbol{\nabla}\boldsymbol{\cdot}\bar{\boldsymbol{F}}_b, \end{equation}

in which $\zeta ={\boldsymbol {k}}\boldsymbol {\cdot }\boldsymbol {\omega }_3=\partial {v}/\partial {x}-\partial {u}/\partial {y}$ is the vertical vorticity and $J({\cdot },{\cdot })$ is the Jacobian operator in Cartesian geometry (see § 3.1 in Dritschel & Jalali 2020). The important point is that $\varXi$ depends only on $h$, ${\boldsymbol {u}}$ and ${\boldsymbol {B}}$. However notice that $D\delta -\delta ^2$ also appears in the derived form of $\bar {p}_n$ in (2.23). Hence, replacing this term by (2.24), we arrive at an explicit, linear, elliptic equation for $\bar {p}_n$:

(2.26)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}(h^{{-}1}\,\boldsymbol{\nabla}\bar{p}_n)-3h^{{-}3}\bar{p}_n=\varXi-{\boldsymbol{B}}\boldsymbol{\cdot}\boldsymbol{\nabla}{\tau}=\tilde{\gamma}_u+\tilde{\gamma}_b, \end{equation}

where

(2.27)\begin{equation} \tilde{\gamma}_u \equiv f\zeta-g\,\nabla^2 h + 2(J(u,v)-\delta^2) \end{equation}

consists of purely hydrodynamic terms (as in Dritschel & Jalali 2020), while

(2.28)\begin{equation} \tilde{\gamma}_b \equiv \boldsymbol{\nabla}\boldsymbol{\cdot}\bar{\boldsymbol{F}}_b-{\boldsymbol{B}}\boldsymbol{\cdot}\boldsymbol{\nabla}{\tau} \end{equation}

consists of purely magnetic terms (apart from $h$).

In summary, given $h$, ${\boldsymbol {u}}$ and ${\boldsymbol {B}}$, we can compute $\bar {p}_n$ from (2.26), then use this in (2.16) to evolve ${\boldsymbol {u}}$. The complete set of equations therefore takes the form

(2.29)$$\begin{gather} D{\boldsymbol{B}}-{\boldsymbol{B}}\boldsymbol{\cdot}\boldsymbol{\nabla}{\boldsymbol{u}} = \eta\,\nabla^2{\boldsymbol{B}}, \end{gather}$$

(2.30)$$\begin{gather}\frac{\partial{h}}{\partial{t}}+\boldsymbol{\nabla}\boldsymbol{\cdot}(h{\boldsymbol{u}})=0, \end{gather}$$

(2.31)$$\begin{gather}D{\boldsymbol{u}}+f{\boldsymbol{k}}\times{\boldsymbol{u}} ={-}g\,\boldsymbol{\nabla}{h}-\frac{\boldsymbol{\nabla}\bar{p}_n}{h} + \bar{\boldsymbol{F}}_b, \end{gather}$$

(2.32)$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}(h^{{-}1}\,\boldsymbol{\nabla}\bar{p}_n)-3h^{{-}3}\bar{p}_n=\tilde{\gamma}_u+\tilde{\gamma}_b, \end{gather}$$

with $\bar {\boldsymbol {F}}_b$ defined in (2.19), $\tilde {\gamma }_u$ defined in (2.27), and $\tilde {\gamma }_b$ defined in (2.28). This is a closed, explicit system of equations governing a thin layer of a conducting magnetised fluid.

2.2. Energy

The total energy $E$ consists of kinetic, magnetic and potential parts. In the parent three-dimensional system,

(2.33)\begin{equation} E=\iint_{\mathcal{D}}\int_0^h\left(\frac{1}{2}(|{\boldsymbol{u}}_3|^2+|{\boldsymbol{B}}_3|^2)+gz\right) {\rm d}\kern 0.06em {x}\,{\rm d}{y}\,{\rm d}{z} \end{equation}

(see (2) in Dellar 2003). Here, ${\mathcal {D}}$ is the horizontal domain of integration. For the vertically averaged (Green–Naghdi) model, the horizontal parts of ${\boldsymbol {u}}_3$ and ${\boldsymbol {B}}_3$ are independent of $z$, while the vertical parts $w$ and $B_z$ are proportional to $z$. Performing the $z$ integration, we obtain

(2.34) \begin{equation} E=\frac{1}{2}\iint_{\mathcal{D}} h\left(|{\boldsymbol{u}}|^2+\frac{1}{3} h^2\delta^2 +|{\boldsymbol{B}}|^2+\frac{1}{3} h^2\tau^2+gh\right){\rm d}\kern 0.06em {x}\,{\rm d}{y}. \end{equation}

This is conserved only when there is no magnetic diffusivity $\eta$, and no normal component of ${\boldsymbol {B}}_3$ on $z=h$. In general, there is a magnetic energy flux through the free surface, physically owing to Alfvén waves propagating along field lines. In principle, the total energy in the extended domain $z\ge 0$ is conserved for the ideal case $\eta =0$. The system (2.29)–(2.32) always conserves total mass, proportional to the domain mean depth $H$, due to the flux form of (2.30). As only the parameter combination $gH\equiv c^2$ (a characteristic squared gravity wave speed) appears in the hydrostatic limit (see below) when using the dimensionless height anomaly

(2.35)\begin{equation} \tilde{h}=\frac{h-H}{H} \end{equation}

instead of $h$, it proves convenient to scale energy $E$ by $H$, and to distinguish the kinetic, magnetic and potential components as

(2.36) \begin{align} E/H &={\mathcal{E}}_u+{\mathcal{E}}_b+{\mathcal{E}}_h, \nonumber\\ {\rm where}\nonumber\\ {\mathcal{E}}_u &=\frac{1}{2}\iint_{\mathcal{D}} (1+\tilde h) \left(|{\boldsymbol{u}}|^2+\frac{1}{3} h^2\delta^2\right) {\rm d}\kern 0.06em {x}\,{\rm d}{y}, \nonumber\\ {\mathcal{E}}_b &=\frac{1}{2}\iint_{\mathcal{D}} (1+\tilde h) \left(|{\boldsymbol{B}}|^2+\frac{1}{3} h^2\tau^2\right) {\rm d}\kern 0.06em {x}\,{\rm d}{y}, \nonumber\\ {\mathcal{E}}_h &=\frac{1}{2} c^2\iint_{\mathcal{D}} \tilde h^2 \,{\rm d}\kern 0.06em {x}\,{\rm d}{y}. \end{align}

Note that the constant background potential energy has been removed to define ${\mathcal {E}}_h$.

2.3. Alternative variables

In geophysical fluid dynamics, it has proven advantageous to derive equations for the potential vorticity (PV) $q$ and for variables representing the leading-order departure from hydrostatic and geostrophic balance (Mohebalhojeh & Dritschel 2000, 2001; Smith & Dritschel 2006). The most convenient variables representing this departure are the velocity divergence $\delta =\boldsymbol {\nabla }\boldsymbol {\cdot }{\boldsymbol {u}}$ and the linearised acceleration divergence

(2.37)\begin{equation} \gamma_l=f\zeta-g\,\nabla^2 h=f\zeta-c^2\,\nabla^2\tilde{h} \end{equation}

(sometimes called the ageostrophic vorticity), due to their linear dependence on $\tilde {h}$ and ${\boldsymbol {u}}$ (recall $\zeta =\partial {v}/\partial {x}-\partial {u}/\partial {y}$). As for the PV, since it has a more complicated form in the vertically averaged (Green–Naghdi) equations, Dritschel & Jalali (2020) advocate using instead the linearised PV, $q_l=\zeta -f\tilde {h}$. The linear dependence of $q_l$, $\delta$ and $\gamma _l$ on $\tilde h$, $u$ and $v$ ensures that the latter variables may be recovered from the former by a simple inversion (for details, see § 3.2 in Dritschel & Jalali 2020).

The evolution equations for $q_l$, $\delta$ and $\gamma _l$ are derived readily from (2.30) and (2.31), and are given by (3.15)–(3.17) in Dritschel & Jalali (2020), apart from new terms arising from the magnetic field. These new terms come from the Lorentz force $\bar {\boldsymbol {F}}_b$ in (2.31). Both $\partial {q_l}/\partial {t}$ and $\partial {\gamma _l}/\partial {t}$ involve the vorticity tendency $\partial {\zeta }/\partial {t}$, which acquires the additional term ${\boldsymbol {k}}\boldsymbol {\cdot }(\boldsymbol {\nabla }\times \bar {\boldsymbol {F}}_b)$. By straightforward manipulation, one can show that

(2.38)\begin{equation} {\boldsymbol{k}}\boldsymbol{\cdot}(\boldsymbol{\nabla}\times\bar{\boldsymbol{F}}_b) = \boldsymbol{\nabla}\boldsymbol{\cdot}(\,j_z{\boldsymbol{B}}) + \tfrac{1}{6}\,J(\tau^2,h^2). \end{equation}

This term is added to (3.15) in Dritschel & Jalali (2020), whereas $f{\boldsymbol {k}}\boldsymbol {\cdot }(\boldsymbol {\nabla }\times \bar {\boldsymbol {F}}_b)$ is added to (3.17) in that paper. The $\delta$ equation was derived above in (2.24), but an equivalent and simpler equation can be found by exploiting (2.26). Then the evolution equations for $q_l$, $\delta$ and $\gamma _l$ take the form

(2.39)$$\begin{gather} \frac{\partial{q_l}}{\partial{t}} = \boldsymbol{\nabla}\boldsymbol{\cdot}(\,j_z{\boldsymbol{B}}-q_l{\boldsymbol{u}}) + J(\bar{p}_n,h^{{-}1}) + \frac{1}{6}\,J(\tau^2,h^2) \equiv S_q, \end{gather}$$

(2.40)$$\begin{gather}\frac{\partial{\delta}}{\partial{t}} = {\boldsymbol{B}}\boldsymbol{\cdot}\boldsymbol{\nabla}{\tau}-3h^{{-}3}\bar{p}_n+2\delta^2-\boldsymbol{\nabla}\boldsymbol{\cdot}(\delta{\boldsymbol{u}}), \end{gather}$$

(2.41)$$\begin{gather}\frac{\partial{\gamma_l}}{\partial{t}} = \mathbb{G}[\delta+\boldsymbol{\nabla}\boldsymbol{\cdot}(\tilde{h}{\boldsymbol{u}})]+fS_q, \end{gather}$$

where

(2.42)\begin{equation} \mathbb{G}\equiv c^2\,\nabla^2-f^2 \end{equation}

is the hydrodynamic ‘gravity-wave operator’. These are supplemented by the induction equation (2.29), as well as by the elliptic equation (2.32) for $\bar {p}_n$.

2.4. The magneto-hydrostatic limit

We now consider the limit where the mean depth satisfies $H\to 0$ while $g\to \infty$, yet keeping the product $gH=c^2$ finite. Using the dimensionless depth anomaly $\tilde h$ (see (2.35)), the reduced, magneto-hydrostatic shallow-water equations follow simply by dropping all terms proportional to positive powers of $H$, or more accurately $H/L$, where $L$ is a characteristic horizontal length. In particular, the non-hydrostatic pressure term in (2.31) is ${O}((H/L)^2)$ smaller than $g\,\boldsymbol {\nabla }{h}$ and is therefore negligible, so the horizontal momentum equation simplifies to

(2.43)\begin{equation} D{\boldsymbol{u}}+f{\boldsymbol{k}}\times{\boldsymbol{u}} ={-}c^2\,\boldsymbol{\nabla}{\tilde h}+\bar{\boldsymbol{F}}_b. \end{equation}

Moreover, the horizontal Lorentz force $\bar {\boldsymbol {F}}_b$ in (2.19) simplifies to

(2.44)\begin{equation} \bar{\boldsymbol{F}}_b=j_z{\boldsymbol{k}}\times{\boldsymbol{B}}, \end{equation}

after dropping a term ${O}((H/L)^2)$ smaller. The mass continuity equation (2.30) rewritten in terms of $\tilde h$ becomes

(2.45)\begin{equation} \frac{\partial{\tilde h}}{\partial{t}}+\boldsymbol{\nabla}\boldsymbol{\cdot}((1+\tilde h){\boldsymbol{u}})=0, \end{equation}

while the equation (2.29) for ${\boldsymbol {B}}$ remains unchanged. Taken together, (2.43), (2.45) and (2.29) constitute a closed set of equations for ${\boldsymbol {u}}$, $\tilde h$ and ${\boldsymbol {B}}$. These are the magneto-hydrostatic shallow-water equations, valid for a three-dimensional magnetic field. Note that the kinetic ${\mathcal {E}}_u$ and magnetic ${\mathcal {E}}_b$ energy components in (2.36) also simplify by dropping the terms proportional to $h^2$ in the integrals.

In alternative variables $q_l$, $\delta$ and $\gamma _l$, the Jacobian term in (2.39) is absent, while (2.41) is formally unchanged. The shallow-water equation for $\delta$ follows from (2.24) by dropping the non-hydrostatic pressure gradient term, which is again ${O}((H/L)^2)$ smaller. Upon rearranging (2.24), the final set of shallow-water equations becomes

(2.46)$$\begin{gather} \frac{\partial{q_l}}{\partial{t}} = \boldsymbol{\nabla}\boldsymbol{\cdot}(\,j_z{\boldsymbol{B}}-q_l{\boldsymbol{u}}) \equiv S_q, \end{gather}$$

(2.47)$$\begin{gather}\frac{\partial{\delta}}{\partial{t}} = \gamma_l+2\,J(u,v)+\boldsymbol{\nabla}\boldsymbol{\cdot}(\bar{\boldsymbol{F}}_b-\delta{\boldsymbol{u}}), \end{gather}$$

(2.48)$$\begin{gather}\frac{\partial{\gamma_l}}{\partial{t}} = \mathbb{G}[\delta+\boldsymbol{\nabla}\boldsymbol{\cdot}(\tilde{h}{\boldsymbol{u}})]+fS_q, \end{gather}$$

where $\mathbb {G}$ is the gravity-wave operator defined in (2.42). These are supplemented by the induction equation (2.29).

2.5. Model with alternative magnetic diffusion

The derivation above is for a model with regular magnetic diffusion controlled by ohmic dissipation. As noted in Appendix C, this leads to the magnetic field ${\boldsymbol {B}}_3$ not remaining tangential to the free surface of the fiducial layer, and a normal component developing. An alternative approach is to consider the evolution of a perfectly conducting fluid $\eta =0$, but then add in a dissipation operator to the two-dimensional equations that forces the free surface to remain a flux surface (Gilbert et al. 2022) – this requires $\boldsymbol {\nabla }\boldsymbol {\cdot }({h{\boldsymbol {B}}})=0$ at all times. The required dissipation operator is given by

(2.49)\begin{equation} \boldsymbol{d} ={-} h^{{-}1}\,\boldsymbol{\nabla} \times\left( \eta h\,\boldsymbol{\nabla} \times {\boldsymbol{B}} \right). \end{equation}

This term replaces the diffusive term in (2.29). However, it is redundant to evolve both components of ${\boldsymbol {B}}$. A more efficient, accurate approach is to instead evolve a scalar potential $A$, defined through the relation

(2.50)\begin{equation} {\boldsymbol{B}}=\frac{1}{1+\tilde{h}}\left(\frac{\partial{A}}{\partial{y}},-\frac{\partial{A}}{\partial{x}}\right), \end{equation}

which automatically satisfies $\boldsymbol {\nabla }\boldsymbol {\cdot }({h{\boldsymbol {B}}})=0$. Then the scalar potential $A$ satisfies the evolution equation

(2.51)\begin{equation} \frac{\partial{A}}{\partial{t}}+({\boldsymbol{u}}+{\boldsymbol{u}}_d)\boldsymbol{\cdot}\boldsymbol{\nabla}{A} = \eta\,\nabla^2{A}, \end{equation}

where ${\boldsymbol {u}}_d\equiv \eta h^{-1}\,\boldsymbol {\nabla }{h}$ is a diffusion velocity (Gilbert et al. 2022). This equation replaces (2.29); all remaining equations are unchanged in this alternative model.

3. Results

In this section, we present results first for the magneto-hydrostatic shallow-water equations derived in § 2.4, then for their non-hydrostatic generalisation, to reveal what new features arise. As a model problem, we discuss the impact of an initially unidirectional magnetic field on an isolated vortex. The set-up is identical to that discussed previously in Dritschel et al. (2018), where the two-dimensional limit ($c=\sqrt {gH}\to \infty$) was studied. In that limit, $\delta =\gamma _l=0$, the free-surface remains flat ($\tilde {h}=0$) and there are no inertia–gravity waves. Here, by contrast, we consider the impact of moderate free-surface variations $\tilde {h}={{O}}(1)$, a finite Rossby number ${Ro}=|\zeta |_{max}/f$, and a finite Rossby deformation length $L_D$ defined by

(3.1)\begin{equation} L_D = \frac{\sqrt{gH}}{f} = \frac{c}{f}, \end{equation}

which controls the deformation of the free surface in a rotating shallow-water system (see e.g. Vallis 2017). The two-dimensional model corresponds to the limit $L_D/L\to \infty$, where $L$ is a characteristic horizontal length scale of the flow (the ratio $(L_D/L)^2$ is known as the Burger number).

3.1. Flow initialisation

The initial flow consists of a Gaussian vortex having (vertical) vorticity

(3.2)\begin{equation} \zeta=\frac{\varepsilon f \exp({-(r/R)^2/2})}{2(1-{\rm e}^{{-}1/2})}-C, \end{equation}

where $r$ is the cylindrical radius (taking the domain centre at the origin), $\varepsilon$ is a nominal Rossby number, $f$ is the Coriolis frequency, and $C$ is a constant required to ensure that the domain average vorticity is zero (a consequence of Stokes’ theorem). We take $f=4{\rm \pi}$ without loss of generality so that a unit of time is one ‘day’. The vortex radius is $R=5{\rm \pi} /32$ as in Dritschel et al. (2018), and the simulation domain is the periodic box $[-{\rm \pi},{\rm \pi} )^2$. The slight discontinuity in the normal derivative of the initial $\zeta$ at the domain boundaries has negligible impact on the results. This choice for $\zeta$ corresponds to a maximum tangential velocity $U_0=\omega _0 R/2$ (where $\omega _0\equiv \varepsilon f$) in an infinite domain.

We initialise all flow fields from the relative vorticity $\zeta$ alone, by enforcing hydro-cyclo-geostrophic balance. This is tantamount to requiring $\bar {p}_n=\delta =\partial \delta /\partial {t}=0$ in the absence of a magnetic field (${\boldsymbol {B}}={\boldsymbol {0}}$). From (2.47), this implies that the ageostrophic vorticity is $\gamma _l=-2J(u,v)$, while $u=-\partial \psi /\partial {y}$ and $v=\partial \psi /\partial {x}$ are found by inverting $\nabla ^2\psi =\zeta$ for the streamfunction $\psi$. Then the definition of $\gamma _l$ in (2.37) provides a Poisson equation for the dimensionless height anomaly, $\nabla ^2\tilde {h}=(f\zeta -\gamma _l)/c^2$, which is easily inverted for $\tilde {h}$. Finally, from $\zeta$ and $\tilde {h}$, the linearised PV $q_l$ is obtained simply from $q_l=\zeta -f\tilde {h}$.

The initial magnetic field has $B_x=B_0/(1+\tilde {h})$ and $B_y=0$. In this way, the initial field ${\boldsymbol {B}}_3$ is tangent to the free surface at $z=H(1+\tilde {h})$ (note that $\tilde {h}\to 0$ in the limit $L_D\to \infty$). The value of $B_0$ is chosen together with the magnetic diffusivity $\eta$ so that the maximum magnetic field ${\boldsymbol {B}}$ becomes comparable in magnitude with $U_0$ when ${\boldsymbol {B}}$ is fully intensified by the vortex (see § 2.2 in Dritschel et al. 2018). In short, we take $\eta =\omega _0({\rm \Delta} x)^2$, where ${\rm \Delta} x=2{\rm \pi} /n_g$, and $n_g$ is the grid resolution in both $x$ and $y$ (most simulations below use $n_g=512$). We then specify the ‘gain’

(3.3)\begin{equation} \mathfrak{g}\equiv\frac{B_0/\mathfrak{d}}{U_0}, \end{equation}

where $\mathfrak {d}={\rm \Delta} x/R$ is the dimensionless diffusion length (note $\sqrt {\eta /\omega _0}={\rm \Delta} x$). The gain is the ratio of the maximum expected magnetic field strength $B_0/\mathfrak {d}$ to the flow speed $U_0$, at least in the two-dimensional limit $L_D\to \infty$. Notably, the magnetic Reynolds number is

(3.4)\begin{equation} {Rm} = \frac{U_0 R}{\eta} = \frac{1}{2\mathfrak{d}^2}. \end{equation}

This is equal to $800$ at the default grid resolution $n_g=512$.

The numerical codes have been adapted directly from those used in Dritschel et al. (2018) for the two-dimensional limit $L_D\to \infty$, and from the hydrodynamic shallow-water and vertically averaged codes used in Dritschel & Jalali (2020) and Jalali & Dritschel (2021). Without a magnetic field, the codes reproduce the results in Jalali & Dritschel (2021) with very minor differences associated with using the linearised PV anomaly $q_l$ in place of full PV, a variable time step, and variable hyperviscosity in the present codes. (The numerical damping rate on the highest resolved wavenumber is $10(f+\zeta _{rms})$ rather than $10f$; see Appendix C in Dritschel & Jalali 2020.) Also, for large $L_D$, the shallow-water code reproduces closely the two-dimensional magnetic results of Dritschel et al. (2018). This is a difficult limit to simulate with a shallow-water code since short-scale gravity waves are fast ($c=fL_D\gg U_0$), requiring a small time step for accuracy. At the same time, the Rossby number $\varepsilon$ must be small to avoid ageostrophic effects. Nonetheless, with $L_D=4$ and ${\varepsilon =0.25}$, the evolution of the vertical current density $j_z$ shown in figure 1 is strikingly similar to that shown in figure 2 (three middle panels, right column) in Dritschel et al. (2018) – note that the field values are a factor $\varepsilon$ smaller here. Only by the latest time are differences apparent, with the vortex rotated slightly less here, an effect due mainly to the finite value of $L_D$. The agreement is remarkable considering that $L_D$ is not particularly large, nor is $\varepsilon$ particularly small. Moreover, the vorticity-based Rossby number ${Ro}=|\zeta |_{max}/f$ rises from approximately $0.305$ to $0.884$ at late times. However, the free surface varies by no more than $0.614\,\%$, and this decays to $0.335\,\%$ by late times.

Figure 1. Current density field $j_z({\boldsymbol {x}},t)$ at three scaled times $\tilde {t}=\varepsilon t$ for an MHD shallow-water simulation with $\mathfrak {g}=2$, ${\varepsilon =0.25}$ and $L_D=4$.

Figure 2. Current density field $j_z({\boldsymbol {x}},t)$ at the final time $\tilde {t}=\varepsilon t=25$ for $\varepsilon =0.1$, ${L_D=0.25}$ and $\mathfrak {g}=2$; (a) regular magnetic diffusion; (b) alternative magnetic diffusion (which keeps $\boldsymbol {\nabla }\boldsymbol {\cdot }(h{\boldsymbol {B}})=0$); and (c) the difference ((b) minus (a)).

3.2. Dependence on the form of the diffusion

We start by comparing the evolution of the models with different forms of magnetic diffusion, holding all other parameters fixed. (Note that in the numerical code, the scalar potential $A$ in (2.51) is split into a mean part $B_0 y$ and a residual $\tilde {A}$ so that the latter can be represented as a periodic field.) To produce a strong interplay between the magnetic field and the fluid motion, we choose the Rossby number $\varepsilon =0.1$, the radius of deformation ${L_D=0.25}$, and the gain $\mathfrak {g}=2$. The current density field $j_z$ at the final time ($t=250$) is shown in figure 2, where the simulation using regular magnetic diffusion in figure 2(a) is compared to that using the alternative model of Gilbert et al. (2022) in figure 2(b); the difference field is shown in figure 2(c). Qualitatively, the alternative diffusion that keeps $\boldsymbol {\nabla }\boldsymbol {\cdot }(h{\boldsymbol {B}})=0$ has little impact, even by this late time in the evolution. The overall structure and amplitude of the field are closely comparable to those found in the regular diffusion case. There are quantitative differences (see figure 2c), primarily at small scales; however, in such a complex flow, the small scales are the least predictable, so this is to be expected. Comparing the energy components in figure 3, there are hardly any visible differences between the two simulations. The kinetic energy differs the most, but still only slightly. The detailed variations of the energy components compare closely despite using different magnetic diffusion mechanisms.

Figure 3. (a–d) Time evolutions of the scaled total, potential, kinetic and magnetic energy, respectively, for a simulation with regular magnetic diffusion (solid curves) and a simulation with alternative magnetic diffusion (dashed curves). The flow parameters are the same as in figure 2. Note that $\mathcal {E}_h$ is defined relative to its initial value to facilitate comparison after scaling by $\varepsilon L_D^2$. This is why this component and the total energy have negative values.

Hence we believe that the precise form of the diffusive operator has little effect on the form of the solution for these problems where the magnetic field is imposed a priori (rather than generated self-consistently via a dynamo process, where the diffusive process is key). For this reason, henceforth in this paper, we will consider the evolution of the model using regular magnetic diffusion, where the magnetic field is allowed to leave (and enter) the upper free surface.

3.3. Dependence on the Rossby deformation length

We next investigate flows for which the Rossby deformation length $L_D$ is comparable with or smaller than the vortex radius $R$ ($\approx 0.5$), and for which the free surface exhibits strong variations. Here, we consider flows in magneto-hydrostatic balance, and defer discussion of non-hydrostatic effects to § 3.5. We vary $L_D$ and $\varepsilon$ together so that the initial height anomaly field $\tilde h$ is similar in each simulation. This is done by choosing $\varepsilon =1.6L_D^2$ for three values of $L_D$: $0.5$, $0.25$ and $0.125$. This gives $\varepsilon =0.4$, $0.1$ and $0.025$, respectively. The corresponding initial minimum and maximum values of $\tilde h$ are $(-0.652,0.164)$, $(-0.604,0.161)$ and $(-0.592,0.161)$, and these values change little as the flow evolves. The initial horizontal magnetic field is specified through the gain parameter $\mathfrak {g}$, and as in figure 1, we take $\mathfrak {g}=2$.

The three simulations were carried out to the same final scaled time $\tilde {t}=\varepsilon t=25$. The vorticity evolution is compared in figure 4 at a few characteristic times. At early times, the evolution is very similar, with a wind-up of alternating positive and negative bands of vorticity that strengthen up to approximately $t=14$ (the Rossby number ${Ro}$ climbs from just under $0.5$ to over $2.1$ when $L_D=1/2$). These bands are created largely by the effect of the Lorentz force; the initial vorticity is only slightly negative outside the vortex. The tightly wound bands then undergo instability and mix, leaving a region of near zero vorticity surrounding a central core (see also Gilbert, Mason & Tobias 2016). In the periphery, strong bands of vorticity stretch in the periodic flow field and decay gradually at late times.

Figure 4. Vertical vorticity field $\zeta ({\boldsymbol {x}},t)$ at three scaled times $\tilde {t}=\varepsilon t$ for three simulations having different $L_D$ as labelled, and for $\mathfrak {g}=2$. Note that $\varepsilon =1.6L_D^2$.

As $L_D$ decreases, the instability is less vigorous and smaller scale, yet leaves a smaller vortex core at late times. The flow is also more axisymmetric. This behaviour is consistent with the fact that flow interactions generally weaken as $L_D$ decreases. This is quantified in figure 5, which shows the Rossby number $|\zeta |_{max}/f$, and the scaled maximum horizontal magnetic field (or actual gain) $\|{\boldsymbol {B}}\|_{max}/U_0$, as functions of scaled time $\tilde {t}=\varepsilon t$ for the three simulations. The Rossby number evolves similarly in all cases, the main difference being the magnitude. There is always an initial drop as the tension in the twisting magnetic field slows the vortex rotation (elaborated in Dritschel et al. 2018), but then there is a significant increase as current sheets form in the periphery of the vortex core. At late times, these current sheets diffuse and the vorticity weakens slightly. Regarding $\|{\boldsymbol {B}}\|_{max}/U_0$, theoretically this should reach $\mathfrak {g}=2$ when the flow is mature, but the actual gain is seen to be $0.4\mathfrak {g}$–$0.5\mathfrak {g}$. In fact, the gain parameter $\mathfrak {g}$ is only a qualitative estimate of the amplification of the magnetic field. Perhaps surprisingly, the actual gain is larger when $L_D$ is smaller, despite the weaker flow interactions occurring in this case. This may explain the reduction in the size of the vortex core as $L_D$ decreases.

Figure 5. Time evolutions of (a) the maximum vorticity scaled by the Coriolis frequency (i.e. Rossby number), and (b) maximum horizontal magnetic field scaled by the characteristic flow speed $U_0$, for the three simulations illustrated in figure 4.

The current density evolution (not shown) is broadly similar to that exhibited by the vorticity in figure 4, with the exception that there is no core; the magnetic field is largely expelled from the centre, as seen previously in figure 1 and first theorised by Weiss (1966) and elucidated by Moffatt & Kamkar (1983). On the other hand, the horizontal divergence field, shown in figure 6, exhibits very different behaviour. While it also winds up at early times, it fragments into small-scale structures that subsequently grow in scale and weaken. These structures are likely to be imbalanced inertia–gravity waves, and an especially novel feature is that these waves appear to be confined to the region of relatively weak magnetic field, and so may be suppressed by the presence of a strong magnetic field. Notably, the amplitude of $\delta$ is much smaller than that of $\zeta$, implying that the associated divergent flow contributes relatively little to the velocity field. Moreover, it decreases in proportion to $L_D^2$ or faster at intermediate and late times.

Figure 6. As in figure 4 but for the horizontal divergence field $\delta ({\boldsymbol {x}},t)$.

The energy evolution in these three flows is shown in figure 7 versus (rescaled) time $\tilde {t}$. In the figure, the energy components are scaled by $\varepsilon L_D^2$ to enable comparison, and moreover, the initial value of the potential energy $\mathcal {E}_h(0)$ is subtracted from $\mathcal {E}_h(\tilde {t})$ (the reference value of potential energy is not important, but the construction here facilitates scaling). The total energy $\mathcal {E}$ in figure 7(a) is the sum of the (redefined) potential $\mathcal {E}_h$, kinetic $\mathcal {E}_u$ and magnetic $\mathcal {E}_b$ components; see (2.36). The initial strong growth in $\mathcal {E}_b$ is offset primarily by the decrease in $\mathcal {E}_h$, especially at small $L_D$. As $L_D$ increases, the kinetic energy $\mathcal {E}_u$ accounts for an increasing proportion of this initial energy change (notice the sharper decrease of the solid curve for $L_D=1/2$ in figure 7(c), together with the slower decrease in $\mathcal {E}_h$ in figure 7(b)). Interestingly, the greatest growth in magnetic energy $\mathcal {E}_b$ occurs for the smallest $L_D$; evidently, small $L_D$ favours a greater intensification of the magnetic field ${\boldsymbol {B}}$, which in turn leads to greater destruction of the vortex core, as seen in figure 4. At early times, the total energy $\mathcal {E}$ is conserved approximately; thereafter, $\mathcal {E}$ decays owing to both magnetic and numerical dissipation. The numerical dissipation is associated indirectly with current sheet formation, which tends to generate collocated intense vorticity sheets. These sheets subsequently stretch and thin, eventually dissipating when they reach the scale of ‘surgery’ in the numerical method (here $1/16$th of the grid scale). This dissipation is inevitable in any numerical method owing to the forward energy cascade (in spectral space) induced by the magnetic field, in particular the Lorentz force (2.44).

Figure 7. (a–d) Time evolutions of the scaled total, potential, kinetic and magnetic energy, respectively, for the three simulations illustrated in figures 4 and 6. Note that $\mathcal {E}_h$ is defined relative to its initial value to facilitate comparison after scaling by $\varepsilon L_D^2$.

3.4. Dependence on the magnetic diffusivity

We next examine the dependence of the flow evolution on the magnetic diffusivity $\eta$. We keep all other parameters the same as in the simulation featuring in figure 4(d–f) ($L_D=1/4$, $\varepsilon =0.1$, $\mathfrak {g}=2$), and vary $\eta$, equivalent to varying the grid resolution $n_g$ in our set-up. The vorticity field $\zeta$ is compared in figure 8. Note that the initial horizontal magnetic field amplitude is varied so that $B_0\propto {Rm}^{-1/2}$ to keep the gain fixed at $\mathfrak {g}=2$, on the understanding that when the flow reaches a mature stage, $\|{\boldsymbol {B}}\|\sim \mathfrak {g}U_0$. At early times, all simulations show the characteristic spiral windup of the vorticity field (echoed by $j_z$, $\delta$ and $\gamma$, not shown). The outer core then destabilises, leading to a mixed region of weak vorticity surrounding a compact inner core. At low ${Rm}$, there is very little vorticity growth before decay at late times (mainly outside the core region). With increasing ${Rm}$, the vorticity growth becomes pronounced, especially on the periphery of the outer core. Filaments thin and intensify owing to interactions with the magnetic field (the current density $j_z$ looks closely similar apart from the absence of an inner core, not shown).

Figure 8. Vertical vorticity field $\zeta ({\boldsymbol {x}},t)$ at three scaled times $\tilde {t}=\varepsilon t$ for three simulations differing in their magnetic Reynolds number ${Rm}$, as indicated. The corresponding grid resolutions are (a–c) $n_g=256$, (d–f) $n_g=512$, and (g–i) $n_g=1024$. Otherwise, the simulations are identical, with ${L_D=0.25}$, $\varepsilon =0.1$ and $\mathfrak {g}=2$.

The growth in both vorticity and maximum horizontal magnetic field strength is quantified in figure 9 for these three simulations. As expected, high ${Rm}$ favours a strong and persistent growth in vorticity, which is caused by the emergence of current sheets (not shown). The amplification of vorticity is approximately proportional to ${Rm}$, which is the expected amplification of $j_z$ based on simple scaling arguments (Dritschel et al. 2018). On the other hand, the actual gain in the horizontal magnetic field decreases with ${Rm}$, though it is sustained over a longer period of time before diffusive effects weaken the field. While three simulations cannot be conclusive, it appears that $\|{\boldsymbol {B}}\|_{max}/U_0$ remains finite as ${Rm}\to \infty$. However, recall that the initial magnetic field is $B_0\propto {Rm}^{-1/2}$; this goes to zero in this limit.

Figure 9. Time evolutions of (a) the maximum vorticity scaled by the Coriolis frequency (i.e. Rossby number, and (b) maximum horizontal magnetic field scaled by the initial maximum flow speed $U_0$, for the three simulations illustrated in figure 8.

The total energy $\mathcal {E}$, together with its potential $\mathcal {E}_h$, kinetic $\mathcal {E}_u$ and magnetic $\mathcal {E}_b$ components, is shown as a function of (rescaled) time $\tilde {t}$ for these three flows in figure 10. The total energy $\mathcal {E}$ is better conserved especially at early times for large ${Rm}$, as expected. The potential energy (or its difference from its initial value) exhibits a similar variation, with a slower decrease occurring for large ${Rm}$. Some of this decrease is dissipative, especially for low ${Rm}$, but some also compensates for the gain in magnetic energy $\mathcal {E}_b$. There is a similar decrease in kinetic energy $\mathcal {E}_u$, though this is smaller due to the small value of $L_D$ (here $1/4$ or approximately half of the vortex radius). The magnetic energy $\mathcal {E}_b$ initially grows from low levels as the horizontal field ${\boldsymbol {B}}$ amplifies, reaching a peak that is delayed and reduced with increasing ${Rm}$. The gain in $\mathcal {E}_b$ is largest for low ${Rm}$. Recall that the scaling of mean magnetic field applied is chosen so that the peak magnetic field in each case is comparable. However, at higher ${Rm}$, the magnetic field is confined to thinner regions (with length scale ${\propto }{Rm}^{1/2}$), so it is reasonable that the integrated energy reduces as ${Rm}$ increases.

Figure 10. (a–d) Time evolutions of the scaled total, potential, kinetic and magnetic energy, respectively, for the three simulations illustrated in figure 8. Note that $\mathcal {E}_h$ is defined relative to its initial value.

3.5. Dependence on finite depth: non-hydrostatic effects

Finally, we turn to the non-hydrostatic magnetic shallow-water or ‘vertically averaged’ model derived in § 2.1 (see also § 2.3). This model extends the traditional hydrostatic shallow-water model by relaxing the long-wave approximation (necessary for the hydrostatic approximation), thereby capturing explicitly horizontal scales comparable to the mean fluid depth $H$. We compare three simulations, varying a dimensionless frequency ratio $f/N$, where $f$ is the Coriolis frequency as before, while $N\equiv \sqrt {3g/H}$ can be thought of as a buoyancy frequency (Dritschel & Jalali 2020). The justification for this appellation is that the dispersion relation for linear waves in the hydrodynamical context then closely resembles that occurring in a three-dimensional linearly stratified fluid. Small values of $f/N$ correspond to strongly stratified flows, characteristic of large to intermediate scale atmospheric and oceanic flows, while moderate to large values of $f/N$ correspond to weakly stratified flows (for observational evidence, see Dritschel & Jalali 2020). When $f/N=1$ identically, linear inertia–gravity waves are entirely trapped; their group velocity vanishes for all wavelengths. In general, inertia–gravity wave frequencies $\omega$ lie between $f$ and $N$, with $\omega \approx f$ for the longest wavelengths, and $\omega \approx N$ for the shortest wavelengths.

Figure 11 compares, for the final time only, the vorticity $\zeta$, divergence $\delta$, and scaled non-hydrostatic pressure $P_n=\bar {p}_n/H$ for three values of $f/N$: $1/4$, $1$ and $4$. The structure of the vorticity field is weakly dependent on $f/N$, with some variation in shape of the central core, and with the lowest extreme field values occurring for $f/N=1$. The root-mean-square vorticity varies even less, taking the values $1.027$, $1.062$ and $1.059$, respectively, for $f/N=1/4$, $1$ and $4$. This weak variation with $f/N$ is consistent with previous findings in the hydrodynamical context (Dritschel & Jalali 2020). Vorticity is dominated by the underlying balance, controlled by the relatively slow advection of PV. By contrast, both divergence and non-hydrostatic pressure vary strongly with $f/N$. Both the amplitude and spatial pattern vary considerably, especially for $P_n$, which appears to scale like $(f/N)^2$. One may also observe that $\delta$ and $P_n$ are approximately $45^\circ$ out of phase but otherwise exhibit a similar spatial pattern. The impact of $f/N$ on the evolution of the magnetic field (e.g. $j_z$, not shown) is comparably weak to that seen above for the vorticity $\zeta$.

Figure 11. Relative vorticity $\zeta$, divergence $\delta$ and scaled vertically integrated non-hydrostatic pressure ${P_n=\bar {p}_n/H}$ at the final (scaled) time $\tilde {t}=\varepsilon t=25$ for three non-hydrostatic shallow-water simulations with varying $f/N$ as indicated. Here, $\mathfrak {g}=2$, ${\varepsilon =0.2}5$ and $L_D=1/2$.

Regarding energy, the impact of $f/N$ is barely visible in comparison to that found above when varying $L_D$ in figure 7 – the weakest energy component, the magnetic energy $\mathcal {E}_b$, differs by no more than 3.5 % across $f/N$ (not shown), while the other components differ by much less than 1 %. This reflects the fact that the balanced part of the flow dominates energetically. In more turbulent flows than those studied here, the impact of $f/N$ is expected to be larger, as found previously in purely hydrodynamical flows (Dritschel & Jalali 2020), owing to greater activity (i.e. more energy) at small scales.

4. Discussion

Following the approach of Green & Naghdi (1976), we have derived a non-hydrostatic shallow-water model for a layer of magnetised fluid that consistently includes finite conductivity (non-zero magnetic diffusivity). Provided that the gas pressure is negligible in the layer above, one may consider two separate models of diffusion. The first allows the magnetic field to exit the domain, and the field in the layer above transits rapidly through a series of force-free equilibria. The application of such a model may be more natural to the outer layers of astrophysical objects (see e.g. Hindle, Bushby & Rogers 2019) having a tenuous atmosphere above. The second model of diffusion confines the field to the primary fluid layer so that the upper free surface remains a magnetic flux surface, which may be a good modelling assumption for the solar tachocline (Tobias & Weiss 2007). In fact, the dynamical difference between these models is exceptionally small at high $Rm$. We note that such models (and even simpler models) have been used previously for studying the stably stratified interiors of stars.

The non-hydrostatic model discussed here gives a representation of small horizontal scales (comparable to the free-surface depth) that is more accurate than traditional shallow-water models. This may be particularly important for situations in which magnetic fields play a key role. It is well known that the presence of magnetic fields can lead to the promotion of a forward cascade, rather than the inverse cascade of quasi-two-dimensional hydrodynamics (see e.g. Seshasayanan, Benavides & Alexakis (2014) for a discussion of the transition between the two types of behaviour). Dellar (2003) has already demonstrated that the linearised form of the model presented here better replicates the dispersion properties of the full three-dimensional system for waves on a smaller scale, and that nonlinear cnoidal wave solutions exist.

The derivation of the non-hydrostatic model presented here follows Green & Naghdi (1976), who pointed out that for the hydrodynamic case, all that is required is a single assumption on the horizontal velocity field – namely, that it is independent of height $z$ in the layer; the equations then follow by vertical averaging (see also Miles & Salmon 1985; Jalali & Dritschel 2021). In the magnetised case, the analogous assumption is that both the velocity and the horizontal magnetic field are independent of height. The non-hydrostatic magnetised shallow-water equations, together with all of their associated conservation laws, then follow directly by vertically averaging the (inviscid) Euler and magnetic induction equations.

The total energy contained within the fluid layer is not conserved owing to magnetic diffusion. Notably, without diffusion, energy is conserved only for a magnetic field confined entirely to the fluid layer (tangent to its upper and lower surfaces). As diffusion causes the magnetic field to develop a non-zero normal component at the surfaces, inevitably energy will pass through the surfaces, as well as diffuse. This situation is common in models of the solar atmosphere (Priest 2014).

Numerical simulations were conducted for the magneto-hydrostatic version of the model (derived assuming finite gravity wave speeds $\sqrt {gh}$ and taking the limit $h/L\to 0$, where $L$ is a characteristic horizontal scale), as well as for the fully non-hydrostatic model. We revisited the problem of a Gaussian vortex initially within a (nearly) uniform horizontal magnetic field, studied previously in Dritschel et al. (2018) in the two-dimensional limit of a flat free surface.

In the hydrostatically balanced cases, finite free-surface variations in the shallow-water model and a finite Rossby number nonetheless produce results that are qualitatively similar to those in Dritschel et al. (2018). The most pronounced differences occur when the Rossby deformation length $L_D=\sqrt {gH}/f$ (where $H$ is the mean fluid depth, and $f$ is the Coriolis frequency) is small compared to the initial vortex radius $R$. Then the core region around the vortex is smaller and more circular, and contains a broader region of well-mixed nearly irrotational flow. The initial magnetic field is expelled to the periphery of this region, where it exhibits current sheets and collocated vorticity sheets. In the parameter regime investigated, the divergent component of the flow is weak, and weakens in proportion to $(L_D/R)^2$.

For the non-hydrostatic model, we investigated the effects of varying the Coriolis buoyancy frequency ratio $f/N$, where $N\equiv \sqrt {3g/H}$ depends on the mean depth $H$. Perhaps as expected, whilst the vorticity (dominated by the underlying balance) is relatively insensitive to the frequency ratio, the divergence and non-hydrostatic pressure are strongly affected by changes in this ratio, exhibiting energy on smaller scales as $f/N$ increases, which may be the signature of magneto-inertia–gravity waves. We note that a greater impact of $f/N$ is likely to occur in fully turbulent flows due to heightened small-scale activity; this is worth investigating further.

We believe that the model may be used to study general features of the solar tachocline, including free-surface and non-hydrostatic effects, exploring the properties and behaviour of linear and nonlinear waves (Rossby, gravity, Alfvén). Finally, an extension to spherical geometry will enable the investigation of global instabilities and MHD turbulence in thin spherical shells, which will provide better models of the solar tachocline and the stability of jets in exoplanetary atmospheres.