1. Introduction
The shallow water system has been a crucial workhorse model for our understanding of geophysical fluid dynamics (GFD) (e.g. Zeitlin Reference Zeitlin2018). It simplifies the more complex GFD equations (e.g. the Boussinesq system) by assuming that one can represent the vertical variation using only one (or a few) stacked stratified layers. This simplification allows for easier numerical and theoretical investigations of its properties, and in many cases shallow water captures the realistic properties of geophysical phenomena in the atmosphere and ocean of Earth and other planets (e.g. Dowling & Ingersoll Reference Dowling and Ingersoll1989; Polvani et al. Reference Polvani, McWilliams, Spall and Ford1994; Cho & Polvani Reference Cho and Polvani1996; Lambaerts, Lapeyre & Zeitlin Reference Lambaerts, Lapeyre and Zeitlin2012; Bembenek, Straub & Merlis Reference Bembenek, Straub and Merlis2020, for just a few).
A further simplification of the shallow water system is the shallow water quasi-geostrophic equation (SWQG) (Vallis Reference Vallis2017; Zeitlin Reference Zeitlin2018). It assumes rotation dominates and thus the geostrophic balance approximation is valid and the Rossby number is small. The shallow water quasi-geostrophic equation is then an asymptotic approximation of the shallow water model in the limit of vanishing Rossby number. It is a prototypical example of a balanced model, in that it can be written with only a single prognostic variable, potential vorticity (PV), from which all other physical variables can be diagnosed. This has the consequence that inertial-gravity waves are not permitted in its dynamics. The SWQG (and its special case, the two-dimensional Euler equations) can capture the dynamics of shallow water in a surprisingly wide parameter regime (Sugimoto, Ishioka & Yoden Reference Sugimoto, Ishioka and Yoden2007), and it has rich turbulence phenomena that closely resemble the balanced dynamics of rotationally dominated atmospheres and oceans (e.g. Gill, Green & Simmons Reference Gill, Green and Simmons1974; Salmon Reference Salmon1980; Maltrud & Vallis Reference Maltrud and Vallis1991; Smith et al. Reference Smith, Boccaletti, Henning, Marinov, Tam, Held and Vallis2002; Dritschel & McIntyre Reference Dritschel and McIntyre2008; Early, Samelson & Chelton Reference Early, Samelson and Chelton2011). The study of SWQG has led to progress in many practical geophysical applications, including the parametrisation of eddies (Held & Larichev Reference Held and Larichev1996; Thompson & Young Reference Thompson and Young2007; Srinivasan & Young Reference Srinivasan and Young2014; Gallet & Ferrari Reference Gallet and Ferrari2020, Reference Gallet and Ferrari2021), the interpolation and interpretation of satellite altimetry measurement of sea surface height (Stammer Reference Stammer1997; Chelton et al. Reference Chelton, Schlax, Samelson and de Szoeke2007; Ubelmann, Klein & Fu Reference Ubelmann, Klein and Fu2015), the analysis of ocean coherent vortices over topographic features (Bretherton & Haidvogel Reference Bretherton and Haidvogel1976; Carnevale & Frederiksen Reference Carnevale and Frederiksen1987; Siegelman & Young Reference Siegelman and Young2023; LaCasce, Palóczy & Trodahl Reference LaCasce, Palóczy and Trodahl2024), the statistics of atmospheric blocking (Marshall & Molteni Reference Marshall and Molteni1993; Lucarini & Gritsun Reference Lucarini and Gritsun2020) and the mechanics of jets and spots on gas giants (Turkington et al. Reference Turkington, Majda, Haven and DiBattista2001; Majda & Wang Reference Majda and Wang2006; Siegelman, Young & Ingersoll Reference Siegelman, Young and Ingersoll2022; Pizzo & Salmon Reference Pizzo and Salmon2024) to name but a few.
Despite the range of phenomenology SWQG captures, it still misses many features of the full shallow water model, even those that might be categorised as `balanced’. The SWQG misses the asymmetry between cyclonic and anticyclonic vortices. Polvani et al. (Reference Polvani, McWilliams, Spall and Ford1994) note that the vorticity skews negative in a freely decaying simulation of the shallow water system. Arai & Yamagata (Reference Arai and Yamagata1994) and Graves, McWilliams & Montgomery (Reference Graves, McWilliams and Montgomery2006) further analyse the detailed free evolution of vortices in the shallow water model and conclude that anticyclones are more robust to small-scale perturbation and to the influence of other vortices, providing a plausible mechanics of the negative skewness of vorticity. For an ocean gyre forced by wind, Thiry et al. (Reference Thiry, Li, Mémin and Roullet2024) found the shallow water model version exhibits asymmetric drift of the strong eastward jet, while the SWQG version lacks this asymmetry. For multi-layer shallow water, Lambaerts et al. (Reference Lambaerts, Lapeyre and Zeitlin2012) find that the evolution of baroclinically unstable jets in two-layer shallow water exhibits complex vorticity asymmetry evolution, further complicated by their inclusion of moisture.
In this paper, we aim at a balanced model of the shallow water system that could capture the evolution of asymmetry between cyclonic and anticyclonic vortices. We focus on the parameter space where SWQG is formally valid but lacking in captured phenomena. As a by-product, we also wish to capture the finite divergence of shallow water flows, missed by SWQG. Some balanced models already exist in the literature. One line of models includes the frontal geostrophic model by Cushman-Roisin (Reference Cushman-Roisin1986) and the generalised geostrophic equation by Cushman-Roisin & Tang (Reference Cushman-Roisin and Tang1989, Reference Cushman-Roisin and Tang1990). These models capture a more extended parameter space of the shallow water model by allowing for large height deviation and can capture vorticity asymmetry. However, they are not higher order in Rossby number than SWQG, and they all insist that the velocities are in geostrophic balance with the height field and thus do not allow for divergent flow. Another model is the Balance Equation (Allen, Barth & Newberger Reference Allen, Barth and Newberger1990; Spall & Mcwilliams Reference Spall and Mcwilliams1992). It is a capable model of shallow water in the quasi-geostrophic (QG) regime that includes higher-order effects and is useful as a way to initialise a shallow water simulation with minimal generation of gravity waves. We use it to generate the initial conditions for our freely decaying simulation. However, evolving the balance equation is challenging. The system has two time derivatives but is still a balanced model in the sense that there is only one prognostic variable. This is because the two time evolution equations are not independent from each other but coupled in complex nonlinear ways. The time evolution usually needs iterative solvers per time step (Barth, Allen & Newberger Reference Barth, Allen and Newberger1990). This poses challenges when using it to study the evolution of the shallow water model and fully turbulent simulations of it are rare. Yet another model is the semi-geostrophic model of Hoskins (Reference Hoskins1975) adapted to the shallow water context (Cloke & Cullen Reference Cloke and Cullen1994; Roulstone & Sewell Reference Roulstone and Sewell1996). It has features and deficiencies similar to those of the Boussinesq semi-geostrophic model. It has energy and PV conservation and can capture vorticity asymmetry. However, it is not of higher order in Rossby number. The semi-geostrophic model is less accurate in modelling curved fronts. The kinetic energy of the semi-geostrophic model includes only the geostrophically balanced component, making it hard to interpret. The inversion from PV involves a change of coordinate and solving a nonlinear elliptic Monge–Ampère problem, both are hard to do numerically. Practically, Fletcher (Reference Fletcher2004) explored its utility as a model for data assimilation in numerical weather prediction. Theoretically, the mathematical challenge of analysing a Monge–Ampère problem coupled to a transport equation has inspired deep mathematical analysis of the semi-geostrophic model for shallow water (see Cullen & Gangbo Reference Cullen and Gangbo2001; De Philippis & Figalli Reference De Philippis and Figalli2014, and references therein). Finally, the class of optimal balance algorithms deals more broadly with the issue of optimal separation of balanced and unbalanced components of a snapshot of the flow, going beyond asymptotics (Viúdez & Dritschel Reference Viúdez and Dritschel2004; Masur & Oliver Reference Masur and Oliver2020; Chouksey et al. Reference Chouksey, Eden, Masur and Oliver2023). Their focus is not on writing a prognostic system that is balanced for all time.
We explore and expand on a balanced model,
$\textrm {SWQG}^{+1}$
, that extends the SWQG into one higher order in Rossby number in an asymptotically consistent manner. Warn et al. (Reference Warn, Bokhove, Shepherd and Vallis1995) first derived
$\textrm {SWQG}^{+1}$
for the one-layer shallow water system, and Vallis (Reference Vallis1996, hereafter V96) extended it to
$\beta$
-plane. The higher-order-in-Rossby terms model the asymmetric correlation between height deviation and vorticity and the cyclostrophic balance. Recently Chouksey et al. (Reference Chouksey, Eden, Masur and Oliver2023) has studied a variant of the Warn et al. (Reference Warn, Bokhove, Shepherd and Vallis1995) model’s ability for balanced initialisation of the shallow water model and has compared it with the optimal balance algorithm of Masur & Oliver (Reference Masur and Oliver2020). In this paper, we rewrite the Warn et al. (Reference Warn, Bokhove, Shepherd and Vallis1995) model using `potentials’, inspired by the Muraki, Snyder & Rotunno (Reference Muraki, Snyder and Rotunno1999) version of
$\textrm {QG}^{+1}$
of the Boussinesq system. This new potential form allows us to extend the model to multi-layer for the first time. We show through nonlinear simulation of the one and two-layer
$\textrm {SWQG}^{+1}$
system that the model captures the asymmetric evolution of vorticity in the one and two-layer shallow model, as well as flow fields with finite divergence.
The SWQG+1model, like
$\textrm {QG}^{+1}$
in the Boussinesq context, is a PV-based balanced model, where the prognostic variable is PV and all other variables are diagnosed from it. Warn et al. (Reference Warn, Bokhove, Shepherd and Vallis1995) first considered abstractly the problem of constructing higher-order balanced models and determined that a stable model should not expand the prognostic variable. For
$\textrm {SWQG}^{+1}$
this means PV is not expanded in the small Rossby number. All other variables are expanded and `slaved’ to PV by the principle of `PV inversion’ of Hoskins, McIntyre & Robertson (Reference Hoskins, McIntyre and Robertson1985). Choosing PV as the central variable is also based on the fact that PV has been crucial in understanding many features of the shallow water model and SWQG (Bretherton & Haidvogel Reference Bretherton and Haidvogel1976; Haynes & McIntyre Reference Haynes and McIntyre1987; Hoskins Reference Hoskins1991; Dritschel & McIntyre Reference Dritschel and McIntyre2008). In addition to theoretical appeal, having PV as the only prognostic variable guarantees that the model is free of inertial-gravity waves, just like SWQG. The additional diagnostic relations of ageostrophic but balanced velocities can be useful in their own right, allowing for inference of the ageostrophic dynamics from sparse observations. In particular, this can be useful in the new era of high-resolution satellite altimetry observation where geostrophic balance fails in the submesoscale (Penven et al. Reference Penven, Halo, Pous and Marié2014).
The rest of the paper is organised as follows. In § 2, we rederive the one-layer
$\textrm {SWQG}^{+1}$
model using the potentials, which is an alternative form of the model of Warn et al. (Reference Warn, Bokhove, Shepherd and Vallis1995) and V96. We explore its turbulent free decay, compared with the corresponding shallow water model. It is shown that
$\textrm {SWQG}^{+1}$
captures the emergent negative vorticity skewness in the shallow water model. We also track the energy evolution. It is shown that, while theoretical energy conservation law is lacking, the energy of
$\textrm {SWQG}^{+1}$
models the energy evolution of the full shallow water model well. In § 3, we extend the
$\textrm {SWQG}^{+1}$
model to many layers. A simulation of a baroclinically unstable jet demonstrates
$\textrm {SWQG}^{+1}$
can capture the vorticity asymmetry in two layers. Contrary to the results of the one-layer simulations, the two-layer simulation shows the initial baroclinic growth stage has a cyclonic bias for vorticity. It also shows patterns of vorticity and divergence similar to strain-driven fronts. The vortex stretching mechanism leads to the cyclonic vorticity bias. We conclude in § 4 and discuss possible applications and extensions to the model.
2. Single-layer shallow water
This section re-derives the
$\textrm {SWQG}^{+1}$
model, providing an alternative to the Warn et al. (Reference Warn, Bokhove, Shepherd and Vallis1995) approach. By using `potentials’ inspired by the Muraki et al. (Reference Muraki, Snyder and Rotunno1999) model for Boussinesq, our form of
$\textrm {SWQG}^{+1}$
involves three elliptic inversions where the elliptic operator to be inverted is the same as the one in SWQG, the screened Poisson operator. We simulate the free decay of random balanced initial conditions of the
$\textrm {SWQG}^{+1}$
model as well as the shallow water model and compare their turbulent statistics. The
$\textrm {SWQG}^{+1}$
model can capture the vorticity asymmetry in the shallow water model.
2.1. Derivation of the
$\textrm {SWQG}^{+1}$
model
2.1.1. The shallow water system and its properties
We start with the one-layer shallow water model on a
$\beta$
-plane centered at the Coriolis parameter
$f$
, mean layer depth of
$H$
, and gravity acceleration
$g$
. The equations for horizontal velocities
$\boldsymbol{u}=(u,v)$
and water layer height perturbation
$h$
are
\begin{align} & \{\varepsilon \} \left(\frac{\mathrm{D}u}{\mathrm{D}t}\right) - (f+\left\{\gamma\varepsilon\right\} \beta y)v = -g h_x,\end{align}
\begin{align} & {\{\varepsilon \}}\left ( \frac {{\mathrm{D}} v}{{\mathrm{D}} t} \right) + (f+\left \{\gamma \varepsilon \right\} \beta y)u = -g h_y,\end{align}
\begin{align} & {\{\varepsilon \}}\left (\frac {{\mathrm{D}} h}{{\mathrm{D}} t}+h\boldsymbol{\nabla }\boldsymbol{\cdot }u\right) + {\left \{\textit {Bu}\right\}} H\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} = 0. \end{align}
To facilitate future asymptotic analysis but also keep the physical constants for interpretability, we write the equation using the notation of V96. All non-dimensional constants are enveloped in curly brackets, and one can ignore them if one wants only the physical equations. One can instead obtain the non-dimensional equations by removing all physical constants, here
$f$
,
$\beta$
,
$g$
and
$H$
. The non-dimensional numbers in this model are the Rossby number
\begin{align} \varepsilon := \frac {U}{fL}, \end{align}
the Burger number
\begin{align} \textit {Bu} := \frac {gH}{f^2L^2}, \end{align}
and the Charney–Green number
\begin{align} \gamma := \frac {\beta L^2}{U}, \end{align}
which capture the effects of
$\beta$
. In this paper,
$\beta$
is retained in the derivation, but its effect is not explored in the simulations. Here, we scale
$h$
using geostrophic balance
\begin{align} h \sim \frac {fUL}{g} \quad \mbox{and} \quad \frac {h}{H} \sim \frac {\varepsilon }{\textit {Bu}}. \end{align}
The shallow water system conserves the PV
\begin{align} &\frac {{\mathrm{D}} Q}{{\mathrm{D}} t} = 0 \quad \mbox{with} \quad Q = \frac {f+{\{\varepsilon \}}\zeta + \left \{\gamma \varepsilon \right\} \beta y}{H+\{\varepsilon /\textit {Bu}\}h}. \end{align}
As a consequence, it conserves the total potential enstrophy density
\begin{align} \text{Potential Enstrophy} = \frac {1}{2}\big\langle{ (H+\{\varepsilon /\textit {Bu}\}h) Q^2} \big\rangle, \end{align}
where
$\left \langle {\boldsymbol{\cdot }} \right\rangle$
is the area average
\begin{align} \left \langle {\boldsymbol{\cdot }} \right\rangle = \frac {1}{A} \iint _A \boldsymbol{\cdot }\;{\mathrm{d}} x{\mathrm{d}} y. \end{align}
It also conserves the total energy made up of eddy kinetic energy (EKE) and available potential energy (APE) density
\begin{align} \text{Total energy} &= \text{EKE}+\text{APE},\end{align}
where
\begin{align} \quad \mbox{} \quad \text{EKE} &= \frac {1}{2}\big\langle {(H+{\{\varepsilon \}} h)(u^2+v^2)} \big\rangle,\end{align}
\begin{align} \text{APE} &= {\big\{\textit {Bu}^{-1}\big\}}\frac {1}{2}\big\langle {gh^2} \big\rangle . \end{align}
To prepare for the derivation of the
$\textrm {SWQG}^{+1}$
model, we propose a form of representing the three variables of shallow water using three `potentials’
\begin{align} u &= -\varPhi _y-F, \end{align}
\begin{align} v &= \varPhi _x-G, \end{align}
\begin{align} h &= \frac {f}{g}\varPhi -{\left \{\textit {Bu}\right\}}\frac {H}{f}G_x+{\left \{\textit {Bu}\right\}}\frac {H}{f}F_y, \end{align}
where we assume the scaling
\begin{align} \varPhi \sim UL \quad \mbox{and} \quad G = F \sim U. \end{align}
This procedure and naming are inspired by
$\textrm {QG}^{+1}$
for the Boussinesq model where the three potentials are the vector potential form of the three-dimensional incompressible velocity field (Muraki et al. Reference Muraki, Snyder and Rotunno1999; Dù et al. Reference Dù, Smith and Bühler2025). Here, for the shallow water case, the natural counterpart is the two-dimensional Helmholtz decomposition as adopted by Warn et al. (Reference Warn, Bokhove, Shepherd and Vallis1995) and V96. Although our potentials might look unnatural at first glance, it is useful for the derivation of
$\textrm {SWQG}^{+1}$
, where it turns out the elliptic inversion problems for all three potentials are the same as SWQG’s screened Poisson problems. This uniformity guides the derivation of
$\textrm {SWQG}^{+1}$
in the more complex multi-layer case.
2.1.2. The shallow water QG and
$\textrm {QG}^{+1}$
model
The SWQG+1 model is a model that makes the same asymptotic assumption as SWQG but extends it to the next level in Rossby number. That is, we assume the small parameter
\begin{align} \varepsilon \ll O(1). \end{align}
and
\begin{align} \textit {Bu} = O(1) \quad \mbox{and} \quad \gamma \leqslant O(1). \end{align}
Using the small parameter, we can asymptotically expand the PV
\begin{align} HQ &= \frac {f+{\{\varepsilon \}}\zeta +\left \{\gamma \varepsilon \right\} \beta y}{1+\{\varepsilon /\textit {Bu}\}h/H} \end{align}
\begin{align} &= [f+{\{\varepsilon \}}\zeta +\left \{\gamma \varepsilon \right\} \beta y]\left (1-\left \{\frac {\varepsilon }{\textit {Bu}}\right\}\frac {h}{H}+\left \{\frac {\varepsilon ^2}{\textit {Bu}^2}\right\}\frac {h^2}{H^2}\right)+O({\{\varepsilon \}}^3) \end{align}
\begin{align} &= f+{\{\varepsilon \}}\left [(\zeta +\left \{\gamma \right\} \beta y)-{\big\{\textit {Bu}^{-1}\big\}}\frac {fh}{H}\right] \nonumber \\ &\qquad +{\{\varepsilon \}}^2\left [{\big\{\textit {Bu}^{-1}\big\}}^2\frac {fh^2}{H^2}-{\big\{\textit {Bu}^{-1}\big\}}\frac {(\zeta +\left \{\gamma \right\} \beta y)h}{H}\right]+O({\{\varepsilon \}}^3) \end{align}
\begin{align} &= f+{\{\varepsilon \}} q+O({\{\varepsilon \}}^3),\end{align}
where we define
$q$
in the last line. We also expand the potential forms under this QG scaling regime
\begin{align} u &= -\varPhi ^0_y+{\{\varepsilon \}}\big(-\varPhi ^1_y-F^1\big), \end{align}
\begin{align} v &= \varPhi ^0_x+{\{\varepsilon \}}\big(\varPhi ^1_x-G^1\big), \end{align}
\begin{align} h &= \frac {f}{g}\varPhi ^0+{\{\varepsilon \}}\left (\frac {f}{g}\varPhi ^1-{\left \{\textit {Bu}\right\}}\frac {H}{f}G^1_x+{\left \{\textit {Bu}\right\}}\frac {H}{f}F^1_y\right),\end{align}
where the zeroth level is just SWQG where all variables are based on the horizontal streamfunction
$\varPhi ^0$
.
The SWQG emerges naturally from the above asymptotic expansions. It is a PV-based balanced model where the PV is advected by the zeroth-order velocities. That is, we have
\begin{align} \frac {{\mathrm{D}}^0 q}{{\mathrm{D}} t} = 0, \end{align}
where
\begin{align} u^0 &= -\varPhi ^0_y, \quad \mbox{and} \quad v^0 = \varPhi ^0_x. \end{align}
To follow the principle of PV inversion, we diagnose
$\varPhi ^0$
from the PV, working with only the zeroth QG level of (2.14)
\begin{align} q &= (\zeta ^0+\left \{\gamma \right\} \beta y)-{\big\{\textit {Bu}^{-1}\big\}}\frac {fh^0}{H} \end{align}
\begin{align} &= {\nabla} ^2\varPhi ^0-{\big\{\textit {Bu}^{-1}\big\}}\frac {f^2}{gH}\varPhi ^0+\left \{\gamma \right\} \beta y. \end{align}
This is the SWQG. It is more common to separate out the
$\beta$
term from the QGPV and include it instead in the advection equation. That is, we write
\begin{align} &\frac {{\mathrm{D}}^0 q}{{\mathrm{D}} t} + \left \{\gamma \right\}\beta v = 0,\end{align}
where
\begin{align} & q = {\nabla} ^2\varPhi ^0-{\big\{\textit {Bu}^{-1}\big\}}\frac {f^2}{gH}\varPhi ^0 = \mathcal{S}(\varPhi ^0). \end{align}
We define the screened Poisson operator
$\mathcal{S}$
for future notational savings.
To extend the model to the next order in Rossby number, we only need to extend the order in the PV inversion. It is worth repeating that, following the theoretical argument of Warn et al. (Reference Warn, Bokhove, Shepherd and Vallis1995), PV undergoes no asymptotic expansion. The elliptic inversion problem for
$\varPhi ^1$
is obtained from the next order of (2.14)
\begin{align} \mathcal{S}(\varPhi ^1) &= C_q-{\big\{\textit {Bu}^{-1}\big\}}^2\frac {f^3}{g^2H^2}(\varPhi ^0)^2+{\big\{\textit {Bu}^{-1}\big\}}\frac {f}{gH}[{\nabla} ^2\varPhi ^0\varPhi ^0 + \left \{\gamma \right\}\beta y\varPhi ^0]. \end{align}
Here,
$C_q$
is the appropriate constant that makes
$\varPhi ^1$
have zero mean
\begin{align} C_q &= {\big\{\textit {Bu}^{-1}\big\}} \frac {f}{gH} \left\langle {{\big\{\textit {Bu}^{-1}\big\}}\frac {f^2}{gH}|\varPhi ^0|^2-{\nabla} ^2\varPhi ^0\varPhi ^0 - \left \{\gamma \right\}\beta y\varPhi ^0} \right\rangle . \end{align}
The parameters
$\varPhi ^0$
and
$\varPhi ^1$
should have zero mean for all time in the simulation to maintain mass conservation (2.15c
) (assuming without loss of generality that
$\left \langle {h} \right\rangle =0$
). Therefore, for
$\textrm {SWQG}^{+1}$
, the
$\varPhi ^0$
inversion (2.18b
) should use
$q-\left \langle {q} \right\rangle$
. Note that, with
$\beta =0$
, the integral of the right-hand side is proportional to the total energy at the QG level after integration by parts
\begin{align} \text{Total energy}^0 &= \frac {H}{2} \left \langle { {\big\{\textit {Bu}^{-1}\big\}} \frac {f^2}{gH} |\varPhi ^0|^2+|\boldsymbol{\nabla }\varPhi ^0|^2} \right\rangle . \end{align}
The inversion for
$F$
and
$G$
is inspired by the derivation of the `
$\omega$
-equation’ of shallow water, where we form the imbalance equation (Hoskins, Draghici & Davies Reference Hoskins, Draghici and Davies1978). We first write the
$O(\epsilon )$
level of the shallow water system.
\begin{align} & \frac {{\mathrm{D}}^0u^0}{{\mathrm{D}} t} - \left \{\gamma \right\} \beta yv^0 = fv^1 -g h^1_x,\\[-10pt]\nonumber\end{align}
\begin{align} & \frac {{\mathrm{D}}^0v^0}{{\mathrm{D}} t} + \left \{\gamma \right\} \beta yu^0 = -fu^1-g h^1_y,\\[-10pt]\nonumber\end{align}
\begin{align} & \frac {{\mathrm{D}}^0h^0}{{\mathrm{D}} t} = -{\left \{\textit {Bu}\right\}} H\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}^1. \end{align}
We cancel the time derivative using the thermal wind balance between
$v^0$
and
$h^0$
\begin{align} &\left (\frac {{\mathrm{D}}^0 h^0}{{\mathrm{D}} t}\right)_x-\frac {f}{g}\frac {{\mathrm{D}} ^0v^0}{{\mathrm{D}} t}-\frac {f}{g}\left \{\gamma \right\} \beta y u^0 \nonumber \\ =& -{\left \{\textit {Bu}\right\}} H(u^1_{xx}+v^1_{xy})-\frac {f}{g}\big(-g h^1_y - fu^1\big). \end{align}
The left-hand side just involves the zeroth-order terms represented by
$\varPhi ^0$
\begin{align} &{\big\{\textit {Bu}^{-1}\big\}}\! \frac{1}{H}\!\left[\!\left (\frac {{\mathrm{D}}^0 h^0}{{\mathrm{D}} t}\right)_x\!-\!\frac {f}{g}\frac {{\mathrm{D}} ^0v^0}{{\mathrm{D}} t}-\frac {f}{g}\left \{\gamma \right\} \beta y u^0\right]\! =\! {\big\{\textit {Bu}^{-1}\big\}}\frac {f}{gH}\left [J(\varPhi _x^0,\varPhi ^0)\!+\!\left \{\gamma \right\} \beta y \varPhi _y^0\right] ,\end{align}
and the right-hand side is in fact
$\mathcal{S}(F^1)$
\begin{align} &-(u^1_{xx}+v^1_{xy})-{\big\{\textit {Bu}^{-1}\big\}}\frac {f}{gH}\big(-g h^1_y- fu^1\big) = \mathcal{S}(F^1). \end{align}
Together, we have the elliptic inversion for
$F^1$
\begin{align} \mathcal{S}(F^1) = {\big\{\textit {Bu}^{-1}\big\}}\frac {f}{gH}\left [J\big(\varPhi _x^0,\varPhi ^0\big)+\left \{\gamma \right\} \beta y \varPhi _y^0\right]. \end{align}
A similar derivation gives
\begin{align} \mathcal{S}(G^1) = {\big\{\textit {Bu}^{-1}\big\}}\frac {f}{gH}\left [J\big(\varPhi _y^0,\varPhi ^0\big)-\left \{\gamma \right\} \beta y \varPhi _x^0\right]. \end{align}
where
$J(A,B) = A_xB_y-A_yB_x$
.
In summary, the
$\textrm {SWQG}^{+1}$
inversions provides a way to obtain the physical variables (2.15) from knowledge of the PV field
$q$
. It consists of four linear elliptic inversions (2.19b
), (2.20), (2.27) and (2.28). These inversions are all the same screened Poisson problem. This is typical for a model based on asymptotic expansions. The above derivation of
$\textrm {SWQG}^{+1}$
uses the same ingredient as the model in Warn et al. (Reference Warn, Bokhove, Shepherd and Vallis1995) and V96, and they are in fact equivalent. The equivalence of the Boussinesq version of the
$\textrm {QG}^{+1}$
model of the model in the Boussinesq version of the V96 model is shown in Dù et al. (Reference Dù, Smith and Bühler2025). We show the equivalency for the shallow water version explicitly in Appendix A. Here, we comment on two points of interest. The
$\textrm {SWQG}^{+1}$
model is equivalent to the Bolin–Charney balance equation (Bolin Reference Bolin1955; Charney Reference Charney1955)
\begin{align} f \zeta ^1-g{\nabla} ^2 h ^1 &= -2J\big(\varPhi ^0_x,\varPhi ^0_y\big) ,\end{align}
where we have set
$\beta =0$
for simplicity and
$\zeta^1=v_x-u_y$
. This implies that
$\textrm {SWQG}^{+1}$
captures the cyclogeostrophic balance, which has been useful in diagnosing ageostrophic velocity from sea surface height observations (Penven et al. Reference Penven, Halo, Pous and Marié2014). Additionally, we can combine the elliptic inversion for
$F^1$
and
$G^1$
to form an elliptic problem for divergence
\begin{align} \mathcal{S}(u_x+v_y) &= {\big\{\textit {Bu}^{-1}\big\}}\frac {f}{gH}\left [J\big(\varPhi ^0,{\nabla} ^2\varPhi ^0\big)+\left \{\gamma \right\} \beta \varPhi _x^0\right]. \end{align}
This is more commonly called the shallow water
$\omega$
-equation.
The SWQG+1 model uses PV as the prognostic variable. Therefore, it retains the PV conservation of (2.6). However, we have not been able to show that it conserved any definition of energy. Instead, we will demonstrate numerically that the
$\textrm {SWQG}^{+1}$
model generally captures the evolution of the shallow water energy (2.9), compared with the full model. However, it is occasionally not monotonic for the largest Rossby number explored. Since the simulations use dissipation, this does not prove but indicates that
$\textrm {SWQG}^{+1}$
does not conserve the usual shallow water energy. The
$\textrm {SWQG}^{+1}$
model does not conserve potential enstrophy (2.7) either, which additionally requires the height to evolve following the height (2.1c
). The
$\textrm {SWQG}^{+1}$
model’s height is inverted from PV and does not necessarily follow (2.1c
). We also show from simulations that the potential enstrophy (2.7) evolved similarly to a full shallow water simulation. Overall, even though
$\textrm {SWQG}^{+1}$
does not explicitly conserve the energy and potential enstrophy, they mirror the behaviour of the underlying shallow water model for long-time simulations. We argue that theoretical thinking based on the conserved quantities can be applied to the
$\textrm {SWQG}^{+1}$
simulations.
2.2. Freely decaying simulations of the one-layer model
The
$\textrm {SWQG}^{+1}$
model was numerically investigated in V96 as a diagnostic tool. That is, given the true PV field of a shallow water simulation, could
$\textrm {SWQG}^{+1}$
improve the inverted height and velocities, over QG? The conclusion in V96 is only in certain cases. Here, we explore the modelling power of
$\textrm {SWQG}^{+1}$
to capture the emergent vorticity asymmetry in shallow water. We evolve an ensemble of shallow water and
$\textrm {SWQG}^{+1}$
simulations from the same random balanced initial conditions and compare their turbulent statistics.
2.2.1. Set-ups of the simulations
The simulations are set in an idealised mid-latitude ocean mixed layer modelled as a one-and-a-half-layer shallow water system. It is translated to our dimensional parameters of height
$H=100 \text{ m}$
, Coriolis parameter
$f=10^{-4}\text{ s}^{-1}$
and a reduced gravity constant
$g = 1{ \rm m\,s}^{- 2}$
(with the usual prime omitted to match the notation of the rest of the section). These parameters imply a deformation radius
$\sqrt {gH}/f$
of 100 km. The domain is doubly periodic with
$Lx=Ly=12\pi \times 100$
km. All the freely decaying simulations have the Burger number (2.3) equal to one, that is, we non-dimensionalise with
$L=100$
km and the non-dimensional domain is
$12\pi \times 12\pi$
. A typical flow speed at energetic regions of the upper ocean (e.g. the Gulf Stream) is around
$U=1.2$
m s−1. This gives a Rossby number (2.2) of
$\varepsilon =0.12$
while a more quiescent region has
$U=0.1$
m s−1 and
$\varepsilon =0.01$
. Our simulation will explore this range of realistic Rossby numbers. Of course, the power of non-dimensionalisation allows our simulation to apply to the atmosphere as well. We do not explore higher Rossby numbers since too high a Rossby number while keeping the Burger number fixed will result in
$\{\varepsilon /Bu\}h/H\gt 1$
, a drying of part of the domain. This is unrealistic for large-scale ocean and atmospheric applications and cannot be simulated using our pseudospectral code.
The initial conditions are carefully constructed to benefit the study of vortex asymmetry. They should have zero vorticity skewness initially and be well balanced to avoid excessive gravity-wave emission. For this, we use the procedure based on the balance equations used in Polvani et al. (Reference Polvani, McWilliams, Spall and Ford1994, (2.5)). The vorticity is set to be isotropic Gaussian random fields (thus no skewness) centred around non-dimensional size
$4$
and non-dimensional wavenumber
$1.6$
, crudely mimicking energy injected by baroclinic instability. The magnitude is determined so that the QG EKE of the vorticity field is one half
\begin{align} -\frac {1}{2}\big\langle {{\nabla} ^{-2}\zeta \boldsymbol{\cdot }\zeta } \big\rangle = \frac {1}{2}. \end{align}
The fixed-point algorithm described in Polvani et al. (Reference Polvani, McWilliams, Spall and Ford1994) is then performed to solve for the well-balanced divergence and height fields to match the vorticity fields.
Vorticity asymmetry emerges from the evolution of the shallow water model. We evolve the simulation pseudospectrally using Dedalus (Burns et al. Reference Burns, Vasil, Oishi, Lecoanet and Brown2020) for 500 non-dimensional time units. The horizontal resolution is
$1024\times 1024$
Fourier modes. We dissipate the small scale for numerical stability using a fourth-order hyper-dissipation just large enough to absorb small-scale noise, as diagnosed by the vorticity spectrum
$\nu {\nabla} ^4 (\boldsymbol{\cdot })$
with non-dimensional value
$\nu =5.5\times 10^{-6}$
on
$u,v,h$
for the shallow water simulation and on
$q$
for the
$\textrm {SWQG}^{+1}$
simulation. Time stepping uses the third-order 4-stage diagonally implicit+explicit Runge–Kutta (DIRK + ERK) scheme (Ascher, Ruuth & Spiteri Reference Ascher, Ruuth and Spiteri1997). The time step size is determined to have a (CFL) Courant-Friedrichs-Lewy number of 0.5.
2.2.2. Simulation results
Figure 1 shows the comparison of snapshots of vorticity and height for the two models at
$t/T=200$
with
$\varepsilon =0.1$
. The domains are filled with roaming vortices with a bias towards anticyclonic ones. The dominant length scale becomes large compared with the initial condition, which is associated with EKE conversion to APE (cf. figure 4) and the inverse cascade. For the one-layer simulations, balanced evolution of vortices does not generate motions with strong divergence, and its effect on the dynamics is minimal. Therefore, we do not further investigate divergence for the one-layer simulations.
Vorticity (
${\{\varepsilon \}}\zeta /f$
) fields (a) and height (
$\{\varepsilon /Bu\}h/H$
) fields (b) for
$\varepsilon =0.1$
at time
$t/T=200$
from the shallow water simulation (left) and
$\textrm {SWQG}^{+1}$
(right).
Ensemble statistics are appropriate for studying the emergent vorticity asymmetry of the shallow water model and whether the
$\textrm {SWQG}^{+1}$
model can capture it. The left of figure 2 shows the time evolution of vorticity skewness of the size ten ensemble of simulations of the shallow water model and
$\textrm {SWQG}^{+1}$
, for
$\varepsilon =0.1$
. The initial conditions are constructed to have symmetric vorticity, and indeed, the vorticity skewness starts at zero for all simulations. All simulations’ vorticity skewness trends towards negative. The spread between members of the ensembles becomes larger over time, as expected for the chaotic dynamics of shallow water turbulence. However, the ensemble means of vorticity skewness of the two models agree remarkably well during the entire time series. They lie within
$1/\sqrt {10}$
times the ensemble standard deviation from each other for most of the time series. We use this range since it is the standard scaling of error of a Monte Carlo estimate of the ensemble mean. The close agreement indicates strongly that
$\textrm {SWQG}^{+1}$
model can capture the mechanism that leads to vorticity asymmetry in the shallow water model (Arai & Yamagata Reference Arai and Yamagata1994; Graves et al. Reference Graves, McWilliams and Montgomery2006).
(a) time series of vorticity (
$\zeta$
) skewness for the
$\varepsilon =0.1$
simulations from the shallow water model as well as
$\textrm {SWQG}^{+1}$
. The lighter lines are the individual ensemble members. The darker lines are the ensemble mean, and the
$1/\sqrt {10}$
of the ensemble standard deviation is the filled colour around the mean. (b) the vorticity (
$\zeta$
) skewness at
$t/T=200$
for
$\varepsilon =0.01,0.03,0.05,0.07,0.1,0.12$
. The error bar is
$1/\sqrt {10}$
of the ensemble standard deviation.
The agreement extends to all Rossby numbers explored. The right of figure 2 shows the ensemble mean and
$1/\sqrt {10}$
times the ensemble standard deviation of vorticity skewness for the shallow water model and
$\textrm {SWQG}^{+1}$
at non-dimensional time
$t/T=200$
. We see that vorticity skewness becomes more negative as the Rossby number gets larger. Polvani et al. (Reference Polvani, McWilliams, Spall and Ford1994, figure 10) shows a similar trend except with the Froude number, which is equal to the Rossby number for all our simulations where the Burger number is one.
Since
$\textrm {SWQG}^{+1}$
is PV-based, we also explore the PV skewness of the evolution. The left of figure 3 shows that the evolution of PV skewness again matches well between the shallow water and the
$\textrm {SWQG}^{+1}$
model, and the right shows the agreement is uniform overall for all Rossby numbers explored. However, the PV skews positive instead. Therefore, the main contributor of the negative vorticity skewness is the nonlinearity in the full shallow water PV (2.6), where positive vortex correlated with negative height deviation and vice versa. In particular the
$\zeta$
-
$h$
correlation term in the asymptotic expansion of
$q$
(2.14) scales as
\begin{align} \frac {\varepsilon }{Bu} = \frac {U^2}{gH} ,\end{align}
which is equal to the Froude number squared. This is consistent with the fact that vorticity skewness is proportional to the Froude number (Polvani et al. Reference Polvani, McWilliams, Spall and Ford1994).
The same as figure 2 but for PV (
$q$
).
The freely evolving simulation in the shallow water model is a simple setting where we can investigate the properties of the
$\textrm {SWQG}^{+1}$
model. We have not been able to show that the
$\textrm {SWQG}^{+1}$
model conserves energy. Instead, we diagnose the energy and show that it mimics the behaviour of the energy of the shallow water model. The left of figure 4 shows the time series of the total energy as well as the EKE and APE components (2.9) for the evolution of the
$\varepsilon =0.1$
ensemble. The EKE decreases as it is converted to the APE, a signature of the inverse cascade. The total energy is well conserved such that more than 90 % of it is in the final state. The
$\textrm {SWQG}^{+1}$
total energy evolution closely follows the one for the full shallow water model, and this is true for all Rossby numbers explored (not shown). However, it is noticeable that the
$\textrm {SWQG}^{+1}$
loses more energy, compared with the shallow water model. The shallow water model dissipates the physical variables separately, while the
$\textrm {SWQG}^{+1}$
model dissipates the PV. Dissipation in the shallow water model does not translate to dissipation of the PV, because of the nonlinearity in the shallow water PV (2.6). This is well known for the Ertel PV of the Boussinesq system, where the dissipation of the physical variables translates to a source term in the PV equations (Haynes & McIntyre Reference Haynes and McIntyre1987). We will not explore the possibility of including source terms in the PV equation in this work any further. The total energy is minimally dissipated, and the dissipation effects are not significant.
It is important to remark that the total energy of the
$\textrm {SWQG}^{+1}$
model is not monotonic with time, although this is hard to notice in the timeseries in figure 4. This implies that the
$\textrm {SWQG}^{+1}$
model does not conserve the shallow water energy, at least with this choice of dissipation of PV. An alternative reasonable measurement of energy is the total energy at the QG level (2.22), but it does not even decrease overall and has a much larger ensemble spread, as shown by the time series in figure 5. This further shows our simulated regime is beyond the simple QG dynamics.
Potential enstrophy is effectively dwindled away by the small-scale dissipation, as shown by the right of figure 4. The evolutions match well between the shallow water model and
$\textrm {SWQG}^{+1}$
.
We have commented in § 2.1.2 that for
$\textrm {SWQG}^{+1}$
total energy at the QG level (2.22) is proportional to
$C_q$
(2.21). It is the
$\textrm {SWQG}^{+1}$
representation of the mean PV (
$\left \langle {q} \right\rangle$
) of the shallow water model. Figure 5 compares the two. Overall, the two values are similar. This confirms the assumption made in the derivation that the difference will be of
$O(\varepsilon ^2)$
. Mean PV evolves in shallow water due to the correlation between PV and divergence
\begin{align} \partial _t \left \langle {q} \right\rangle = \left \langle {(u_x+v_y)q} \right\rangle . \end{align}
The weak time tendency implies that this correlation term is weak, mostly likely due to the fact that divergence is weak. In fact, for
$\textrm {SWQG}^{+1}$
evolution, the mean of the advected PV (
$\left \langle {q} \right\rangle$
) does not change. Indeed
\begin{align} \left \langle { (u_x+v_y) q } \right\rangle &= \big\langle { (F^1_x+G^1_y)(\mathcal{S}\varPhi ^0 + \left \langle {q} \right\rangle ) } \big\rangle \\[-10pt]\nonumber\end{align}
\begin{align} &= \big\langle { \mathcal{S}(F^1_x+G^1_y)\varPhi ^0 } \big\rangle \\[-10pt]\nonumber\end{align}
\begin{align} &= \big\langle { J({\nabla} ^2\varPhi ^0,\varPhi ^0)\varPhi ^0} \big\rangle = 0 . \end{align}
Therefore,
$C_q$
is the only representation of the evolution of the mean of PV (
$\left \langle {q} \right\rangle$
) in the
$\textrm {SWQG}^{+1}$
model.
Time series of the total energy at the QG level (2.22) in the
$\textrm {SWQG}^{+1}$
simulations, compared with the mean PV (
$\langle {q} \rangle$
) of the shallow water model.
3. Multi-layer shallow water
We extend the
$\textrm {QG}^{+1}$
model to a multi-layer configuration. Here, we only present the two-layer version, but the extension to more layers is straightforward, following the steps laid out here. A simulation of the two-layer version of
$\textrm {SWQG}^{+1}$
under baroclinic instability shows that it can capture the vortex asymmetry in nonlinear baroclinic waves.
3.1. The two-layer shallow water system and the derivation of the two-layer
$\textrm {SWQG}^{+1}$
model
We start with the two-layer shallow water system on the
$\beta$
-plane with two equal layers of height
$H$
and inter-layer reduced gravity
$g'$
\begin{align} &{\{\varepsilon \}}\frac {{\mathrm{D}}_1 \boldsymbol{u}_1}{{\mathrm{D}} t}+( f+\left \{\gamma \varepsilon \right\} \beta y)\boldsymbol{e}_z\times \boldsymbol{u}_1 = -g\boldsymbol{\nabla }(h_1+h_2),\end{align}
\begin{align} &{\{\varepsilon \}}\frac {{\mathrm{D}}_2 \boldsymbol{u}_2}{{\mathrm{D}} t}+( f+\left \{\gamma \varepsilon \right\} \beta y)\boldsymbol{e}_z\times \boldsymbol{u}_2 = -g\boldsymbol{\nabla }(h_1+h_2)-g'\boldsymbol{\nabla }h_2, \end{align}
\begin{align} &{\{\varepsilon \}}\left [\frac {{\mathrm{D}}_1 h_1}{{\mathrm{D}} t}+h_1\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}_1\right]+{\left \{\textit {Bu}\right\}} H\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}_1 = 0, \end{align}
\begin{align} &{\{\varepsilon \}}\left [\frac {{\mathrm{D}}_2 h_2}{{\mathrm{D}} t}+h_2\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}_2\right]+{\left \{\textit {Bu}\right\}} H\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}_2 = 0. \end{align}
where we scale the perturbation height geostrophically
\begin{align} &h_1+h_2\sim \frac {fUL}{g},\end{align}
\begin{align} &h_1\sim h_2\sim \frac {fUL}{g'}. \end{align}
All non-dimensional numbers are the same as those for the one-layer shallow water, except we use the baroclinic Burger number
\begin{align} Bu := \frac {g'H}{f^2L^2}. \end{align}
The system conserves the shallow water PV in each layer (
$i=1,2$
)
\begin{align} HQ_i =& \frac {f+{\{\varepsilon \}} (\zeta _i+\left \{\gamma \right\} \beta y)}{1+\{\varepsilon /\textit {Bu}\} h_i/H} \end{align}
\begin{align} =& f+{\{\varepsilon \}}\left [(\zeta _i+\left \{\gamma \right\} \beta y)-{\big\{\textit {Bu}^{-1}\big\}}\frac {fh_i}{H}\right] \nonumber \\ &\quad\quad+{\{\varepsilon \}}^2\left [{\big\{\textit {Bu}^{-1}\big\}}^2\frac {fh_i^2}{H^2}-{\big\{\textit {Bu}^{-1}\big\}}\frac {(\zeta _i+\left \{\gamma \right\} \beta y)h_i}{H}\right]+O \big ({\{\varepsilon \}}^3\big),\end{align}
where we expanded in the small Rossby number.
The potential form for the two-layer shallow water velocity follows from the one-layer version. The
$\varPhi$
components of the height are modified to the familiar QG form
\begin{align} u_1 &= -\varPhi _{1,y}-F_1,\end{align}
\begin{align} u_2 &= -\varPhi _{2,y}-F_2,\end{align}
\begin{align} v_1 &= \varPhi _{1,x}-G_1,\end{align}
\begin{align} v_2 &= \varPhi _{2,x}-G_2,\end{align}
\begin{align} h_1 &= \frac {f}{g'}\big(\varPhi _1-\varPhi _2\big) + \frac {f}{g}\varPhi _1-{\left \{\textit {Bu}\right\}}\frac {H}{f}G_{1,x}+{\left \{\textit {Bu}\right\}}\frac {H}{f}F_{1,y},\end{align}
\begin{align} h_2 &= \frac {f}{g'}\big(\varPhi _2-\varPhi _1\big)-{\left \{\textit {Bu}\right\}}\frac {H}{f}G_{2,x}+{\left \{\textit {Bu}\right\}}\frac {H}{f}F_{2,y}. \end{align}
To derive SWQG and
$\textrm {SWQG}^{+1}$
, we expand the potential form (3.6), similar to (2.15). The zeroth order of the PV expression (3.5) gives the PV inversion of SWQG
\begin{align} \mathcal{L}\begin{pmatrix} \varPhi _1^0\\[6pt] \varPhi _2^0 \end{pmatrix} = \left(\begin{array}{lllll} \!\!{\nabla} ^2\varPhi _1^0+{\big\{\textit {Bu}^{-1}\big\}}\frac {f^2}{g'H}\big(\varPhi ^0_2-\varPhi ^0_1\big)-{\big\{\textit {Bu}^{-1}\big\}}\frac {f^2}{gH}\varPhi _1^0\\[10pt] {\nabla} ^2\varPhi _2^0+{\big\{\textit {Bu}^{-1}\big\}}\frac {f^2}{g'H}\big(\varPhi ^0_1-\varPhi ^0_2\big)\!\!\!\! \end{array}\right) = \begin{pmatrix} q_1\\ q_2 \end{pmatrix}-\left \{\gamma \right\}\beta y. \end{align}
The operator on the left-hand side will appear in all higher-order inversions, and we name it
$\mathcal{L}$
for convenience.
We follow the same derivation steps as in § 2.1.2 but with more algebra manipulations. Inversions for
$\varPhi ^1$
’s come from the first order of the PV expression (3.5)
\begin{align} \mathcal{L}\begin{pmatrix} \varPhi _1^1\\[6pt] \varPhi _2^1 \end{pmatrix} = \left(\begin{array}{llllll}\!\!\! C_{q,1}-\left [{\big\{\textit {Bu}^{-1}\big\}}^2\frac {f}{H^2}\left (\frac {f}{g'}\big(\varPhi ^0_1-\varPhi ^0_2\big)+\frac {f}{g}\varPhi ^0_1\right)^2\right.\\[8pt] \left .\qquad -{\big\{\textit {Bu}^{-1}\big\}}\frac {1}{H}({\nabla} ^2\varPhi ^0_1+\left \{\gamma \right\} \beta y)\left (\frac {f}{g'}\big(\varPhi ^0_1-\varPhi ^0_2\big)+\frac {f}{g}\varPhi ^0_1\right)\right]\\[8pt] C_{q,2}-\left [{\big\{\textit {Bu}^{-1}\big\}}^2\frac {f}{H^2}\left (\frac {f}{g'}\big(\varPhi ^0_2-\varPhi ^0_1\big)\right)^2\right.\\[8pt] \left .\qquad -{\big\{\textit {Bu}^{-1}\big\}}\frac {1}{H}({\nabla} ^2\varPhi ^0_2+\left \{\gamma \right\} \beta y)\left (\frac {f}{g'}\big(\varPhi ^0_2-\varPhi ^0_1\big)\right)\right]\!\!\!\! \end{array}\right)\!. \end{align}
The
$C_{q,i}$
’s are again to make sure
$\left \langle {\varPhi ^1_i} \right\rangle =0$
.
The elliptic problems for the rest of the potentials
$F$
and
$G$
come from appealing to the thermal wind balance. We form the so-called `imbalance’ equation using the
$O(\varepsilon )$
equations. Then one side becomes elliptic operators applied to the potentials, while the other side only depends on
$\varPhi ^0_{1,2}$
. We leave the details of the derivation to Appendix B. The inversions for
$F^1_1,F^1_2$
are
\begin{align} \mathcal{L}\begin{pmatrix} F^1_1\\[6pt] F^1_2 \end{pmatrix} &= {\big\{\textit {Bu}^{-1}\big\}}\frac {f}{H} \left(\begin{array}{lllll} \! \left .\left [J\big(\varPhi ^0_{1,x}+\varPhi ^0_{2,x},\varPhi ^0_{1}-\varPhi ^0_{2}\big)+\left \{\gamma \right\}\beta y \big(\varPhi _{1}^0-\varPhi _{2}^0\big)_{,y}\right]\right/g'\\[8pt] +\left [J\big(\varPhi ^0_{1,x},\varPhi ^0_{1}\big)+\left \{\gamma \right\}\beta y \varPhi _{1,y}^0\right]/g \\[8pt] \left .\left [J\big(\varPhi _{1,x}^0+\varPhi _{2,x}^0,\varPhi ^0_2-\varPhi ^0_1\big)+\left \{\gamma \right\}\beta y \big(\varPhi _{2}^0-\varPhi _{1}^0\big)_{,y}\right]\right/g' \!\!\!\!\! \end{array}\right)\!, \end{align}
and for
$G_1^1,G_2^1$
are
\begin{align} \mathcal{L}\begin{pmatrix} G^1_1 \\[6pt] G^1_2 \end{pmatrix} &= {\big\{\textit {Bu}^{-1}\big\}}\frac {f}{H} \left(\begin{array}{lllll} \! \left [J\big(\varPhi ^0_{1,y}+\varPhi ^0_{2,y},\varPhi ^0_1-\varPhi ^0_2\big)+\left \{\gamma \right\}\beta y \big(\varPhi _1^0-\varPhi _2^0\big)_{,x}\right]/g'\\[8pt] +\left [J\big(\varPhi ^0_{1,y},\varPhi ^0_1 \big)-\left \{\gamma \right\}\beta y\varPhi _{1,x}^0\right]/g \\[8pt] \left [J\big(\varPhi _{1,y}^0+\varPhi _{2,y}^0,\varPhi ^0_2-\varPhi ^0_1\big)+\left \{\gamma \right\}\beta y\big(\varPhi ^0_2-\varPhi ^0_1\big)_{,x}\right]/g' \!\!\!\!\! \end{array}\right)\!. \end{align}
The multi-layer
$\textrm {SWQG}^{+1}$
inversions again have the same elliptic problems for the potentials.
3.2. Baroclinic unstable jets simulations
To explore the modelling power of multi-layer
$\textrm {SWQG}^{+1}$
, we simulate the evolution of a baroclinically unstable jet. Baroclinic instability is only possible in shallow water with at least two layers. It spins up vorticity features from small noise, at the expense of large-scale potential energy. This is in contrast and in complement to the freely decaying set-up in § 2.2.
3.2.1. Set-ups of the simulations
We set our simulation in the atmospheric context, inspired by the work of Lambaerts et al. (Reference Lambaerts, Lapeyre and Zeitlin2012). The idealised two-layer atmosphere has
$f=10^4\text{ s}^{-1}$
,
$g'=2{ \,\rm m\,s}^{- 2}$
,
$H=10$
km mean thickness and a typical flow speed of
$U=5{\, \rm m\,s}^{- 1}$
. If we non-dimensionalise with a length scale of
$L=500$
km, the non-dimensional number will be
$Bu = 8$
and
$\varepsilon =0.1$
. The domain is doubly periodic with
$Ly=32\times L$
and
$Lx = 64\times L$
.
Initially, the domain is filled with two large baroclinic jets of meridional size
$4000$
km (
$8\times L$
) of maximum speed 13 m s−1, as shown by figure 6. The northern one has a westerly jet in the top layer and an easterly jet in the bottom layer, a rough sketch of the polar jet stream. The southern jet is the opposite and is present only to restore periodicity for the height perturbation. We will focus on the northern jet. The initial heights are set to match the initial jets in geostrophic balance. Since this initial condition only varies in one direction, it is already a solution of the two-layer shallow water model and in perfect balance. No divergence field correction is needed.
The initial jets’ non-dimensional velocity in the upper (
$u_1/U$
) and lower (
$u_2/U$
) layers.
The instability is triggered by the seeding of small magnitude noise, to height in the shallow water simulation and to PV in the
$\textrm {SWQG}^{+1}$
simulation. We again use Dedalus to simulate the nonlinear evolution. The resolution is
$Nx\times Ny = 512\times 256$
modes. Everything else about the numerics is the same as the simulation in § 2.2 except that the fourth-order hyper-dissipation constant is larger with non-dimensional value
\begin{align} \nu =7.3\times 10^{-4}. \end{align}
3.2.2. Simulation results
The
$\textrm {SWQG}^{+1}$
captures the vorticity evolution of the baroclinic life cycle to remarkable accuracy. Figure 7 displays snapshots of the barotropic (
$(\zeta _1+\zeta _2)/2$
) and baroclinic (
$(\zeta _1-\zeta _2)/2$
) vorticity (in colour) and divergence (in contour) during the growth of the baroclinic instability. Note that the time of the snapshot is offset. This is because the growth of the instability from small-scale noise happens at different times in the model, due to the chaotic nature of the instability. However, the growth rate is remarkably similar (see figure 8). The match is indeed impressive. Not only does the
$\textrm {SWQG}^{+1}$
model capture the asymmetric evolution of the cyclonic and anticyclonic part of the baroclinic wave, it is also able to recover the associated divergence fields. Note that one could infer the same divergence, using the
$\omega$
-equation. The difference of
$\textrm {SWQG}^{+1}$
is that the full velocity participates in the advection, while in SWQG, the higher-order divergence is implicit and is not used in the advection of PV. The divergence is of crucial atmospheric interest as a convergent region in the lower layer corresponds to precipitation. For our example, baroclinic convergence correlates with baroclinic cyclonic vorticity, and baroclinic divergence is large near regions of barotropic strain. On the other hand, barotropic divergence is uniformly small. This is reminiscent of the strain-induced fronts configuration studied in Hoskins & Bretherton (Reference Hoskins and Bretherton1972). We will explore further strain-induced fronts in § 3.3. The similarity of the evolution for the two models go beyond this early snapshot. In the supplementary movies are available at https://doi.org/10.1017/jfm.2025.10979, we provide videos of the
${\{\varepsilon \}}=0.1$
simulations.
The vorticity field of the evolution of the jets for the shallow water model (a) and
$\textrm {SWQG}^{+1}$
model (b). The contour is divergence
${\{\varepsilon \}}\delta /f$
at
$[-0.1,-0.005,0.005,0.1]$
, with the solid being positive. The shallow water snapshot is taken at
$t/T=30$
while the
$\textrm {SWQG}^{+1}$
snapshot is taken at
$t/T=24$
.
The match of the vorticity snapshots translates to the ability of the
$\textrm {SWQG}^{+1}$
to capture the vorticity skewness of the simulation. Figure 8 shows the agreement of the overall trend over a large time (
$t/T=100$
), for both layers. The vorticity is now skewed towards positive, in comparison with the freely decaying result. This change in vorticity skewness is due to the new effects of the baroclinic mode now included in the two-layer model. The vortex stretching due to the correlation of baroclinic vorticity and divergence drives the initial growth of positive vorticity skewness shown in figure 8. This effect is entirely missing in the one-layer simulation. The vorticity skewness is an ageostrophic phenomenon. We diagnose the first local maximum of the vorticity skewness, as it is an indicator of the strength of the ageostrophic frontogenesis. Figure 9 shows clearly that this measurement of vorticity skewness is linearly proportional to the Rossby number, and
$\textrm {SWQG}^{+1}$
is effective over a large parameter space.
Our results agree with those of Lambaerts et al. (Reference Lambaerts, Lapeyre and Zeitlin2012) in the upper layer for the onset of instability. However, their lower layer does not have a mean jet while ours does, thus the results differ.
Vorticity skewness of two-layer shallow water (a) and
$\textrm {SWQG}^{+1}$
(b) during the evolution of unstable jets simulations. The circles and crosses mark the first local maximum of vorticity skewness, which is further explored in figure 9.
The first local maximum of vorticity skewness of a set of simulations with varying Rossby numbers.
3.3. Baroclinic strain-driven fronts modelled by
$\textrm {SWQG}^{+1}$
The two-layer simulation reveals that strain-driven fronts are important in the dynamics and contribute to the positive vorticity skewness, which is opposite to the results of our one-layer simulations. This is a new feature of the two-layer shallow water simulations and worth further investigation. Here, we study a minimal model of strain-driven fronts by adapting the set-up of Hoskins & Bretherton (Reference Hoskins and Bretherton1972) to the two-layer model.
We study a strain-driven cold filament across the
$y$
-direction. That is, we ignore the
$x$
-variation for all perturbation fields. The filament has a centre with a large
$h_2$
perturbation. We prescribe the perturbation PV field to be Gaussian
\begin{align} q_1(y) = -q_2(y) = e^{-y^2}, \end{align}
and solve for the physical fields under the response of a barotropic incompressible, irrotational strain field
\begin{align} {U}^M=\alpha x, \qquad {V}^M = -\alpha y, \end{align}
which can be captured by a horizontal streamfunction
\begin{align} {\varPhi }^M = -\alpha xy. \end{align}
The strain field has no imprint on the PV if we take the rigid-lid limit of
$\textrm {SWQG}^{+1}$
. That is, we ignore all terms that are divided by
$g$
. All inversions are not affected by the strain field except for
$G$
(and
$F=0$
), where we have
\begin{align} \mathcal{L}\begin{pmatrix} G^1_1 \\[6pt] G^1_2 \end{pmatrix} &= {\big\{\textit {Bu}^{-1}\big\}}\frac {2\alpha f}{g'H}\begin{pmatrix} \big(\varPhi ^0_2-\varPhi ^0_1\big)_{,y} \\[6pt] \big(\varPhi ^0_1-\varPhi ^0_2\big)_{,y} \end{pmatrix}. \end{align}
Note that the operator
$\mathcal{L}$
now only has
$y$
-derivatives. Without the strain field (
$\alpha =0$
), there would be no response since
$G=0$
. This is expected since the filament is a geostrophically balanced solution to the shallow water equations.
The Hoskins & Bretherton (Reference Hoskins and Bretherton1972) model is posed in an infinite domain in
$y$
, where perturbation velocities decay at infinity. This implies boundary conditions for the potentials
\begin{align} &\left .G^{1}\right|_{y=\pm \infty } = 0,\end{align}
\begin{align} &\left .\varPhi ^{0,1}_{,y}\right|_{y=\pm \infty } = 0. \end{align}
Numerically, we solve for the potentials in a finite domain
$y/L\in [-5,5]$
, which is large enough so that the solution is not changed with a large domain. It is represented numerically using
$N=256$
Chebyshev modes, using the Dedalus solver.
We solve the system with
$Bu=1$
and
$\varepsilon =0.3$
. Figure 10 shows the divergence
\begin{align} \delta _1 = v_{1,y} &= -{\{\varepsilon \}} G_{1,y} \quad \mbox{and} \quad \delta _2 = v_{2,y} = -{\{\varepsilon \}} G_{2,y} \end{align}
overlaid on the height perturbations. In the top layer, a cold filament has cyclonic vorticity, which is correlated with convergence. This correlation is a key feature of fronts and explains the cyclonic bias of vorticity in the two-layer simulations by appealing to vortex stretching: baroclinic cyclonic vorticity is correlated with baroclinic negative divergence. Note that the above equation is equivalent to the ‘
$\omega$
-equation’. Therefore, it is implicit in the SWQG model. However, SWQG does not use this ageostrophic velocity in its advective evolution. It does not capture the vorticity skewness evolution.
Finally, a comparison with the surface buoyancy gradient driven Boussinesq frontogenesis of Hoskins & Bretherton (Reference Hoskins and Bretherton1972) is warranted. In the layered shallow water set-up, the convergence will smooth away the height perturbation and release the stored APE. Shallow water fronts would be less strong with no finite time singularity, compared with the Boussinesq case.
Divergence of a cold filament driven by a barotropic strain, modelled by a two-layer rigid-lid
$\textrm {SWQG}^{+1}$
model. The domain is divided vertically into two layers by the height perturbation. The colour shows divergence values. This is a zoomed-in view showing
$y/L = [-3,3]$
.
4. Conclusion and outlook
In this paper, we have recapitulated and derived the PV-based balanced model,
$\textrm {SWQG}^{+1}$
, for one- and multiple-layer shallow water systems. Through simulations, we have shown that
$\textrm {SWQG}^{+1}$
can capture the balanced ageostrophic effects of the shallow water model in the finite Rossby number regime. In particular, it can model the negative skewness of vorticity in freely decaying one-layer shallow water, as well as the positive skewness in the baroclinic instability of two-layer shallow water. It does this by capturing the correlation between vorticity, divergence and height for coherent vortices and baroclinic fronts. The time evolution of the energy and potential enstrophy is close to that of the full shallow water model, even though we cannot prove exact conservation. This shows that
$\textrm {SWQG}^{+1}$
has potential as a tool for the study of many GFD problems in the atmosphere and ocean on Earth and other planets. Its dynamical evolution is free of inertial-gravity waves by design, focusing the model on transport active balanced motion. The diagnostic relations between ageostrophic physics can be useful in inverse problems when observed physics is limited. There is much work to be done to adapt and extend the
$\textrm {SWQG}^{+1}$
model. The rest of this section describes some possible applications and the desirable extensions to the model.
One central challenge in the parametrisation of ocean eddies is to model correctly the mean location and variability of the ocean jets (Chassignet & Marshall Reference Chassignet and Marshall2008). Thiry et al. (Reference Thiry, Li, Mémin and Roullet2024) show that SWQG is inaccurate in modelling the jet location of a barotropic wind-driven gyre, perhaps the simplest model of the general ocean circulation. Therefore, QG does not capture enough phenomenology to be sufficient for eddy parameterisation. However, PV-based thinking and balanced thinking have enabled much progress in eddy parametrisation and should not be abandoned. The
$\textrm {SWQG}^{+1}$
model offers a nice intermediate option. It is still simple like SWQG, where one prognostic PV variable governs the entire model evolution. However, it is able to capture the important balanced ageostrophic effects. Future work should aim to replicate the asymmetric lean of the jet in a wind-driven gyre in the
$\textrm {SWQG}^{+1}$
model. Boundary conditions for the
$\textrm {SWQG}^{+1}$
model can be adapted from the boundary condition specifications in Muraki et al. (Reference Muraki, Snyder and Rotunno1999).
The shallow water model in one layer (Lahaye & Zeitlin Reference Lahaye and Zeitlin2016) and two layers (Lambaerts et al. Reference Lambaerts, Lapeyre and Zeitlin2012; Bembenek et al. Reference Bembenek, Straub and Merlis2020) has been a useful simple model to study the nonlinear and non-smooth effects of moisture in the mid-latitude atmosphere. The SWQG model has also been adapted to include moisture effects. In particular, definitions of moist PV have been useful to interpret model results (Lapeyre & Held Reference Lapeyre and Held2004; Smith & Stechmann Reference Smith and Stechmann2017; Lutsko & Hell Reference Lutsko and Hell2021; Brown, Pauluis & Gerber Reference Brown, Pauluis and Gerber2023). The
$\textrm {SWQG}^{+1}$
model, with an appropriate definition of moist PV, can be a useful bridge between these two models. The removal of gravity waves removes an undesirable complication in the study of the large-scale evolution of storm tracks. More broadly, it would be interesting to see if moist effects can be sufficiently captured by inverting from the moist PV. That is, are moist effects fundamentally `balanced’ or not?
The research on the effect of ocean topography on the coherent structure of the ocean eddy has enjoyed a resurgence (Solodoch, Stewart & McWilliams Reference Solodoch, Stewart and McWilliams2021; Siegelman & Young Reference Siegelman and Young2023; LaCasce et al. Reference LaCasce, Palóczy and Trodahl2024). Theoretical understanding of this problem typically comes from SWQG, or the even simpler two-dimensional Euler equation (e.g. Bretherton & Haidvogel Reference Bretherton and Haidvogel1976; Siegelman & Young Reference Siegelman and Young2023). However, the necessary symmetry between differently signed vortices in the SWQG model is unrealistic. The study of the realistic formation and destruction of coherent vortices requires numerical simulation using complex ocean models.
$\textrm {SWQG}^{+1}$
and Boussinesq
$\textrm {QG}^{+1}$
can bridge this gap. The
$\textrm {SWQG}^{+1}$
model improves upon SWQG and allows for the modelling of realistic evolution of the coherent vortices (although it still assumes asymptotically small topography height). In the meantime, there is still only one variable to predict, PV, and it is possible that PV-based theoretical thinking can be applied to
$\textrm {SWQG}^{+1}$
. Of course, the study of coherent structure on other planets can also benefit from upgrading to the
$\textrm {SWQG}^{+1}$
model (Siegelman et al. Reference Siegelman, Young and Ingersoll2022). For these applications, the
$\textrm {SWQG}^{+1}$
model needs to be extended and tested on the cases with non-constant
$f$
, topography and on spherical domains (Cho & Polvani Reference Cho and Polvani1996).
From these examples, it is clear that there are many applications of
$\textrm {SWQG}^{+1}$
and much work to be done. The results in this paper serve as the bare minimum showcase of its ability to capture phenomena of interest of one fundamental model of GFD, the shallow water model.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2025.10979.
Acknowledgements
We thank O. Bühler, G. Lapeyre, and O. Pauluis for helpful discussions, without implying their endorsement. This work was supported in part through the NYU IT High Performance Computing resources, services, and staff expertise. We thank Shenglong Wang for his work on configuring Dedalus on the NYU Greene cluster. We thank three anonymous reviewers for their helpful suggestions.
Funding
K.S.S. and R.D. acknowledge support from NASA contract 80NSSC20K1142.
Declaration of interests
The authors report no conflict of interest.
Data availability statements
The code that generates the data and figures in this paper is available on GitHub: https://github.com/Empyreal092/SWQGp1_freedecay_BCIjets_Public.
Appendix A. Equivalency of the
$\textrm {SWQG}^{+1}$
potential form to the model of V96
We show that our model is equivalent to the model in Warn et al. (Reference Warn, Bokhove, Shepherd and Vallis1995) and the extended version under the
$\beta$
-approximation in § 2(a) of V96 up to the constant terms in the PV inversions.
First, the two models are the same at the QG level. Then the V96 model inverts for the next-order quantities using their (V15–V17). We show that the two models are equivalent by recovering (V15–V17) from the
$\textrm {SWQG}^{+1}$
inversions. V96’s (V15) is the next-order PV inversion by noting its left-hand side is
$\mathcal{S}(\varPhi ^1)$
while the right-hand side is equal to the right-hand side of (2.20). For (V17), the shallow water
$\omega$
-equation, we have
\begin{align} \mathcal{S}(u_x+v_y) &= -\mathcal{S}(F^1_x+G^1_y) \nonumber \\ &= {\big\{\textit {Bu}^{-1}\big\}}\frac {f}{gH}\left [J\big(\varPhi ^0,{\nabla} ^2\varPhi ^0\big)+\left \{\gamma \right\} \beta \varPhi _x^0\right]. \end{align}
To recover (V16), the Bolin–Charney balance equation (Bolin Reference Bolin1955; Charney Reference Charney1955), we have
\begin{align} f \zeta ^1-g{\nabla} ^2 h ^1 &= f ({\nabla} ^2\varPhi ^1-G^1_x+F^1_y) - f{\nabla} ^2\varPhi ^1 - {\left \{\textit {Bu}\right\}}\frac {gH}{f}(-{\nabla} ^2 G^1_x+{\nabla} ^2 F^1_y) \end{align}
\begin{align} &= {\left \{\textit {Bu}\right\}}\frac {gH}{f}\left [\mathcal{S}(G^1_x)-\mathcal{S}(F^1_y)\right] \end{align}
\begin{align} &= -2J\big(\varPhi ^0_x,\varPhi ^0_y\big)-\left \{\gamma \right\} \beta \varPhi _y^0-\left \{\gamma \right\} \beta y {\nabla} ^2\varPhi ^0,\end{align}
where we used the elliptic problem for
$F^1$
and
$G^1$
. The last term on the right-hand side is not included in (V16). This is a mistake in V96 which stems from the same term missing in the ‘divergence equation’ (V6). Forming the divergence of the shallow water momentum (2.1) reveals (V6) should additionally include a
$-\beta y \zeta$
term.
Appendix B. The derivations for the elliptic inversions of
$\boldsymbol{F_1^1}$
,
$\boldsymbol{F_2^1}$
,
$\boldsymbol{G_1^1}$
and
$\boldsymbol{G_2^1}$
for the two-layer
$\boldsymbol{\textrm {QG}^{+1}}$
model
We start with the two-layer shallow water system (3.1) and the potential form for the physical variables in the two-layer case (3.6). For the QG-level variables, they are in thermal wind balance. For
$h_{1,x}$
, this is
\begin{align} h^0_{1,x} &= \frac {f}{g'}\big(\varPhi ^0_{1,x}-\varPhi ^0_{2,x}\big)+\frac {f}{g}\varPhi ^0_{1,x}= \frac {f}{g'}(v^0_1-v^0_2)+\frac {f}{g}v^0_1. \end{align}
The ‘imbalance’ equation is the time evolution equation of the difference between the thermally balanced variables, that is, the ‘imbalance’. For the QG-level geostrophically balanced variables, the difference should be zero for all time. This allows us to form a diagnostic equation free of time partial derivatives
\begin{align} &\!\!\left(\frac{{\mathrm{D}}^0_1 h^0_1}{{\mathrm{D}} t}\right)_x\! -\! \frac {f}{g'}\left (\frac {{\mathrm{D}}^0_1 v_1^0}{{\mathrm{D}} t}-\frac {{\mathrm{D}}^0_2 v_2^0}{{\mathrm{D}} t}+\left \{\gamma \right\}\beta y u_1^0-\left \{\gamma \right\}\beta y u_2^0\right)-\frac {f}{g}\left (\frac {{\mathrm{D}}^0_1 v_1^0}{{\mathrm{D}} t}+\left \{\gamma \right\}\beta y u_1^0 \right)\nonumber \\ &= \!-\!{\left \{\textit {Bu}\right\}} H(u^1_{1,xx}+v^1_{1,xy})-\!\frac{f}{g'}\!\left(-fu_1^1+fu_2^1 \right)\!-\!\frac{f}{g}\left ( -fu_1^1-gh^1_{1,y} \right) \end{align}
\begin{align} &\Rightarrow {\nabla} ^2 F^1_{1}+ {\big\{\textit {Bu}^{-1}\big\}}\frac {f^2}{g'H}\left ( F_2^1-F_1^1 \right)-{\big\{\textit {Bu}^{-1}\big\}}\frac {f^2}{gH} F_1^1 \end{align}
\begin{align} &= {\big\{\textit {Bu}^{-1}\big\}}\frac {f}{H}\left (\left .\left [J\big(\varPhi ^0_{1,x}+\varPhi ^0_{2,x},\varPhi ^0_{1}-\varPhi ^0_{2}\big)+\left \{\gamma \right\}\beta y \big(\varPhi _{1}^0-\varPhi _{2}^0\big)_{,y}\right]\right/g'\right. \end{align}
\begin{align} &\qquad \qquad \quad \left .\left .+\left [J\big(\varPhi ^0_{1,x},\varPhi ^0_{1}\big)+\left \{\gamma \right\}\beta y \varPhi _{1,y}^0\right]\right/g\right). \end{align}
We see that the right-hand side becomes elliptic operators applied to the potentials, while the left-hand side only depends on
$\varPhi ^0_{1,2}$
.
For
$F_2^1$
, we form the `imbalance’ equation
\begin{align} &\left (\frac {{\mathrm{D}}^0_2 h^0_2}{{\mathrm{D}} t}\right)_x - \frac {f}{g'}\left (\frac {{\mathrm{D}}_2^0 v_2^0}{{\mathrm{D}} t}-\frac {{\mathrm{D}}_1^0 v_1^0}{{\mathrm{D}} t}+\left \{\gamma \right\}\beta y u_2^0-\left \{\gamma \right\}\beta y u_1^0\right) \nonumber \\=& -{\left \{\textit {Bu}\right\}} H(u^1_{2,xx}+v^1_{2,xy})- \frac {f}{g'}\left ( -fu_2^1-g'h^1_{2,y}+fu_1^1 \right)\!. \end{align}
Similar steps gives
\begin{align} &{\nabla} ^2 F^1_{2}+ {\big\{\textit {Bu}^{-1}\big\}}\frac {f^2}{g'H}\left ( F_1^1-F_2^1 \right) \nonumber \\=& {\big\{\textit {Bu}^{-1}\big\}}\frac {f}{g'H}\left [J\left (\varPhi _{1,x}^0+\varPhi _{2,x}^0,\varPhi ^0_2-\varPhi ^0_1\right)+\left \{\gamma \right\}\beta y \big(\varPhi _{2}^0-\varPhi _{1}^0\big)_{,y}\right]\! . \end{align}
Together with (B5) gives the compact form (3.9).
For
$G_1^1$
and
$G_2^1$
\begin{align} &\left (\frac {{\mathrm{D}}^0_1 h^0_1}{{\mathrm{D}} t}\right)_y - \frac {f}{g'}\left (\frac {{\mathrm{D}}^0_2 u_2^0}{{\mathrm{D}} t}-\frac {{\mathrm{D}}^0_1 u_1^0}{{\mathrm{D}} t}-\left \{\gamma \right\}\beta yv_2^0+\left \{\gamma \right\}\beta yv_1^0\right)+\frac {f}{g}\left (\frac {{\mathrm{D}}^0_1 u_1^0}{{\mathrm{D}} t}-\left \{\gamma \right\}\beta yv_1^0 \right) \nonumber \\ =& -{\left \{\textit {Bu}\right\}} H(u^1_{1,xy}+v^1_{1,yy})- \frac {f}{g'}\left ( fv_2^1-fv_1^1 \right)+\frac {f}{g}\left (fv_1^1-gh^1_{1,x}\right) \end{align}
\begin{align} &\left (\frac {{\mathrm{D}}^0_2 h^0_2}{{\mathrm{D}} t}\right)_y - \frac {f}{g'}\left (\frac {{\mathrm{D}}^0_1 u_1^0}{{\mathrm{D}} t}-\frac {{\mathrm{D}}^0_2 u_2^0}{{\mathrm{D}} t}-\left \{\gamma \right\}\beta y v_1^0+\left \{\gamma \right\}\beta y v_2^0\right) \nonumber \\ =& -{\left \{\textit {Bu}\right\}} H(u^1_{2,xy}+v^1_{2,yy})- \frac {f}{g'}\left ( fv_1^1-fv_2^1+g'h^1_{2,x} \right). \end{align}
Standard but somewhat tedious algebraic manipulation, using the shallow water version of the thermal wind balance gives (3.10).
Open access
