On the hydrostatic approximation in rotating stratified flow

Hydrostatic models were and still are the workhorses for realistic simulations of ocean dynamics, especially for climate applications. Introducing a Fourier space projection method and using the Heisenberg–Gabor limit, a formalism is developed to systematically evaluate the role of flatness, stratification, rotation and friction for the fidelity of the hydrostatic approximation. The hydrostatic approximation is formally first order in $\mathit{\gamma }=H/L$, where H is the vertical and L the horizontal scale of the phenomenon considered. For linear (low-amplitude) and unforced stratified rotating flow, the dynamics can be separated into balanced flow and wave motion. It is shown that for the linear balanced motion the hydrostatic approximation is exact and for wave motion it is second order, obtaining the leading prefactors. The fidelity of the hydrostatic approximation therefore also relies on the ratio of the amplitude of wave motion to balanced motion. This ratio adds considerably to the quality of the hydrostatic approximation for larger-scale flows in the atmosphere and the ocean.

Imposing the divergenceless condition is a linear projection of the dynamical variables into the subspace of divergenceless vector fields, for both the Navier–Stokes and the hydrostatic formalism. Both projections are local in Fourier space. The former is well known, while the latter, developed here, asks for an extension of the dynamical space to four dimensions. The projection is followed by a time-evolution operator, which differs in the wave frequencies only. Combining the projection and the linear evolution operators in both formalisms leads to the linear projection-evolution operator.

Calculating the difference of the two projection-evolution operators, the expression of the error, scaling and prefactors done by the hydrostatic approximation is obtained. Analyzing the eigenspace of the projector-evolution operators, it is shown that for rotating buoyant vortical flow, the hydrostatic approximation is of third order for buoyant forcing, second order for horizontal and first order for vertical dynamical forcing. Balanced dynamics is in the kernel of the linear projection-evolution operator, and conservation of potential vorticity is expressed by the kernel of its adjoint.

Using the Heisenberg–Gabor limit, it is shown that for large-scale ocean dynamics, the difference of the dynamics of the projection-evolution operator between the two formalisms is insignificant. It is shown that the hydrostatic approximation is appropriate for realistic ocean simulations with vertical viscosities larger than $\approx {\mathrm{10}}^{-\mathrm{2}}$ m² s⁻¹. A special emphasis is on unveiling the physical interpretation of the calculations.

1 Introduction

At large scales, ocean velocities point predominantly in the horizontal direction, due to the flatness of the ocean basin, the action of gravity, the rotation of the earth and the feeble variations of the density of ocean waters. These properties allow for the derivation of idealized mathematical models of ocean dynamics from the fully three-dimensional Navier–Stokes equations (see, e.g., Vallis, 2017). Idealized models are not only implemented numerically to efficiently integrate atmosphere and ocean dynamics, but are at the basis of human understanding (Wirth, 2010) and therefore of the development of storylines explaining physical processes and their interactions. Most of these models are based on the hydrostatic approximation. In the quasi-geostrophic model, for example, the entire dynamics is projected on the direction of the balanced eigenvector (to be defined below), and the influence of dynamics along the other vectors due to their linear and nonlinear interaction is estimated from the balanced flow. At small scales, the turbulent motion is fully three-dimensional, and the three-dimensional Navier–Stokes equations have to be solved. As numerical simulations, especially when coastal processes are considered, move towards finer horizontal resolution, it is no longer clear to what extent the hydrostatic approximation is valid. This has sparked an increased interest in quasi-hydrostatic or fully three-dimensional Navier–Stokes models for ocean modeling (see, e.g., Marshall et al., 1997a; Guillaume et al., 2017; Auclair et al., 2018; Popinet, 2020; Calandrini et al., 2024).

To accelerate the iteratively solved elliptic problem of finding the non-hydrostatic pressure, it is often preconditioned by the more easily obtained hydrostatic pressure. At the same time, hydrostatic models are still widely used and developed for larger-scale simulations of the atmosphere (see, e.g., Janjic et al., 2001; Wan et al., 2013; Milewski and Tabak, 2015; Snodin and Wood, 2022; Spensberger et al., 2022; Bouvier et al., 2024), the ocean (see, e.g., Kärnä et al., 2018; Kevlahan and Dubos, 2019) and the climate (see, e.g., Boucher et al., 2020). Hydrostatic models are also employed to increase physical insight (Atoufi et al., 2023; Winn et al., 2023). Comparisons of the two approaches using numerical simulations are reported (Tseng et al., 2005; Zeman et al., 2021) as well as the coupling of the two model types (Blayo and Rousseau, 2016; Qu et al., 2019). For slow and large-scale dynamics, the hydrostatic approximation is valid; for fast small-scale dynamics, it is not. The present paper discusses the gray zone between the two limits and establishes where it is situated using a purely analytical approach.

The justification of the hydrostatic approximation relies on the flatness of the dynamics – that is, the smallness of the ratio of the vertical scale to the horizontal scale, $\mathit{\gamma }=H/L$. For the basin-wide circulation $\mathit{\gamma }\approx {\mathrm{10}}^{-\mathrm{4}}$ for a mesoscale eddy $\mathit{\gamma }\approx {\mathrm{10}}^{-\mathrm{2}}$ and their dynamics is well described by the hydrostatic models. They acceptably model gravity currents, coastal currents and internal inertial waves of large wavelength, with $\mathit{\gamma }\approx \mathrm{0.01}-\mathrm{0.1}$, but are inappropriate to integrate shorter internal waves $\mathit{\gamma }\approx \mathrm{0.1}-\mathrm{1}$ or the mixed-layer dynamics and three-dimensional mixing, where γ≈1. Scaling arguments on the validity of the hydrostatic approximation are often based on inertial gravity wave dynamics. In Marshall et al. (1997b) it was estimated, comparing the horizontal advection time to the wave frequency, that the motion is hydrostatic if ${\mathit{\gamma }}^{\mathrm{2}}/Ri\ll \mathrm{1}$, where $Ri=(NH/U{)}^{\mathrm{2}}$ is the Richardson number, N the buoyancy frequency and U a typical horizontal velocity scale. These scaling arguments do not include the Coriolis parameter f, which is known to further suppress vertical accelerations; this can be seen in the prominent Taylor–Proudman–Poincaré theorem (see, e.g., Wirth and Barnier, 2008). Furthermore, the effect of resonances are not included. Note also that Gill (1982) states that the hydrostatic approximation is appropriate if the frequency of the wave motion is smaller than the buoyancy frequency ω≪N. For inertial-gravity waves, $N>\mathit{\omega }>f$, and for balanced motion, $\mathit{\omega }\ll N,f$. In the ocean and atmosphere we typically have N>10f, but there are also areas where the ocean and atmosphere are weakly stratified or unstratified, or where the stratification is unstable. The action of gravity waves on the vortical flow is, however, supposed to be small, rendering such arguments insignificant for most oceanic phenomena (see, e.g., Vanneste, 2013; Wirth, 2013). A counter example is, however, trapped lee waves generated by topography, of which the dynamics differs qualitatively, when the hydrostatic approximation is employed, as discussed by Yu and Teixeira (2015). In the present work I develop a framework that allows for discussing the fidelity of the hydrostatic approximation as a function of the flatness γ, stratification N, rotation f and friction ν. The influence of these parameters revealed by the calculations given in the text unveil the physics of the hydrostatic approximation. The calculations are performed in Fourier space, as projections into the space of incompressible vector fields are local, but results are by no means restricted to spectral models.

The action of pressure, which ensures the vanishing divergence of the velocity field, is discussed in Sect. 2. The present work combines three techniques. First (T1), I use the fact that imposing the divergenceless condition is a linear and local projection in Fourier space, for both the Navier–Stokes and the hydrostatic formalism; it is discussed in Sect. 3. Second (T2), I extend the dynamical space to four dimensions considering the dynamical-vertical and buoyant acceleration as separate variables; it is discussed in Sect. 4. Third (T3) I use the linearized equations of the dynamics in both formalisms and compare their eigenvectors and eigenvalues, to bring out the similarities and differences, it is discussed in Sect. 5. This leads to a formal framework to investigate the validity of the hydrostatic approximation using linear algebra in low (≤5) dimensional spaces.

To evaluate the hydrostatic approximation, I consider the response of the linearized Navier–Stokes and the hydrostatic scheme to a forcing. Formally, the linear terms of the equations can be written on the left-hand side of the equality sign and the other terms on the right-hand side; the latter then appear, formally, as forcings to the linear system. The origin of this generalized forcing can be external or internal, by the boundary or by nonlinear terms and other processes, which are explicitly resolved or parameterized. Two main steps are necessary to solve the two linearity schemes (see Fig. 3). The first step (S1) considers the effect of a forcing on the eigenvectors of the corresponding linearity projection-evolution operator. This step is explained in Sect. 6. Section 7 details the different components of the forcing vector. The second step (S2) considers the response of the system to the forcing in the eigenspaces of the linear projection-evolution operator, which is given in Sect. 7.1 (S2a) for stationary forcing and in Sect. 7.2 (S2b) for time-dependent forcing. To evaluate the hydrostatic approximation, for each of the two steps, two sub-steps are necessary: one has to determine what is the error due to the hydrostatic approximation in each eigenspace, the absolute error, and what is the amplitude of the motion in each eigenspace, to obtain the relative error. The sub-steps are given in Sect. 7.1 for a stationary forcing and in Sect. 7.3 for a periodic forcing. In Sect. 7.4, I briefly discuss the inclusion of the nonlinear terms in the formalism. The discussion in Sect. 8 compares the error due to the hydrostatic approximation to other causes of error in numerical models such as the increased friction drag coefficient and the Doppler shift due to unresolved advection. This section also discusses whether the finite time observations of waves allow us to discriminate between the two formalisms, based on the Heisenberg–Gabor uncertainty principle. The conclusions, Sect. 9, give a summary of the scaling due to different forcings and compare the deficiencies of the hydrostatic equations to unresolved processes, internal variability, finite resolution in space and uncertainties in the observations.

1.1 The scales of ocean dynamics

Before evaluating the hydrostatic formalism by comparing its dynamics to the Navier–Stokes model, it is important to introduce the scales in space and time of the different regimes of ocean dynamics, as shown in Fig. 1. For small scales and short times, the dynamics is three-dimensional and the hydrostatic approximation deficient. We therefore consider the dynamics at scales larger than 10 m and slower than 1 h. Four main regions can be identified. The first is the viscous region, below the blue line associated with the eddy viscosity employed in the model. The actual kinematic viscosity of sea water (O(10⁻⁶ m² s⁻¹)) is many orders of magnitudes smaller, and the dynamics in the viscous region is misrepresented in both formalisms. It is important to note that the vertical dynamics with typical vertical velocities of $w={\mathrm{10}}^{-\mathrm{4}}$ m s⁻¹ are in the viscous regime. The vertical velocities are crucial for the fidelity of the hydrostatic approximation. The second regime is where balanced vortical flow dominates and the hydrostatic and the Navier–Stokes formalism are close, as confirmed by our findings below. This regime is operating at larger scales than the viscous regime and at scales faster than the Coriolis parameter. At timescales faster than the Coriolis parameter, the third regime representing inertia-gravity waves is dominant at larger scales and can be separated into barotropic and baroclinic modes. The former are much faster than the latter; both interact through topographic variations of the ocean floor. The former have negligible direct influence on the ocean circulation averaged over a few days. The latter can break in the ocean interior, leading to vertical mixing. The differences in the wave dynamics between the two formalisms are of second order in the aspect ratio O(γ²). The fourth regime is fully three-dimensional; it is called the mixing regime in Fig. 1. Due to its three-dimensional dynamics, the horizontal scale is constrained by the thickness of the ocean layer considered. In numerical simulations it is squeezed between the viscous and the wave regime, and it disappears for larger viscosities and hyper-viscosities. The validity of the hydrostatic approximation is questioned in this regime. It is, however, in this range of scales that a large number of specific processes of ocean dynamics occur. Examples are convective plumes (Marshall and Schott, 1999; Wirth and Barnier, 2006; Griffies and Treguier, 2013), down-slope gravity currents (Wirth, 2009; Manucharyan et al., 2014; Gačić et al., 2021), interactions with topography (Jayne et al., 2004), breaking of internal waves (Lamb, 2014), nonlinears (top, bottom and interfacial) (Laanaia et al., 2010; Elipot and Gille, 2009), Langmuir cells (McWilliams et al., 1997), ocean fronts (McWilliams, 2021), small-scale instabilities (Dewar et al., 2015) and others. Dedicated numerical simulations are necessary to validate the hydrostatic approximation for these processes. In the present work I do not want to rival these efforts, but rather present a theoretical framework to consider its validity.

Figure 1 Schematic view of the ocean dynamics in (logarithmic) timescales and space scales. The timescale is normalized by the Coriolis period $\mathrm{1}/f$, with $f={\mathrm{10}}^{-\mathrm{4}}$ s⁻¹. The buoyancy frequency is given by the red dashed line, a typical value is N=10f. Inertia-gravity waves have periods shorter than $\mathrm{1}/f$ and their group velocity cg decreases with the wave length. At horizontal scales smaller than the depth of a homogeneous layer and not subject to substantial influence of the Coriolis force (red area), the dynamics can lead to overturning and mixing for which the hydrostatic approximation fails. Slow large-scale motion is close to balance (green area). Thin red lines (slope = 1) present constant velocity. Balanced motion in the ocean is rarely faster than 1 m s⁻¹, Barotropic gravity (bt.) waves (yellow area) have velocities of the order of 100 m s⁻¹. Internal gravity (bc.) waves (orange) are slower than 10 m s⁻¹ and typical vertical velocities are less than 10⁻⁴ m s⁻¹. The blue lines (slope = $\mathrm{1}/\mathrm{2}$) represent the effect of viscous damping for different values of viscosity and also for bi-Laplace damping (blue dashed line, slope $=\mathrm{1}/\mathrm{4}$). The area below is dominated by viscous friction, not correctly represented in the hydrostatic or the Navier–Stokes formalism.

In the ocean, the dynamics also takes place at scales faster than 1 h and smaller than 10 m. For this dynamics, however, neither rotation nor stratification plays a dominant role; it is close to three-dimensional isotropic turbulence, and standard large eddy simulation schemes developed in other fields of turbulent fluid dynamics can be employed to parameterize their effect on the larger-scale dynamics. For even smaller-scale processes, for example wave breaking and other surface processes, dedicated characterizations have to be employed even in a far future.

2 Pressure

The starting points of our investigation are the three-dimensional, incompressible Navier–Stokes equations in a rotating frame, subject to the Boussinesq approximation and the traditional approximation (see, e.g., Wirth and Barnier, 2008), which supposes the rotation vector to be aligned with gravity:

$$\begin{array}{}\text{(1)}& {\displaystyle }{\partial }_t\mathit{u}=-(\mathit{u}\cdot \mathrm{\nabla })\mathit{u}-f\mathit{z}\times \mathit{u}-g\mathit{z}{\displaystyle \frac{\mathit{\rho }}{{\mathit{\rho }}_{\mathrm{0}}}}+\mathrm{\nabla }\left(\mathit{\nu }\mathrm{\nabla }\mathit{u}\right)-\mathrm{\nabla }P,\text{(2)}& {\displaystyle }\mathrm{\nabla }\cdot \mathit{u}=\mathrm{0},\end{array}$$

where u is the three-dimensional velocity vector, z is the unit vector in the vertical upward direction, g is gravity, ν viscosity and P is the pressure divided by the average density (see below). When the velocity vanishes, the above equations simplify to the hydrostatic equation:

$$\begin{array}{}\text{(3)}& {\partial }_{z}P=-g{\displaystyle \frac{\mathit{\rho }}{{\mathit{\rho }}_{\mathrm{0}}}}=\stackrel{\mathrm{̃}}{b},\end{array}$$

showing that the vertical pressure gradient equals the buoyancy.

The density can be separated into three parts:

$$\begin{array}{}\text{(4)}& \mathit{\rho }(x,z,t)={\mathit{\rho }}_{\mathrm{0}}\left(t\right)+{\mathit{\rho }}_z(z,t)+{\mathit{\rho }}^{\prime }(x,y,z,t),\end{array}$$

where the terms from left to right on the right-hand side of Eq. 4 represent (1) the average density in the domain that can vary in time, (2) a part of the density variation that is horizontally averaged and only varies in the vertical direction and in time, and (3) the deviations from the sum of the two. In the incompressible Navier–Stokes equations, only the last term has to be considered for the acceleration, as the other terms lead to a vertical force which is independent of the horizontal direction and is therefore countered by a horizontally independent hydrostatic pressure force. We can replace ρ by ρ^′ and write the buoyancy $b=-g{\mathit{\rho }}^{\prime }/{\mathit{\rho }}_{\mathrm{0}}$ in Eq. (1) and the pressure will change from P to P^′, where the latter does not include the hydrostatic pressure force due to ρ₀(t) and ρ_z(z,t). The evolution equation of the buoyancy anomaly is

$$\begin{array}{}\text{(5)}& {\partial }_{t}b=-\left(\right(\mathit{u}\cdot \mathrm{\nabla })b-u_z{N}^{\mathrm{2}}+\mathit{\kappa }{\mathrm{\nabla }}^{\mathrm{2}}b.\end{array}$$

If, considering the terms on the right-hand side of Eq. (1) from left to right, we see that the total acceleration is the sum of the advection of inertia a_u (first term on the right-hand side of Eq. 1), the Coriolis acceleration a_f (second term), buoyancy acceleration $a_{\mathrm{b}}=(\mathrm{0},\mathrm{0},b{)}^t$ (third term, the symbol ^t stands for the transpose of a vector) and the viscous acceleration a_ν (fourth term). The pressure term on the right-hand side ensures the vanishing divergence of the total acceleration, $\mathrm{\nabla }\cdot {\partial }_t\mathit{u}=\mathrm{0}$. Taking the divergence of Eq. (1), we can calculate the pressure by solving the elliptic equation

$$\begin{array}{}\text{(6)}& {\mathrm{\nabla }}^{\mathrm{2}}P=-\mathrm{\nabla }\cdot (a_u+a_{\mathrm{b}}+a_f+a_{\mathit{\nu }})\end{array}$$

subject to boundary conditions. Note that a_ν does not appear in Eq. (6) if the viscosity is isotropic and constant in space, as the correct implementation of the viscous term is $a_{\mathit{\nu }}=({\mathrm{\nabla }}^t\cdot \mathit{\nu }\mathrm{\nabla })\mathit{u}$ and $\mathrm{\nabla }\cdot \mathit{u}=\mathrm{0}$ implies $\mathrm{\nabla }\cdot a_{\mathit{\nu }}=\mathrm{0}$.

If the flow is in geostrophic equilibrium (or more precisely, satisfies the thermal wind relation) $a_{\mathrm{b}}=-a_f$, $a_u=a_{\mathit{\nu }}=\mathrm{0}$ and there is no time evolution of the velocity field. The quasi-geostrophic dynamics is based on the hydrostatic balance. If a three-dimensional perturbation is added to the dynamics, it will approach quasi-geostrophy if rotation is a dominant feature of the dynamics. How fast and how close this tendency towards quasi-geostrophy happens, at different scales, is today an open question, and it depends on the problem considered (see, e.g., Vanneste, 2013 and Wirth, 2013). Quasi-geostrophic flow is always close to balanced flow, which relies on the hydrostatic balance in the vertical and geostrophic balance in the horizontal. Viscosity perturbs the latter balance but not the former.

We can write the pressure as a sum of four parts, each ensuring the zero divergence of the corresponding term on the right-hand side of Eq. (12):

$$\begin{array}{}\text{(7)}& \begin{array}{rl}P(x,z,t)& =-P_u(x,z,t)+P_{\mathrm{f}}(x,z,t)+P_{\mathrm{b}}(x,z,t)\\ & +P_{\mathit{\nu }}(x,z,t).\end{array}\end{array}$$

The first is related to the advection:

$$\begin{array}{}\text{(8)}& {\mathrm{\nabla }}^{\mathrm{2}}{P}_u=-\mathrm{\nabla }\cdot \left(\right(\mathit{u}\cdot \mathrm{\nabla }\left)\mathit{u}\right),\end{array}$$

and it is nonlinear. The second is due to the Coriolis force:

$$\begin{array}{}\text{(9)}& {\mathrm{\nabla }}^{\mathrm{2}}{P}_{\mathrm{f}}=-\mathrm{\nabla }\cdot (f\mathit{z}\mathrm{\nabla }\times \mathit{u}).\end{array}$$

The third term is due to buoyancy:

$$\begin{array}{}\text{(10)}& {\mathrm{\nabla }}^{\mathrm{2}}{P}_{\mathrm{b}}=\mathrm{\nabla }\cdot b.\end{array}$$

Calculating the pressures above requires solving an elliptic (non-local) problem. It is instructive to distinguish the different terms, as the vertical acceleration of the of the right-hand side of Eqs. (8) and (9), only, are neglected in the hydrostatic approximation.

If the viscosity is homogeneous and isotropic, the fourth term,

$$\begin{array}{}\text{(11)}& {\mathrm{\nabla }}^{\mathrm{2}}{P}_{\mathit{\nu }}=\mathrm{\nabla }\cdot (\mathrm{\nabla }\cdot \mathit{\nu }\mathrm{\nabla }\mathit{u}),\end{array}$$

vanishes.

The linear operator ℙ projects the acceleration in the subspace of zero divergence (see Fig. 2). Equations (1) and (2) can therefore be written as

$$\begin{array}{}\text{(12)}& {\displaystyle }{\partial }_t\mathit{u}=-\mathbb{P}\left(\right(\mathit{u}\cdot \mathrm{\nabla }\left)\mathit{u}\right)-\mathbb{P}(f\mathit{z}\times \mathit{u})+\mathbb{P}\left(\mathit{z}b\right)-\mathbb{P}\left(\mathit{\nu }{\mathrm{\nabla }}^{\mathrm{2}}\mathit{u}\right)\text{(13)}& {\displaystyle }\mathrm{\nabla }\cdot \mathbb{P}\left(\mathit{u}\right)=\mathrm{0}.\end{array}$$

Note that although the projection is a linear operation, the pressure term is a nonlinear term as it depends nonlinearly on the velocity field through the nonlinear advection term. Separating each step of a fluid dynamics integration in an evolution and a projection sub-step is a classical procedure and is widely applied, also in ocean simulations (see discussion on projection in Marshall et al., 1997a; Guillaume et al., 2017; Popinet, 2020). In other numerical schemes (Auclair et al., 2018), the evolution stays close to the subspace of divergenceless vector fields by introducing a spurious compressibility.

Linear, elliptic, non-local problems with constant coefficients become local algebraic problems in Fourier space and their solution is reduced to solving linear algebraic systems. The orthogonal projection on the subspace of divergence-free functions is a local operation in Fourier space, and it is therefore natural to present our formalism in Fourier space. The dilemma with the Fourier space is that boundary conditions, which are local in physical space, become non-local in Fourier space, see Wirth (2005) for a detailed discussion and application to ocean modeling.

3 Fourier space projections (T1)

In this section I discuss the well-known technique (labeled T1 in the introduction) of projections of a vector field into the subspace of divergenceless vector fields for the Navier–Stokes formalism and develop the equivalent for the hydrostatic approximation (see Fig. 2). The Fourier transform of any L³-periodic function f is

$$\begin{array}{}\text{(14)}& \widehat{f}\left(\mathit{k}\right)={\displaystyle \frac{\mathrm{1}}{L^{\mathrm{3}}}}\underset{L^{\mathrm{3}}}{\int }{e}^{i(k_{x}x+k_{y}y+k_{z}z)}f(x,y,z)\mathrm{d}x\mathrm{d}y\mathrm{d}z,\end{array}$$

where the wave numbers are given by $k_x,k_y,k_z=\mathrm{2}\mathit{\pi }n/L$ and n∈ℤ.

Figure 2 Schematic view of the Fourier subspace spanned by the real part of the Fourier coefficients of the velocity (${\widehat{u}}_x(k_x,k_z),{\widehat{u}}_z(k_x,k_z)$) for the mode $\mathit{k}=(k_x,k_z\ne \mathrm{0})$, given by the thin red vector. The corresponding perpendicular subspace of zero divergence is shown by the double red line; it is orthogonal to the k vector. Points are velocities and vectors are accelerations multiplied by the time step Δt. One time step in the hydrostatic (blue) and Navier–Stokes (olive) formalism is shown. It consists in an evolution followed by a projection step. The evolution step (thin full olive vector) ${\partial }_t\widehat{\mathit{u}}$ starts from the initial velocity ${\widehat{u}}_i$. The buoyancy acceleration in the Navier–Stokes formalism is given by the dashed olive vector. When added to the dynamic vertical acceleration a_z it gives the total vertical acceleration. The Navier–Stokes projection (thick olive vector) is orthogonal onto the subspace of zero divergence, ending at $\widehat{\mathit{u}}$. In the hydrostatic formalism, the vertical acceleration is neglected during the projection step. This corresponds to a projection in the ${\widehat{u}}_z$ direction (thick dotted blue vector) leading to the thin full blue vector. The buoyancy correction in the hydrostatic formalism is given by the dashed blue line; it is in the horizontal and added to the horizontal acceleration. It is obtained by projecting the buoyancy acceleration along the wave vector (dashed red line, parallel to the wave vector) on the line defined by the horizontal acceleration. The third part of the hydrostatic projection (thick blue vector) is in the ${\widehat{u}}_z$ direction, down on the subspace of zero divergence. The hydrostatic projection is composed of three vectors (the thick dotted, dashed and full blue vectors.) The projection step in both formalisms end on the double red line, representing vanishing divergence, but at different locations. The difference between the hydrostatic and the Navier–Stokes formalism, the error, is given by the thick red vector.

A geometric illustration of the projections, for both formalisms is given in Fig. 2. The actual operators are presented below. For both formalisms, the integration starts from a divergence-free initial condition ${\widehat{\mathit{u}}}_i$ using Eq. (1) (see Fig. 2). For the integration step, we write the evolution equation, Eq. (12), without the projector and the buoyancy equation, Eq. (5), in Fourier space:

$$\begin{array}{}\text{(15)}& {\displaystyle }{\widehat{a}}_x=-i{k}_x\widehat{u_x^{\mathrm{2}}}-i{k}_y\widehat{u_x{u}_y}-i{k}_z\widehat{u_x{u}_z}+f\widehat{u_y}-\stackrel{\mathrm{˘}}{\mathit{\nu }}\widehat{u_x}\text{(16)}& {\displaystyle }{\widehat{a}}_y=-i{k}_x\widehat{u_x{u}_y}-i{k}_y\widehat{u_y^{\mathrm{2}}}-i{k}_z\widehat{u_y{u}_z}-f\widehat{u_x}-\stackrel{\mathrm{˘}}{\mathit{\nu }}\widehat{u_y}\text{(17)}& {\displaystyle }{\widehat{a}}_z=-i{k}_x\widehat{u_x{u}_z}-i{k}_y\widehat{u_y{u}_z}-i{k}_z\widehat{u_z^{\mathrm{2}}}-\stackrel{\mathrm{˘}}{\mathit{\nu }}\widehat{u_z},\text{(18)}& {\displaystyle }\dot{\widehat{b}}=-i{k}_x\widehat{u_{x}b}-i{k}_y\widehat{u_{y}b}-i{k}_z\widehat{u_{z}b}-\widehat{u_z}{N}^{\mathrm{2}}-\mathit{\kappa }{k}^{\mathrm{2}}\widehat{b},\end{array}$$

with the wavelength $k=\sqrt{k_x^{\mathrm{2}}+k_y^{\mathrm{2}}+k_z^{\mathrm{2}}}$, the inverse friction time $\stackrel{\mathrm{˘}}{\mathit{\nu }}=\mathit{\nu }{k}^{\mathrm{2}}$ and the Brunt–Väisälä frequency $N=\sqrt{{\partial }_z{b}_{\mathrm{BG}}}$ of the background stratification. The three-dimensional acceleration vector is

$$\begin{array}{}\text{(19)}& {\widehat{\mathit{a}}}^{\mathrm{N}}=\left(\begin{array}{c}{\widehat{a}}_x\\ {\widehat{a}}_y\\ {\widehat{a}}_z+\widehat{b}\end{array}\right).\end{array}$$

The hydrostatic pressure (Eq. 10) in Fourier space is $i{k}_z\widehat{P}=\widehat{b}$ and therefore $i{k}_x\widehat{P}=\frac{k_x}{k_z}\widehat{b}$, and the horizontal hydrostatic acceleration becomes $-\frac{(k_x,k_y{)}^t}{k_z}\widehat{b}$. The two-dimensional acceleration in the hydrostatic formalism is

$$\begin{array}{}\text{(20)}& {\widehat{\mathit{a}}}^{\mathrm{H}}=\left(\begin{array}{c}{\widehat{a}}_x-\frac{k_x}{k_z}\widehat{b}\\ {\widehat{a}}_y-\frac{k_y}{k_z}\widehat{b}\\ \mathrm{0}\end{array}\right).\end{array}$$

That is, the hydrostatic pressure is added to the horizontal acceleration and the vertical acceleration is omitted.

The integration step is followed by a projection on the divergence-free subspace, which is formed by vector fields perpendicular to the wave vector, i.e., $\mathbb{P}\left(\widehat{\mathit{a}}\right)\cdot \mathit{k}=\mathrm{0}$. It is orthogonal to the divergence-free subspace in the Navier–Stokes formalism:

$$\begin{array}{}\text{(21)}& {\stackrel{\mathrm{‾}}{\mathbb{P}}}^{\mathrm{N}}\left({\widehat{\mathit{a}}}^{\mathrm{N}}\right)={\widehat{\mathit{a}}}^{\mathrm{N}}-{\displaystyle \frac{(\mathit{k}\cdot {\widehat{\mathit{a}}}^{\mathrm{N}})\mathit{k}}{{\mathit{k}}^{\mathrm{2}}}}.\end{array}$$

This projection is local in Fourier space. It is performed independently for every wave vector k and makes the power of the pseudo-spectral method based on the Fourier series expansion of the velocity field. In Fourier space it is a matrix multiplication, whereas in real space it is a (non-local) integration over the whole domain. It was applied to ocean modeling in Wirth (2005).

For k_z=0, the hydrostatic projection equals the Navier–Stokes projection applied to the first two components. For k_z≠0, the projection is in the vertical direction in the hydrostatic formalism leaving the horizontal components unchanged (see Fig. 2):

$$\begin{array}{}\text{(22)}& {\displaystyle }{k}_z=\mathrm{0}:{\stackrel{\mathrm{‾}}{\mathbb{P}}}^{\mathrm{H}}\left({\widehat{\mathit{a}}}^{\mathrm{H}}\right)=\left(\begin{array}{ccc}\mathrm{1}-\frac{k_x^{\mathrm{2}}}{k_{\mathrm{h}}^{\mathrm{2}}}& -\frac{k_x{k}_y}{k_{\mathrm{h}}^{\mathrm{2}}}& \mathrm{0}\\ -\frac{k_x{k}_y}{k_{\mathrm{h}}^{\mathrm{2}}}& \mathrm{1}-\frac{k_x^{\mathrm{2}}}{k_{\mathrm{h}}^{\mathrm{2}}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right)\left(\begin{array}{c}{\widehat{a}}_x\\ {\widehat{a}}_y\\ \mathrm{0}\end{array}\right),\text{(23)}& {\displaystyle }{k}_z\ne \mathrm{0}:{\stackrel{\mathrm{‾}}{\mathbb{P}}}^{\mathrm{H}}\left({\widehat{\mathit{a}}}^{\mathrm{H}}\right)=\left(\begin{array}{ccc}\mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}\\ -\frac{k_x}{k_z}& -\frac{k_y}{k_z}& \mathrm{0}\end{array}\right){\widehat{\mathit{a}}}_H\end{array}$$

with $k_{\mathrm{h}}^{\mathrm{2}}=\sqrt{k_x^{\mathrm{2}}+k_y^{\mathrm{2}}}$. It is clearly seen that the vertical acceleration prior to the projection is ignored; the entries of the third column vanish. The vertical acceleration is obtained, through the projection, entirely based on the horizontal components, such that the velocity field is divergenceless. Note that the correction for k_z≠0 are of order

$$\begin{array}{}\text{(24)}& \mathit{\gamma }=\sqrt{k_{\mathrm{h}}^{\mathrm{2}}/k_z^{\mathrm{2}}},\end{array}$$

which is small for dynamical processes which have a horizontal extension larger than the vertical. Both projectors lead to a divergenceless vector field, which is readily verified by multiplying the projectors by the wave vector, k, on their left.

The major difference between the hydrostatic and the Navier–Stokes projection is that in the latter the vertical acceleration from the dynamical and gravitational acceleration are treated equally, while in the former the gravitational acceleration is added to the horizontal acceleration using the hydrostatic equation and the vertical dynamical acceleration is calculated from its horizontal counterparts. A geometrical interpretation of the two projections and their differences is shown in Fig. 2.

In numerical models based on Fourier representation, the projection is exact at the machine precision employed. In other models, the projection in the Navier–Stokes formalism is done via an elliptic solver of which the precision is specified and usually far less than the machine precision. In this case, the double line in Fig. 2, representing the subspace of zero divergence, has a finite thickness corresponding to the precision prescribed to the elliptic solver. The question of how to employ Fourier methods in ocean modeling, where boundary conditions have to be satisfied, is discussed in detail in Wirth (2005).

4 Fourier space extension (T2)

In this section, I introduce the extension of the Fourier space to four dimensions; this technique is labeled T2 in the introduction and in Fig. 3. It allows to use a common mathematical formalism, based on four-dimensional linear algebra, for the Navier–Stokes and the hydrostatic projection. In the four-dimensional space, both projections apply to the same vector ($\stackrel{\mathrm{̃}}{\mathit{a}}$, introduced below) and it is then straightforward to obtain the difference between the two projections as four-dimensional matrices (𝔸, given at the end of this section). Only the entries of the matrices differ between the two formalisms. Increasing the dimension from three to four is necessary as the two parts of the vertical acceleration, dynamical and buoyancy, are treated differently between the formalisms.

Figure 3 Schematic of the projection-evolution formalism. Step S1 gives the difference in the projection of a forcing in the two formalisms and step S2 the difference in the linear evolution. Nonlinear terms (NL) can be represented as (time-dependent) forcings and the scheme has to be iterated. The different forcings (vertical F_z, divergent F_d, rotational F_r and buoyancy flux B) are projected on wave forcing F_w and balanced forcing F_b, which span the subspace of divergenceless vector fields. The rotational and the buoyancy flux forcings are identical in both formalisms. The vertical forcing F_z is in the kernel of the hydrostatic projection (missing arrow) and the error of the divergent forcing F_d due to the hydrostatic projection is second order O(γ²) (orange arrow). The projection-evolution operator acts identically on balanced forcing in both formalisms, whereas it differs to second order O(γ²) (orange arrow) on the wave forcing (the symbol □ is a place holder for “N” for the Navier–Stokes or “H” for the hydrostatic formalism). The nonlinear terms formed by the wave and balanced amplitudes appear as retro-acting forces (NL) in the formalism.

To emphasize the similarities and differences of both projections, I split the vertical acceleration into two parts: dynamical (${\widehat{a}}_z$) and buoyant (${\widehat{a}}_b=\widehat{b}$) part:

$$\begin{array}{}\text{(25)}& \stackrel{\mathrm{̃}}{\mathit{a}}=\left(\begin{array}{c}{\widehat{a}}_x\\ {\widehat{a}}_y\\ {\widehat{a}}_z\\ {\widehat{a}}_b\end{array}\right).\end{array}$$

The buoyancy adds to the dynamical vertical acceleration. In the Navier–Stokes formalism, there is no difference in the projection of the two, while in the hydrostatic projection, the dynamical vertical acceleration is ignored, while the buoyancy acceleration appears in the horizontal component. More precisely, the hydrostatic projection becomes

$$\begin{array}{}\text{(26)}& \begin{array}{rl}& {\mathbb{P}}^{\mathrm{H}}=\left(\begin{array}{cccc}\mathrm{1}& \mathrm{0}& \mathrm{0}& -\frac{k_x}{k_z}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}& -\frac{k_y}{k_z}\\ -\frac{k_x}{k_z}& -\frac{k_y}{k_z}& \mathrm{0}& \frac{k_{\mathrm{h}}^{\mathrm{2}}}{k_z^{\mathrm{2}}}\end{array}\right)\\ \\ & =\mathbb{1}-{\displaystyle \frac{\mathrm{1}}{k_z}}\left[\left(\begin{array}{c}\mathrm{0}\\ \mathrm{0}\\ \mathrm{1}\end{array}\right)(k_x,k_y,\mathrm{0},-{\displaystyle \frac{k_{\mathrm{h}}^{\mathrm{2}}}{k_z}})+\left(\begin{array}{c}{k}_x\\ k_y\\ \mathrm{0}\end{array}\right)(\mathrm{0},\mathrm{0},\mathrm{0},\mathrm{1})\right].\end{array}\end{array}$$

The third column is zero as the dynamical vertical acceleration is ignored, while it can be seen in the fourth column that the gravitational acceleration acts on the horizontal and vertical. The Navier–Stokes projection is

$$\begin{array}{}\text{(27)}& \begin{array}{rl}{\mathbb{P}}^{\mathrm{N}}& =\left(\begin{array}{ccc}\mathrm{1}-\frac{k_x^{\mathrm{2}}}{k^{\mathrm{2}}}& -\frac{k_x{k}_y}{k^{\mathrm{2}}}-\frac{k_x{k}_z}{k^{\mathrm{2}}}& -\frac{k_x{k}_z}{k^{\mathrm{2}}}\\ -\frac{k_x{k}_y}{k^{\mathrm{2}}}& \mathrm{1}-\frac{k_y^{\mathrm{2}}}{k^{\mathrm{2}}}& -\frac{k_y{k}_z}{k^{\mathrm{2}}}& -\frac{k_y{k}_z}{k^{\mathrm{2}}}\\ -\frac{k_x{k}_z}{k^{\mathrm{2}}}& -\frac{k_y{k}_z}{k^{\mathrm{2}}}& \mathrm{1}-\frac{k_z^{\mathrm{2}}}{k^{\mathrm{2}}}& \mathrm{1}-\frac{k_z^{\mathrm{2}}}{k^{\mathrm{2}}}\end{array}\right)\\ & =\mathbb{1}-{\displaystyle \frac{\mathrm{1}}{k^{\mathrm{2}}}}\left(\begin{array}{c}{k}_x\\ k_y\\ k_z\end{array}\right)(k_x,k_y,k_z,k_z).\end{array}\end{array}$$

The third and fourth columns are identical, as the vertical acceleration due to the dynamical and gravitational parts are treated equally. If buoyancy vanishes, the divergence-free subspace is spanned by the orthogonal vectors ${\mathit{e}}_{\mathrm{1}}^{\mathrm{P}}=(k_y,-k_x,\mathrm{0},\mathrm{0}{)}^t$ and ${\mathit{e}}_{\mathrm{2}}^{\mathrm{P}}=(k_x,k_y,-\frac{k_{\mathrm{h}}^{\mathrm{2}}}{k_z},\mathrm{0}{)}^t$. There is no difference between the Navier–Stokes and the hydrostatic projections if applied to these two vectors. The first vector represents a divergenceless horizontal rotational motion with no vertical component. In two dimensions, this is the only way to have a divergenceless mode. The second vector is a compensation between the horizontal and vertical divergence; it represents vertical stretching and compression of the water column. The vertical component of the curl of the second vector is vanishing.

To emphasize the projection on ${\mathit{e}}_{\mathrm{1}}^{\mathrm{P}}$ and ${\mathit{e}}_{\mathrm{2}}^{\mathrm{P}}$, the projection operators can also be written as

$$\begin{array}{}\text{(28)}& {\displaystyle }\begin{array}{rl}{\mathbb{P}}^{\mathrm{H}}& ={\displaystyle \frac{\mathrm{1}}{k_{\mathrm{h}}^{\mathrm{2}}}}\left[\left(\begin{array}{c}{k}_y\\ -k_x\\ \mathrm{0}\end{array}\right)(k_y,-k_x,\mathrm{0},\mathrm{0})\right.\\ & \left.+\left(\begin{array}{c}{k}_x\\ k_y\\ -\frac{k_{\mathrm{h}}^{\mathrm{2}}}{k_z}\end{array}\right)(k_x,k_y,\mathrm{0},-{\displaystyle \frac{k_{\mathrm{h}}^{\mathrm{2}}}{k_z}})\right],\end{array}\text{(29)}& {\displaystyle }\begin{array}{rl}{\mathbb{P}}^{\mathrm{N}}& ={\displaystyle \frac{\mathrm{1}}{k_{\mathrm{h}}^{\mathrm{2}}}}\left[\left(\begin{array}{c}{k}_y\\ -k_x\\ \mathrm{0}\end{array}\right)(k_y,-k_x,\mathrm{0},\mathrm{0})\right.\\ & \left.+{\displaystyle \frac{k_z^{\mathrm{2}}}{k^{\mathrm{2}}}}\left(\begin{array}{c}{k}_x\\ k_y\\ -\frac{k_{\mathrm{h}}^{\mathrm{2}}}{k_z}\end{array}\right)(k_x,k_y,-{\displaystyle \frac{k_{\mathrm{h}}^{\mathrm{2}}}{k_z}},-{\displaystyle \frac{k_{\mathrm{h}}^{\mathrm{2}}}{k_z}})\right].\end{array}\end{array}$$

There is no difference in the projection of the purely horizontal divergenceless motion, which is the first term in both projections, whereas, if horizontally divergent, buoyant and vertical motion are considered, there are two differences: the horizontally divergent and buoyant accelerations differ by a factor $k_z^{\mathrm{2}}/k^{\mathrm{2}}=(\mathrm{1}+{\mathit{\gamma }}^{\mathrm{2}}{)}^{-\mathrm{1}}$, and the vertical dynamic acceleration is neglected in the hydrostatic approximation. The former difference is second order O(γ²), while the latter is first order O(γ); a higher order can be assigned if the dynamics constrain the order of the vertical dynamical acceleration ${\widehat{a}}_z$; this will be discussed in Sect. 5 based on the linearized, low-amplitude dynamics. Note that the projections are closely related to the well-known Craya decomposition or wave–vortex compositions; see Deriaz et al. (2010) and references therein for a detailed discussion.

The vector ${\mathit{e}}_{\mathrm{3}}^{\mathrm{P}}=(k_x,k_y,\mathrm{0},k_z{)}^t$ is in the kernel of both projections and represents balanced flow: the hydrostatic balance in the vertical, ${\widehat{a}}_z=\mathrm{0}$, saying that all the acceleration in the vertical is due to buoyancy $\widehat{b}$, and the geostrophic balance (thermal wind balance, see Vallis, 2017) in the horizontal: ${\widehat{a}}_x=f{\widehat{u}}_y=\frac{k_x}{k_z}\widehat{b}$ and ${\widehat{a}}_y=-f{\widehat{u}}_x=\frac{k_y}{k_z}\widehat{b}$. These compensations are the essence of geophysical fluid dynamics (GFD) – that is, the dynamics of a buoyant fluid in a rotating frame. It includes geostrophy, thermal wind balance and cyclostrophic balance.

The difference between the two projections lies in the second vector in the kernel; for the hydrostatic operator it is formally ${\mathit{e}}_{\mathrm{0}}^{\mathrm{H}}=(\mathrm{0},\mathrm{0},\mathrm{1},\mathrm{0})$, but as in the hydrostatic formalism, ${\widehat{a}}_z=\mathrm{0}$, it has zero amplitude. For the Navier–Stokes operator, the second vector in the kernel is ${\mathit{e}}_{\mathrm{0}}^{\mathrm{N}}=(\mathrm{0},\mathrm{0},\mathrm{1},-\mathrm{1})$; it expresses the trivial fact that compensations between the dynamical and gravitational accelerations are in the kernel of the operator. In the hydrostatic approximation (${\widehat{a}}_z=\mathrm{0}$), no such compensation exists. Furthermore, ${\mathit{e}}_{\mathrm{00}}^{\mathrm{N}}=k_z{\mathit{e}}_{\mathrm{0}}^{\mathrm{N}}+{\mathit{e}}_{\mathrm{3}}^{\mathrm{P}}=(k_x,k_y,k_z,\mathrm{0}{)}^t$ is also in the kernel of the Navier–Stokes operator; if multiplied by ${\widehat{P^{\prime }}}_{\mathrm{d}}$ it is nothing other than the dynamic pressure correction $\mathrm{\nabla }{P}_{\mathrm{d}}^{\prime }$. Note that ${\mathit{e}}_{\mathrm{1}}^{\mathrm{P}},{\mathit{e}}_{\mathrm{2}}^{\mathrm{P}},{\mathit{e}}_{\mathrm{00}}^{\mathrm{N}}$ are an orthogonal subset.

For k_z≠0, the difference between the operators is

$$\begin{array}{}\text{(30)}& \begin{array}{rl}& \mathbb{A}\left(\stackrel{\mathrm{̃}}{\mathit{a}}\right)={\mathbb{P}}^{\mathrm{H}}\left(\stackrel{\mathrm{̃}}{\mathit{a}}\right)-{\mathbb{P}}^{\mathrm{N}}\left(\stackrel{\mathrm{̃}}{\mathit{a}}\right)=\\ & =\left[\left(\begin{array}{cccc}\mathrm{1}& \mathrm{0}& \mathrm{0}& -\frac{k_x}{k_z}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}& -\frac{k_y}{k_z}\\ -\frac{k_x}{k_z}& -\frac{k_y}{k_z}& \mathrm{0}& \frac{k_{\mathrm{h}}^{\mathrm{2}}}{k_z^{\mathrm{2}}}\end{array}\right)\right.\\ & \left.-\left(\begin{array}{cccc}\mathrm{1}-\frac{k_x^{\mathrm{2}}}{k^{\mathrm{2}}}& -\frac{k_x{k}_y}{k^{\mathrm{2}}}& -\frac{k_x{k}_z}{k^{\mathrm{2}}}& -\frac{k_x{k}_z}{k^{\mathrm{2}}}\\ -\frac{k_x{k}_y}{k^{\mathrm{2}}}& \mathrm{1}-\frac{k_y^{\mathrm{2}}}{k^{\mathrm{2}}}& -\frac{k_y{k}_z}{k^{\mathrm{2}}}& -\frac{k_y{k}_z}{k^{\mathrm{2}}}\\ -\frac{k_x{k}_z}{k^{\mathrm{2}}}& -\frac{k_y{k}_z}{k^{\mathrm{2}}}& \mathrm{1}-\frac{k_z^{\mathrm{2}}}{k^{\mathrm{2}}}& \mathrm{1}-\frac{k_z^{\mathrm{2}}}{k^{\mathrm{2}}}\end{array}\right)\right]\stackrel{\mathrm{̃}}{\mathit{a}}\\ & =\underset{\mathbb{A}}{\underbrace{{\displaystyle \frac{\mathrm{1}}{k_z{k}^{\mathrm{2}}}}\left(\begin{array}{cccc}{k}_x^{\mathrm{2}}{k}_z& k_x{k}_y{k}_z& k_x{k}_z^{\mathrm{2}}& -k_x{k}_{\mathrm{h}}^{\mathrm{2}}\\ k_x{k}_y{k}_z& k_y^{\mathrm{2}}{k}_z& k_y{k}_z^{\mathrm{2}}& -k_y{k}_{\mathrm{h}}^{\mathrm{2}}\\ -k_x{k}_{\mathrm{h}}^{\mathrm{2}}& -k_y{k}_{\mathrm{h}}^{\mathrm{2}}& -k_z{k}_h^{\mathrm{2}}& (k_{\mathrm{h}}^{\mathrm{2}}{)}^{\mathrm{2}}/k_z\end{array}\right)}}\stackrel{\mathrm{̃}}{\mathit{a}}\\ & ={\displaystyle \frac{\mathrm{1}}{k^{\mathrm{2}}}}\left(\begin{array}{c}{k}_x\\ k_y\\ -\frac{k_{\mathrm{h}}^{\mathrm{2}}}{k_z}\end{array}\right)(k_x,k_y,k_z,-{\displaystyle \frac{k_{\mathrm{h}}^{\mathrm{2}}}{k_z}})\stackrel{\mathrm{̃}}{\mathit{a}}.\end{array}\end{array}$$

By construction, the difference of the two operators is in the space of divergenceless functions spanned by ${\mathit{e}}_{\mathrm{1}}^{\mathrm{P}}$ and ${\mathit{e}}_{\mathrm{2}}^{\mathrm{P}}$ (with the last component omitted), but the above shows that it is aligned with ${\mathit{e}}_{\mathrm{2}}^{\mathrm{P}}$. This illustrates that the difference between the hydrostatic and the Navier–Stokes operator is never in the vertical vorticity, but in the compensation of horizontal divergence through vertical divergence, the water column is stretched differently if $(k_x,k_y,k_z,-\frac{k_{\mathrm{h}}^{\mathrm{2}}}{k_z})\stackrel{\mathrm{̃}}{\mathit{a}}\ne \mathrm{0}$. The Navier–Stokes projector does not lead to a change in vorticity, while a non-vanishing hydrostatic correction alters the horizontal components of the vorticity vector only.

If there is no constraint on the acceleration vector, $\sqrt{{\left[\mathbb{A}\left(\stackrel{\mathrm{̃}}{\mathit{a}}\right)\right]}^{\mathrm{2}}}<\mathit{\gamma }\sqrt{{\stackrel{\mathrm{̃}}{\mathit{a}}}^{\mathrm{2}}}$ if γ<1. Furthermore, the error in the vertical is by a factor γ smaller than the horizontal error. If the vertical dynamic acceleration is γ-times smaller than the horizontal or buoyant accelerations, $\sqrt{{\left[\mathbb{A}\left(\stackrel{\mathrm{̃}}{\mathit{a}}\right)\right]}^{\mathrm{2}}}\le {\mathit{\gamma }}^{\mathrm{2}}\sqrt{{\stackrel{\mathrm{̃}}{\mathit{a}}}^{\mathrm{2}}}$. Note that for the hydrostatic approximation to be of second order, O(γ²), only the vertical dynamic acceleration has to be small, not the buoyant. The relative error is $\mathbb{A}\left(\stackrel{\mathrm{̃}}{\mathit{a}}\right)/{\mathbb{P}}^{\mathrm{N}}\left(\stackrel{\mathrm{̃}}{\mathit{a}}\right)$; for $\stackrel{\mathrm{̃}}{\mathit{a}}={\mathit{e}}_{\mathrm{0}}^{\mathrm{N}}$ it is ∞.

So far, I have only considered the error in the hydrostatic approximation for an acceleration field, without any reference to the underlying dynamics. In the next section, I will discuss the error for low-amplitude motion if the dynamics is well approximated by the linearized equations.

5 Linearized equations (T3)

In this section, I combine the projections in the extended Fourier space with the linearized evolution equation. If the velocity is small, the first term on the right-hand side of Eq. (6) is negligible and the nonlinear terms in Eqs. (15), (16), (17) and (18) can be neglected. When we further introduce a buoyancy velocity $u_{\mathrm{b}}=b/N$ and write the four-dimensional velocity vector $\widehat{\mathit{u}}=(u_x,u_y,u_z,u_{\mathrm{b}}{)}^t$, the linearized dynamics is given by

$$\begin{array}{}\text{(31)}& \left(\begin{array}{c}\stackrel{\mathrm{̃}}{\mathit{a}}\\ {\dot{\widehat{u}}}_{\mathrm{b}}\end{array}\right)={\mathbb{L}}^{\square }\widehat{\mathit{u}},\end{array}$$

where the symbol □ is a placeholder for “N” for the Navier–Stokes or “H” for the hydrostatic formalism. Note that the vector on the right-hand side has five dimensions as it consists of the accelerations was well as the buoyancy and the buoyancy flux divided by the Brunt–Väisälä frequency. This normalization renders the vectors and operators, introduced below, dimensionally homogeneous. Note that the extension to five dimensions is necessary as, firstly, the two formalisms treat the buoyancy and the vertical dynamic acceleration differently and I want to present both in the same formalism, and, secondly, the vertical velocity acts on the time derivative of the buoyancy through the background stratification.

The extension to five dimensions is necessary to calculate the vertical acceleration from the dynamical and buoyant part and to keep track of the buoyancy anomaly, which is unchanged by the projection operator. The explicit form of the operators are

$$\begin{array}{}\text{(32)}& \begin{array}{rl}& {\mathbb{L}}^{\mathrm{N}}=\left(\begin{array}{cccc}-\stackrel{\mathrm{˘}}{\mathit{\nu }}& f& \mathrm{0}& \mathrm{0}\\ -f& -\stackrel{\mathrm{˘}}{\mathit{\nu }}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& -\stackrel{\mathrm{˘}}{\mathit{\nu }}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& N\\ \mathrm{0}& \mathrm{0}& -N& -\mathit{\kappa }{k}^{\mathrm{2}}\end{array}\right)\mathrm{and}\\ & {\mathbb{L}}^{\mathrm{H}}=\left(\begin{array}{cccc}-\stackrel{\mathrm{˘}}{\mathit{\nu }}& f& \mathrm{0}& \mathrm{0}\\ -f& -\stackrel{\mathrm{˘}}{\mathit{\nu }}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& -\stackrel{\mathrm{˘}}{\mathit{\nu }}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& N\\ N\frac{k_x}{k_z}& N\frac{k_y}{k_z}& \mathrm{0}& -\mathit{\kappa }{k}^{\mathrm{2}}\end{array}\right).\end{array}\end{array}$$

I used ${\partial }_{t}b=-u_z{N}^{\mathrm{2}}-\mathit{\kappa }{k}^{\mathrm{2}}b$.

We define the extended projection operator:

$$\begin{array}{}\text{(33)}& {\stackrel{\mathrm{̃}}{\mathbb{P}}}^{\square }=\left(\begin{array}{cc}{\mathbb{P}}^{\square }& \mathrm{0}\\ \mathrm{0}& \mathrm{1}\end{array}\right).\end{array}$$

The change in buoyancy is unaffected by the divergenceless condition, but the acceleration has to be projected into the subspace of divergence-free vector fields. The operator ${\stackrel{\mathrm{̃}}{\mathbb{L}}}^{\square }={\stackrel{\mathrm{̃}}{\mathbb{P}}}^{\square }{\mathbb{L}}^{\square }$ governs the time evolution via

$$\begin{array}{}\text{(34)}& {\partial }_t{\widehat{\mathit{u}}}^{\square }\left(t\right)={\stackrel{\mathrm{̃}}{\mathbb{L}}}^{\square }{\widehat{\mathit{u}}}^{\square }\left(t\right)+{\stackrel{\mathrm{̃}}{\mathbb{P}}}^{\square }\widehat{\mathit{F}}\left(t\right).\end{array}$$

The force $\widehat{\mathit{F}}=({\widehat{F}}_x,{\widehat{F}}_y,{\widehat{F}}_z,\mathrm{0},\widehat{B}{)}^t$ is formed of forces in the first three components. The fourth component is vanishing as there is no buoyancy imposed by the exterior but a buoyancy flux, which is in the fifth component, when divided by the Brunt-Väisälä frequency. The forces and fluxes are divided by the reference density and represent accelerations and relative fluxes in the Boussinesq approximation used here. The explicit forms of the evolution operators are

$$\begin{array}{}\text{(35)}& {\displaystyle }{\stackrel{\mathrm{̃}}{\mathbb{L}}}^{\mathrm{N}}=\left(\begin{array}{cccc}-\mathit{\nu }(k^{\mathrm{2}}-k_x^{\mathrm{2}})+f{\displaystyle \frac{k_x{k}_y}{k^{\mathrm{2}}}}& \mathit{\nu }{k}_x{k}_y+f\frac{k^{\mathrm{2}}-k_x^{\mathrm{2}}}{k^{\mathrm{2}}}& \mathit{\nu }{k}_x{k}_z& -N\frac{k_x{k}_z}{k^{\mathrm{2}}}\\ \mathit{\nu }{k}_x{k}_y-f\frac{k^{\mathrm{2}}-k_y^{\mathrm{2}}}{k^{\mathrm{2}}}& -\mathit{\nu }(k^{\mathrm{2}}-k_y^{\mathrm{2}})-f\frac{k_x{k}_y}{k^{\mathrm{2}}}& \mathit{\nu }{k}_y{k}_z& -N\frac{k_y{k}_z}{k^{\mathrm{2}}}\\ \mathit{\nu }{k}_x{k}_z+f\frac{k_y{k}_z}{k^{\mathrm{2}}}& \mathit{\nu }{k}_y{k}_z-f\frac{k_x{k}_z}{k^{\mathrm{2}}}& -\mathit{\nu }{k}_{\mathrm{h}}^{\mathrm{2}}& N\frac{k_{\mathrm{h}}^{\mathrm{2}}}{k^{\mathrm{2}}}\\ \mathrm{0}& \mathrm{0}& -N& -\mathit{\kappa }{k}^{\mathrm{2}}\end{array}\right),\text{(36)}& {\displaystyle }{\stackrel{\mathrm{̃}}{\mathbb{L}}}^{\mathrm{H}}=\left(\begin{array}{cccc}-\stackrel{\mathrm{˘}}{\mathit{\nu }}& f& \mathrm{0}& -N\frac{k_x}{k_z}\\ -f& -\stackrel{\mathrm{˘}}{\mathit{\nu }}& \mathrm{0}& -N\frac{k_y}{k_z}\\ \stackrel{\mathrm{˘}}{\mathit{\nu }}\frac{k_x}{k_z}+f\frac{k_y}{k_z}& \stackrel{\mathrm{˘}}{\mathit{\nu }}\frac{k_y}{k_z}-f\frac{k_x}{k_z}& \mathrm{0}& N\frac{k_{\mathrm{h}}^{\mathrm{2}}}{k_z^{\mathrm{2}}}\\ N\frac{k_x}{k_z}& N\frac{k_y}{k_z}& \mathrm{0}& -\mathit{\kappa }{k}^{\mathrm{2}}\end{array}\right).\end{array}$$

If diffusivity equals viscosity, $\mathit{\kappa }=\mathit{\nu }>\mathrm{0}$, both matrices have two complex conjugated eigenvalues ${\mathit{\lambda }}_{\mathrm{1},\mathrm{2}}^{\square }=\pm i{\mathit{\omega }}^{\square }-\stackrel{\mathrm{˘}}{\mathit{\nu }}$, where

$$\begin{array}{}\text{(37)}& ({\mathit{\omega }}^{\mathrm{N}}{)}^{\mathrm{2}}={\mathit{\omega }}^{\mathrm{2}}={\displaystyle \frac{\mathit{\alpha }}{k^{\mathrm{2}}}},({\mathit{\omega }}^{\mathrm{H}}{)}^{\mathrm{2}}={\displaystyle \frac{\mathit{\alpha }}{k_z^{\mathrm{2}}}}\mathrm{with}\mathit{\alpha }=f^{\mathrm{2}}{k}_z^{\mathrm{2}}+N^{\mathrm{2}}{k}_{\mathrm{h}}^{\mathrm{2}}\end{array}$$

are the dispersion relations of internal waves in the Navier–Stokes and the hydrostatic formalisms, respectively. They are equal if and only if the horizontal wave numbers vanish, k_h=0. Note that for the Navier–Stokes formalism, the frequency of the waves, ω, is between the Coriolis f and the buoyancy frequency N, whereas the frequency for hydrostatic waves, ω^H, is above the Coriolis frequency. For $k_{\mathrm{h}}^{\mathrm{2}}>k_z^{\mathrm{2}}$, the hydrostatic approximation allows for waves with square frequencies faster than the maximum of f² and N². I use the superscript ^□ when either formalism is concerned. The other two eigenvalues are ${\mathit{\lambda }}_{\mathrm{3}}=-\stackrel{\mathrm{˘}}{\mathit{\nu }}$ and λ₄=0. The corresponding eigenvectors are

$$\begin{array}{}\text{(38)}& \begin{array}{rl}& {\mathit{e}}_{\mathrm{1},\mathrm{2}}^{\square }=\left(\begin{array}{c}\pm i{\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }\mathit{\omega }\frac{k_x}{k_{\mathrm{h}}}+f\frac{k_y}{k_{\mathrm{h}}}\\ \pm i{\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }\mathit{\omega }\frac{k_y}{k_{\mathrm{h}}}-f\frac{k_x}{k_{\mathrm{h}}}\\ \mp \frac{i{\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }\mathit{\omega }{k}_{\mathrm{h}}}{k_z}\\ \frac{N{k}_{\mathrm{h}}}{k_z}\end{array}\right),{\mathit{e}}_{\mathrm{3}}={\displaystyle \frac{f}{k_{\mathrm{h}}}}\left(\begin{array}{c}{k}_y\\ -k_x\\ \mathrm{0}\\ -\frac{f}{N}{k}_z\end{array}\right),\\ & {\mathit{e}}_{\mathrm{4}}^{\mathrm{N}}={\displaystyle \frac{\mathrm{1}}{f^{\mathrm{2}}{k}_{\mathrm{h}}}}\left(\begin{array}{c}({\stackrel{\mathrm{˘}}{\mathit{\nu }}}^{\mathrm{2}}+N^{\mathrm{2}})\left[f{k}_y+\stackrel{\mathrm{˘}}{\mathit{\nu }}{k}_x\right]\\ ({\stackrel{\mathrm{˘}}{\mathit{\nu }}}^{\mathrm{2}}+N^{\mathrm{2}})\left[-f{k}_x+\stackrel{\mathrm{˘}}{\mathit{\nu }}{k}_y\right]\\ \stackrel{\mathrm{˘}}{\mathit{\nu }}(f^{\mathrm{2}}+{\stackrel{\mathrm{˘}}{\mathit{\nu }}}^{\mathrm{2}})k_z\\ -N(f^{\mathrm{2}}+{\stackrel{\mathrm{˘}}{\mathit{\nu }}}^{\mathrm{2}})k_z\end{array}\right),{\mathit{e}}_{\mathrm{4}}^{\mathrm{H}}=\left(\begin{array}{c}\mathrm{0}\\ \mathrm{0}\\ f\\ \mathrm{0}\end{array}\right).\end{array}\end{array}$$

The first two vectors present inertial oscillations which are fast motion. We see that the eigenvectors for the fast motion and the corresponding eigenvalues differ by ${\mathit{\omega }}^{\mathrm{H}}={\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\mathrm{H}}\mathit{\omega }$, where we set ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\mathrm{H}}=\sqrt{\mathrm{1}+{\mathit{\gamma }}^{\mathrm{2}}}$ and ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\mathrm{N}}=\mathrm{1}$. This shows that the error is of second order, O(γ²). The first two, complex-valued, vectors can be substituted by the real orthogonal vectors,

$$\begin{array}{}\text{(39)}& {\stackrel{\mathrm{̃}}{\mathit{e}}}_{\mathrm{1}}=\left(\begin{array}{c}f\frac{k_y}{k_{\mathrm{h}}}\\ -f\frac{k_x}{k_{\mathrm{h}}}\\ \mathrm{0}\\ \frac{N{k}_{\mathrm{h}}}{k_z}\end{array}\right),{\stackrel{\mathrm{̃}}{\mathit{e}}}_{\mathrm{2}}=\mathit{\omega }\left(\begin{array}{c}\frac{k_x}{k_{\mathrm{h}}}\\ \frac{k_y}{k_{\mathrm{h}}}\\ -\frac{k_{\mathrm{h}}}{k_z}\\ \mathrm{0}\end{array}\right),\end{array}$$

of which the directions are independent of the projection, Navier–Stokes or hydrostatic, which facilitates the comparison of the two formalisms.

The basis $D_{\mathrm{4}}^{\square }=({\stackrel{\mathrm{̃}}{\mathit{e}}}_{\mathrm{1}},{\stackrel{\mathrm{̃}}{\mathit{e}}}_{\mathrm{2}},{\mathit{e}}_{\mathrm{3}},{\mathit{e}}_{\mathrm{4}}^{\square })$ is dimensionally homogeneous and we have ${\stackrel{\mathrm{̃}}{\mathbb{L}}}^{\square }{\stackrel{\mathrm{̃}}{\mathit{e}}}_{\mathrm{1}}=-{\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }\mathit{\omega }{\stackrel{\mathrm{̃}}{\mathit{e}}}_{\mathrm{2}}$ and ${\stackrel{\mathrm{̃}}{\mathbb{L}}}^{\square }{\stackrel{\mathrm{̃}}{\mathit{e}}}_{\mathrm{2}}=({\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }{)}^{-\mathrm{1}}\mathit{\omega }{\stackrel{\mathrm{̃}}{\mathit{e}}}_{\mathrm{1}}$, where $\stackrel{\mathrm{̃}}{\mathit{\gamma }}=\mathrm{1}$ in the Navier–Stokes formalism. The third vector represents the thermal wind balance, a generalization of geostrophy. Its eigenvalue vanishes with viscosity; it converges to a stationary state. We see that the first three vectors represent divergenceless motion as $({\stackrel{\mathrm{̃}}{\mathit{e}}}_i{)}^t\cdot ({\mathit{k}}^t,\mathrm{0})=\mathrm{0}$, $\forall i=\mathrm{1},\mathrm{2},\mathrm{3}$. They are identical in both formalisms. The fourth vector has a divergence: $({\stackrel{\mathrm{̃}}{\mathit{e}}}_{\mathrm{4}}{)}^t\cdot ({\mathit{k}}^t,\mathrm{0})=\stackrel{\mathrm{˘}}{\mathit{\nu }}({\stackrel{\mathrm{˘}}{\mathit{\nu }}}^{\mathrm{2}}+\mathit{\alpha })$ and $({\stackrel{\mathrm{̃}}{\mathit{e}}}_{\mathrm{4}}^{\mathrm{H}}{)}^t\cdot ({\mathit{k}}^t,\mathrm{0})=k_z$, in both formalisms and therefore does not participate in the dynamics. The error between the projections along the fourth direction, with a vanishing eigenvalue, is first order O(γ¹) if $f,N\ll \stackrel{\mathrm{˘}}{\mathit{\nu }}$ and zeroth order O(γ⁰) if ν→0 as the vectors ${\mathit{e}}_{\mathrm{4}}^{\mathrm{H}}$ and e₄ become orthogonal. A perturbation in the vertical acceleration triggers inertial-gravity waves in the Navier–Stokes formalism, whereas it has no effect in the hydrostatic case. If viscosity is zero, the matrix ${\stackrel{\mathrm{̃}}{\mathbb{L}}}^{\mathrm{N}}$ is not diagonalizable and has a Jordan normal form. I will therefore in the sequel discuss the case ν→0 rather than ν=0, which is also appropriate if numerical models are considered. Note that the first three components of ${\stackrel{\mathrm{̃}}{\mathit{e}}}_{\mathrm{1}}$ and e₃ are aligned with ${\mathit{e}}_{\mathrm{1}}^{\mathrm{P}}$, and those of ${\stackrel{\mathrm{̃}}{\mathit{e}}}_{\mathrm{2}}$ with ${\mathit{e}}_{\mathrm{2}}^{\mathrm{P}}$. The matrix that transforms from the sub-basis $D=({\stackrel{\mathrm{̃}}{\mathit{e}}}_{\mathrm{1}},{\stackrel{\mathrm{̃}}{\mathit{e}}}_{\mathrm{2}},{\mathit{e}}_{\mathrm{3}})$ to the Cartesian coordinates is formed by the three basis vectors

$$\begin{array}{}\text{(40)}& {\mathbb{E}}_D=\left(\begin{array}{ccc}f\frac{k_y}{k_{\mathrm{h}}}& \mathit{\omega }\frac{k_x}{k_{\mathrm{h}}}& f\frac{k_y}{k_{\mathrm{h}}}\\ -f\frac{k_x}{k_{\mathrm{h}}}& \mathit{\omega }\frac{k_y}{k_{\mathrm{h}}}& -f\frac{k_x}{k_{\mathrm{h}}}\\ \mathrm{0}& -\mathit{\omega }\frac{k_{\mathrm{h}}}{k_z}& \mathrm{0}\\ N\frac{k_{\mathrm{h}}}{k_z}& \mathrm{0}& -\frac{f^{\mathrm{2}}}{N}\frac{k_z}{k_{\mathrm{h}}}\end{array}\right).\end{array}$$

For any vector in the subspace formed by D, the vertical acceleration to the horizontal counterpart is always at least first order O(γ), as the first two components of ${\stackrel{\mathrm{̃}}{\mathit{e}}}_{\mathrm{2}}$ are orthogonal to the first two components of the two other basis vectors: a result which can be also seen directly from the divergence-free condition.

The validity of the hydrostatic approximation relies also on how strongly the modes ${\stackrel{\mathrm{̃}}{\mathit{e}}}_{\mathrm{1}},{\stackrel{\mathrm{̃}}{\mathit{e}}}_{\mathrm{2}}$ on one side and e₃ on the other get excited. If all the energy resides in the third mode, the hydrostatic approximation is exact, while the hydrostatic solution can diverge at a rate γ for the energy in the wave modes. For a comparison of the amplitudes of the forcing for each vector, the vectors are normalized by the square root of twice the energy density, kinetic plus potential. The modules are

$$\begin{array}{}\text{(41)}& {\mathbb{e}}_{\mathrm{1}}={\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }{\mathbb{e}}_{\mathrm{2}}=\sqrt{{\displaystyle \frac{\mathit{\alpha }}{k_z^{\mathrm{2}}}}},{\mathbb{e}}_{\mathrm{3}}={\mathit{\gamma }}^{-\mathrm{1}}\sqrt{{\displaystyle \frac{f^{\mathrm{2}}}{N^{\mathrm{2}}}}}{\mathbb{e}}_{\mathrm{1}},\end{array}$$

and the energies of the first and the third mode are equal when the horizontal scale equals the Rossby radius of deformation. Note that in the hydrostatic approximation, the vertical velocity is not a dynamical variable and does not contribute to the kinetic energy and, therefore, the energy density of ${\stackrel{\mathrm{̃}}{\mathit{e}}}_{\mathrm{2}}$ in the Navier–Stokes formalism agrees to the energy density of ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\mathrm{H}}{\stackrel{\mathrm{̃}}{\mathit{e}}}_{\mathrm{2}}$ in the hydrostatic formalism. For small aspect ratios the energy of the balanced eigenvector grows as γ⁻¹ with respect to the wave modes.

A vector that represents an equilibrium, geostrophy in our case, is in the kernel of both operators, as its time evolution vanishes. A conserved scalar quantity that is a linear combination of the components of the time evolution vector, divergence and potential vorticity in our case, is represented by a vector that is in the kernel of the adjoint of the operators, as its dot product with the time evolution vector vanishes. The adjoint operators are given by the transposed matrices $({\stackrel{\mathrm{̃}}{\mathbb{L}}}^{\mathrm{N}}{)}^t$ and $({\stackrel{\mathrm{̃}}{\mathbb{L}}}^{\mathrm{H}}{)}^t$ for Navier–Stokes and the hydrostatic case, respectively. Note that neither of the operators is self-adjoint – that is, the matrices are not symmetric. The kernel of the adjoint operators contains the wave vector ${\mathbb{r}}_{\mathrm{3}}=(\mathit{k},\mathrm{0}{)}^t$, which is therefore a conserved quantity. This shows that the acceleration is divergence free in both formalisms, the divergence of the velocity field is conserved and the vanishing divergence is maintained during the evolution of the velocity field. When viscosity and diffusivity vanish, the vector ${\mathbb{r}}_{\mathrm{4}}=(k_y,-k_x,\mathrm{0},-f{k}_z{N}^{-\mathrm{1}}{)}^t$ is also in the kernel of the adjoint hydrostatic operator; this corresponds to the conservation of linear potential vorticity, whereas in the Navier–Stokes case we have $({\stackrel{\mathrm{̃}}{\mathbb{L}}}^{\mathrm{N}}{)}^t{\mathbb{r}}_{\mathrm{4}}=f{\mathbb{r}}_{\mathrm{3}}$, which corresponds to the conservation of linear potential vorticity when the velocity field is divergence-free. As the velocity field is divergence free, linear potential vorticity is conserved in both formalisms. It is interesting to note that for the Navier–Stokes formalism, the kernel is one-dimensional as the sub-matrix corresponding to the zero eigenvalue has a Jordan normal form. The discussion on the adjoint operators also shows that there are no more than two conserved quantities in the linear approximations.

6 Projection on the basis of eigenvectors (S1)

The projection of a given acceleration $\stackrel{\mathrm{̃}}{\mathit{a}}$ on the modes of the linear dynamics represents step S1 as defined in the introduction and shown in Fig. 3. The identities

$$\begin{array}{}\text{(42)}& {\stackrel{\mathrm{̃}}{\mathit{e}}}_{\mathrm{1}}-{\mathit{e}}_{\mathrm{3}}={\displaystyle \frac{\mathit{\alpha }}{N{k}_z{k}_{\mathrm{h}}}}\left(\begin{array}{c}\mathrm{0}\\ \mathrm{0}\\ \mathrm{0}\\ \mathrm{1}\end{array}\right),f^{\mathrm{2}}{k}_z^{\mathrm{2}}{\stackrel{\mathrm{̃}}{\mathit{e}}}_{\mathrm{1}}+N^{\mathrm{2}}{k}_{\mathrm{h}}^{\mathrm{2}}{\mathit{e}}_{\mathrm{3}}={\displaystyle \frac{f\mathit{\alpha }}{k_{\mathrm{h}}}}\left(\begin{array}{c}{k}_y\\ -k_x\\ \mathrm{0}\\ \mathrm{0}\end{array}\right)\end{array}$$

allow us to write the extended projection operators that convert the forcing $\widehat{\mathit{F}}$ to a divergenceless acceleration in the first three components and an unchanged buoyancy flux in the fourth. The projection operators given in Eqs. (28), (29) that project the generalized acceleration $\stackrel{\mathrm{̃}}{\mathit{a}}$ written in terms of the basis D^□ are

$$\begin{array}{}\text{(43)}& {\displaystyle }\begin{array}{rl}& D{\stackrel{\mathrm{̃}}{\mathbb{P}}}^{\mathrm{N}}={\displaystyle \frac{\mathrm{1}}{k_{\mathrm{h}}}}\left[{\displaystyle \frac{f{k}_z^{\mathrm{2}}}{\mathit{\alpha }}}{\stackrel{\mathrm{̃}}{\mathit{e}}}_{\mathrm{1}}+{\displaystyle \frac{N^{\mathrm{2}}{k}_{\mathrm{h}}^{\mathrm{2}}}{f\mathit{\alpha }}}{\mathit{e}}_{\mathrm{3}}\right](k_y,-k_x,\mathrm{0},\mathrm{0},\mathrm{0})\\ & +{\displaystyle \frac{k_z^{\mathrm{2}}}{\mathit{\omega }{k}^{\mathrm{2}}{k}_{\mathrm{h}}}}{\stackrel{\mathrm{̃}}{\mathit{e}}}_{\mathrm{2}}(k_x,k_y,-{\displaystyle \frac{k_{\mathrm{h}}^{\mathrm{2}}}{k_z}},\mathrm{0},\mathrm{0})\\ & +k_z{k}_{\mathrm{h}}\left[{\displaystyle \frac{N}{\mathit{\alpha }}}({\stackrel{\mathrm{̃}}{\mathit{e}}}_{\mathrm{1}}-{\mathit{e}}_{\mathrm{3}})-{\displaystyle \frac{\mathrm{1}}{\mathit{\omega }{k}^{\mathrm{2}}}}{\stackrel{\mathrm{̃}}{\mathit{e}}}_{\mathrm{2}}\right](\mathrm{0},\mathrm{0},\mathrm{0},\mathrm{1},\mathrm{0})\\ & +{\displaystyle \frac{N{k}_z{k}_{\mathrm{h}}}{\mathit{\alpha }}}({\stackrel{\mathrm{̃}}{\mathit{e}}}_{\mathrm{1}}-{\mathit{e}}_{\mathrm{3}})(\mathrm{0},\mathrm{0},\mathrm{0},\mathrm{0},\mathrm{1}),\end{array}\text{(44)}& {\displaystyle }\begin{array}{rl}& D{\stackrel{\mathrm{̃}}{\mathbb{P}}}^{\mathrm{H}}={\displaystyle \frac{\mathrm{1}}{k_{\mathrm{h}}}}\left[{\displaystyle \frac{f{k}_z^{\mathrm{2}}}{\mathit{\alpha }}}{\stackrel{\mathrm{̃}}{\mathit{e}}}_{\mathrm{1}}+{\displaystyle \frac{N^{\mathrm{2}}{k}_{\mathrm{h}}^{\mathrm{2}}}{f\mathit{\alpha }}}{\mathit{e}}_{\mathrm{3}}\right](k_y,-k_x,\mathrm{0},\mathrm{0},\mathrm{0})\\ & +{\displaystyle \frac{\mathrm{1}}{\mathit{\omega }{k}_{\mathrm{h}}}}{\stackrel{\mathrm{̃}}{\mathit{e}}}_{\mathrm{2}}(k_x,k_y,\mathrm{0},\mathrm{0},\mathrm{0})\\ & +k_z{k}_{\mathrm{h}}\left[{\displaystyle \frac{N}{\mathit{\alpha }}}({\stackrel{\mathrm{̃}}{\mathit{e}}}_{\mathrm{1}}-{\mathit{e}}_{\mathrm{3}})-{\displaystyle \frac{\mathrm{1}}{\mathit{\omega }{k}_z^{\mathrm{2}}}}{\stackrel{\mathrm{̃}}{\mathit{e}}}_{\mathrm{2}}\right](\mathrm{0},\mathrm{0},\mathrm{0},\mathrm{1},\mathrm{0})\\ & +{\displaystyle \frac{N{k}_z{k}_{\mathrm{h}}}{\mathit{\alpha }}}({\stackrel{\mathrm{̃}}{\mathit{e}}}_{\mathrm{1}}-{\mathit{e}}_{\mathrm{3}})(\mathrm{0},\mathrm{0},\mathrm{0},\mathrm{0},\mathrm{1}).\end{array}\end{array}$$

In matrix form this is

$$\begin{array}{}\text{(45)}& {\displaystyle }\begin{array}{rl}& D{\stackrel{\mathrm{̃}}{\mathbb{P}}}^{\mathrm{N}}={\displaystyle \frac{\mathrm{1}}{\mathit{\alpha }{k}_{\mathrm{h}}}}\\ & \left(\begin{array}{ccccc}f{k}_z^{\mathrm{2}}{k}_y& -f{k}_z^{\mathrm{2}}{k}_x& \mathrm{0}& N{k}_z{k}_{\mathrm{h}}^{\mathrm{2}}& N{k}_z{k}_{\mathrm{h}}^{\mathrm{2}}\\ \frac{\mathit{\alpha }{k}_x{k}_z^{\mathrm{2}}}{\mathit{\omega }{k}^{\mathrm{2}}}& \frac{\mathit{\alpha }{k}_y{k}_z^{\mathrm{2}}}{\mathit{\omega }{k}^{\mathrm{2}}}& -\frac{\mathit{\alpha }{k}_{\mathrm{h}}^{\mathrm{2}}{k}_z}{\mathit{\omega }{k}^{\mathrm{2}}}& -\frac{\mathit{\alpha }{k}_{\mathrm{h}}^{\mathrm{2}}{k}_z}{\mathit{\omega }{k}^{\mathrm{2}}}& \mathrm{0}\\ \frac{N^{\mathrm{2}}{k}_{\mathrm{h}}^{\mathrm{2}}{k}_y}{f}& -\frac{N^{\mathrm{2}}{k}_{\mathrm{h}}^{\mathrm{2}}{k}_x}{f}& \mathrm{0}& -N{k}_z{k}_{\mathrm{h}}^{\mathrm{2}}& -N{k}_z{k}_{\mathrm{h}}^{\mathrm{2}}\end{array}\right),\end{array}\text{(46)}& {\displaystyle }\begin{array}{rl}& D{\stackrel{\mathrm{̃}}{\mathbb{P}}}^{\mathrm{H}}={\displaystyle \frac{\mathrm{1}}{\mathit{\alpha }{k}_{\mathrm{h}}}}\\ & \left(\begin{array}{ccccc}f{k}_z^{\mathrm{2}}{k}_y& -f{k}_z^{\mathrm{2}}{k}_x& \mathrm{0}& N{k}_z{k}_{\mathrm{h}}^{\mathrm{2}}& N{k}_z{k}_{\mathrm{h}}^{\mathrm{2}}\\ \frac{\mathit{\alpha }{k}_x}{\mathit{\omega }}& \frac{\mathit{\alpha }{k}_y}{\mathit{\omega }}& \mathrm{0}& -\frac{\mathit{\alpha }{k}_{\mathrm{h}}^{\mathrm{2}}}{\mathit{\omega }{k}_z}& \mathrm{0}\\ \frac{N^{\mathrm{2}}{k}_{\mathrm{h}}^{\mathrm{2}}{k}_y}{f}& -\frac{N^{\mathrm{2}}{k}_{\mathrm{h}}^{\mathrm{2}}{k}_x}{f}& \mathrm{0}& -N{k}_z{k}_{\mathrm{h}}^{\mathrm{2}}& -N{k}_z{k}_{\mathrm{h}}^{\mathrm{2}}\end{array}\right).\end{array}\end{array}$$

The operators only differ in the second line as all the differences are restricted to the second eigenvector. The projection operators are the left pseudo-inverse of 𝔼_D – that is, $D{\stackrel{\mathrm{̃}}{\mathbb{P}}}^{\square }{\mathbb{E}}_D={\mathbb{1}}_D$, which means that the vector space spanned by the basis $D=({\stackrel{\mathrm{̃}}{\mathit{e}}}_{\mathrm{1}},{\stackrel{\mathrm{̃}}{\mathit{e}}}_{\mathrm{2}},{\mathit{e}}_{\mathrm{3}})$ is composed of divergenceless vectors.

6.1 Forcing of the linear system

An exterior force $\widehat{\mathit{F}}=({\widehat{F}}_x,{\widehat{F}}_y,{\widehat{F}}_z,\mathrm{0},\widehat{B}{)}^t$ was introduced in Sect. 5. The exterior forcing is projected on the basis D by the projector $D{\stackrel{\mathrm{̃}}{\mathbb{P}}}^{\square }$. I then determine the amplitude of the three basis vectors subject to different forcings. To this end, I write the forcing vector as

$$\begin{array}{}\text{(47)}& \widehat{\mathit{F}}=\left(\begin{array}{c}\frac{k_x{\widehat{F}}_{\mathrm{d}}-k_y{\widehat{F}}_{\mathrm{r}}}{k_{\mathrm{h}}}\\ \frac{k_y{\widehat{F}}_{\mathrm{d}}+k_x{\widehat{F}}_{\mathrm{r}}}{k_{\mathrm{h}}}\\ {\widehat{F}}_z\\ \mathrm{0}\\ \widehat{B}\end{array}\right),\end{array}$$

where I split the horizontal forcing in a divergent and a rotational part, using the Helmholtz decomposition (${\widehat{F}}_{\mathrm{d}}=(k_x{\widehat{F}}_x+k_y{\widehat{F}}_y)/k_{\mathrm{h}}$ and ${\widehat{F}}_{\mathrm{r}}=(k_x{\widehat{F}}_y-k_y{\widehat{F}}_x)/k_{\mathrm{h}}$), to obtain

$$\begin{array}{}\text{(48)}& {\displaystyle }\begin{array}{rl}& D{\stackrel{\mathrm{̃}}{\mathbb{P}}}^{\mathrm{N}}\widehat{\mathit{F}}=-{\widehat{F}}_{\mathrm{r}}{\displaystyle \frac{\mathrm{1}}{f\mathit{\alpha }}}\left(\begin{array}{c}{f}^{\mathrm{2}}{k}_z^{\mathrm{2}}\\ \mathrm{0}\\ N^{\mathrm{2}}{k}_{\mathrm{h}}^{\mathrm{2}}\end{array}\right)\\ & +\widehat{B}{\displaystyle \frac{N{k}_{\mathrm{h}}{k}_z}{\mathit{\alpha }}}\left(\begin{array}{c}\mathrm{1}\\ \mathrm{0}\\ -\mathrm{1}\end{array}\right)+{\displaystyle \frac{k_z}{\mathit{\omega }{k}^{\mathrm{2}}}}({\widehat{F}}_{\mathrm{d}}{k}_z-{\widehat{F}}_z{k}_{\mathrm{h}})\left(\begin{array}{c}\mathrm{0}\\ \mathrm{1}\\ \mathrm{0}\end{array}\right),\end{array}\text{(49)}& {\displaystyle }\begin{array}{rl}& D{\stackrel{\mathrm{̃}}{\mathbb{P}}}^{\mathrm{H}}\widehat{\mathit{F}}=-{\widehat{F}}_{\mathrm{r}}{\displaystyle \frac{\mathrm{1}}{f\mathit{\alpha }}}\left(\begin{array}{c}{f}^{\mathrm{2}}{k}_z^{\mathrm{2}}\\ \mathrm{0}\\ N^{\mathrm{2}}{k}_{\mathrm{h}}^{\mathrm{2}}\end{array}\right)+\widehat{B}{\displaystyle \frac{N{k}_{\mathrm{h}}{k}_z}{\mathit{\alpha }}}\left(\begin{array}{c}\mathrm{1}\\ \mathrm{0}\\ -\mathrm{1}\end{array}\right)\\ & +{\displaystyle \frac{{\widehat{F}}_{\mathrm{d}}}{\mathit{\omega }}}\left(\begin{array}{c}\mathrm{0}\\ \mathrm{1}\\ \mathrm{0}\end{array}\right)=D{\stackrel{\mathrm{̃}}{\mathbb{P}}}^{\mathrm{N}}\widehat{\mathit{F}}+{\displaystyle \frac{k_{\mathrm{h}}}{\mathit{\omega }{k}^{\mathrm{2}}}}({\widehat{F}}_{\mathrm{d}}{k}_{\mathrm{h}}+{\widehat{F}}_z{k}_z)\left(\begin{array}{c}\mathrm{0}\\ \mathrm{1}\\ \mathrm{0}\end{array}\right).\end{array}\end{array}$$

The difference between the two operators is twofold: ${\widehat{F}}_z$ is ignored in the hydrostatic projection, and for the projection of ${\widehat{F}}_{\mathrm{d}}$ it is of second order, O(γ²). Both projections are identical if rotational and buoyant forcings are considered. All the differences are only in the projections on the second eigenvector. Therefore, the effect of all forcings on the balanced mode, which is the third component in the above equations, is identical in both formalisms. That is, if forcings act differently in the hydrostatic approximation, this difference does not directly affect the balanced mode. If the buoyancy frequency vanishes, N→0, all the energy from a rotational forcing goes to the first wave mode as no geostrophy is possible in this case. Furthermore, a divergent or vertical forcing does not act on the balanced mode. For a vertical forcing, all the work is done on the wave part in the Navier–Stokes formalism, whereas it has no effect in the hydrostatic approximation. If ${\widehat{F}}_{\mathrm{r}}f{k}_z=\widehat{B}N{k}_{\mathrm{h}}$, the forcing is balanced and all the work is done on the balanced mode.

6.2 Linearized equation in the eigenspace

As the fourth vector, in both formalisms, is in the kernel of the corresponding operator, the problem simplifies when it is restricted to the linear operator in the basis $D=({\stackrel{\mathrm{̃}}{\mathit{e}}}_{\mathrm{1}},{\stackrel{\mathrm{̃}}{\mathit{e}}}_{\mathrm{2}},{\mathit{e}}_{\mathrm{3}})$, which is common to both operators,

$$\begin{array}{}\text{(50)}& \begin{array}{rl}& D{\stackrel{\mathrm{̃}}{\mathbb{L}}}_D^{\square }=\left(\begin{array}{ccc}-\stackrel{\mathrm{˘}}{\mathit{\nu }}& \mathit{\omega }& \mathrm{0}\\ -({\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }{)}^{\mathrm{2}}\mathit{\omega }& -\stackrel{\mathrm{˘}}{\mathit{\nu }}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& -\stackrel{\mathrm{˘}}{\mathit{\nu }}\end{array}\right)\mathrm{and}\\ & D({\stackrel{\mathrm{̃}}{\mathbb{L}}}^{\square }{)}_D^{-\mathrm{1}}={\mathit{\beta }}^{\square }\left(\begin{array}{ccc}-\stackrel{\mathrm{˘}}{\mathit{\nu }}& -\mathit{\omega }& \mathrm{0}\\ ({\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }{)}^{\mathrm{2}}\mathit{\omega }& -\stackrel{\mathrm{˘}}{\mathit{\nu }}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& -(\stackrel{\mathrm{˘}}{\mathit{\nu }}{\mathit{\beta }}^{\square }{)}^{-\mathrm{1}}\end{array}\right),\end{array}\end{array}$$

where ${\mathit{\beta }}^{\square }=({\stackrel{\mathrm{˘}}{\mathit{\nu }}}^{\mathrm{2}}+({\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }\mathit{\omega }{)}^{\mathrm{2}}{)}^{-\mathrm{1}}$. The linearized dynamics consists in a damped thermal wind balance along e₃, for which the hydrostatic approximation is exact, and damped inertia-gravity waves around this balanced state along the directions ${\stackrel{\mathrm{̃}}{\mathit{e}}}_{\mathrm{1}}$ and ${\stackrel{\mathrm{̃}}{\mathit{e}}}_{\mathrm{2}}$. The only difference between the Navier–Stokes and the hydrostatic approximation in the subspace spanned by $\mathit{\{}{\stackrel{\mathrm{̃}}{\mathit{e}}}_{\mathrm{1}},{\stackrel{\mathrm{̃}}{\mathit{e}}}_{\mathrm{2}}\mathit{\}}$ lies in the frequencies, given in Eq. (37). Note that the correction ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }$ appears only as a square, and the approximations for the operator and the inverse operator are second order: O(γ²).

7 Forced motion (S2)

7.1 Slow forced motion, Ekman dynamics

This section discusses the special case of the response of the linearized Navier–Stokes and hydrostatic operators subject to an almost time-independent forcing. In geophysical fluid dynamics, the large-scale processes are evolving and interacting on a timescale much slower than the wave frequencies ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }\mathit{\omega }$. In this case, the left-hand side of Eq. (34) can be neglected. In the linear approximation, an external (boundary) forcing can be balanced by a viscous force, the Coriolis force or both as in Ekman dynamics. Boundary conditions are local in physical space and therefore non-local in Fourier space – that is, they are applied to different modes, which also interact through the boundary condition (see, e.g., Wirth, 2005). If the horizontal forcing is divergent, so is the horizontal velocity field, and vertical velocities are induced, which represent Ekman pumping. If we suppose that the forcing $\widehat{\mathit{F}}$ is weak, the linearized evolution in Eq. (34) is suitable. Using the Boussinesq approximation, the forcing is represented as an acceleration. If the forcing is time-independent, the stationary velocity field is restricted to the basis D and can be obtained by performing steps S1 and S2a through in Eqs. (34), (45), (46) and (50):

$$\begin{array}{}\text{(51)}& {\widehat{\mathit{x}}}^{\square }=-D({\stackrel{\mathrm{̃}}{\mathbb{L}}}^{\square }{)}_D^{-\mathrm{1}}D{\stackrel{\mathrm{̃}}{\mathbb{P}}}^{\square }\widehat{\mathit{F}}.\end{array}$$

The image of the operator $D{\stackrel{\mathrm{̃}}{\mathbb{P}}}^{\square }$ is spanned by the basis D for which $D({\stackrel{\mathrm{̃}}{\mathbb{L}}}^{\square }{)}_D^{-\mathrm{1}}$ exists if viscosity is non-vanishing.

7.1.1 Absolute error

After having determined the amplitude of each mode due to the forcing (S1 + S2a), the difference between the two formalisms in the basis D is determined. The total error is

$$\begin{array}{}\text{(52)}& \begin{array}{rl}& {\widehat{E}}_{\mathrm{0}}={\widehat{\mathit{x}}}^{\mathrm{H}}-{\widehat{\mathit{x}}}^{\mathrm{N}}=\left[-D({\stackrel{\mathrm{̃}}{\mathbb{L}}}^{\mathrm{H}}{)}_D^{-\mathrm{1}}D{\stackrel{\mathrm{̃}}{\mathbb{P}}}^{\mathrm{H}}+D({\stackrel{\mathrm{̃}}{\mathbb{L}}}^{\mathrm{N}}{)}_D^{-\mathrm{1}}D{\stackrel{\mathrm{̃}}{\mathbb{P}}}^{\mathrm{N}}\right]\\ & \widehat{\mathit{F}}=\left[-D({\stackrel{\mathrm{̃}}{\mathbb{L}}}^{\mathrm{H}}{)}_D^{-\mathrm{1}}+D({\stackrel{\mathrm{̃}}{\mathbb{L}}}^{\mathrm{N}}{)}_D^{-\mathrm{1}}\right]D{\stackrel{\mathrm{̃}}{\mathbb{P}}}^{\mathrm{H}}\widehat{\mathit{F}}-D({\stackrel{\mathrm{̃}}{\mathbb{L}}}^{\mathrm{N}}{)}_D^{-\mathrm{1}}\\ & \left[D{\stackrel{\mathrm{̃}}{\mathbb{P}}}^{\mathrm{H}}-D{\stackrel{\mathrm{̃}}{\mathbb{P}}}^{\mathrm{N}}\right]\widehat{\mathit{F}}\\ & ={\displaystyle \frac{\mathit{\omega }}{\mathit{\alpha }}}\left[(\mathit{\beta }-{\mathit{\beta }}^{\mathrm{H}})\left(\begin{array}{c}-\stackrel{\mathrm{˘}}{\mathit{\nu }}/\mathit{\omega }\\ \mathrm{1}\\ \mathrm{0}\end{array}\right)-{\mathit{\gamma }}^{\mathrm{2}}{\mathit{\beta }}^{\mathrm{H}}\left(\begin{array}{c}\mathrm{0}\\ \mathrm{1}\\ \mathrm{0}\end{array}\right)\right]\\ & \left[f{k}_z^{\mathrm{2}}{\widehat{F}}_{\mathrm{r}}+N{k}_{\mathrm{h}}{k}_z\widehat{B}\right]+\left(\begin{array}{c}\mathrm{1}\\ \stackrel{\mathrm{˘}}{\mathit{\nu }}/\mathit{\omega }\\ \mathrm{0}\end{array}\right)\\ & \left[({\mathit{\beta }}^{\mathrm{H}}-({\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\mathrm{H}}{)}^{-\mathrm{2}}\mathit{\beta }){\widehat{F}}_{\mathrm{d}}+\mathit{\beta }{\displaystyle \frac{k_z{k}_{\mathrm{h}}}{k^{\mathrm{2}}}}{\widehat{F}}_z\right].\end{array}\end{array}$$

We immediately see that the response to a stationary forcing of the balanced dynamics (e₃) is identical in both formalisms. Furthermore, a balanced forcing with $-f{k}_z{\widehat{F}}_{\mathrm{r}}+N{k}_{\mathrm{h}}\widehat{B}={\widehat{F}}_{\mathrm{d}}={\widehat{F}}_z=\mathrm{0}$ leads to a vanishing error. The last equality in the first line shows that the error can be separated into an evolution error, due to the differences in the wave frequencies (S2), and the error done by the hydrostatic projection (S1). The forcing F_z is in the kernel of the hydrostatic projection (S1) only. This and the subsequent evolution (S2) of the projection error gives the error shown by the second part of the last term of the second line. The order can be written as $O(\mathit{\gamma }{\widehat{F}}_z+{\mathit{\gamma }}^{\mathrm{2}}{\widehat{F}}_{\mathrm{d}}+{\mathit{\gamma }}^{\mathrm{2}}{\widehat{F}}_{\mathrm{r}}+{\mathit{\gamma }}^{\mathrm{3}}\widehat{B})$. If the vertical component of the forcing ${\widehat{F}}_z$ itself is first order, due to a boundary with a slope ≈γ, the approximation becomes second order, which shows that the hydrostatic approximation has two ingredients: the horizontal scales are larger than the vertical and the vertical forcing is smaller than the horizontal.

For a vanishing viscosity $\stackrel{\mathrm{˘}}{\mathit{\nu }}\to \mathrm{0}$, the error associated with a divergent forcing ${\widehat{F}}_{\mathrm{d}}$ is of second order $O\left({\stackrel{\mathrm{˘}}{\mathit{\nu }}}^{\mathrm{2}}\right)$, while it is of zero order for the other forcings. Furthermore, a vertical forcing leads to a velocity of ${\widehat{F}}_z\frac{k_{\mathrm{h}}{k}_z}{{\mathit{\omega }}^{\mathrm{2}}{k}^{\mathrm{2}}}{\stackrel{\mathrm{̃}}{\mathit{e}}}_{\mathrm{1}}$ in the Navier–Stokes formalism, while it has no effect in the hydrostatic approximation. As the vertical component of ${\stackrel{\mathrm{̃}}{\mathit{e}}}_{\mathrm{1}}$ is vanishing, so is the error in the vertical component. For a strong viscous damping $\stackrel{\mathrm{˘}}{\mathit{\nu }}\to \mathrm{\infty }$, present in numerical models just above the grid scale, where non-hydrostatic effects are also important, we have $\mathit{\beta }=O\left({\stackrel{\mathrm{˘}}{\mathit{\nu }}}^{-\mathrm{2}}\right)$, and the error is of order $O({\stackrel{\mathrm{˘}}{\mathit{\nu }}}^{-\mathrm{1}}{\widehat{F}}_z+{\stackrel{\mathrm{˘}}{\mathit{\nu }}}^{-\mathrm{1}}{\widehat{F}}_{\mathrm{d}}+{\stackrel{\mathrm{˘}}{\mathit{\nu }}}^{-\mathrm{2}}{\widehat{F}}_{\mathrm{r}}+{\stackrel{\mathrm{˘}}{\mathit{\nu }}}^{-\mathrm{2}}\widehat{B})$. For scales $f\approx \stackrel{\mathrm{˘}}{\mathit{\nu }}$, Ekman layer dynamics occurs.

7.1.2 Relative error

We have seen above that the dynamics can be split into two parts: wave dynamics along the first two eigenvectors and balanced flow along the third eigenvector. As the latter is free of error, the fidelity of the hydrostatic approximation depends also on how the energy is distributed among the eigenmodes. This is achieved by considering the relative error,

$$\begin{array}{}\text{(53)}& E_{\mathrm{R}}={\displaystyle \frac{\left|\widehat{\mathit{E}}\right|}{\left|\widehat{\mathit{x}}\right|}},\end{array}$$

rather than the absolute error, with

$$\begin{array}{}\text{(54)}& \begin{array}{rl}& \widehat{\mathit{x}}=-D{{\stackrel{\mathrm{̃}}{\mathbb{L}}}^{\mathrm{N}}}_D^{-\mathrm{1}}D{\stackrel{\mathrm{̃}}{\mathbb{P}}}^{\mathrm{N}}\widehat{\mathit{F}}=\left(\begin{array}{c}-f^{\mathrm{2}}{k}_z^{\mathrm{2}}\stackrel{\mathrm{˘}}{\mathit{\nu }}\mathit{\beta }\\ f^{\mathrm{2}}{k}_z^{\mathrm{2}}\mathit{\omega }\mathit{\beta }\\ -N^{\mathrm{2}}{k}_{\mathrm{h}}^{\mathrm{2}}{\stackrel{\mathrm{˘}}{\mathit{\nu }}}^{-\mathrm{1}}\end{array}\right){\displaystyle \frac{\mathrm{1}}{f\mathit{\alpha }}}{\widehat{F}}_{\mathrm{r}}\\ & +\left(\begin{array}{c}-\stackrel{\mathrm{˘}}{\mathit{\nu }}\mathit{\beta }\\ \mathit{\omega }\mathit{\beta }\\ -{\stackrel{\mathrm{˘}}{\mathit{\nu }}}^{-\mathrm{1}}\end{array}\right){\displaystyle \frac{N{k}_{\mathrm{h}}{k}_z}{\mathit{\alpha }}}\widehat{B}+\left(\begin{array}{c}\mathrm{1}\\ \stackrel{\mathrm{˘}}{\mathit{\nu }}/\mathit{\omega }\\ \mathrm{0}\end{array}\right){\displaystyle \frac{k_z\mathit{\beta }}{k^{\mathrm{2}}}}({\widehat{F}}_{\mathrm{d}}{k}_z-{\widehat{F}}_z{k}_{\mathrm{h}}),\end{array}\end{array}$$

for the different types of forcing. Note that the absolute error is a vector while the relative error is a scalar. A constant forcing on the wave mode leads to a finite-amplitude motion, while in the balanced mode it is only countered by viscosity.

The balanced mode is excited by rotational and buoyancy forcing only, and its amplitude compared to the amplitude of the wave motion is proportional to ${\stackrel{\mathrm{˘}}{\mathit{\nu }}}^{-\mathrm{1}}$, as it is damped by viscosity only. Furthermore, Eq. (41) shows that the square root of the energy of the third mode, as compared to the wave modes, grows with decreasing aspect ratio as γ⁻¹. The hydrostatic approximation for the relative error is therefore $O\left(E_{\mathrm{R}}\right)=\stackrel{\mathrm{˘}}{\mathit{\nu }}\mathit{\gamma }O\left(\widehat{\mathit{E}}\right)$, smaller, if the third component in Eq. (54) does not vanish – that is if $N^{\mathrm{2}}{k}_{\mathrm{h}}{\widehat{F}}_{\mathrm{r}}\ne f{k}_z\widehat{B}$. This adds considerably to the validity of the hydrostatic approximation for scales much larger than the viscous cut-off. If the forcing is purely vertical, ${\widehat{F}}_z\ne \mathrm{0}$ and ${\widehat{F}}_{\mathrm{r}}=\widehat{B}={\widehat{F}}_{\mathrm{d}}=\mathrm{0}$, the relative error is unity.

7.2 Projection-evolution operator (S2b)

After having discussed the case of a stationary force, I now consider the general case of a force varying in time. The linearized equation governing the dynamics subject to a low-amplitude forcing $\widehat{\mathit{F}}\left(t\right)$ is given by Eq. (34). A forcing can only lead to divergenceless acceleration and has therefore to be projected into the subspace of divergenceless functions by the corresponding projector given in Eqs. (48) and (49). The dynamics takes place in the subspace spanned by D. Its solution, starting from the initial condition $\left(\widehat{\mathit{u}}\right(\mathrm{0}),b(\mathrm{0}){)}^t$, is

$$\begin{array}{}\text{(55)}& \begin{array}{rl}& \left(\begin{array}{c}{\widehat{\mathit{u}}}^{\square }\left(t\right)\\ b^{\square }\left(t\right)\end{array}\right)={\mathbb{E}}_D\exp \left(D{\stackrel{\mathrm{̃}}{\mathbb{L}}}_D^{\square }\right(t\left)\right)\\ & \left[\underset{\mathrm{0}}{\overset{t}{\int }}\exp \left(D{\stackrel{\mathrm{̃}}{\mathbb{L}}}_D^{\square }\right(-s\left)\right)D{\stackrel{\mathrm{̃}}{\mathbb{P}}}^{\square }\widehat{\mathit{F}}\left(s\right)\mathrm{d}s+D{\stackrel{\mathrm{̃}}{\mathbb{P}}}^{\square }\left(\begin{array}{c}\widehat{\mathit{u}}\left(\mathrm{0}\right)\\ b\left(\mathrm{0}\right)\end{array}\right)\right]\\ & ={\mathbb{E}}_D\underset{\mathrm{0}}{\overset{t}{\int }}\exp \left(D{\stackrel{\mathrm{̃}}{\mathbb{L}}}_D^{\square }\right(t-s\left)\right)D{\stackrel{\mathrm{̃}}{\mathbb{P}}}^{\square }\widehat{\mathit{F}}\left(s\right)\mathrm{d}s\\ & +{\mathbb{E}}_D\exp \left(D{\stackrel{\mathrm{̃}}{\mathbb{L}}}_D^{\square }\right(t\left)\right)D{\stackrel{\mathrm{̃}}{\mathbb{P}}}^{\square }\left(\begin{array}{c}\widehat{\mathit{u}}\left(\mathrm{0}\right)\\ b\left(\mathrm{0}\right)\end{array}\right),\end{array}\end{array}$$

with

$$\begin{array}{}\text{(56)}& \begin{array}{rl}& \exp \left(D{\stackrel{\mathrm{̃}}{\mathbb{L}}}_D^{\square }\right(t\left)\right)=\\ & e^{-\stackrel{\mathrm{˘}}{\mathit{\nu }}t}\left(\begin{array}{ccc}\cos \left({\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }\mathit{\omega }t\right)& \left({\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }{)}^{-\mathrm{1}}\sin \right({\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }\mathit{\omega }t)& \mathrm{0}\\ -{\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }\sin \left({\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }\mathit{\omega }t\right)& \cos \left({\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }\mathit{\omega }t\right)& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right),\end{array}\end{array}$$

where I used $\exp \left(D{\stackrel{\mathrm{̃}}{\mathbb{L}}}_D^{\square }\right(t\left)\right)\exp \left(D{\stackrel{\mathrm{̃}}{\mathbb{L}}}_D^{\square }\right(-s\left)\right)=\exp \left(D{\stackrel{\mathrm{̃}}{\mathbb{L}}}_D^{\square }\right(t-s\left)\right)$. Once the dynamics is projected on the basis D, the evolution is identical in both formalisms along e₃, and the oscillations between ${\stackrel{\mathrm{̃}}{\mathit{e}}}_{\mathrm{1}}$ and ${\stackrel{\mathrm{̃}}{\mathit{e}}}_{\mathrm{2}}$ have a difference in the frequencies, and the hydrostatic solution diverges as $\sin \left(\frac{\mathit{\omega }{\mathit{\gamma }}^{\mathrm{2}}}{\mathrm{2}}t\right)$, for a slowly varying forcing, if the forcing frequency ζ≪ω. If the forcing has frequencies equal or close to ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }\mathit{\omega }$, resonances can occur and the divergence can take place at a faster rate.

The unforced linear equations allow two types of motion: (decaying, if ν≠0) internal waves with a frequency faster than f and (decaying, if ν≠0) thermal wind equilibrium. This dynamics is described by the second term on the right-hand side of Eq. (55). The error done by the hydrostatic approximation in the basis D, starting from rest and subject to a time-dependent force $\widehat{\mathit{F}}\left(t\right)$ is

$$\begin{array}{}\text{(57)}& \begin{array}{rl}& \textcolor{red}{}\widehat{\mathit{F}}=\underset{\mathrm{0}}{\overset{t}{\int }}\left[\exp \left(D{\stackrel{\mathrm{̃}}{\mathbb{L}}}_D^{\mathrm{H}}\right(t-s\left)\right)D{\stackrel{\mathrm{̃}}{\mathbb{P}}}^{\mathrm{H}}\right.\\ & \left.-\exp \left(D{{\stackrel{\mathrm{̃}}{\mathbb{L}}}^{\mathrm{N}}}_D\right(t-s\left)\right)D{\stackrel{\mathrm{̃}}{\mathbb{P}}}^{\mathrm{N}}\right]\widehat{\mathit{F}}\left(s\right)\mathrm{d}s.\end{array}\end{array}$$

7.3 Periodic forcing

When solving linear equations, solutions can be considered independently for each wave vector k and frequency ζ. For a periodic forcing with frequency ζ, the evolution in Eq. (34) in the basis D is

$$\begin{array}{}\text{(58)}& {\partial }_t{\widehat{\mathit{x}}}^{\square }={\stackrel{\mathrm{̃}}{\mathbb{L}}}_D^{\square }{\widehat{\mathit{x}}}^{\square }+{\widehat{\mathit{F}}}_D\cos \left(\mathit{\zeta }t\right).\end{array}$$

The general solution is

$$\begin{array}{}\text{(59)}& \begin{array}{rl}& {\widehat{\mathit{x}}}^{\square }=\exp \left(D{\stackrel{\mathrm{̃}}{\mathbb{L}}}_D^{\square ,\mathit{\zeta }}\right(t\left)\right){\widehat{\mathit{F}}}_D\\ & =\underset{\mathrm{0}}{\overset{t}{\int }}\exp \left(D{\stackrel{\mathrm{̃}}{\mathbb{L}}}_D^{\square }\right(t-s\left)\right)\cos \left(\mathit{\zeta }s\right)\mathrm{d}s{\widehat{\mathit{F}}}_D.\end{array}\end{array}$$

Before discussing the general case, I examine the case of a very small viscosity (${\stackrel{\mathrm{˘}}{\mathit{\nu }}}^{\mathrm{2}}\ll ({\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }\mathit{\omega }{)}^{\mathrm{2}},{\mathit{\zeta }}^{\mathrm{2}}$). The solution is

$$\begin{array}{}\text{(60)}& \begin{array}{rl}& {\widehat{\mathit{x}}}^{\square }={\displaystyle \frac{\mathrm{1}}{{\mathit{\zeta }}^{\mathrm{2}}-({\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }\mathit{\omega }{)}^{\mathrm{2}}}}\\ & \left(\begin{array}{ccc}\mathit{\zeta }\sin \left(\mathit{\zeta }t\right)-a{\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }\mathit{\omega }\sin \left({\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }\mathit{\omega }t\right)& \mathit{\omega }(-\cos (\mathit{\zeta }t)+a\cos ({\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }\mathit{\omega }t\left)\right)& \mathrm{0}\\ \left({\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }{)}^{\mathrm{2}}\mathit{\omega }\right(\cos \left(\mathit{\zeta }t\right)-a\cos \left({\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }\mathit{\omega }t\right))& \mathit{\zeta }\sin \left(\mathit{\zeta }t\right)-a{\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }\mathit{\omega }\sin \left({\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }\mathit{\omega }t\right)& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \frac{{\mathit{\zeta }}^{\mathrm{2}}-({\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }\mathit{\omega }{)}^{\mathrm{2}}}{\mathit{\zeta }}\sin \left(\mathit{\zeta }t\right)\end{array}\right)\\ & {\widehat{\mathit{F}}}_D.\end{array}\end{array}$$

The constant a is there to satisfy the initial condition. For high-frequency forcing ${\mathit{\zeta }}^{\mathrm{2}}\gg ({\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }\mathit{\omega }{)}^{\mathrm{2}}$, there are no resonances and the part independent of the initial condition, which is “forgotten” after a long time $t\gg {\stackrel{\mathrm{˘}}{\mathit{\nu }}}^{-\mathrm{1}}$, is

$$\begin{array}{}\text{(61)}& {\widehat{\mathit{x}}}^{\square }=\exp \left(D{\stackrel{\mathrm{̃}}{\mathbb{L}}}_D^{\square ,\mathit{\zeta }}\right(t\left)\right){\widehat{\mathit{F}}}_D=\mathit{\zeta }\sin \left(\mathit{\zeta }t\right){\widehat{\mathit{F}}}_D.\end{array}$$

In this case, the operator ${\stackrel{\mathrm{̃}}{\mathbb{L}}}_D^{\square }$ has no time to act on the dynamics.

The general solution for $t\gg {\stackrel{\mathrm{˘}}{\mathit{\nu }}}^{-\mathrm{1}}$, when the initial condition is forgotten, is

$$\begin{array}{}\text{(62)}& {\widehat{\mathit{x}}}^{\square }=\left(\begin{array}{ccc}\mathit{\zeta }{\stackrel{\mathrm{̃}}{s}}^{\square }& -\mathit{\omega }{\stackrel{\mathrm{̃}}{c}}^{\square }& \mathrm{0}\\ ({\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }{)}^{\mathrm{2}}\mathit{\omega }{\stackrel{\mathrm{̃}}{c}}^{\square }& \mathit{\zeta }{\stackrel{\mathrm{̃}}{s}}^{\square }& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \stackrel{\mathrm{̃}}{b}\end{array}\right){\mathit{F}}_D,\end{array}$$

with

$$\begin{array}{}\text{(63)}& {\displaystyle }\begin{array}{rl}{\stackrel{\mathrm{̃}}{s}}^{\square }& ={\stackrel{\mathrm{̃}}{q}}^{\square }\left[{\displaystyle \frac{\stackrel{\mathrm{˘}}{\mathit{\nu }}}{\mathit{\zeta }}}({\stackrel{\mathrm{˘}}{\mathit{\nu }}}^{\mathrm{2}}+{\mathit{\zeta }}^{\mathrm{2}}+({\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }\mathit{\omega }{)}^{\mathrm{2}}\left)\cos \right(\mathit{\zeta }t)\right.\\ & \left.+({\stackrel{\mathrm{˘}}{\mathit{\nu }}}^{\mathrm{2}}+{\mathit{\zeta }}^{\mathrm{2}}-({\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }\mathit{\omega }{)}^{\mathrm{2}}\left)\sin \right(\mathit{\zeta }t)\right],\end{array}\text{(64)}& {\displaystyle }{\stackrel{\mathrm{̃}}{c}}^{\square }={\stackrel{\mathrm{̃}}{q}}^{\square }\left[-({\stackrel{\mathrm{˘}}{\mathit{\nu }}}^{\mathrm{2}}-{\mathit{\zeta }}^{\mathrm{2}}+({\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }\mathit{\omega }{)}^{\mathrm{2}}\left)\cos \right(\mathit{\zeta }t)-\mathrm{2}\stackrel{\mathrm{˘}}{\mathit{\nu }}\mathit{\zeta }\sin (\mathit{\zeta }t)\right],\text{(65)}& {\displaystyle }\begin{array}{rl}({\stackrel{\mathrm{̃}}{q}}^{\square }{)}^{-\mathrm{1}}& =({\stackrel{\mathrm{˘}}{\mathit{\nu }}}^{\mathrm{2}}+(\mathit{\zeta }+{\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }\mathit{\omega }{)}^{\mathrm{2}}\left)\right({\stackrel{\mathrm{˘}}{\mathit{\nu }}}^{\mathrm{2}}+(\mathit{\zeta }-{\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }\mathit{\omega }{)}^{\mathrm{2}})\\ & ={\stackrel{\mathrm{˘}}{\mathit{\nu }}}^{\mathrm{4}}+\mathrm{2}{\stackrel{\mathrm{˘}}{\mathit{\nu }}}^{\mathrm{2}}({\mathit{\zeta }}^{\mathrm{2}}+({\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }\mathit{\omega }{)}^{\mathrm{2}})+({\mathit{\zeta }}^{\mathrm{2}}-({\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }\mathit{\omega }{)}^{\mathrm{2}}{)}^{\mathrm{2}},\end{array}\text{(66)}& {\displaystyle }\stackrel{\mathrm{̃}}{b}={\displaystyle \frac{\stackrel{\mathrm{˘}}{\mathit{\nu }}\cos \left(\mathit{\zeta }t\right)+\mathit{\zeta }\sin \left(\mathit{\zeta }t\right)}{{\stackrel{\mathrm{˘}}{\mathit{\nu }}}^{\mathrm{2}}+{\mathit{\zeta }}^{\mathrm{2}}}}.\end{array}$$

For a stationary forcing, ζ=0, the results from Sect. 7.1 are obtained, and if ${\stackrel{\mathrm{˘}}{\mathit{\nu }}}^{\mathrm{2}}\ll ({\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }\mathit{\omega }{)}^{\mathrm{2}},{\mathit{\zeta }}^{\mathrm{2}}$, the solution given in Eq. (61) is recovered. The wave frequency only appears as a square, which assures that the error of the projection-evolution operator is of second order. The eigenvalues of the evolution matrix are

$$\begin{array}{}\text{(67)}& {\mathit{\lambda }}_{W\mathrm{1}/W\mathrm{2}}^{\square }=\mathit{\zeta }{\stackrel{\mathrm{̃}}{s}}^{\square }\pm i{\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }\mathit{\omega }{\stackrel{\mathrm{̃}}{c}}^{\square }\mathrm{and}{\mathit{\lambda }}_b=\stackrel{\mathrm{̃}}{b}.\end{array}$$

The corresponding eigenvectors are ${\mathit{e}}_{\mathrm{1},\mathrm{2}}^{\square }$ and e₃. Resonances are efficiently damped by viscosity if for a forcing at the resonant frequency, ζ²=ω², the viscous timescale is larger than the resonance period, ${\stackrel{\mathrm{˘}}{\mathit{\nu }}}^{\mathrm{2}}>{\mathit{\omega }}^{\mathrm{2}}$.

7.3.1 Absolute error

The difference between the two evolution operators in the case of a long viscous damping time ${\stackrel{\mathrm{˘}}{\mathit{\nu }}}^{\mathrm{2}}\ll {\mathit{\zeta }}^{\mathrm{2}},({\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }\mathit{\omega }{)}^{\mathrm{2}}$ is

$$\begin{array}{}\text{(68)}& \begin{array}{rl}{\mathbb{G}}_{\mathit{\zeta }}^{\mathrm{0}}& ={\displaystyle \frac{{\mathit{\gamma }}^{\mathrm{2}}}{({\mathit{\zeta }}^{\mathrm{2}}-({\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }{)}^{\mathrm{2}}{\mathit{\omega }}^{\mathrm{2}}\left)\right({\mathit{\zeta }}^{\mathrm{2}}-{\mathit{\omega }}^{\mathrm{2}})}}\\ & \left(\begin{array}{ccc}\mathit{\zeta }{\mathit{\omega }}^{\mathrm{2}}\sin \left(\mathit{\zeta }t\right)& -{\mathit{\omega }}^{\mathrm{3}}\cos \left(\mathit{\zeta }t\right)& \mathrm{0}\\ {\mathit{\zeta }}^{\mathrm{2}}\mathit{\omega }\cos \left(\mathit{\zeta }t\right)& \mathit{\zeta }{\mathit{\omega }}^{\mathrm{2}}\sin \left(\mathit{\zeta }t\right)& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right).\end{array}\end{array}$$

The hydrostatic approximation is of second order. In the general case, the difference between the two evolution operators is

$$\begin{array}{}\text{(69)}& {\mathbb{G}}_{\mathit{\zeta }}^{\mathit{\nu }}=\left(\begin{array}{ccc}\mathit{\zeta }({\stackrel{\mathrm{̃}}{s}}^{\mathrm{H}}-\stackrel{\mathrm{̃}}{s})& -\mathit{\omega }({\stackrel{\mathrm{̃}}{c}}^{\mathrm{H}}-\stackrel{\mathrm{̃}}{c})& \mathrm{0}\\ {\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\mathrm{2}}\mathit{\omega }{\stackrel{\mathrm{̃}}{c}}^{\mathrm{H}}-\mathit{\omega }\stackrel{\mathrm{̃}}{c}& \mathit{\zeta }({\stackrel{\mathrm{̃}}{s}}^{\mathrm{H}}-\stackrel{\mathrm{̃}}{s})& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right).\end{array}$$

All the components of the matrix are of at least second order, O(γ²).

The error done by the hydrostatic approximation with a periodic forcing of frequency ζ is

$$\begin{array}{}\text{(70)}& {\widehat{\mathit{E}}}_{\mathit{\zeta }}^{\mathit{\nu }}=\left[{\mathbb{G}}_{\mathit{\zeta }}^{\mathit{\nu }}D{\stackrel{\mathrm{̃}}{\mathbb{P}}}^{\mathrm{H}}+\exp \left(D{{\stackrel{\mathrm{̃}}{\mathbb{L}}}^{\mathrm{N}}}_D^{\mathit{\zeta }}\right(t\left)\right)(D{\stackrel{\mathrm{̃}}{\mathbb{P}}}^{\mathrm{H}}-D{\stackrel{\mathrm{̃}}{\mathbb{P}}}^{\mathrm{N}})\right]\widehat{\mathit{F}}.\end{array}$$

It has two terms: the first is due to the different evolution in the two formalisms (see Eq. 69) and the second to the difference in the projections (see Eq. 49). The explicit form is

$$\begin{array}{}\text{(71)}& \begin{array}{rl}& {\widehat{\mathit{E}}}_{\mathit{\zeta }}^{\mathit{\nu }}={\displaystyle \frac{k_z}{\mathit{\alpha }}}(-f{k}_z{\widehat{F}}_{\mathrm{r}}+k_{\mathrm{h}}\widehat{B})\left(\begin{array}{c}\mathit{\zeta }({\stackrel{\mathrm{̃}}{s}}^{\mathrm{H}}-\stackrel{\mathrm{̃}}{s})\\ ({\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }{)}^{\mathrm{2}}\mathit{\omega }{\stackrel{\mathrm{̃}}{c}}^{\mathrm{H}}-\mathit{\omega }\stackrel{\mathrm{̃}}{c}\\ \mathrm{0}\end{array}\right)\\ & +{\widehat{F}}_{\mathrm{d}}\left[\left(\begin{array}{c}-({\stackrel{\mathrm{̃}}{c}}^{\mathrm{H}}-\stackrel{\mathrm{̃}}{c})\\ \frac{\mathit{\zeta }}{\mathit{\omega }}({\stackrel{\mathrm{̃}}{s}}^{\mathrm{H}}-\stackrel{\mathrm{̃}}{s})\\ \mathrm{0}\end{array}\right)+{\displaystyle \frac{k_{\mathrm{h}}^{\mathrm{2}}}{k^{\mathrm{2}}}}\left(\begin{array}{c}-\stackrel{\mathrm{̃}}{c}\\ \frac{\mathit{\zeta }}{\mathit{\omega }}\stackrel{\mathrm{̃}}{s}\\ \mathrm{0}\end{array}\right)\right]-{\widehat{F}}_z{\displaystyle \frac{k_{\mathrm{h}}{k}_z}{k^{\mathrm{2}}}}\left(\begin{array}{c}-\stackrel{\mathrm{̃}}{c}\\ \frac{\mathit{\zeta }}{\mathit{\omega }}\stackrel{\mathrm{̃}}{s}\\ \mathrm{0}\end{array}\right).\end{array}\end{array}$$

For the case of a long viscous damping time I obtain

$$\begin{array}{}\text{(72)}& \begin{array}{rl}& {\widehat{\mathit{E}}}_{\mathit{\zeta }}^{\mathrm{0}}={\displaystyle \frac{{\mathit{\gamma }}^{\mathrm{2}}\mathit{\omega }\mathit{\zeta }}{({\mathit{\zeta }}^{\mathrm{2}}-({\mathit{\omega }}^{\mathrm{H}}{)}^{\mathrm{2}}\left)\right({\mathit{\zeta }}^{\mathrm{2}}-{\mathit{\omega }}^{\mathrm{2}})}}\left[{\displaystyle \frac{k_z}{\mathit{\alpha }}}(-f{k}_z{\widehat{F}}_{\mathrm{r}}+k_{\mathrm{h}}\widehat{B})\right]\\ & \left(\begin{array}{c}\mathit{\omega }\sin \left(\mathit{\zeta }t\right)\\ \mathit{\zeta }\cos \left(\mathit{\zeta }t\right)\\ \mathrm{0}\end{array}\right)+{\displaystyle \frac{\mathrm{1}}{{\mathit{\zeta }}^{\mathrm{2}}-{\mathit{\omega }}^{\mathrm{2}}}}\left[{\widehat{F}}_{\mathrm{d}}({\displaystyle \frac{{\mathit{\gamma }}^{\mathrm{2}}\mathit{\omega }}{{\mathit{\zeta }}^{\mathrm{2}}-({\mathit{\omega }}^{\mathrm{H}}{)}^{\mathrm{2}}}}+{\displaystyle \frac{k_{\mathrm{h}}^{\mathrm{2}}}{\mathit{\omega }{k}^{\mathrm{2}}}})\right.\\ & \left.-{\widehat{F}}_z{\displaystyle \frac{k_{\mathrm{h}}{k}_z}{\mathit{\omega }{k}^{\mathrm{2}}}}\right]\left(\begin{array}{c}-\mathit{\omega }\cos \left(\mathit{\zeta }t\right)\\ \mathit{\zeta }\sin \left(\mathit{\zeta }t\right)\\ \mathrm{0}\end{array}\right).\end{array}\end{array}$$

This shows, once more, that the third component, the balanced dynamics, is free of error. The scaling of the error with the aspect ratio is

$$\begin{array}{}\text{(73)}& O\left(\widehat{\mathit{E}}\right)=O({\mathit{\gamma }}^{\mathrm{2}}{\widehat{F}}_{\mathrm{r}}+{\mathit{\gamma }}^{\mathrm{3}}\widehat{B}+{\mathit{\gamma }}^{\mathrm{2}}{\widehat{F}}_{\mathrm{d}}+\mathit{\gamma }{\widehat{F}}_z).\end{array}$$

If the forcing evolves much slower than the buoyancy frequency, the dynamics is close to the stationary state calculated in Sect. 7.1. If the forcing occurs at the wave frequency of either formalism, resonances or near resonances dominate the linear dynamics. The hydrostatic approximation breaks down as ω^H≠ω if k_h≠0 and a resonance occurs in one formalism but not in the other. During resonances or close resonances, the scaling is still given by Eq. (73) but the prefactor diverges. Note that if ${\stackrel{\mathrm{˘}}{\mathit{\nu }}}^{\mathrm{2}}\gg ({\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }\mathit{\omega }{)}^{\mathrm{2}}$, resonances are efficiently damped.

7.3.2 Relative error

To obtain the relative error the procedure is identical to the case of stationary forcing considered in Sect. 7.1.2. We compare the absolute error calculated above to the solution of the Navier–Stokes model, using Eq. (62):

$$\begin{array}{}\text{(74)}& \begin{array}{rl}{\widehat{\mathit{x}}}^{\mathrm{N}}& ={\widehat{F}}_{\mathrm{r}}{\displaystyle \frac{\mathrm{1}}{f\mathit{\alpha }}}\left(\begin{array}{c}-f^{\mathrm{2}}{k}_z^{\mathrm{2}}\mathit{\zeta }{\stackrel{\mathrm{̃}}{s}}^{\mathrm{N}}\\ -f^{\mathrm{2}}{k}_z^{\mathrm{2}}{\mathit{\omega }}^{\mathrm{2}}{\stackrel{\mathrm{̃}}{c}}^{\mathrm{N}}\\ N^{\mathrm{2}}{k}_{\mathrm{h}}^{\mathrm{2}}\stackrel{\mathrm{̃}}{b}\end{array}\right)+\widehat{B}{\displaystyle \frac{k_{\mathrm{h}}{k}_z}{\mathit{\alpha }}}\left(\begin{array}{c}\mathit{\zeta }{\stackrel{\mathrm{̃}}{s}}^{\mathrm{N}}\\ {\mathit{\omega }}^{\mathrm{2}}{\stackrel{\mathrm{̃}}{c}}^{\mathrm{N}}\\ -\stackrel{\mathrm{̃}}{b}\end{array}\right)\\ & +({\widehat{F}}_{\mathrm{d}}{\displaystyle \frac{k_z^{\mathrm{2}}}{k^{\mathrm{2}}}}-{\widehat{F}}_z{\displaystyle \frac{k_z{k}_{\mathrm{h}}}{k^{\mathrm{2}}}})\left(\begin{array}{c}-{\stackrel{\mathrm{̃}}{c}}^{\mathrm{N}}\\ \mathit{\zeta }{\stackrel{\mathrm{̃}}{s}}^{\mathrm{N}}\\ \mathrm{0}\end{array}\right),\end{array}\end{array}$$

which is taken as the reference. Calculations confirm the results from Sect. 7.1.2 also for the case of periodic forcing, except that prefactors in the scaling can diverge due to resonances.

7.4 Nonlinear terms

Formally, the nonlinear terms can be added to the forcing in Eq. (55) by $\widehat{\mathit{F}}\to \widehat{\mathit{F}}+{\widehat{F}}_{\mathrm{nl}}$. Using Eqs. (15) to (18), I obtain

$$\begin{array}{}\text{(75)}& {\widehat{F}}_{\mathrm{nl}}=-i\left(\begin{array}{c}{k}_x\widehat{u_x^{\mathrm{2}}}+k_y\widehat{u_x{u}_y}+k_z\widehat{u_x{u}_z}\\ k_x\widehat{u_x{u}_y}+k_y\widehat{u_y^{\mathrm{2}}}+k_z\widehat{u_y{u}_z}\\ k_x\widehat{u_x{u}_z}+k_y\widehat{u_y{u}_z}+k_z\widehat{u_z^{\mathrm{2}}}\\ k_x\widehat{u_{x}b}+k_y\widehat{u_{y}b}+k_z\widehat{u_{z}b}\end{array}\right).\end{array}$$

Due to the zero divergence of the velocity field, the horizontal velocity components are at least a factor γ larger than the vertical velocity, $\widehat{u_x}\sim \widehat{u_y}\sim {\mathit{\gamma }}^{-\mathrm{1}}\widehat{u_z}$. If this property is conserved by the nonlinear term, the hydrostatic approximation is of second order O(γ²), following Eq. (73).

8 Discussion

The hydrostatic approximation is exact on the balanced mode, and if this mode receives all the energy by a given forcing then the approximation is exact. If the forcing strongly excites wave motion, then the deficiencies of the hydrostatic approximation in predicting the wave dynamics determines its overall quality. The hydrostatic approximation fails if the forcing is a vertical acceleration which is not a dynamical variable in the hydrostatic approximation and is therefore ignored, while in the Navier–Stokes formalism, it generates horizontal velocities of first order in the scale ratio, O(γ).

In my evaluation, done using the techniques T1, T2 and T3 in Fourier space, the calculation is done separately for every three-dimensional wave vector k. The projection step, S1, is instantaneous; it does not depend on the temporal variation of the forcing, whereas the evolution step, S2, asks for integrating the evolution equations in time. In the linear approximation, forcings at different frequencies superpose without interacting and step S2 can be calculated for each forcing frequency separately. The quality of the hydrostatic approximation depends on the Coriolis parameter f, the buoyancy frequency N and the viscosity ν. It has to be determined for every three-dimensional wave vector k, forcing vector $\widehat{\mathit{F}}$ and frequency ζ (see Fig. 3).

The mathematical formalism to evaluate the hydrostatic approximation was developed and applied in the previous sections. It is now important to compare the deficiencies of the hydrostatic approximation to other uncertainties, to evaluate its fidelity. The previous discussion shows that viscosity and damping can mask the differences in the wave frequencies between the two formalisms. This will be assessed in Sect. 8.1. The wave packets move at the group velocity and are advected by the mean flow, which shows internal variability. In Sect. 8.2, I compare this Doppler effect to the error due to the hydrostatic approximation. Every measurement or observation is limited in time and space. Section 8.3 applies the Heisenberg–Gabor uncertainty principle to discuss for which processes the difference between the two formalisms is actually detectable, given the finite space-time observations.

8.1 Masking of resonance differences by friction

I first note that the amplitude of the balanced mode is bounded by $max(\stackrel{\mathrm{˘}}{\mathit{\nu }},\mathit{\zeta }{)}^{-\mathrm{1}}$ times the forcing amplitude (see Eq. 66). Resonances occur when $\sqrt{({\mathit{\zeta }}^{\mathrm{2}}-({\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }\mathit{\omega }{)}^{\mathrm{2}}{)}^{\mathrm{2}}}\ll ({\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }\mathit{\omega }{)}^{\mathrm{2}}$. For not too small viscous damping, $({\mathit{\zeta }}^{\mathrm{2}}-({\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }\mathit{\omega }{)}^{\mathrm{2}}{)}^{\mathrm{2}}\ll \mathrm{2}({\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }\mathit{\omega }{)}^{\mathrm{2}}{\stackrel{\mathrm{˘}}{\mathit{\nu }}}^{\mathrm{2}}$, Eqs. (63) and (64) approach

$$\begin{array}{}\text{(76)}& {\displaystyle }{\stackrel{\mathrm{̃}}{s}}^{\square }\approx {\displaystyle \frac{\mathrm{1}}{\mathrm{2}{\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }\mathit{\omega }\stackrel{\mathrm{˘}}{\mathit{\nu }}}}\cos \left({\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }\mathit{\omega }t\right),\text{(77)}& {\displaystyle }{\stackrel{\mathrm{̃}}{c}}^{\square }\approx {\displaystyle \frac{-\mathrm{1}}{\mathrm{2}{\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }\mathit{\omega }\stackrel{\mathrm{˘}}{\mathit{\nu }}}}\sin \left({\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }\mathit{\omega }t\right).\end{array}$$

The amplitude of the wave motion is bounded by viscosity and diverges as ${\stackrel{\mathrm{˘}}{\mathit{\nu }}}^{-\mathrm{1}}$. There is no significant difference between the two formalisms. For smaller viscosities $({\mathit{\zeta }}^{\mathrm{2}}-({\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }\mathit{\omega }{)}^{\mathrm{2}}{)}^{\mathrm{2}}\gg \mathrm{2}({\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }\mathit{\omega }{)}^{\mathrm{2}}{\stackrel{\mathrm{˘}}{\mathit{\nu }}}^{\mathrm{2}}$:

$$\begin{array}{}\text{(78)}& {\displaystyle }{\stackrel{\mathrm{̃}}{s}}^{\square }\approx {\displaystyle \frac{\mathrm{2}{\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }\mathit{\omega }\stackrel{\mathrm{˘}}{\mathit{\nu }}}{({\mathit{\zeta }}^{\mathrm{2}}-({\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }\mathit{\omega }{)}^{\mathrm{2}}{)}^{\mathrm{2}}}}\cos \left({\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }\mathit{\omega }t\right),\text{(79)}& {\displaystyle }{\stackrel{\mathrm{̃}}{c}}^{\square }\approx {\displaystyle \frac{\mathrm{2}{\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }\mathit{\omega }\stackrel{\mathrm{˘}}{\mathit{\nu }}}{({\mathit{\zeta }}^{\mathrm{2}}-({\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }\mathit{\omega }{)}^{\mathrm{2}}{)}^{\mathrm{2}}}}\sin \left({\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }\mathit{\omega }t\right).\end{array}$$

The amplitude of the wave motion shows a quadratic divergence: $O\left(\mathit{\zeta }\pm \left({\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{\square }\mathit{\omega }\right){)}^{-\mathrm{2}}\right)$ and resonances can occur in one formalism and not in the other when

$$\begin{array}{}\text{(80)}& {\mathit{\gamma }}^{\mathrm{2}}\ll {\displaystyle \frac{{\stackrel{\mathrm{˘}}{\mathit{\nu }}}^{\mathrm{2}}}{{\mathit{\omega }}^{\mathrm{2}}}},\end{array}$$

which occurs at scales $k\gg k_{\mathrm{r}}=(\mathit{\gamma }\mathit{\omega }/\mathit{\nu }{)}^{\mathrm{1}/\mathrm{2}}$. If hyper- or hypo-dissipation operators ν_α∇^2α are used (hyper-viscosity, α=1 corresponding to Laplacian dissipation and α=0 to scale-independent bottom friction), this occurs at scales $k\gg k_{\mathrm{r}}=(\mathit{\gamma }\mathit{\omega }/{\mathit{\nu }}_{\mathit{\alpha }}{)}^{\mathrm{1}/\left(\mathrm{2}\mathit{\alpha }\right)}$.

In simulations of the ocean dynamics, we either resolve the turbulent mixed-layer dynamics at the surface and the bottom, which asks for a resolution on the meter scale in the vertical and the horizontal, or we parameterize it, which asks for vertical viscosity above ${\mathit{\nu }}_z>{\mathrm{10}}^{-\mathrm{2}}$ m² s⁻¹ (supposing an Ekman layer thickness of a few tens of meters). In numerical models of the ocean dynamics, resolution and viscous damping are very different in the horizontal and the vertical. If the model is fine resolution, an isotropic eddy viscosity is chosen based on the Smagorinsky large eddy simulation model. Its viscous damping time is $\stackrel{\mathrm{˘}}{\mathit{\nu }}=C_{\mathrm{s}}^{\mathrm{2}}u/\mathrm{\Delta }x$. A Smagorinsky constant C_s=0.16, a current speed of 10⁻¹ m s⁻¹ and a grid scale Δx lead to $\stackrel{\mathrm{˘}}{\mathit{\nu }}\approx \mathrm{3}\cdot {\mathrm{10}}^{-\mathrm{3}}$ m s${}^{-\mathrm{1}}/\left(\mathrm{\Delta }x\right)$. If we suppose that $\mathit{\omega }>{\mathrm{10}}^{-\mathrm{4}}$ s⁻¹, Eq. (80) gives $\mathit{\gamma }\ll \frac{\mathrm{30}\mathrm{m}}{\mathrm{\Delta }x}$, meaning that the difference between the resonant frequencies in both formalisms is hidden by viscosity if the vertical scale is below 30 m, which is a typical scale for the oceanic Ekman layer. This shows that viscosity masks differences between the two formalisms in the bottom and surface Ekman layer dynamics. At larger horizontal scales friction is parameterized through turbulent bottom friction, which has an inverse friction time of $c_{\mathrm{D}}u/H$ A bottom friction drag coefficient of $c_{\mathrm{D}}\approx {\mathrm{10}}^{-\mathrm{3}}$, a current speed of $u={\mathrm{10}}^{-\mathrm{1}}$ m s⁻¹ and a resonance frequency $\mathit{\omega }\approx {\mathrm{10}}^{-\mathrm{4}}$ s⁻¹ show that bottom friction is sufficient to hide the difference in resonance frequencies (γω) only if the ocean layer thickness in meters is below γ⁻¹. The above shows that neither viscous nor bottom friction can mask the difference in the resonance frequencies for large-scale dynamical processes.

8.2 Comparison to Doppler shift and effective Coriolis parameter

Horizontal waves evolve in a medium that is moving with a velocity, u, which leads to a Doppler shift in the frequency at an Eulerian point. If the velocity, u, has a turbulent variability Δu, due to the turbulent dynamics of the flow field, also the Doppler shift in the frequency is modified and the modification can be compared to the error done by the hydrostatic approximation. The phase velocities are given by $c^{\mathrm{N}}={\mathit{\omega }}^{\mathrm{N}}/k$ and $c^{\mathrm{H}}={\mathit{\omega }}^{\mathrm{H}}/k$, where the wave frequencies are defined by the dispersion relation in Eq. (37). The square root of difference in the square of the horizontal wave velocity between the two formalisms, $\mathrm{\Delta }c=\sqrt{(c^{\mathrm{H}}{)}^{\mathrm{2}}-(c^{\mathrm{N}}{)}^{\mathrm{2}}}$, is $\mathrm{\Delta }c=\mathit{\gamma }c=\mathit{\gamma }{\mathit{\omega }}^{\mathrm{N}}/k_{\mathrm{h}}={\mathit{\omega }}^{\mathrm{N}}/k_z$ (see Eq. 37), which equals the vertical phase speed. The difference between the formalisms is significant if Δc>Δu – that is, if $\mathit{\gamma }>\mathrm{\Delta }u/c$. In the Strait of Gibraltar, the uncertainty in the horizontal advection due to natural variability is on the order of Δu=0.1 and c=1 m s⁻¹, so that for γ<0.1, the hydrostatic approximation can not be distinguished from the Navier–Stokes solution. With a depth of around 300 m, the horizontal wave length has to be below 3000 m, which asks for a resolution of 300 m or finer to distinguish between the two formalisms. The same can be applied to the group velocity of a moving wave packet, but as it is lower than the phase speed, it is even more affected by the Doppler shift. We see that internal variability of the horizontal velocity efficiently masks the differences between the two formalisms.

Horizontal gradients of the horizontal velocities lead to vorticity, ζ, and change the effective Coriolis parameter $f_{\mathrm{eff}}=\sqrt{\mathrm{1}+\mathit{\zeta }/f}f$ (for f>0) (see Wirth, 2013). If the Rossby number $\mathit{\zeta }/f>\mathrm{2}{\mathit{\gamma }}^{\mathrm{2}}{\mathit{\omega }}^{\mathrm{2}}/f^{\mathrm{2}}$, the change in the effective Coriolis parameter is larger than the difference in the formalisms. For oceanic values of $N^{\mathrm{2}}/f^{\mathrm{2}}=\mathrm{100}$ and γ<0.01, changes in the formalisms for Rossby numbers larger than $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{2}}$ are masked. A horizontal velocity difference across the Strait of Gibraltar of 0.3 m s⁻¹, leading to a Rossby number of 0.1, masks the changes in the frequencies of the formalisms for a horizontal wave length above 3 km.

8.3 Uncertainty estimation

This section discusses the estimation of uncertainty in time and space based on the theory of the Heisenberg–Gabor limit (see Folland and Sitaram, 1997 for a general discussion and Mordant, 2010 for an application similar to the present). Considering the time-frequency formalism, it shows that the resolution of the frequency (described by the standard deviation, σ_f, in the frequency domain of the normalized signal) and the length of the signal (described by the standard deviation of the temporal extent of the signal σ_t) are restrained by

$$\begin{array}{}\text{(81)}& {\mathit{\sigma }}_{\mathrm{f}}{\mathit{\sigma }}_{\mathrm{t}}\ge {\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}},\end{array}$$

equality being obtained only by a normal distribution. When the space-wavenumber formalism is considered, the standard deviations are replaced by σ_x and σ_k. Note that the inequality (81) applies to perfect data not disturbed by superposed signals or noise, so it is independent of the quality of the observation or the model; it is purely theoretical. The difference in the square frequencies between the two formalisms, obtained from Eq. (37), is (ω^Nγ)². To allow the observations to discriminate between the two frequencies, we need σ_f≤ω^Nγ; the Heisenberg–Gabor limit leads to an observation time ${\mathit{\sigma }}_{\mathrm{t}}\ge (\mathrm{2}\mathit{\pi }\mathit{\omega }\mathit{\gamma }{)}^{-\mathrm{1}}$. For short waves, for example wave trains generated by tidal forcing in the Strait of Gibraltar, we have ω≈10f. A resolved traveling wave asks for a grid scale which is at least $\mathrm{1}/\mathrm{10}$ of the wave length, and if the horizontal resolution is about 10 times the water depth, we obtain γ>0.01. This leads to an observation time of the phenomenon longer than the Coriolis period, which exceeds the lifetime of waves generated by tidal motion in, for example, the Strait of Gibraltar. This means that even in the absence of noise, currents and other phenomena, it is mathematically impossible to distinguish between the hydrostatic and the Navier–Stokes wave frequency. Furthermore, during this time, the advective and the group velocity will have transported the wave packet several tens of kilometers to regions where the depth and stratification largely differ, further altering the wave frequency and the phase and group velocity.

If waves with larger horizontal extensions are considered, γ≫0.01 with ω≈f, the observation time to distinguish between the hydrostatic and the Navier–Stokes wave-frequencies increases to several weeks. This means the differences between the two formalisms cannot be detected in the temporal domain.

If the spatial resolution is fine enough so that vertical convection can be explicitly resolved, $\mathrm{\Delta }x=\mathrm{\Delta }z<\mathrm{10}$ m, we have γ≈1, the domain is essentially unstratified, so that ω≈f, and we obtain $t_s\ge {\mathit{\sigma }}_{\mathrm{t}}=\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }f}=\mathrm{2}$ d, which is about the time an oceanic deep convection plume takes to go from the surface to the bottom. So also in the case of oceanic deep convection, the difference between the two projection-evolution operators cannot be detected. The difference in the projection operator in the two formalisms might be significant due to the advective transport of vertical velocity, which is present in the Navier–Stokes formalism but neglected in the hydrostatic approximation.

In the present paragraph, I will compare the difference in the frequencies between the two formalisms to the observation error in the horizontal and the vertical wave number. The calculations are based on the dispersion relations given in Eq. (37). In the spatial domain, the wave lengths can be considered in the vertical and the horizontal domain. This is done separately in the horizontal (supposing $k_z^{\mathrm{H}}=k_z$) and the vertical (supposing $k_{\mathrm{h}}^{\mathrm{H}}=k_{\mathrm{h}}$). In the horizontal, the relation between the hydrostatic and the Navier–Stokes horizontal wave number is

$$\begin{array}{}\text{(82)}& (k_{\mathrm{h}}^{\mathrm{H}}{)}^{\mathrm{2}}-(\mathrm{1}-{\displaystyle \frac{f^{\mathrm{2}}}{N^{\mathrm{2}}}}){\displaystyle \frac{k_z^{\mathrm{2}}}{k^{\mathrm{2}}}}{k}_{\mathrm{h}}^{\mathrm{2}}=\mathrm{0}.\end{array}$$

For a strong stratification and small γ this approaches $(k_{\mathrm{h}}^{\mathrm{H}}{)}^{\mathrm{2}}=(\mathrm{1}-{\mathit{\gamma }}^{\mathrm{2}})k_{\mathrm{h}}^{\mathrm{2}}$. Using the Heisenberg–Gabor limit (Eq. 81) means that the observation length has to exceed ${\mathit{\gamma }}^{-\mathrm{2}}/\left(\mathrm{4}\mathit{\pi }\right)$ the wave length, which means that for a rather large γ≈0.1, we need at least 10 wave lengths to distinguish between the two formalisms, and γ≈0.01 asks for an observation of around 1000 wave periods; no wave packets of such extension are observed in the ocean.

In the vertical, the solution is

$$\begin{array}{}\text{(83)}& (\mathrm{1}-{\displaystyle \frac{f^{\mathrm{2}}}{N^{\mathrm{2}}}})(k_z^{\mathrm{H}}{)}^{\mathrm{2}}-(\mathrm{1}+{\mathit{\gamma }}^{\mathrm{2}})k_z^{\mathrm{2}}=\mathrm{0}.\end{array}$$

For the leading vertical modes, there are only from one to a few half periods imposed by the depth of the domain. The question I answer is by how much the depth has to vary to perfectly compensate the differences of the two formalisms. If N=10f, a typical value for the ocean with a buoyancy frequency of a few hours, and γ=0.01, the difference is 1 %. With a vertical extension of 500 m of the phenomena, for example in the Strait of Gibraltar, this asks for a vertical resolution of the topography below 5 m. Such precision is not possible with today's models; for example, the topography or the thermocline depth largely exceeds 5 m from one horizontal grid point to the next. So also in this case, the difference in the projection-evolution operator between the two formalisms is insignificant.

Finally, I consider by how much the buoyancy frequency has to differ to compensate the error done by the hydrostatic approximation. To make the wave frequency of both formalisms agree at the same wave vector, the difference in the buoyancy frequency is $(N^{\mathrm{H}}{)}^{\mathrm{2}}=(N^{\mathrm{2}}-f^{\mathrm{2}})/(\mathrm{1}+{\mathit{\gamma }}^{\mathrm{2}})$; for N²=100f² and γ=0.1 this leads to a difference in buoyancy frequency of less than 2 %, far below the precision of today's models or observations.

The above discussions on the uncertainty of the horizontal, vertical wave number and buoyancy frequency demonstrate that the difference in the projection-evolution operators of the formalisms is not significant in realistic simulations of ocean dynamics, for small-amplitude linearized motion. Only nonlinear or externally forced processes can lead to configurations that project differently into the subspace of divergenceless fields in the two formalisms.

9 Conclusions

The present paper displays a theoretical framework to investigate the differences between the Navier–Stokes model and the hydrostatic approximation. To investigate the fidelity of the hydrostatic approximation, I constructed a four-dimensional local projection of the dynamics in Fourier space by separating the buoyant and dynamic vertical acceleration (T1 + T2) and combined it with the linear projection-evolution operator (T3). The linearized dynamics in the hydrostatic and Navier–Stokes formalism is compared. It shows that the hydrostatic approximation relies on two properties: the difference in the eigenvectors of the linear operator between the two formalisms and the difference in the corresponding eigenvalues. The former shows in the projection on the eigenvectors (S1), the latter in the evolution in the eigenspaces (S2). The projection depends on the Coriolis frequency f, the buoyancy frequency N and the three-dimensional wave vector k. The projection-evolution operator depends on the same parameters plus the forcing frequency ζ.

The eigenspaces of the wave dynamics (2D) and the balanced dynamics (1D) are identical between both formalisms. The eigenvalue for the balanced flow agrees between the two formalisms, whereas the eigenvalues in the wave dynamics differ due to the different wave frequencies. Therefore, for the projection-evolution operator (S2), the hydrostatic approximation is exact if the balance dynamics is considered, whereas it differs to second order O(γ²) for the wave dynamics. The fourth eigenvector is in the kernel of the corresponding projection. It is different in both formalisms to zeroth order O(γ⁰), but as the corresponding eigenvalue is vanishing, it does not participate in the dynamics.

The difference in the projection (S1) between the two formalisms is always along the second vector of the eigenspace of the wave dynamics (${\stackrel{\mathrm{̃}}{\mathit{e}}}_{\mathrm{2}}$). It corresponds to horizontally divergent flow, with a vanishing vertical component of the vorticity vector and is caused either by a purely vertical force, with error O(γ), or a horizontally divergent force, with error O(γ²). The projection of a rotational forcing and a buoyant forcing on the eigenvectors is identical in both formalisms. This makes the hydrostatic approximation first order, except if the vertical dynamical acceleration is first order, which would make the hydrostatic approximation second order O(γ²). An expansion of the nonlinear terms in powers of γ suggests that its vertical acceleration is indeed first order.

The error in the projection-evolution operator is second order O(γ²) for the wave dynamics and exact for the balanced flow, as the dynamics of the latter is perfectly described by the hydrostatic approximation. The validity of the hydrostatic approximation relies not only on the order of approximation in the scale ratio, γ, but also on the amount of energy that is found in the subspaces spanned by the wave modes and the balanced mode. For fast motion, if the frequency of the dynamics is larger than the minimum of either f and N, resonances can occur and the hydrostatic approximation breaks down even if the projection-evolution operator is formally second order. But for slower and slower motion, of frequency ζ, relatively more and more energy goes into the balanced mode, for which the hydrostatic approximation is exact. I showed that this depends also on the viscous damping. Another important result is that I manage to write the linear projection-evolution operators so that the lowest power of the wave frequency ω is quadratic, which ensures that the approximation for the linear projection-evolution operator is second order.

There are two important points one has to keep in mind. First, the hydrostatic approximation is a singular perturbation of the full Navier–Stokes equations, as the highest (temporal) derivative in the equation for the vertical velocity is omitted and expansions in small parameters are no guarantee for the smallness of the error (see, e.g., Bender et al., 1999). Second, the hydrostatic approximation breaks down where neither rotation nor stratification nor scale ratio constrain the fluid – that is, in fast-moving unstratified fluid dynamics. Where, when and how such areas are spontaneously created by the turbulent dynamics or external forcing depends on the situation considered. For these reasons, dedicated numerical experiments of oceanic phenomena have to be employed to determine the fidelity of the hydrostatic approximation for the scales and phenomena studied.

Considering a non-vanishing viscosity is of importance as in numerical models, the non-hydrostatic effects are present at the smallest scales of the model at which the numerical viscosity parameter has a dominant effect. It is straightforward to extend the calculations to an anisotropic viscosity, which has vertical values much smaller than the horizontal counterparts, usually employed due to the disparate resolution in the vertical and horizontal. Such anisotropy, however, hinders the commencement of the fully three-dimensional dynamics. Comparing the impact of anisotropic friction to the hydrostatic approximation has not been done here, as in fully non-hydrostatic simulations, which resolve the dynamics down to the meter scale, an anisotropic viscous damping lacks justification.

It is shown that if the vertical turbulent viscosity parameterizes the turbulent Ekman layer dynamics, the viscosity is high enough to mask the differences between the two formalisms in the wave frequencies at scales below the Ekman layer thickness (≈30 m). If it is lower ${\mathit{\nu }}_z<{\mathrm{10}}^{-\mathrm{2}}$ m² s⁻¹, the vertical turbulent exchange of momentum has to be resolved explicitly, which asks for a grid scale in the vertical and horizontal of 1 m or finer. At these scales, the free surface variations become important, leading to Langmuir cells (see McWilliams et al., 1997). For an explicit representation of these processes, the fully three-dimensional Navier–Stokes equations with a free surface should be used. Such models will not be available for basin-wide simulations in the near future, and parameterizations for surface processes have to be used.

At larger scales, the Doppler shift, the change in the effective Coriolis parameter due to unresolved processes and the internal variability of the dynamics, and the variation of the topography do not allow us to discriminate between the two formalisms in realistic ocean simulations. Furthermore, the uncertainty obtained through the Heisenberg–Gabor limit shows that the observation time of propagating inertia-gravity waves does not allow us to discriminate between the two formalisms. The calculations in the previous section also show that the Heisenberg–Gabor limit allows us to determine the gray zone for an oceanic process, between the large scales where the hydrostatic approximation is valid, down to the scale where it clearly fails. In conclusion: if the vertical viscosity is larger than ${\mathit{\nu }}_z>{\mathrm{10}}^{-\mathrm{2}}$ m² s⁻¹, the hydrostatic approximation is appropriate to represent motion slower then the Coriolis frequency in realistic simulations of the ocean dynamics. Differences between the two formalisms can be apparent in simulations of specific oceanic processes in an idealized setting.