1.4. Non-Hydrostatic, Non-Boussinesq Equations

The full set of Navier-Stokes equations for a free-surface ocean is explicitly integrated in the non-hydrostatic, non-Boussinesq version of CROCO (#define NBQ). In this approach, acoustic waves are solved explicitly to avoid Boussinesq-degeneracy, which inevitably leads to a 3D Poisson system in non-hydrostatic Boussinesq methods – detrimental to computational costs and challenging to implement within a split-explicit barotropic/baroclinic model.

NBQ equations include the momentum and continuity equations, the surface kinematic relation (for free surface), temperature, salinity – or other tracer \(C\) – and the equation of state, which reads in Cartesian coordinates:

\[ \begin{align}\begin{aligned}\frac{\partial \rho u}{\partial t} + \vec{\bf \nabla} . \left ( \rho \vec{\textbf v} u \right ) - \rho f v + \rho \tilde{f} w & = - \frac{\partial P}{\partial x} + \lambda \frac{\partial \vec{\nabla} . \vec{\textbf v}}{\partial x} + \mathcal{F}_u + \mathcal{D}_u\\\frac{\partial \rho v}{\partial t} + \vec{\bf \nabla} . \left ( \rho \vec{\textbf v} v \right ) + \rho f u & = - \frac{\partial P}{\partial y} + \lambda \frac{\partial \vec{\nabla} . \vec{\textbf v}}{\partial y} + \mathcal{F}_v + \mathcal{D}_v\\\frac{\partial \rho w}{\partial t} + \vec{\bf \nabla} . \left ( \rho \vec{\textbf v} w \right ) - \rho \tilde{f} u & = - \frac{\partial P}{\partial z} -\rho g + \lambda \frac{\partial (\vec{\nabla} . \vec{\textbf v} )}{\partial z} + \mathcal{F}_w + \mathcal{D}_w\\\frac{\partial \rho}{\partial t} & = - \vec{\nabla} . (\rho \vec{\textbf v})\\\frac{\partial \xi}{\partial t} & = ~w_f \rvert_{z=\xi} -\vec{\textbf v} \rvert_{z=\xi} . \vec{\nabla} \xi &\\\frac{\partial \rho C }{\partial t} & = - \vec{\nabla} . (\rho \vec{\textbf v} C) + \mathcal{F}_C + \mathcal{D}_C\end{aligned}\end{align} \]

\(\lambda\) is the second (bulk) viscosity associated with compressibility (it can be used to damp acoustic waves).

A relation between \(\rho\) and \(P\) is now required. To that end, and as part of a time-splitting approach, density is decomposed into slow and fast components based on a first-order decomposition concerning total pressure. In the following, \(s\) and \(f\) subscripts refer to slow and fast-mode components, respectively:

\[ \begin{align}\begin{aligned}\rho & = \rho_s (T,S,P) + \overbrace{\frac{\partial \rho }{\partial P} \bigg\rvert_{T,S}\delta P}^{\rho_f= c_s^{-2} P_f} + O(\delta P^2 )\\P & = \underbrace{ P_{atm} + \int_{z}^{\xi} (\rho_s - \rho_0) g ~dz' }_{{\color{red} {SLOW}}} +\underbrace{ \rho_0 g (\xi -z) + \overbrace{\delta P}^{P_f} }_{{\color{red} {FAST}}}\end{aligned}\end{align} \]

\(c_s\) is the speed of sound and \(\delta P=P_f\) is the nonhydrostatic pressure.

The Navier-Stokes equations are then integrated with two different time steps within the time-splitting approach. The slow mode is identical to ROMS, whereas the fast mode (in the NBQ equations) is 3D and the fast time step includes the integration of the compressible terms of the momentum and continuity equations. In vector form:

\[ \begin{align}\begin{aligned}\begin{split}\frac{\partial \rho \vec{\textbf v}}{\partial t} = &\underbrace{ - \vec{\nabla} . (\rho \vec{\textbf v} \otimes \vec{\textbf v}) - 2 \rho \vec{\Omega} \times \vec{\textbf v} - \vec{\nabla} (\int_{z}^{\xi_f} (\rho_s - \rho_0) g ~dz') + \vec{ \mathcal{F}}_{\vec{\textbf v} }+ \vec{ \mathcal{D}}_{\vec{\textbf v} } }_{{\color{red} {SLOW}}}\\ &\underbrace{ -\rho_0 g \vec{\nabla} \xi_f - \vec{\nabla} P_f + \rho \vec{\textbf g} + \lambda \vec{\nabla} (\vec{\nabla} . \vec{\textbf v} )}_{{\color{red} {FAST}}}\\ \frac{\partial \rho_f}{\partial t} = & - \frac{\partial \rho_s}{\partial t} - \vec{\nabla} . (\rho \vec{\textbf v}) &\\ P_f = & ~c_s^{2} ~\rho_f &\end{split}\\\frac{\partial \xi_f}{\partial t} = & ~w_f \rvert_{z=\xi} -\vec{\textbf v}_f \rvert_{z=\xi} . \vec{\nabla} \xi_f &\\\frac{\partial \rho C_s }{\partial t} = & - \vec{\nabla} . (\rho \vec{\textbf v} C_s) + \mathcal{F}_C + \mathcal{D}_C &\\\rho_s = & ~\rho (T_s,S_s,\xi_f) &\\\rho = & ~\rho_s + \rho_f\end{aligned}\end{align} \]

The momentum is integrated both in slow and fast modes but the right-hand-side of the equation is split in two parts: a slow part, made of slowly varying terms (advection, Coriolis force, baroclinic pressure force and viscous dissipation), and a fast part, made of fast-varying terms (the surface-induced and compressible pressure force, the weight and the dissipation associated with bulk-viscosity). This momentum equation is numerically integrated twice, once with a large time-step keeping the fast part constant, and once with a smaller time-step keeping the slow part constant. This is much more computationally efficient than integrating the whole set of equations at the same fast time step. More details can be found in Auclair et al. (2018).

Note that the solved acoustic waves can become pseudo-acoustic if their phase speed \(c_s\) is artificially slowed down (it is a model input). In this case, high-frequency processes associated with bulk compressibility may be unphysical, but a coherent solution for slow non-hydrostatic dynamics is preserved, while the CFL constraint is relaxed.