# Usual Approximations to the Equations of Atmospheric Motion: A Variational Perspective

Marine Tort, Thomas Dubos

► To cite this version:

Marine Tort, Thomas Dubos. Usual Approximations to the Equations of Atmospheric Motion: A Variational Perspective. Journal of the Atmospheric Sciences, 2014, 71 (7), pp.2452-2466. 10.1175/jas-d-13-0339.1. hal-01088889

HAL Id: hal-01088889

https://hal.science/hal-01088889

Submitted on 12 Nov 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Usual Approximations to the Equations of Atmospheric Motion: A Variational Perspective

MARINE TORT AND THOMAS DUBOS

IPSL/Laboratoire de Météorologie Dynamique, Ecole Polytechnique, Palaiseau, France

(Manuscript received 25 October 2013, in final form 25 February 2014)

## ABSTRACT

The usual geophysical approximations are reframed within a variational framework. Starting from the Lagrangian of the fully compressible Euler equations expressed in a general curvilinear coordinates system, Hamilton's principle of least action yields Euler-Lagrange equations of motion. Instead of directly making approximations in these equations, the approach followed is that of Hamilton's principle asymptotics; that is, all approximations are performed in the Lagrangian. Using a coordinate system where the geopotential is the third coordinate, diverse approximations are considered. The assumptions and approximations covered are 1) particular shapes of the geopotential; 2) shallowness of the atmosphere, which allows for the approximation of the relative and planetary kinetic energy; 3) small vertical velocities, implying quasi-hydrostatic systems; and 4) pseudoincompressibility, enforced by introducing a Lagrangian multiplier.

This variational approach greatly facilitates the derivation of the equations and systematically ensures their dynamical consistency. Indeed, the symmetry properties of the approximated Lagrangian imply the conservation of energy, potential vorticity, and momentum. Justification of the equations then relies, as usual, on a proper order-of-magnitude analysis. As an illustrative example, the asymptotic consistency of recently introduced shallow-atmosphere equations with a complete Coriolis force is discussed, suggesting additional corrections to the pressure gradient and gravity.

## 1. Introduction

Numerical models for weather prediction and global climate seek to simulate the behavior of the atmosphere by using accurate representations of the governing equations of motion, thermodynamics, and continuity. The governing equations of motion can be approximated using geometrical or dynamical order-of-magnitude arguments but the retained equation set has also to be dynamically consistent in the sense that it possesses conservation principles for mass, energy, absolute angular momentum (AAM), and potential vorticity. For instance, the widely used hydrostatic primitive equations (HPE) make use of the following approximations:

* the spherical geopotential approximation, whereby the small angle between the radial direction and the local vertical is neglected;

* the shallow-atmosphere approximation, whereby the distance to the center of the earth is assumed constant, simplifying many metric terms arising when expressing the equations of motion in spherical coordinates;

* the traditional approximation, which neglects those components of the Coriolis force that vary as the cosine of the latitude; and

* the hydrostatic approximation, which neglects some terms in the vertical momentum budget, turning vertical velocity into a diagnostic quantity.

The HPE describe quite accurately large-scale atmospheric and oceanic motions. Furthermore, they filter out the acoustic waves supported by the fully compressible Euler equations, which avoids certain numerical difficulties. For certain applications like high-resolution global weather forecasting, the use of the hydrostatic approximation becomes inappropriate. Hence, less drastic approximations have been sought to filter out the acoustic waves: the soundproof approximations share the feature that the relationship between density and pressure is suppressed, while a more or less accurate representation of the relationship between density and entropy/potential

<footnote>Corresponding author address: Marine Tort, Laboratoire Météorologique, Ecole Polytechnique, Route de Saclay, 91128 Palaiseau CEDEX, France.
E-mail: marine.tort@lmd.polytechnique.fr</footnote>
temperature is retained (see Ogura and Phillips 1962; Lipp and Hemler 1982; Durran 1989, 2008; Klein and Pauluis 2012; Cotter and Holm 2013).

As more accurate equations sets are sought, it becomes desirable to also relax the spherical-geopotential approximation in order to take into account the flattening of the planet, which also implies a latitudinal variation of the gravity acceleration *g* between the poles and the equator. The flattening at the pole of the giant gas planets Saturn and Jupiter could have important dynamical effects on the large-scale atmospheric motion because of their high speed rotation rate. It may be worthwhile, then, to include this effect by allowing a nonspherical geopotential. Gates (2004) first derived such equations of motion using oblate spheroidal coordinates. Unfortunately, this coordinate system leads to the wrong sign for the variation of *g* between the poles and the equator. Richer coordinate systems were suggested to overcome this problem. White et al. (2008) have introduced a similar oblate spheroid geometry that allows qualitatively correct, but quantitatively incorrect variations of *g* between the poles and the equator. White and Inverarity (2012) have proposed a quasi-spheroidal geometry for which the resulting ratio of *g* between the poles and the equator is unity and, in that sense, will not be useful to model meridional gravity variation. Nevertheless, it could be relevant to quantify geometric differences comparing to purely spherical geometry. Very recently, Bénard (2014) presented a “fitted oblate spheroid” coordinate system relevant for global numerical weather prediction. This coordinate system has the advantage of being defined analytically and allowing a realistic horizontal variation of *g*. Last, White and Wood (2012) have derived the equations of motion using a general orthogonal coordinate system, subject only to the assumption of zonal symmetry, extending their previous work (White et al. 2005) to zonally symmetric (i.e., axisymmetric) geopotential.

Moreover, while the dynamical effects of the non-traditional Coriolis force are not fully understood, several studies have demonstrated its important role for certain geophysical and astrophysical applications (Gerkema et al. 2008). Particularly, oceanic equatorial flows are subjected to nontraditional dynamical effects (Hua et al. 1997; Gerkema and Shrira 2005). Closely related, the large depth of the atmosphere should be taken into account to specifically model other planets such as Jupiter or Saturn (or Saturn's planet-like moon Titan) (Gerkema et al. 2008).

Thus, for certain applications some of the usual approximations may not be satisfactory, which raises the question of whether and how they can be relaxed, fully or partially, and combined together, without compromising the model consistency.

The dynamical consistency of a model can be checked by explicitly deriving the relevant budgets. Within the approximation of a spherical geopotential, four dynamically consistent approximated models correspond to whether the shallow-atmosphere and hydrostatic approximations are individually made or not made (Phillips 1966; White and Bromley 1995; White et al. 2005). These authors use a combination of intuition and ingenuity to identify the terms that need to cancel each other in the various budgets.

However, it can be more straightforward to derive the approximated equations following the approach of Hamilton's principle asymptotics (Holm et al. 2002): all approximations are performed in the Lagrangian, then Hamilton's principle of least action produces the equations of motion following standard variational calculus (Morrison 1998). The desired conservation properties are ensured by the symmetry properties of the approximated Lagrangian. This approach was used recently to derive nontraditional shallow-atmosphere equations, that is, shallow-atmosphere equations with a complete Coriolis force representation (Tort and Dubos 2013). In addition to the nontraditional $\cos\phi$ Coriolis force part, extra terms need to be taken into account in the equations of motion to restore the angular momentum budget. The physical origin of those terms is not trivial and, in fact, arises from the vertical dependence of planetary angular momentum, of which the $\cos\phi$ Coriolis force is only one aspect.

Although many known approximate systems have been shown to derive from Hamilton's principle, this has typically been done in hindsight (Müller 1989; Roulstone and Brice 1995; Cotter and Holm 2013). For example, the variational formulation of the anelastic and pseudoincompressible approximations has only recently been obtained (Cotter and Holm 2013).

The overarching goal of the present work is to frame the above-mentioned approximations within a systematic framework starting from the unapproximated compressible Euler equations, with very mild assumptions regarding the geopotential field. Hamilton's principle of least action, despite its perceived technicality, is the ideal tool for this. Fortunately, for the purpose just stated, it is possible to invoke Hamilton's principle just once with a simple, but sufficiently general, form of the Lagrangian. This leads to the Euler-Lagrange equations of motion in (13). This sufficiently general form relies on general curvilinear coordinates, in order to be able to use the geopotential later as a vertical coordinate.

The necessary notations are introduced in section 2, and the conservation laws are obtained from the Euler-Lagrange equations in (13) without further variational calculus. The next step is to actually construct a curvilinear
system where the geopotential is a vertical coordinate. This problem is addressed in section 3. Then, the dominant force—gravity—acts only in an accurately defined vertical direction, and it becomes possible to simplify the equations of motion without jeopardizing the conservation laws by approximating directly the Lagrangian itself. This is done in section 4. Many well-known approximate systems of equations are “rediscovered” this way, a number of which had already been formulated from a variational principle. Nevertheless, we still obtain new variational formulations for recently derived approximate systems (White and Wood 2012; Klein and Pauluis 2012). Furthermore, a new set of shallow-atmosphere nontraditional equations in a zonally symmetric, nonspherical geopotential is derived, combining White and Wood (2012) and Tort and Dubos (2013). As in Tort and Dubos (2013), the derivation is based on an asymptotic expansion of kinetic energy and planetary terms. In section 5, a more general discussion addresses the asymptotic consistency of the complete Lagrangian, especially between the terms retained—neglected in the kinetic and Coriolis terms, and those retained—neglected in potential and internal energy. The main results are then summarized in section 6.

## 2. Euler–Lagrange equations of motion in general curvilinear coordinates

### a. The action functional

Hamilton’s principle of least action states that flows satisfying the equations of motion render the action stationary; that is, $\delta \int \mathcal{L} dt = 0$, where the Lagrangian $\mathcal{L}$ is defined as the mass-weighted integral of a Lagrangian density $L(\mathbf{r}, \rho, s, \dot{\mathbf{r}})$:

$$ \mathcal{L} = \int_V L(\mathbf{r}, \rho, s, \dot{\mathbf{r}}) dm, \quad dm = \rho d^3\mathbf{r}, \quad (1) $$

where $\mathbf{r}$ is the position within a Cartesian frame attached to the planet, $V$ is the spatial domain containing the fluid of density $\rho$, and $[t_0, t_1]$ is the time domain. Notice that $L$ is a function of $\mathbf{r}, \rho, s$, and $\dot{\mathbf{r}}$ only. This is a restriction to the family of equations that can be considered. As will become apparent, a wide-ranging family of approximated equations can be derived from this restricted form of the action.

We follow Morrison (1998) and adopt the Lagrangian point of view. Fluid parcels are identified by their Lagrangian labels $\mathbf{a} = (a_1, a_2, a_3)$. We use $\mathbf{r}(\mathbf{a}, \tau)$ as the position of a fluid parcel and $\dot{\mathbf{r}} = \partial\mathbf{r}(\mathbf{a}, \tau)/\partial\tau$ as its three-dimensional velocity; they are both functions of labels $\mathbf{a}$ and time $\tau$. The variable $\tau$ is used to emphasize that partial time derivatives $\partial/\partial\tau$ are taken at fixed particle labels $\mathbf{a}$, not at fixed spatial coordinates, so that

$\partial/\partial\tau = D/Dt$ is in fact the Lagrangian time derivative. Furthermore, the mass of an infinitesimal volume surrounding a fluid parcel is $dm = \mu d^3\mathbf{a} = \rho d^3\mathbf{r}$, where $\mu = \rho\det(\partial\mathbf{r}/\partial\mathbf{a})$ does not depend on time and is therefore determined by the initial value of $\rho$ and $\mathbf{r}$. When invoking Hamilton’s principle, $\int \mathcal{L} d\tau$ is considered to be a functional of the label-time field $\mathbf{r}(\mathbf{a}, \tau)$. Variations $\delta u^k$ and $\delta\rho$ can be expressed in terms of variations $\delta\mathbf{r}$ taken at fixed Lagrangian labels. Variations $\delta\mathbf{r}$ vanish at $\tau = t_0, t_1$.

Letting $e(\alpha, s)$ be the specific internal energy with $\alpha = 1/\rho$ (specific volume), $s$ is the specific entropy, $p = -\partial e/\partial\alpha$ (pressure), and $T = \partial e/\partial s$ (temperature), the compressible Euler equations with Coriolis force result from the Lagrangian density $L(\rho, \dot{\mathbf{r}}, s; \mathbf{r}, t)$:

$$ L(\rho, \dot{\mathbf{r}}, s; \mathbf{r}, t) = \frac{1}{2} \dot{\mathbf{r}}^2 + \dot{\mathbf{r}} \cdot \mathbf{R}(\mathbf{r}) - e(\rho, s) - \Phi(\mathbf{r}), \quad (2) $$

where $\mathbf{R}(\mathbf{r}) = \mathbf{\Omega} \times \mathbf{r}$ is the solid-body velocity due to the planetary rotation $\mathbf{\Omega}$ and the geopotential $\Phi(\mathbf{r})$ is the sum of the gravitational and centrifugal potentials. In this section, (2) is expressed in a general curvilinear coordinate system. Hamilton’s principle of least action then yields the equations of motion.

### b. Motion and transport in general curvilinear coordinates

We now consider general curvilinear coordinates $\xi^k$, that is, a mapping $(\xi^k) \mapsto \mathbf{r}(\xi^k)$. The chain rule shows that motion in the curvilinear system $\xi^k$ is described by the Lagrangian derivatives:

$$ u^k = \frac{D\xi^k}{Dt}, \quad (3) $$

$$ \dot{\mathbf{r}} = u^k \partial_k \mathbf{r}, \quad \text{and} \quad (4) $$

$$ \frac{Ds}{Dt} = \frac{\partial s}{\partial t} + u^k \partial_k s, \quad (5) $$

where $s$ is some scalar field, $\partial_k$ is the partial derivative of a space-time field with respect to $\xi^k$, and we use Einstein’s summation convention (unless explicitly stated otherwise), with indices $k, l = 1, 2, 3$. Later, we will need to distinguish between horizontal ($k = 1, 2$) and vertical directions ($k = 3$), and use indices $i, j = 1, 2$ instead of $k, l$. Equation (4) shows that $(u^k)$ are the contravariant components of velocity $\dot{\mathbf{r}}$. Squaring (4) yields

$$ \dot{\mathbf{r}} \cdot \dot{\mathbf{r}} = G_{kl} u^k u^l \quad \text{with} \quad G_{kl} = \partial_k \mathbf{r} \cdot \partial_l \mathbf{r} \quad (k, l = 1, 2, 3). $$

The metric tensor associated with coordinates $\xi^k$ is $G_{kl}$. At this point no orthogonality is assumed. For a vector field $w^k \partial_k \mathbf{r}$, the Lagrangian derivative is
$$ \frac{D}{Dt}(w^k \partial_k \mathbf{r}) = \left[ \left( \frac{\partial}{\partial t} + u^m D_m \right) w^k \right] \partial_k \mathbf{r}, $$

where the covariant derivative $D_m$ is defined via the Christoffel symbol $\Gamma^k_{ml}$:

$$ D_m w^k = \partial_m w^k + \Gamma^k_{ml} w^l, $$

$$ 2G_{kl} \Gamma^k_{ml} = \partial_m G_{kl} + \partial_l G_{km} - \partial_k G_{lm}. $$

Noting $J = \sqrt{\det G_{kl}}$ the Jacobian such that $d^3\mathbf{r} = J d^3\xi$, the divergence operator is

$$ \operatorname{div}(u^k \partial_k \mathbf{r}) = \frac{1}{J} \partial_k (J u^k). $$

Hence, the budget for the pseudodensity $\hat{\rho} = J\rho$, where $\rho$ is the mass per unit volume, is

$$ \frac{\partial \hat{\rho}}{\partial t} + \partial_k (\hat{\rho} u^k) = 0, \quad \frac{D \hat{\rho}}{Dt} + \hat{\rho} \partial_k u^k = 0. \qquad (6) $$

Finally, we will note $R^k$ and $R_k$, the contravariant and covariant components of $\mathbf{R}$, respectively, as

$$ \Omega \times \mathbf{r} = R^k \partial_k \mathbf{r}, \quad R_k = (\Omega \times \mathbf{r}) \cdot \partial_k \mathbf{r} = G_{kl} R^l. $$

### c. Euler–Lagrange equations in a general curvilinear coordinate system

With the above definitions the Lagrangian can be rewritten as

$$ \mathcal{L} = \int \hat{L}(\hat{\rho}, \xi^k, u^k, s) dm, \quad dm = \hat{\rho} d^3 \xi, \qquad (7) $$

$$ \hat{L} = K + C - \Phi(\xi^k) - e \left[ \frac{J(\xi^k)}{\hat{\rho}}, s \right], \qquad (8) $$

$$ K = \frac{1}{2} G_{kl} u^k u^l, \quad \text{and} \qquad (9) $$

$$ C = G_{kl} u^k R^l = u^k R_k. \qquad (10) $$

In what follows, we will need to distinguish between $\partial_k \hat{L}$ and $\partial \hat{L}/\partial \xi^k$. The latter retains only the explicit dependence of $\hat{L}$ on $\xi^k$, and not its indirect dependence through the fields $u^k$, $\hat{\rho}$, and $s$. For instance, $\partial K/\partial \xi^k = \partial_k G_{lm} u^l u^m / 2$, and we have the chain rule:

$$ \partial_k \hat{L} = \frac{\partial \hat{L}}{\partial \xi^k} \partial_k \xi^l + \frac{\partial \hat{L}}{\partial u^l} \partial_k u^l + \frac{\partial \hat{L}}{\partial \hat{\rho}} \partial_k \hat{\rho} + \frac{\partial \hat{L}}{\partial s} \partial_k s. \qquad (11) $$

The action $\int \mathcal{L} d\tau$ is now considered as a functional of the label-time field $\xi^k(\mathbf{a}, \tau)$. By requiring the stationarity of the action $\int \mathcal{L} dt = 0$, we obtain

$$ \int_{t_0}^{t_1} \left( \int_{V} \frac{\partial \hat{L}}{\partial \hat{u}^k} \cdot \delta \hat{u}^k + \frac{\partial \hat{L}}{\partial \xi^k} \cdot \delta \xi^k + \frac{\partial \hat{L}}{\partial \hat{\rho}} \cdot \delta \hat{\rho} \right) dm \, d\tau = 0, $$

where $\delta s = 0$ due to Lagrangian conservation of specific entropy $s$. The integral involving $\delta \hat{\rho}$ can be expressed as

$$ \int_V \frac{\partial \hat{L}}{\partial \hat{\rho}} \delta \hat{\rho} \, dm = \int_A \frac{1}{\hat{\rho}} \partial_k \left( \frac{\partial \hat{L}}{\partial \hat{\rho}} \hat{\rho}^2 \right) \delta \xi^k \, dm, \qquad (12) $$

by using $\delta \hat{\rho}/\hat{\rho}^2 dm = -(\partial_k \delta \xi^k) d^3 \xi$ and integrating by parts with respect to $\xi^k$ (see Morrison 1998). In (12), we have omitted boundary terms that vanish with appropriate boundary conditions (see Morrison 1998). Using (12), expressing $\delta u^k$ as $\delta u^k = (\partial/\partial\tau) \delta \xi^k$ and integrating by parts with respect to $\tau$ yields

$$ \int_{t_0}^{t_1} \int_V \left[ -\frac{\partial}{\partial \tau} \frac{\partial \hat{L}}{\partial u^k} + \frac{\partial \hat{L}}{\partial \xi^k} + \frac{1}{\hat{\rho}} \partial_k \left( \frac{\partial \hat{L}}{\partial \hat{\rho}} \hat{\rho}^2 \right) \right] \delta \xi^k \, dm \, d\tau = 0. $$

Requiring that $\int \mathcal{L} dt = 0$ for arbitrary variations $\delta \xi^k$ yields the Euler–Lagrange equations of motion:

$$ \frac{D}{Dt} \frac{\partial \hat{L}}{\partial u^k} - \frac{\partial \hat{L}}{\partial \xi^k} = \frac{1}{\hat{\rho}} \partial_k \left( \frac{\partial \hat{L}}{\partial \hat{\rho}} \hat{\rho}^2 \right). \qquad (13) $$

### d. Interpretation of Euler–Lagrange equations

To decipher (13), we first note that the terms $K$ and $C$ produce the covariant components of acceleration $(D/Dt)\dot{\mathbf{r}}$ and the Coriolis force curl $\mathbf{R} \times \dot{\mathbf{r}}$, respectively:

$$ \left( \frac{D}{Dt} \frac{\partial}{\partial u^k} - \frac{\partial}{\partial \xi^k} \right) K = G_{kl} \left( \frac{\partial}{\partial t} + u^m D_m \right) u^l, $$

$$ \left( \frac{D}{Dt} \frac{\partial}{\partial u^k} - \frac{\partial}{\partial \xi^k} \right) C = (\partial_m R_k - \partial_k R_m) u^m. $$

Moreover, $\partial \Phi/\partial \xi^k$ are the covariant components of $\nabla \Phi$ and $\hat{\rho}^2 \partial \hat{L}/\partial \hat{\rho} = -pJ$.

$$ \frac{\partial}{\partial \xi^k} e(J/\hat{\rho}, s) = -\frac{p}{\hat{\rho}} \partial_k J, $$

which partially cancel with $\partial_k (Jp)/\hat{\rho}$ to leave only the covariant components of $-\nabla p$ on the right-hand side:
$$ G_{kl} \left( \frac{\partial}{\partial t} + u^m D_m \right) u^l + (\partial_m R_k - \partial_k R_m) u^m \\ = -\partial_k \Phi - \frac{1}{\rho} \partial_k p. \qquad (14) $$

Therefore, (13) is, as expected, nothing other than the covariant components of the Euler equation

$$ \frac{D}{Dt} \dot{\mathbf{r}} + \mathrm{curl} \, \mathbf{R} \times \dot{\mathbf{r}} = -\nabla \Phi - \frac{1}{\rho} \nabla p. \qquad (15) $$

### e. Vector-invariant form

Expanding $D/Dt = \partial_t + u^l \partial_l$ in (13), we obtain

$$ \frac{\partial v_k}{\partial t} + u^l \partial_l v_k - \frac{\partial \hat{L}}{\partial \xi_k} = \partial_k \left( \frac{\partial \hat{L}}{\partial \hat{\rho}} \hat{\rho} \right) + \frac{\partial \hat{L}}{\partial \hat{\rho}} \partial_k \hat{\rho}, \qquad (16) $$

where $v_k = \partial \hat{L}/\partial u^k$. Introducing the Bernoulli function,

$$ \hat{B} = v_l u^l + s \frac{\partial \hat{L}}{\partial s} - \frac{\partial \hat{L}}{\partial \hat{\rho}} \hat{\rho} - \hat{L}, \qquad (17) $$

and using the chain rule (11) in (16), the vector-invariant form of (13) is finally obtained:

$$ \frac{\partial v_k}{\partial t} + u^l (\partial_l v_k - \partial_k v_l) + \partial_k \hat{B} - s \partial_k \left( \frac{\partial \hat{L}}{\partial s} \right) = 0. \qquad (18) $$

Notice that $v_k$ are the covariant components of absolute velocity $\mathbf{R} + \dot{\mathbf{r}}$, and $\partial \hat{L}/\partial s = -T$. The thermodynamic contribution to the Bernoulli function (17) is Gibbs' free energy: $s \partial \hat{L}/\partial s - \hat{\rho} \partial \hat{L}/\partial \hat{\rho} + e = e + \alpha p - Ts$. For an ideal perfect gas this simplifies if one uses potential entropy $\theta(s)$ instead of $s$ as a prognostic variable. Indeed, (18) and (17) become

$$ \frac{\partial v_k}{\partial t} + u^l (\partial_l v_k - \partial_k v_l) + \partial_k \hat{B} - \theta \partial_k \left( \frac{\partial \hat{L}}{\partial \theta} \right) = 0, $$

where now

$$ \hat{B} = v_l u^l + \theta \frac{\partial \hat{L}}{\partial \theta} - \frac{\partial \hat{L}}{\partial \hat{\rho}} \hat{\rho} - \hat{L}. $$

The thermodynamic contribution to $\hat{B}$ is then $e + \alpha p - \theta \pi$, which vanishes in the particular case of an ideal perfect gas, $e + \alpha p = c_p T = \theta \pi$, where $c_p$ is the specific heat at constant pressure and $\pi = \partial e/\partial \theta$ is the Exner function. Then, $\hat{B} = v_l u^l - K - C + \Phi = (1/2) G_{kl} u^k u^l + \Phi$. One recovers therefore the well-known vector-invariant form of (15):

$$ \partial_t \dot{\mathbf{r}} + \mathrm{curl}(\mathbf{R} + \dot{\mathbf{r}}) \times \dot{\mathbf{r}} + \nabla \left( \frac{\dot{\mathbf{r}}^2}{2} + \Phi \right) + \theta \nabla \pi = 0. $$

This form is often derived from the advective form (15) by algebraic manipulations and using $\alpha dp = \theta d\pi$. However, it is important to stress that $\alpha dp = \theta d\pi$ holds for an ideal perfect gas only and that, for a general equation of state, thermodynamics will contribute to the Bernoulli function. There is then no obvious advantage to using potential temperature instead of entropy. The vector-invariant form can also be obtained directly and naturally from Hamilton's principle of least action using the Lie derivative formulation of Holm et al. (2002).

### f. Conservation laws

We now briefly state the conservation properties of (13). Due to Noether's theorem and the invariance of the Lagrangian $\hat{L}$ with respect to time, conservation of energy is expected. Due to the restricted form (2) that we consider for $\hat{L}$, the action is invariant under parcel relabeling, which implies the conservation of Ertel's potential vorticity (Newcomb 1967; Salmon 1988; Müller 1995; Padhye and Morrison 1996). If the geopotential is zonally symmetric, conservation of AAM should hold. We now provide explicit derivations of these expected results.

The absolute and potential vorticities are defined as

$$ J \omega^k = \epsilon^{klm} \partial_l v_k, \quad \eta = \frac{J \omega^k \partial_k s}{\hat{\rho}}, $$

where $\epsilon^{klm}$ is totally antisymmetric. Using the vector-invariant form (18), an expression for $\partial_t(J\omega^k)$ is derived. Combining this expression with the evolution equation for $\partial_k s$ obtained by differentiating (4) and with the mass budget (6) yields the Lagrangian conservation of potential vorticity:

$$ D\eta/Dt = 0. $$

These algebraic manipulations are strictly identical to the Cartesian case (Vallis 2006).

The local energy budget is

$$ \frac{\partial \hat{E}}{\partial t} + \frac{\partial}{\partial \xi^k} [(\hat{E} + Jp)u^k] = 0, \qquad (19) $$

where $\hat{E} = \hat{\rho} E$,

$$ E = u^k v_k - \hat{L} = K + \Phi(\xi^k) + e \left( \frac{J}{\hat{\rho}}, s \right). $$
Indeed, using the chain rule,

$$ \frac{DE}{Dt} = u^k \frac{Dv_k}{Dt} + v_k \frac{Du^k}{Dt} - \frac{\partial \hat{L}}{\partial \xi^k} u^k - \frac{\partial \hat{L}}{\partial u^k} \frac{Du^k}{Dt} - \frac{\partial \hat{L}}{\partial \hat{\rho}} \frac{D\hat{\rho}}{Dt}, $$

and using the Euler–Lagrange equations of motion, (19) follows.

To derive the AAM budget, we multiply (13) by $\hat{\rho}$ to obtain

$$ \frac{\partial}{\partial t}(\hat{\rho}v_k) + \partial_l(\hat{\rho}u^l v_k) + \partial_k(Jp) = \hat{\rho} \frac{\partial \hat{L}}{\partial \xi^k}. \quad (20) $$

Therefore, apart from boundary terms, $\int v_k dm$ is conserved provided the source term on the r.h.s of (20) vanishes. If the coordinate system is zonally symmetric, that is, $\xi^1$ is longitude, $\partial_1 G_{kl} = 0$, and $\partial_1 R_k = 0$, the source term for $\int v_1 dm$ reduces to $\hat{\rho} \partial_1 \Phi$. Hence, $\int v_1 dm$ is conserved if the geopotential is zonally symmetric.

Not much seems to have been achieved at this point, since (13) only restates the well-known compressible Euler equations, together with their conservation properties. However, we are now in a position to make approximations without jeopardizing the conservation properties. Indeed, the vector-invariant form (18) and the conservation laws for energy, potential vorticity, and AAM depend only on the equations of motion taking the Euler–Lagrange form (13), but not on the details of the Lagrangian $\hat{L}$. We can therefore approximate $\hat{L}$ as we wish. In particular, the metric tensor $G_{kl}$, the covariant components of planetary velocity $R_k = G_{kl}R^l$, and the Jacobian $J$ can be approximated, and these approximations can be made independently.

Useful and accurate approximations will take place in a coordinate system adapted to the dominance of the gravitational force in geophysical flows. Such a coordinate system, where the geopotential depends only on $\xi^3$, is constructed in the next section.

### 3. Geopotential-based curvilinear coordinates

With a nonspherical geopotential one must distinguish between the radial direction parallel to $\mathbf{r}$ and the vertical direction along $\nabla\Phi$. Similarly, one distinguishes between the tangential directions, orthogonal to $\mathbf{r}$, and the horizontal directions, orthogonal to $\nabla\Phi$. In this section we examine the construction of curvilinear coordinates $\xi^1, \xi^2$, and $\xi^3$, where

* the geopotential depends only on $\xi^3$ and, therefore, the third direction is vertical and

* furthermore the directions $\xi^1, \xi^2$ are horizontal; hence, $G_{13} = G_{23} = 0$.

We now need to distinguish between the horizontal ($k=1,2$) and vertical directions ($k=3$), and we use indices $i, j=1,2$ instead of $k, l$. The problem boils down to finding a mapping $(\xi^1, \xi^2, \xi^3) \mapsto \mathbf{r}$ such that $\partial_3 \mathbf{r} \parallel \nabla\Phi$ and $\partial_i \mathbf{r} \cdot \partial_j \mathbf{r} = 0$. We first show how a construction can be found in principle with a general geopotential $\Phi(\mathbf{r})$. Then, an approximate but explicit construction is sketched, and implemented for a specific, zonally symmetric geopotential, taking into account the leading aspherical corrections of the earth's geopotential.

#### a. General geopotential field

Let $\bar{\Phi}(\xi^3)$ be the desired dependence of $\Phi$ on $\xi^3$, and $\Phi_{\text{ref}} = \bar{\Phi}(\xi_{\text{ref}}^3)$ be a reference geopotential. It is generally possible, although not necessarily simple, in practice to find a system of curvilinear coordinates $\xi^1, \xi^2$ on the geoid $\Phi = \Phi_{\text{ref}}$, that is, a mapping $(\xi^1, \xi^2) \mapsto \mathbf{r}_{\text{ref}}$ such that $\Phi[\mathbf{r}_{\text{ref}}(\xi^1, \xi^2)] = \Phi_{\text{ref}}$. Notice that such a coordinate system must have singularities, like the pole for standard latitude–longitude coordinates. Rigorously, one must consider several such curvilinear systems and patch them together to cover the whole sphere/spheroid. This procedure is unambiguous provided one manipulates only expressions that transform properly under a change of curvilinear coordinates. This is what we do in sections 4 and 5. In fact, although we do not do it here, it is possible to adopt an intrinsic formulation of all that follows by replacing the coordinates $\xi^1$ and $\xi^2$ by a vector $\mathbf{n}$ belonging to the unit sphere. Now, let us follow the vertical curve passing through $\mathbf{r}_{\text{ref}}(\xi^1, \xi^2)$; that is, we integrate

$$ \partial_3 \mathbf{r} = \frac{\nabla \Phi}{\|\nabla \Phi\|^2} \frac{d\bar{\Phi}}{d\xi^3}, \quad \mathbf{r}(\xi^1, \xi^2, \xi_{\text{ref}}^3) = \mathbf{r}_{\text{ref}}(\xi^1, \xi^2). \quad (21) $$

Then,

$$ \frac{\partial}{\partial \xi^3} \Phi[\mathbf{r}(\xi^k)] = \frac{d\bar{\Phi}}{d\xi^3} \Rightarrow \Phi[\mathbf{r}(\xi^k)] = \bar{\Phi}(\xi^3), $$

which implies $\partial_i \mathbf{r} \cdot \partial_3 \mathbf{r} = 0$.

Notice that even if $(\xi^1, \xi^2)$ is orthogonal on the reference surface, nothing can be said when $\xi^3 \neq \xi_{\text{ref}}^3$. In the sequel we do not assume $G_{12} = 0$, although it is possible to obtain $G_{12} = 0$ when the geopotential is zonally symmetric.

Furthermore, (21) does not guarantee that the mapping $(\xi^k) \mapsto \mathbf{r}(\xi^k)$ is invertible. Since it clearly is for a spherical geopotential, and the actual geopotential is close to spherical, we assume that a breakdown does not occur in the spatial domain of interest.

#### b. A perturbative approach for nearly spherical geopotential

In section 3a, a general shape of the geopotential was considered. However, since the earth and more generally
telluric planets are quite well described by a sphere, it can be sufficient and more explicit to construct **r**(**ξ**<sup>k</sup>) by a perturbative procedure starting from a spherical geometry.

Let us recall that geopotential Φ(**r**) is defined as the sum of the gravitational and centrifugal potentials:

$$ \Phi(\mathbf{r}) = V(\mathbf{r}) - \frac{1}{2} \|\mathbf{\Omega} \times \mathbf{r}\|^2. \quad (22) $$

Assuming r = ||**r**|| is of order O(a) with a suitably defined planetary radius, V is of order g₀a, where ||∇V|| = O(g₀) at r = O(a). The nondimensional parameter

$$ \gamma = \frac{a\Omega^2}{g_0} \quad (23) $$

is typically small (γ ~ 1/300 for the earth). Since γ measures the relative strength of centrifugal and gravitational accelerations, it also defines the order of magnitude of the planetary ellipticity and, therefore, the deviation of V(**r**) from spherical symmetry. Therefore, one can decompose Φ/(ag₀) as

$$ \frac{\Phi}{ag_0} = \Phi_0(r/a) + \gamma\Phi_1(\mathbf{r}/a) + \dots, \quad (24) $$

where Φ₀(r/a) = (r/a)⁻¹ and Φ₁ collects the nonspherical part of the gravitational potential and the centrifugal potential. We can now explicitly construct a corresponding expansion of **r**(**ξ**<sup>k</sup>) in powers of γ,

$$ \mathbf{r}(\xi^k) = R(\xi^3)\mathbf{r}_0(\xi^i) + \gamma\mathbf{r}_1(\xi^k) + \dots, \quad (25) $$

satisfying, order by order,

$$ \Phi[\mathbf{r}(\xi^k)] = \bar{\Phi}(\xi^3), \quad \partial_3 \mathbf{r} \cdot \partial_i \mathbf{r} = 0. $$

The leading order is satisfied if **r**₀ defines curvilinear coordinates on the unit sphere and ag₀Φ₀[R(**ξ**<sup>3</sup>)] = Φ̄(**ξ**<sup>3</sup>); that is, R(**ξ**<sup>3</sup>) = g₀a²/Φ(**ξ**<sup>3</sup>). At order γ,

$$ \frac{\Phi}{ag_0} = \frac{a}{R} - \gamma \mathbf{r}_0 \cdot \mathbf{r}_1 \frac{a}{R^2} + \gamma \Phi_1(R \mathbf{r}_0), $$

hence, the condition Φ[**r**(**ξ**<sup>k</sup>)] = Φ̄(**ξ**<sup>3</sup>) determines the radial part of the correction **r**₁,

$$ \mathbf{r}_0 \cdot \mathbf{r}_1 = \frac{R^2}{a} \Phi_1(R \mathbf{r}_0), \quad (26) $$

while a tangential correction is required to maintain orthogonality ∂₃**r** · ∂ᵢ**r** = 0:

$$ R \partial_3 \mathbf{r}_1 \cdot \partial_i \mathbf{r}_0 = - \frac{dR}{d\xi^3} \mathbf{r}_0 \cdot \partial_i \mathbf{r}_1 = \frac{dR}{d\xi^3} [\mathbf{r}_1 \cdot \partial_i \mathbf{r}_0 - \partial_i (\mathbf{r}_0 \cdot \mathbf{r}_1)]. $$

Differentiating (26),

$$ \mathbf{r}_0 \cdot \partial_3 \mathbf{r}_1 = \partial_3 \left[ \frac{R^2}{a} \Phi_1(R \mathbf{r}_0) \right], $$

so that all of the covariant components of ∂₃**r**₁ in the basis (∂₁**r**₀, ∂₂**r**₀, **r**₀) are known as a function of **r**₁ and **ξ**<sup>3</sup>. Therefore, at fixed **ξ**<sup>1</sup>, **ξ**<sup>2</sup>, we face a simple ordinary differential equation for **r**₁(**ξ**<sup>3</sup>). If Φ₁ is given as a sum of spherical harmonics, each with a power-law dependence on **r**, an explicit solution can be found. We provide an example in the next subsection.

### c. A simple set of nearly spherical coordinates

We now apply the procedure outlined in section 3b to the dominant terms considered by White et al. (2008):

$$ \frac{\Phi}{g_0 a} = \frac{a}{r} + \gamma \left[ \alpha_1 \left(\frac{a}{r}\right)^3 \left( \sin^2 \chi - \frac{1}{3} \right) + \alpha_2 \left(\frac{r}{a}\right)^2 \cos^2 \chi \right], \quad (27) $$

where α₁, α₂ are O(1) constants and χ is the geocentric latitude such that

$$ \mathbf{r} = r(\cos\lambda \cos\chi \mathbf{e}_x + \sin\lambda \cos\chi \mathbf{e}_y + \sin\chi \mathbf{e}_z). $$

This geopotential is zonally symmetric. To define **r**₀ and express **r**₁, we use longitude–latitude coordinates; that is, **ξ**<sup>1</sup> = λ, **ξ**<sup>2</sup> = φ, **ξ**<sup>3</sup> = ξ and

$$ \begin{aligned} \mathbf{r}_0 &= \mathbf{e}_R = \cos\lambda \cos\phi \mathbf{e}_x + \sin\lambda \cos\phi \mathbf{e}_y + \sin\phi \mathbf{e}_z, \\ \partial_\phi \mathbf{r}_0 &= \mathbf{e}_\phi = -\cos\lambda \sin\phi \mathbf{e}_x - \sin\lambda \sin\phi \mathbf{e}_y + \cos\phi \mathbf{e}_z, \\ \mathbf{r}_1 &= a_\phi(\phi, \xi)\mathbf{e}_\phi + a_R(\phi, \xi)\mathbf{e}_R. \end{aligned} \quad (28) $$

Now, neglecting terms O(γ²) and using expansion (25),

$$ \begin{aligned} \frac{\partial \mathbf{r}}{\partial \xi} &= \frac{dR}{d\xi} \mathbf{r}_0 + \gamma \left( \frac{\partial a_\phi}{\partial \xi} \mathbf{e}_\phi + \frac{\partial a_R}{\partial \xi} \mathbf{e}_R \right), \\ \frac{\partial \mathbf{r}}{\partial \phi} &= R \mathbf{e}_\phi + \gamma \left( \frac{\partial a_\phi}{\partial \phi} \mathbf{e}_\phi + \frac{\partial a_R}{\partial \phi} \mathbf{e}_R - a_\phi \mathbf{e}_R + a_R \mathbf{e}_\phi \right), \\ \frac{\partial \mathbf{r}}{\partial \xi} \cdot \frac{\partial \mathbf{r}}{\partial \phi} &= \gamma \frac{dR}{d\xi} \left( \frac{\partial a_R}{\partial \phi} - a_\phi \right) + \gamma R \frac{\partial a_\phi}{\partial \xi}, \\ &= \gamma R^2 \left[ \frac{\partial}{\partial \xi} \left( \frac{a_\phi}{R} \right) - \frac{\partial}{\partial \phi} \left( R^{-2} \frac{dR}{d\xi} a_R \right) \right], \end{aligned} \quad (29) $$
and a way to satisfy $(\partial\mathbf{r}/\partial\xi) \cdot (\partial\mathbf{r}/\partial\phi) = 0$ is to introduce the nondimensional potential $\psi(\phi, R)$ such that

$$a_{\phi} = R \frac{\partial \psi}{\partial \phi}, \quad a_{R} = R^{2} \frac{\partial \psi}{\partial R}. \qquad (30)$$

Finally, $\psi$ is determined by the condition that $\Phi[\mathbf{r}(\lambda, \phi, \xi)] = \bar{\Phi}(\xi)$. At leading order this implies $\bar{\Phi}(\xi) = a g_0 \Phi_0[R(\xi)]$, while at order $\gamma$ we obtain

$$\frac{\partial \psi}{\partial R} = -R^{-2} a_R = -R^{-2} \left( \frac{d\Phi_0}{dR} \right)^{-1} \Phi_1, \\ a \frac{\partial \psi}{\partial R} = \alpha_1 \left(\frac{a}{R}\right)^3 \left(\sin^2\phi - \frac{1}{3}\right) + \alpha_2 \left(\frac{R}{a}\right)^2 \cos^2\phi, \\ \psi = -\frac{\alpha_1}{2} \left(\frac{a}{R}\right)^2 \left(\sin^2\phi - \frac{1}{3}\right) + \frac{\alpha_2}{3} \left(\frac{R}{a}\right)^3 \cos^2\phi. \quad (31)$$

Notice that the coordinate system defined by (25), (28), (30), and (31) is horizontally orthogonal; that is, $(\partial\mathbf{r}/\partial\phi) \cdot (\partial\mathbf{r}/\partial\lambda) = 0$. As noted above, this seems to be allowed only by a zonally symmetric geopotential.

## 4. Approximations

In (13) no approximation has been made to the fully compressible Euler equations. However, if we use a geopotential-based coordinate system as defined and constructed in section 3 (i.e., $\partial_i \Phi = 0$ and $G_{i3} = 0$), the kinetic energy and Euler-Lagrange equations in (14) become

$$K = \frac{1}{2} G_{ij} u^i u^j + \frac{1}{2} G_{33} u^3 u^3, \\ \left( \frac{D}{Dt} \frac{\partial}{\partial u^i} - \frac{\partial}{\partial \xi^i} \right) K + (\partial_m R_i - \partial_i R_m) u^m = - \frac{1}{\rho} \partial_i p, \quad (32)$$

$$\left(\frac{D}{Dt} \frac{\partial}{\partial u^3} - \frac{\partial}{\partial \xi^3}\right) K + (\partial_m R_3 - \partial_3 R_m) u^m = -\partial_3 \Phi - \frac{1}{\rho} \partial_3 p, \quad (33)$$

where we remind that $i, j = 1, 2$ while $m = 1, 2, 3$. Notice that, for the sake of completeness, we keep $R_3 \neq 0$. However, $R_3 = 0$ as soon as the geopotential is zonally symmetric, which seems a good enough approximation for the vast majority of applications. Equations of motion (32) and (33) are written in terms of relative kinetic energy $K$, Coriolis force, and the pressure gradient. In what follows, we emphasize for each kind of approximation how they approximate each of these three terms.

The main improvement is that gravity, $\partial_3 \Phi$, enters only the third equation of motion. This will simplify the derivation of the usual approximations, and allow the

derivation of new ones, while preserving dynamical consistency. We first show how the introduction of a hydrostatic switch $\delta_{NH}$ into the exact Lagrangian yields quasi-hydrostatic equations in a general, nonaxisymmetric geopotential. Turning then to the shallow-atmosphere approximation, we recover and generalize previously obtained equation sets (White and Wood 2012; Tort and Dubos 2013).

### a. Quasi-hydrostatic approximation

A defining feature of quasi-hydrostatic systems (White and Bromley 1995; White and Wood 2012) is that the vertical balance loses its prognostic character and becomes a diagnostic equation. Equation (13) shows that this will be the case if $\partial L/\partial u^3 = 0$. In the Lagrangian density (7), only $K$ and $C$ depend on $u^3$. We therefore introduce a hydrostatic switch $\delta_{NH}$ and redefine the $K$ and $C$ as

$$K = \frac{1}{2} G_{ij} u^i u^j + \frac{\delta_{NH}}{2} G_{33} u^3 u^3, \quad \text{and} \\ C = u^i R_j + \delta_{NH} u^3 R_3.$$

Setting $\delta_{NH} = 1$ gives the full equation set while setting $\delta_{NH} = 0$ modifies the vertical momentum balance. From the energetic point of view, the total energy is now

$$E = \frac{1}{2} G_{ij} u^i u^j + \frac{\delta_{NH}}{2} G_{33} u^3 u^3 + \Phi + e.$$

Hence, by setting $\delta_{NH} = 0$, the vertical kinetic energy is neglected from the energy budget, which is a feature of hydrostatic systems (Holm et al. 2002). From a physical point of view, neglecting vertical kinetic energy is equivalent to setting the inertia of the fluid to zero for vertical motion, which imposes that vertical forces balance each other. As is well known, the ratio of vertical to horizontal velocity scales like the ratio of the vertical to horizontal characteristic scales of the flow, so that those terms should be retained to model small-scale flows.

To compare with White and Wood (2012), we now obtain the evolution equations for the physical components of velocity. For this the coordinates $\xi^k$ need to be orthogonal, that is, $G_{12} = 0$, which seems to require a zonally symmetric geopotential. We therefore assume for the remainder of this section that the geopotential is zonally symmetric; hence, $R_3 = 0$. Then, it makes sense to define the metric factors $h_k = \sqrt{G^{kk}}$. The physical components of velocity are then $u_k = h_k u^k$ (the reader will note the absence of summation in the expressions of $h_k$ and $u_k$ and also that the notation $u_k$ does not refer to the covariant components of relative velocity) and
$$ 
\begin{aligned}
\frac{D}{Dt} \frac{\partial K}{\partial u^i} - \frac{\partial K}{\partial \xi^i} &= h_i \frac{Du_i}{Dt} + u^i (u^i h_i \partial_j h_j - u^j h_i \partial_j h_j) \\
&\quad + u^i u^j h_i \partial_j h_j - \delta_{NH} u^3 u^3 h_3 \partial_i h_3, \\
\frac{D}{Dt} \frac{\partial K}{\partial u^3} - \frac{\partial K}{\partial \xi^3} &= \delta_{NH} h_3 \left( \frac{Du_3}{Dt} + u^3 u^j \partial_j h_3 \right) - u^i u^j h_j \partial_i h_3.
\end{aligned}
$$

Assuming now zonally symmetric coordinates (i. e., $R_1 = \Omega h_1^2$, $R_2 = R_3 = 0$, $\partial_1 h_k = 0$), the Euler–Lagrange equations simplify to

$$ 
\begin{aligned}
h_1 \frac{Du_1}{Dt} &+ u^2 u^1 h_1 \partial_2 h_1 + u^1 u^3 h_1 \partial_3 h_1 \\
&+ 2\Omega h_1 (u^2 \partial_2 h_1 + u^3 \partial_3 h_1) = -\frac{1}{\rho} \partial_1 p, \\
h_2 \frac{Du_2}{Dt} &- u^1 u^1 h_1 \partial_2 h_1 + u^2 u^3 h_2 \partial_3 h_2 - \delta_{NH} u^3 u^3 h_3 \partial_2 h_3 \\
&- 2\Omega h_1 u^1 \partial_2 h_1 = -\frac{1}{\rho} \partial_2 p, \\
\delta_{NH} h_3 \left( \frac{Du_3}{Dt} + u^3 u^2 \partial_2 h_3 \right) &- u^i u^j h_j \partial_3 h_j - 2\Omega h_1 \partial_3 h_1 u^1 \\
&= -\frac{1}{\rho} \partial_3 p - \partial_3 \Phi.
\end{aligned}
\quad (34) $$

When $\delta_{NH} = 1$, (34) are precisely the nonhydrostatic (A.10)–(A.12) from White and Wood (2012), while when $\delta_{NH} = 0$, the quasi-hydrostatic (A.13)–(A.15) are recovered. Notice that we have also checked that equations from Gates (2004) are recovered with oblate spheroidal coordinates with $(\xi^1, \xi^2, \xi^3) = (\lambda, \phi, \xi)$:

$$ \mathbf{r} = c(\cosh\xi \cos\phi \cos\lambda \mathbf{e}_x + \cosh\xi \cos\phi \sin\lambda \mathbf{e}_y + \sinh\xi \sin\phi \mathbf{e}_z). \quad (35) $$

Compared to the derivation by White and Wood (2012), we arrive here straightforwardly at several nontrivial results:

* the necessity to neglect $u^3 u^2 \partial_2 h_3$ in the quasi-hydrostatic equations follows from the neglect of vertical kinetic energy in the Lagrangian, while White and Wood (2012) needed to utilize the energy budget to justify it;

* the expression of the Coriolis force in terms of $\partial_2 h_1$ and $\partial_3 h_1$ derives naturally from the expression of the covariant component of planetary velocity $R_1 = \Omega h_1^2$, while a geometric reasoning was used in White and Wood (2012); and

* the nontraditional Coriolis term exists because $\partial_3 R_1 \neq 0$; an approximate system neglecting the vertical variations of $R_1$ necessarily makes the traditional approximation.

### b. Soundproof approximations

Generally speaking, acoustic waves are suppressed if the feedback of pressure on density is suppressed. This can be achieved by constraining the value of $\rho$. The simplest soundproof approximation is the Boussinesq approximation, whereby $\rho = \rho_r = \text{cst}$, but density modifications $\delta\rho$ due to entropy $s$ are taken into account only in the potential energy [i.e., $\Phi(\delta\rho)$]. In this section we omit for brevity the kinetic (K), Coriolis (C), and geopotential (Φ) terms of the Lagrangian as they are left untouched. One should bring these terms back into the Lagrangian in order to obtain the complete equations of motion. While the Boussinesq approximation is adequate for oceanic applications, it is important for atmospheric applications to allow for large variations of $\rho$. This can be achieved with $\rho \approx \rho^*(s, \xi^k) = \rho[s, p^*(\xi^k)]$, where $p^*(\xi^k)$ is a background pressure, often taken to be hydrostatically balanced. Within a variational principle, such a constraint is enforced by augmenting the Lagrangian through the introduction of a Lagrangian multiplier $\lambda$. The corresponding augmented Lagrangian is

$$ \hat{L}(\hat{\rho}, s, \lambda, \xi^k) = -e\left(\frac{J}{\hat{\rho}}, s\right) + \lambda \left[ \frac{J}{\hat{\rho}} - \frac{1}{\rho^*(s, \xi^k)} \right], \quad (36) $$

$$ \hat{\rho}^2 \frac{\partial \hat{L}}{\partial \hat{\rho}} = -J(p^* + \lambda), \quad \text{and} \quad (37) $$

$$ \frac{\partial \hat{L}}{\partial \xi^k} = \frac{1}{\hat{\rho}} \left[ (p^* + \lambda) \partial_k J + \lambda \rho^{*-2} \frac{\partial \rho^*}{\partial \xi^k} \right]. \quad (38) $$

In (36), $\lambda$ enforces the condition that the expression that it multiplies vanishes; that is, $\rho = \hat{\rho}/J = \rho^*$. The specific form chosen here gives $\lambda$ the dimension of a pressure. Inserting the above into (13), one obtains the adiabatic equations of motion. It turns out that they coincide with those derived in Cartesian coordinates by Klein and Pauluis (2012) by letting $p = p^* + \lambda$ with $\lambda \ll p^*$ and expanding up to first order in $\lambda$:

$$ \rho^{-1} \partial_k p \approx \left( \rho^{*-1} - \lambda \rho^{*-2} \frac{\partial \rho}{\partial p} \right) \partial_k p^* + \rho^{*-1} \partial_k \lambda. $$

Hence, the Lagrange multiplier has the physical interpretation of a deviation of total pressure from $p^*$. The variational derivation of the equations obtained by Klein and Pauluis (2012) directly shows that they conserve potential vorticity. Conservation of energy holds if the background state $p^*$ is stationary and conservation of angular momentum holds for a zonally symmetric background state.
Equation (36) simplifies for an ideal perfect gas because $e/\theta = \kappa\pi$ and $\alpha/\theta = \pi/p$ depend only on pressure. Using $\theta$ as a prognostic variable and taking into account the constraint $p = p^*, \pi = \pi^*$ in $e = \kappa\theta\pi$ yields

$$ \hat{L}(\hat{\rho}, \theta, \lambda, \xi^k) = -\kappa\theta\pi^* + \lambda\theta \left( \frac{J}{\theta\hat{\rho}} - \frac{\pi^*}{p^*} \right). \quad (39) $$

Variations of density with potential temperature are neglected (i.e., $\rho^* \approx p^*\theta/\pi^*$) if

$$ \hat{L}(\hat{\rho}, \theta, \lambda, \xi^k) = -\kappa\theta\pi^* + \lambda \left[ \frac{J}{\hat{\rho}} - \frac{1}{\rho^*(\xi^k)} \right]. \quad (40) $$

Cotter and Holm (2013) have shown that the Lagrangian in (39) generates the pseudoincompressible equations where $\lambda' = \lambda\theta$ is their Lagrangian multiplier and (40) generates the Lipps-Hemler anelastic equations, respectively. The more general Lagrangian in (36) is, to the best of our knowledge, new.

### c. Shallow-atmosphere approximation

To analyze the shallow-atmosphere approximation, we first need to define quantitatively the shallowness of the atmosphere. For this, we let $\Delta\Phi$ be the order of magnitude of the geopotential difference between the top and bottom of the atmosphere, both supposed to be close to a geopotential surface $\Phi = \Phi_{\text{ref}} = \Phi(\xi^3 = 0)$. Since $\nabla\Phi = O(g_0)$ where the reference gravity $g_0$ has been introduced in section 3, an order of magnitude of the atmospheric thickness is $H = \Delta\Phi/g_0$ and a measure of its shallowness is

$$ \varepsilon = \frac{H}{a} = \frac{\Delta\Phi}{g_0 a}. \quad (41) $$

Then, $\mathbf{r}(\xi^k)$ can be expanded in powers of $\varepsilon$:

$$ \mathbf{r}(\xi^k) = \mathbf{r}(\xi^i, \xi_{\text{ref}}^3) + \xi^3 \partial_3 \mathbf{r}(\xi^i, \xi_{\text{ref}}^3) + O(a\varepsilon^2), \quad (42) $$

where the first two terms are $O(a)$ and $O(H = \varepsilon a)$. Using expansion (42) to approximate the metric tensor $G_{kl}$ implies at leading order that the vertical dependence of $G_{kl}$ is neglected:

$$ \begin{aligned} G_{kl}(\xi^i, \xi^3) &\approx G_{kl}^{\text{ref}}(\xi^i), \\ K &\approx \frac{1}{2} G_{ij}^{\text{ref}} u^i u^j + \delta_{\text{NH}} \left( \frac{1}{2} G_{33}^{\text{ref}} u^3 u^3 \right), \end{aligned} \quad (43) $$

where $G_{ij}^{\text{ref}} = G_{ij}(\xi^1, \xi^2, \xi^3 = 0)$.

Similarly, $\Phi(\xi^3)$ can be expanded as

$$ \Phi = \Phi_{\text{ref}} + \xi^3 \partial_3 \Phi_{\text{ref}} + O(\varepsilon^2 a g_0). \quad (44) $$

Notice that $\partial_3\Phi_{\text{ref}}$ is a constant (independent from $\xi^i$); hence, (44) is equivalent to using an affine function of $\Phi$ as a vertical coordinate. With (43) and (44), gravity $\partial_3\Phi/h_3 \approx \partial_3\Phi_{\text{ref}}/h_3^{\text{ref}}$ becomes independent from $\xi^3$ but can still depend on $\xi^i$.

Regarding the Coriolis term $C = u^i R_i + \delta_{\text{NH}} u^3 R_3$, it is tempting to simply evaluate it at $\xi^3 = 0$. However, this would neglect both the vertical dependence of the metric, which is justified by $\varepsilon \ll 1$, and the vertical dependence of the planetary velocity. As argued by Tort and Dubos (2013), the latter approximation requires more care and may not be justified if the planetary velocity is large compared to the fluid velocity, as measured by the smallness of the planetary Rossby number $\mu$:

$$ \mu = \frac{U}{a\Omega}, \quad (45) $$

where $U$ is a characteristic velocity scale. Indeed, $K = \dot{\mathbf{r}} \cdot \dot{\mathbf{r}} = O(U^2)$, while $C = \dot{\mathbf{r}} \cdot \mathbf{R} = O(U\Omega a)$. If $\mu \sim \varepsilon$, or more generally if $\mu \ll \varepsilon$ does not hold, an approximation to $C$ should retain terms of order $O(U^2)$. Explicitly,

$$ \begin{aligned} C = u^i R_j + \delta_{\text{NH}} u^3 R_3, \\ \approx u^i R_j^{\text{ref}} + \delta_{\text{NT}} \xi^3 u^i \partial_3 R_j^{\text{ref}} \\ + \delta_{\text{NH}} (u^3 R_3^{\text{ref}} + \delta_{\text{NT}} \xi^3 u^3 \partial_3 R_3^{\text{ref}}). \end{aligned} $$

where the switch $\delta_{\text{NT}}$ (nontraditional) is used to tag the terms that retain a dependence on $\xi^3$. At this point we have retained the terms proportional to $R_3$ for completeness. However, since they vanish for a zonally symmetric geopotential, they are likely small in the vast majority of applications. Therefore, for the sake of simplicity, we now drop them to finally retain only

$$ C = u^i (R_j^{\text{ref}} + \delta_{\text{NT}} \xi^3 \partial_3 R_j^{\text{ref}}). \quad (46) $$

The first term is $O(U\Omega a)$ and the second one is $O(U\Omega H)$ and can be safely neglected only if $\Omega H \ll U$ (i.e., $\varepsilon \ll \mu$). This condition is not met by typical oceanic flows and is marginally met by typical atmospheric flows (Tort and Dubos 2013).

To compare the equations resulting from (43), (44), and (46) with White and Wood (2012), we assume again that the coordinate system is orthogonal and zonally symmetric and we derive the evolution equations for the physical velocity components. For this, starting from (34), it suffices to let $\partial_3 h_k = 0$, except for $R_1$:
$$
\begin{align*}
\partial_2 R_1 &= \Omega \partial_2 [(1 + \delta_{\text{NT}} \xi^3 \partial_3) h_1^2], \\
&= 2\Omega \left( h_1 \partial_2 h_1 + \delta_{\text{NT}} \frac{\xi^3}{2} \partial_{23} h_1^2 \right), \\
\partial_3 R_1 &= \Omega \partial_3 [(1 + \delta_{\text{NT}} \xi^3 \partial_3) h_1^2] = 2\delta_{\text{NT}} \Omega h_1 \partial_3 h_1,
\end{align*}
$$

where it is implied that $h_1$ and $\partial_3 h_1$ are evaluated at $\xi^3 = 0$. Hence, (34) become

$$
\begin{equation}
\begin{split}
& h_1 \frac{Du_1}{Dt} + u^2 u^1 h_1 \partial_2 h_1 + 2 \Omega u^2 \left( h_1 \partial_2 h_1 + \delta_{\text{NT}} \frac{\xi^3}{2} \partial_{32} h_1^2 \right) \\
& \quad + \delta_{\text{NT}} 2 \Omega u^3 h_1 \partial_3 h_1 = -\frac{1}{\rho} \partial_1 p,
\end{split}
\tag{47}
\end{equation}
$$

$$
\begin{equation}
\begin{split}
& h_2 \frac{Du_2}{Dt} - u^1 u^1 h_1 \partial_2 h_1 - \delta_{\text{NH}} u^3 u^3 h_3 \partial_2 h_3 \\
& \quad - 2 \Omega u^1 \left( h_1 \partial_2 h_1 + \delta_{\text{NT}} \frac{\xi^3}{2} \partial_{32} h_1^2 \right) = -\frac{1}{\rho} \partial_2 p, \quad \text{and}
\end{split}
\tag{48}
\end{equation}
$$

$$
\begin{align}
& \delta_{\text{NH}} h_3 \left( \frac{Du_3}{Dt} + u^3 u^2 \partial_2 h_3 \right) - \delta_{\text{NT}} 2 \Omega u^1 h_1 \partial_3 h_1 \nonumber \\
& \quad = -\frac{1}{\rho} \partial_3 p - \partial_3 \Phi. \tag{49}
\end{align}
$$

If $\delta_{\text{NT}} = 0$, the traditional shallow-atmosphere (A16)-(A21) from White and Wood (2012) are recovered. On the other hand, if one makes the spherical-geoid approximation and uses as coordinates $(\xi^1, \xi^2, \xi^3) = (\lambda, \phi, z)$ with $z = (\Phi - \Phi_{\text{ref}})/g$, then $h_3 \approx 1$, $h_1 \approx (a+z) \cos\phi$ leading at $z=0$ to $h_1 = a \cos\phi$, $\partial_3 h_1 = \cos\phi$, $\partial_{32} h_1^2 = -4a \sin\phi \cos\phi$, so that

$$
h_1 \partial_2 h_1 + \delta_{\text{NT}} \frac{\xi^3}{2} \partial_{32} h_1^2 = -a^2 \sin\phi \cos\phi \left(1 + \delta_{\text{NT}} \frac{2z}{a}\right),
$$

and (14) from Tort and Dubos (2013) is recovered. As in Tort and Dubos (2013), the reintroduction of the non-traditional Coriolis terms (proportional to $\Omega u^3$ for $Du_1/Dt$ and $\Omega u^1$ for $Du_3/Dt$) must be accompanied by a correction of the traditional Coriolis terms (proportional to $\Omega u^2$ for $Du_1/Dt$ and $\Omega u^1$ for $Du_2/Dt$) in order to retain all conservation laws.

d. Nonspherical geopotential corrections

Although the spherical geometry is relevant to describe
our planet, under specific circumstances and for some
other planets, the flattening at the poles, which induces
a latitudinal variation for g, may have significant effects
and should be taken into account to estimate those effects.
The nonspherical corrections from a spherical model at

the leading order described in section 3b allow us to
consider nearly spherical geometry and geopotential.
In this section, we assume an axisymmetric geopotential
but one slightly flattened at the poles. Flattening is
characterized by the small parameter $\gamma$ and the set of
nearly spherical coordinates defined in section 3c is used.
In addition to (29), we have

$$
\partial_{\lambda} \mathbf{r} = [R \cos\phi + \gamma(-a_{\phi} \sin\phi + a_{R} \cos\phi)] \mathbf{e}_{\lambda}.
$$

The coordinate system is horizontally orthogonal ($G_{12} = 0$) and we neglect coefficients $G_{i3}$ because they are $O(\gamma^2)$. Hence, the metric tensor is diagonal with $G_{ii} = h_i^2$, such as

$$
\begin{align}
h_{\lambda}^{2} &= R \cos\phi (R \cos\phi + 2\gamma H), \notag \\
h_{\phi}^{2} &= R(R + 2\gamma G), \notag \\
h_{\xi}^{2} &= d_{\xi} R (d_{\xi} R + 2\gamma \partial_{\xi} a_{R}), \tag{50}
\end{align}
$$

where $H(\phi, \xi) = -a_{\phi}\sin\phi + a_{R}\cos\phi$ and $G(\phi, \xi) = a_{R} + \partial_{\phi}a_{\phi}$. To obtain the nonhydrostatic deep equations of motion with nonspherical corrections at leading order, the corrections to $K$ and $C$ at order $O(\gamma)$ are being retained:

$$
K = \frac{1}{2} (h_{\lambda}^{2} \dot{\lambda}^{2} + h_{\phi}^{2} \dot{\phi}^{2} + \delta_{\text{NH}} h_{\xi}^{2} \dot{\xi}^{2}), \quad (51)
$$

$$
C = \Omega h_{\lambda}^{2} \dot{\lambda}, \qquad (52)
$$

where the expressions of ($h_{\lambda}^2$, $h_{\phi}^2$, $h_{\xi}^2$) are given in (50). In terms of the scaling used in section 4c, $K$ and $C$ are asymptotically correct up to $O(\mu\gamma)$ and $O(\gamma)$, respectively. Therefore, if $\mu \ll 1$, it is asymptotically consistent to retain only the nonspherical corrections to $C$ and neglect those to $K$. To usefully retain corrections $O(\mu\gamma)$ to $K$, the expansion of $C$ should be more accurate and the expansion sketched in section 3b should be pursued to order $O(\gamma^2)$. It is of course possible to retain the corrections $O(\mu\gamma)$ to $K$ and expand $C$ only to $O(\gamma)$, but the resulting model would not be more accurate than the simpler model where $K$ is not corrected for nonspherical effects.

Notice that if the atmosphere is shallow ($\epsilon \ll 1$), re-
taining the full dependence of $R$ and $\Phi$ on $\xi^3$ in the
Lagrangian amounts to retaining terms of all orders in $\epsilon^n$
in the expansion of the Lagrangian in powers of $\epsilon$. For
large-scale atmospheric motion on the earth, the orders
of magnitude of the dimensionless parameters are typ-
ically equal to

$$
\gamma \sim 3.4 \times 10^{-3}, \quad \mu \sim 2.1 \times 10^{-2}, \quad \epsilon \sim 1.6 \times 10^{-3}.
$$
Because of the fast planetary rotation of Saturn (index s) and Jupiter (index j), the flattening is quite important; in fact, it is larger than ($\mu, \varepsilon$):

$$\gamma_s \sim 1.8 \times 10^{-1}, \quad \mu_s \sim 3.1 \times 10^{-2}, \quad \varepsilon_s \sim 6.7 \times 10^{-4},$$

and

$$\gamma_j \sim 9.6 \times 10^{-2}, \quad \mu_j \sim 4.0 \times 10^{-3}, \quad \varepsilon_j \sim 2.8 \times 10^{-4}.$$

In the above regimes it would be consistent to neglect the $O(\mu\varepsilon)$ and $O(\mu\gamma)$ corrections to K while retaining the $O(\varepsilon)$ and $O(\gamma)$ corrections to C. The leading-order corrections to K + C are therefore those of C. Compared to a spherical-geopotential, shallow-atmosphere model, the corrections $O(\varepsilon)$ are those of Tort and Dubos (2013) and restore a complete Coriolis force, while corrections $O(\gamma)$ modify slightly the traditional Coriolis term.

## 5. Discussion: Full asymptotic consistency

Invoking Hamilton's principle of least action from an approximated Lagrangian will systematically lead to a dynamically consistent equation set with the appropriate conserved quantities, no matter how carefully the Lagrangian is approximated. Hence, dynamical consistency is not dependent on asymptotic consistency, which is understood as the condition that the smallest terms retained are larger than the largest terms neglected. Depending on the dynamical regime considered, it is still desirable to check the asymptotic consistency of the model in order not to overestimate its accuracy.

Tort and Dubos (2013) have derived equations resulting from a consistent asymptotic development of kinetic energy $K + C$ in the approximated Lagrangian. However, no asymptotic expansion was performed for geopotential $\Phi$ nor internal energy $e$, and only the leading-order term was retained for them. Until such an expansion has been completed, it is not known whether the model obtained by Tort and Dubos (2013) is actually more accurate at order $O(\varepsilon)$ than would be a traditional shallow-atmosphere model. We now investigate this point as an illustrative example of the issue of asymptotic consistency. Hence, we expand the potential and internal energy terms to include a $O(\varepsilon)$ correction, and we analyze whether the additional terms appearing in the equations of motion can be consistently neglected compared to all other retained terms. The dependence of internal energy on $\varepsilon$ comes from the ratio $r^2\cos\phi$ between density and pseudodensity.

$$\begin{align*} \alpha &= \frac{r^2 \cos\phi}{\hat{\rho}} \approx \frac{a^2 \cos\phi}{\hat{\rho}} + \frac{2az \cos\phi}{\hat{\rho}} = \alpha_s + \alpha', \\ e(\alpha, s) &\approx e(\alpha_s, s) - p_s \alpha', \\ p(\alpha, s) &\approx p(\alpha_s, s) - \left(\frac{c}{\alpha_s}\right)^2 \alpha' = p_s + p', \end{align*}$$

where $\alpha' \ll \alpha_s$ and $p' \ll p_s$ are $O(\varepsilon)$ corrections to the shallow-atmosphere expressions $\alpha_s$ and $p_s$, respectively, and we have used $(\partial p/\partial \alpha)|_{\alpha_s} = -(c/\alpha_s)^2$ with c the sound speed. Regarding potential energy,

$$\Phi(z) = a^2 g_0 \left( \frac{1}{a} - \frac{1}{a+z} \right) \approx g_0 z \left( 1 - \frac{z}{a} \right), \quad (53)$$

which results in a $O(\varepsilon)$ correction to $g = d\Phi/dz$ in the vertical momentum balance.

To proceed and compare these corrections to the other terms retained in the equations of motion, we need to make an assumption on the order of magnitude of the pressure terms. Assuming a nearly geostrophic regime, the horizontal pressure gradient is of the same order of magnitude as the traditional Coriolis term. Since the latter has been expanded to the next order in $\varepsilon$, the $O(\varepsilon)$ corrections to the pressure gradient should be retained also. Having retained these corrections in the Lagrangian, they will appear in the vertical balance as $O(\varepsilon)$ corrections to $\partial_z p_s$, whose order of magnitude is $\rho g_0$. Hence, $O(\varepsilon)$ corrections to $\Phi(z)$ should be included in the Lagrangian, with the effect of taking into account small vertical variations of gravity. Finally, the dynamically and asymptotically consistent density Lagrangian will be

$$\begin{aligned} \hat{L} = \frac{1}{2} a [a \cos^2 \phi \dot{\lambda} + a \dot{\phi} + 2 \cos^2 \phi \Omega (a + 2z)] \\ - g_0 z \left(1 - \frac{z}{a}\right) - e \left[ \frac{a \cos^2 \phi (a + 2z)}{\hat{\rho}}, s \right]. \end{aligned} \quad (54)$$

The first term in (54) corresponds to the nontraditional shallow-atmosphere kinetic energy of Tort and Dubos (2013) for which only the vertical dependence of the planetary part is retained. The second term is the potential energy where the vertical variation of gravity acceleration is retained at order $O(\varepsilon)$. The last term is the internal energy and takes into account the slightly conical shape of an atmospheric columns at order $O(\varepsilon)$.

A consistent asymptotic development is then obtained at leading order in an expansion in $\varepsilon$, retaining vertical variations of planetary velocity, of the Jacobian and of the gravity acceleration. Note that in a different dynamical regime (e.g., near the equator where geostrophic balance breaks down) it may be asymptotically consistent
<table><caption>TABLE 1. Estimations of large-scale parameters $\mu$, $\epsilon$, and $\gamma$ ($10^{-2}$).</caption><thead><tr><th rowspan="2">Giant planets</th><th colspan="4">Solar</th><th colspan="2">Extrasolar</th></tr><tr><th>Jupiter</th><th>Saturn</th><th>Uranus</th><th>Neptune</th><th>HD 209458b</th><th>HD 189733b</th></tr></thead><tbody><tr><th>μ</th><td>0.4</td><td>3</td><td>11</td><td>12</td><td>19</td><td>15</td></tr><tr><th>ε</th><td>0.03</td><td>0.07</td><td>0.13</td><td>0.12</td><td>2</td><td>0.24</td></tr><tr><th>γ</th><td>9.5</td><td>18</td><td>2.9</td><td>2.2</td><td>0.55</td><td>0.42</td></tr></tbody></table>

to neglect the above corrections to the internal and potential energy. Furthermore,

* using the vertical coordinate Z introduced by Dellar (2011), such as $r^2dr = a^2dZ$ and $Z|_{r=a} = 0$, will give the exact pressure gradient without any approximation, because $J = a^2 \cos\phi$, while on the other hand, the geopotential and planetary velocity will take a non-trivial form as a function of the vertical coordinate Z, and

* using the geopotential as vertical coordinate $\xi^3 = \Phi$ will give the exact geopotential term in the vertical balance, that is, $d_{\Phi}\Phi = 1$, but will give a non-trivial form of the planetary velocity and pressure gradient.

The above discussion may be generalized to all other approximations. If one wants to add a nonzero vertical acceleration, one may check its order of magnitude by introducing a horizontal scale $L \ll a$ and comparing with all the other terms in the equations. Taking into account nonspherical corrections in the geopotential as in the previous section, the Jacobian J has to be developed at $O(\gamma)$ order for asymptotic consistency.

White et al. (2005) pointed out that taking into account a latitudinal variation of gravity acceleration g, within the spherical geopotential approximation, will produce spurious sources of potential vorticity and then will lead to a dynamically inconsistent model. Our analysis confirms this: latitudinal variations of g will arise only if at least $O(\gamma)$ corrections to the Jacobian J are included in the expressions for internal and potential energy (as done above with $O(\epsilon)$ corrections).

## 6. Conclusions

In this paper, we have described a general variational framework that allows a systematic derivation of equations of motion, for a large panel of approximations. The derivation highlights the essence of usual geophysical approximations and provides dynamically consistent systems in the sense that all physical properties of conservation are ensured.

We first considered a general class of Euler–Lagrange equations from which a wide range of dynamically

consistent equations of motion can be obtained without doing any variational calculus. We then identified new Lagrangians corresponding to existing equations of motion, originally derived by manipulating and approximating the exact equations of motion rather than the Lagrangian (Klein and Pauluis 2012; White and Wood 2012). We also extended Tort and Dubos (2013) by

* considering a zonally symmetric (not spherical) geopotential in a general orthogonal coordinate system and

* expanding the geopotential and internal energy at next order in atmospheric shallowness $\epsilon$ to achieve asymptotic consistency in geostrophically balanced flow.

The last extension underlines the difference between the dynamical and asymptotic consistency. All approximations in (32) and (33) can be made independently as soon as expressions for $K$, $R_i$, $\Phi$, and $\partial_i p/\rho$ are kept identical in the equations. This sometimes leads to rather exotic but still dynamically consistent models. As an example, keeping the exact metric terms in K and neglecting vertical variations in C, a dynamically consistent deep-atmosphere model with incomplete Coriolis force is obtained. But it will not be asymptotically consistent because some of the terms that are retained are smaller than some terms that are neglected. The asymptotic consistency typically depends on the dynamical regime that is considered. Close to a geostrophic regime, Tort and Dubos's (2013) equations are not asymptotically consistent. To be consistent, next-order vertical variations of $J = a \cos\phi(a + 2z)$ and $\Phi = g_0z(1 - z/a)$ have also to be retained. The equations derived by Tort and Dubos (2013) should nevertheless correctly capture the full Coriolis force in far-from-geostrophic situations (e.g., near the equator).

Finally, we provided a method of obtaining explicit metric terms corresponding to a nearly spherical geopotential. It should be quite easy, then, to include nonspherical corrections in an existing general circulation model. As for nontraditional regimes, to achieve asymptotic consistency close to a geostrophic regime, corrections at first order in $\gamma$ should be retained:
• in the Coriolis term, $\Omega R \cos\phi(R \cos\phi + 2\gamma H) + O(\gamma^2)$;

• in the Jacobian $J = h_{\lambda}h_{\phi}h_{\xi}$ from (50), where $J = R^2 \cos\phi d_{\xi}R + \gamma R(\partial_{\xi}a_{R}R \cos\phi + G \cos\phi d_{\xi}R + Hd_{\xi}R) + O(\gamma^2)$; and

• in the geopotential (27), where $\Phi = (a^2g_0/r) + \gamma(\dots) + O(\gamma^2)$.

In Table 1, $\mu$, $\epsilon$, and $\gamma$ have been estimated for a few giant planets using large-scale parameters from Cho and Polvani (1996), Cho et al. (2003), and Showman and Polvani (2011). Taking into account nonspherical geopotential corrections could be relevant to model rapidly rotating giant gas planets for which the equatorial bulge is more significant than that of the earth. A few exoplanets have already been modeled (Cho et al. 2003; Showman and Polvani 2011; Mayne et al. 2014). For the exoplanets HD202458b and HD189733b, the aspect ratio $\epsilon$ and the flattening $\gamma$ are of order $O(\mu^2)$ or even smaller. Therefore, the contribution of the centrifugal force should be taken into account if one wishes to include the nontraditional and/or nonspherical effects.

Note that, shortly before submitting this manuscript, the authors became aware of independent work by Andrew Staniforth, sharing a number of goals and results, recently submitted to the *Quarterly Journal of the Royal Meteorological Society* (Staniforth 2014a,b).

*Acknowledgments.* The authors are grateful to the anonymous reviewers for careful examination of the manuscript and constructive remarks, which helped improve its clarity and the discussion about the implications of the results.

## REFERENCES

Bénard, P., 2014: An oblate-spheroid geopotential approximation for global meteorology. Quart. J. Roy. Meteor. Soc., **140**, 170–184, doi:10.1002/qj.2141.

Cho, J. Y. K., and L. M. Polvani, 1996: The morphogenesis of bands and zonal winds in the atmospheres on the giant outer planets. Science, **273**, 335–337, doi:10.1126/science.273.5273.335.

—, K. Menou, B. M. S. Hansen, and S. Seager, 2003: The changing face of the extrasolar giant planet HD 209458b. Astrophys. J. Lett., **587**, L117, doi:10.1086/375016.

Cotter, C. J., and D. Holm, 2013: Variational formulations of sound-proof models. Quart. J. Roy. Meteor. Soc., doi:10.1002/qj.2260, in press.

Dellar, P. J., 2011: Variations on a beta-plane: Derivation of non-traditional beta-plane equations from Hamilton’s principle on a sphere. J. Fluid Mech., **674**, 174–195, doi:10.1017/S0022112010006464.

Durran, D. R., 1989: Improving the anelastic approximation. J. Atmos. Sci., **46**, 1453–1461, doi:10.1175/1520-0469(1989)046<1453: ITAA>2.0.CO;2.

—, 2008: A physically motivated approach for filtering acoustic waves from the equations governing compressible stratified flow. J. Fluid Mech., **601**, 365–379, doi:10.1017/S0022112008000608.

Gates, W. L., 2004: Derivation of the equations of atmospheric motion in oblate spheroidal coordinates. J. Atmos. Sci., **61**, 2478–2487, doi:10.1175/1520-0469(2004)061<2478: DOTEOA>2.0.CO;2.

Gerkema, T., and V. I. Shira, 2005: Near-inertial waves in the ocean: Beyond the traditional approximation. J. Fluid Mech., **529**, 195–219, doi:10.1017/S0022112005003411.

—, J. T. F. Zimmerman, L. R. M. Maas, and H. van Haren, 2008: Geophysical and astrophysical fluid dynamics beyond the traditional approximation. Rev. Geophys., **46**, 1–33, doi:10.1029/2006RG000220.

Holm, D. D., J. E. Marsden, and T. S. Ratiu, 2002: The Euler-Poincaré equations in geophysical fluid dynamics. Large-Scale Atmosphere-Ocean Dynamics, J. Norbury and I. Roulstone, Eds., Cambridge University Press, 251–300.

Hua, B. L., D. W. Moore, and S. Le Gentil, 1997: Inertial nonlinear equilibration of equatorial flows. J. Fluid Mech., **331**, 345–371, doi:10.1017/S0022112096004016.

Klein, R., and O. Pauluis, 2012: Thermodynamic consistency of a pseudo-incompressible approximation for general equations of state. J. Atmos. Sci., **69**, 961–968, doi:10.1175/JAS-D-11-0110.1.

Lipps, F. B., and R. S. Hemler, 1982: A scale analysis of deep moist convection and some related numerical calculations. J. Atmos. Sci., **39**, 2192–2210, doi:10.1175/1520-0469(1982)039<2192: ASAODM>2.0.CO;2.

Mayne, N. J., and Coauthors, 2014: The unified model, a fully compressible, non-hydrostatic, deep atmosphere global circulation model, applied to hot jupiters. Astron. Astrophys., **561**, A1, doi:10.1051/0004-6361/201322174.

Morrison, P. J., 1998: Hamiltonian description of the ideal fluid. Rev. Mod. Phys., **70**, 467–521, doi:10.1103/RevModPhys.70.467.

Müller, P., 1995: Ertel’s potential vorticity theorem in physical oceanography. Rev. Geophys., **33**, 67–87, doi:10.1029/94RG03215.

Müller, R., 1989: A note on the relation between the traditional approximation and the metric of the primitive equations. Tellus, **41A**, 175–178, doi:10.1111/j.1600-0870.1989.tb00374.x.

Newcomb, W. A., 1967: Exchange invariance in fluid systems in magneto-fluid and plasma dynamics. Proc. Symp. Appl. Math., **18**, 152–161.

Ogura, Y., and N. A. Phillips, 1962: Scale analysis of deep and shallow convection in the atmosphere. J. Atmos. Sci., **19**, 173–179, doi:10.1175/1520-0469(1962)019<0173:SAODAS>2.0.CO;2.

Padhye, N., and P. J. Morrison, 1996: Fluid element relabeling symmetry. Phys. Lett., **219A**, 287–292, doi:10.1016/0375-9601(96)00472-0.

Phillips, N. A., 1966: The equations of motion for a shallow rotating atmosphere and the traditional approximation. J. Atmos. Sci., **23**, 626–628, doi:10.1175/1520-0469(1966)023<0626: TEOMFA>2.0.CO;2.

Roulstone, I., and S. J. Brice, 1995: On the Hamiltonian formulation of the quasi-hydrostatic equations. Quart. J. Roy. Meteor. Soc., **121**, 927–936, doi:10.1002/qj.49712152410.

Salmon, R., 1988: Hamiltonian fluid dynamics. Annu. Rev. Fluid Mech., **20**, 225–256, doi:10.1146/annurev.fl.20.010188.001301.

Showman, A. P., and L. M. Polvani, 2011: Equatorial superrotation on tidally locked exoplanets. Astrophys. J., **738**, 71, doi:10.1088/0004-637X/738/1/71.
Staniforth, A., 2014a: Spheroidal and spherical geopotential approximations. *Quart. J. Roy. Meteor. Soc.*, doi:10.1002/qj.2324, in press.

—, 2014b: Deriving consistent approximate models of the global atmosphere using Hamilton’s principle. *Quart. J. Roy. Meteor. Soc.*, doi:10.1002/qj.2273, in press.

Tort, M., and T. Dubos, 2013: Dynamically consistent shallow-atmosphere equations with a complete Coriolis force. *Quart. J. Roy. Meteor. Soc.*, doi:10.1002/qj.2274, in press.

Vallis, G. K., 2006: Atmospheric and Oceanic Fluid Dynamics: Fundamentals and Large-scale Circulation. Cambridge University Press, 745 pp.

White, A. A., and R. A. Bromley, 1995: Dynamically consistent, quasi-hydrostatic equations for global models with a complete representation of the Coriolis force. *Quart. J. Roy. Meteor. Soc.*, 121, 399–418, doi:10.1002/qj.49712152208.

—, and G. W. Inverarity, 2012: A quasi-spheroidal system for modelling global atmospheres: Geodetic coordinates. *Quart. J. Roy. Meteor. Soc.*, 138, 27–33, doi:10.1002/qj.885.

—, and N. Wood, 2012: Consistent approximate models of the global atmosphere in non-spherical geopotential coordinates. *Quart. J. Roy. Meteor. Soc.*, 138, 980–988, doi:10.1002/qj.972.

—, B. J. Hoskins, I. Roulstone, and A. Staniforth, 2005: Consistent approximate models of the global atmosphere: Shallow, deep, hydrostatic, quasi-hydrostatic and non-hydrostatic. *Quart. J. Roy. Meteor. Soc.*, 131, 2081–2107, doi:10.1256/qj.04.49.

—, A. Staniforth, and N. Wood, 2008: Spheroidal coordinate systems for modelling global atmospheres. *Quart. J. Roy. Meteor. Soc.*, 134, 261–270, doi:10.1002/qj.208.