.. _nbq:

*****************************************
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 :math:`C` -- and 
the equation of state, which reads in Cartesian coordinates:

.. math::

     \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


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

A relation between :math:`\rho` and :math:`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, :math:`s` 
and :math:`f` subscripts refer to slow and fast-mode components, respectively:

.. math::

    \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}}}

:math:`c_s` is the speed of sound and :math:`\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:

.. math::

    \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 &

    \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


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 :cite:t:`auclair_non-hydrostatic_2018`. 

Note that the solved acoustic waves can become pseudo-acoustic if their phase 
speed :math:`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.