Equations used in the solvers
=============================



Shallow-water equations
-----------------------

We consider a layered ocean with upward oriented z-axis:

.. image:: nlayers.png
   :scale: 70 %

.. |uu| mathmacro:: \mathbf{u}
.. |p| mathmacro:: \partial
.. |ez| mathmacro:: \mathbf{e_z}
.. |bnabla| mathmacro:: \boldsymbol{\nabla}
.. |D| mathmacro:: \mbox{D}
.. |xx| mathmacro:: \mathbf{x}

The adiabatic and non-dissipating shallow water equations for :math:`N` layers can
be written as

.. math::

   \D_t \uu_n = - \bnabla g m_n,

.. math::

   \p_t h_n(x,y) + \bnabla \cdot (\uu_n h_n ) = 0.

for :math:`n=0,\cdots, N-1` (:math:`n=0` corresponding to the top layer) and where
:math:`\uu_n` and :math:`h_n` are the layer 2d velocity and thickness, respectively,
and :math:`g` is the gravitational acceleration.

The normalized Montgomery potentials :math:`m_n = M_n / (g\rho_n)` are
proportional to the Montgomery potentials :math:`M_n` and can be computed as
:math:`m_0 = z_0` and, for :math:`n>0`, as

.. math::

    m_n = \frac{\rho_{n-1}}{\rho_n} m_{n-1}
    + \left(1 - \frac{\rho_{n-1}}{\rho_n}\right) z_n.


Here :math:`\rho_n` is the density in layer :math:`n` and :math:`z_n` is the
:math:`z`-coordinate of layer :math:`n` upper boundary, given by:

.. math::

   z_n &= z_{n-1} - h_{n-1},

       &= - H(\xx) + \sum_{j=n,N-1} h_j.

For example for :math:`N=2`, :math:`z_0 = - H(\xx) + h_1 + h_0` and :math:`z_1 = -
H(\xx) + h_1`. The total depth can be decomposed as the sum of a constant "full
depth" :math:`\tilde H` and a topography :math:`T(\xx)` (the height of the
submarine mountains), :math:`H(\xx) = \tilde H -  T(\xx)`, so that we find

.. math::

   \bnabla z_n = \bnabla \big(T + \sum_{j=n,N-1} h_j\big).

The equation for normalized Montgomery potential may be rewritten into a useful
equation for interface vertical position:

.. math::

   z_n = \frac{\rho_n m_n -\rho_{n-1} m_{n-1}}{\rho_n - \rho_{n-1}}.

Note that since the velocity fields :math:`\uu_n` are not divergent-free, we have
:math:`\D_t \uu_n = \p_t \uu_n + \bnabla |\uu_n|^2/2 + \zeta_n \ez\times \uu_n`, with
:math:`\zeta_n = \ez \cdot (\bnabla \times \uu_n )` the relative vorticity.


Considering now Earth rotation, forcing, and, dissipating terms, the shallow water
equation for conservation of momentum can be written as:

.. math::

   \p_t \uu_n(x,y) + (\zeta_n + f) \ez\times\mathbf{u}_n =
   -\bnabla \Big \{ g (m_n + \Pi_n) + \frac{1}{2} |\uu_n|^2 \Big \}
   + \mathbf{F}_{un} + \mathbf{D}(\uu_n),

for :math:`n=0,\cdots, N-1` and where:

- :math:`f` is the Coriolis frequency,
- :math:`\Pi_n`  is a forcing potential, presumably tidal,
- :math:`\mathbf{F}_{un}` gathers non potential forcing (typically relaxations),
- :math:`\mathbf{D}(\uu_n)` is dissipation.

The thickness tendency equation is:

.. math::

   \p_t h_n(x,y) + \bnabla \cdot (\mathbf{u}_n h_n ) = F_{hn},

where :math:`F_{hn}` represents thickness forcing and/or diabatic effects.

The pressure can be recomputed as:

.. math::

   p_n(x,y,z) = g \rho_n ( m_n(x,y) -  z).

The one-layer case for N=1 `is implemented in fluidsim
<http://fluidsim.readthedocs.io/en/latest/generated/fluidsim.solvers.sw1l.solver.html>`_.


Conserved quantities
~~~~~~~~~~~~~~~~~~~~

We need to write the formula for the main quantities conserved by the
equations.

- kinetic energy:

.. math::

   E_k = \sum_{0,N-1} \frac{1}{2} \rho_n h_n |\uu_n(x,y)|^2.

- potential energy:

.. math::

   E_p &= \sum_{0,N-1} \rho_n g \int_{z_{n+1}}^{z_n} z dz

       &= \sum_{0,N-1} \frac{1}{2} \rho_n g (2z_{n+1} + h_n) h_n.

- available potential energy:

.. math::

   E_{a} = \sum_{0,N-1} \frac{1}{2} \rho_n g (\eta_n^2 - \eta_{n+1}^2),

where :math:`\eta_n` is the departure of the n'th interface from its position at rest.

- potential enstrophy:

.. math::

   E_e = \sum_{0,N-1} \frac{1}{2} \rho_n \frac{(\zeta_n +f)^2}{h_n}.

Note that the potential vorticity is :math:`q_n=(\zeta_n+f)/h_n`


Vertical normal modes
---------------------

We assume now the layers are weakly departing from a rest state characterized by
layer thicknesses :math:`H_n` and associated interface vertical positions :math:`Z_n`.
We also assume the topography is flat.
Interface displacement from their positions at rest are labelled :math:`\eta_n`, such
that:

.. math::

   z_n = Z_n + \eta_n

The temporal evolution of layer thicknesses is then given by:

.. math::

   \p_t h_n(x,y) + H_n \bnabla \cdot \mathbf{u}_n = F_h.

One can then separate vertical and horizontal dependencies according to:

.. math::

   u_n &= u(x,y,t) \phi_n

   m_n &= m(x,y,t) \phi_n

   h_n &= h(x,y,t) \times \lambda \phi_n H_n

The vertical mode structure functions :math:`\phi_n` satisfy the following
eigenvalue problem:

.. math::

   \lambda \phi_n H_n = -\frac{\rho_{n-1}}{\rho_{n}-\rho_{n-1}} \phi_{n-1}
      + \Big \{ \frac{\rho_{n-1}}{\rho_{n}-\rho_{n-1}} + \frac{\rho_{n}}{\rho_{n+1}-\rho_{n}} \Big \} \phi_n
      - \frac{\rho_{n}}{\rho_{n+1}-\rho_{n}} \phi_{n},

for :math:`0<n<N-1`, and:

.. math::

   \lambda \phi_0 H_0 = \phi_0
      - \frac{\rho_{1} \phi_1 - \rho_{0} \phi_0}{\rho_{1}-\rho_{0}},

   \lambda \phi_{N-1} H_{N-1} =
      \frac{\rho_{N-1} \phi_{N-1} - \rho_{N-2} \phi_{N-2}}{\rho_{N-1}-\rho_{N-2}},

The eigenvalue if related to an equivalent speed :math:`c_n` according :math:`c_n=\sqrt{g/\lambda}`


Quasi-geostrophic equations
---------------------------

We consider now that the flow weakly depart from a reference state given by constant and
uniform layer thicknesses :math:`H_n`.

Each interface is departing from its mean position by :math:`\eta_n`.
The layer thickness is then given by :math:`h_n=H_n+\eta_{n}-\eta_{n+1}`.

Normalized Montgomery potentials anomaly are then given by :math:`\mu_0 = \eta_0` and
for :math:`n>0` as:

.. math::

   \mu_n = \frac{\rho_{n-1}}{\rho_n} \mu_{n-1}
         - \left(1 - \frac{\rho_{n-1}}{\rho_n}\right) \eta_n.

This last expression can be rewritten to express :math:`\eta_n` as a function of :math:`\mu`'s:

.. math::

   \eta_n &= \frac{\rho_n \mu_n - \rho_{n-1} \mu_{n-1}}{\rho_n - \rho_{n-1}},

      &\sim \frac{\rho_n}{\rho_n - \rho_{n-1}} \times (\mu_n - \mu_{n-1}),

where the last approximation will be considered as an equality next.

Under typical quasi-geostrophic conditions (weak Rossby number, weak topography), the flow
is at zeroth order geostrophic:

.. math::

   \uu_n = \ez \cdot \bnabla \big \{ g \mu_n/f_0 \big \}.

The geostrophic streamfunction is thus :math:`\psi_n=g \mu_n/f_0`.

The curl of momentum conservation is at first order given by:

.. math::

   \D_t \big \{ \zeta_n+f' \big \} + f_0 \bnabla\cdot\mathbf{u}_n =
   \bnabla \times \big \{ \mathbf{F}_{un} + \mathbf{D}(\uu_n) \big \},

where :math:`f'=f-f_0`.

The divergence of the horizontal flow is given by the thickness tendency equation:

.. math::

   \bnabla\cdot\mathbf{u}_n &= -\D_t h_n/H_n + F_{hn}/H_n,

      &= \D_t ( \eta_{n}-\eta_{n-1} )/H_n + F_{hn}/H_n,

which leads to the conservation of QG potential vorticity:

.. math::

   \D_t q_n =
   \bnabla \times \big \{ \mathbf{F}_{un} + \mathbf{D}(\uu_n) \big \} - f_0 F_{hn}/H_n .

where :math:`q_n` is the quasi-geostrophic potential vorticity:

.. math::

   q_n = \zeta_n + f' + f_0 ( \eta_{n}-\eta_{n-1} )/H_n,

For interior layers (:math:`0<n<N-1`):

.. math::

      q_n = \nabla^2 \psi_n + f'
       + \frac{f_0^2}{H_n} \Big \{ \frac{\psi_n - \psi_{n-1}}{g_n'}
       - \frac{\psi_{n-1} - \psi_{n-2}}{g_{n-1}'} \Big \}

where :math:`g_n'=g(\rho_n - \rho_{n-1})/\rho_n`.


Forcing: waves
--------------

We want to force internal waves with a given frequency :math:`\omega` and modal vertical
structure corresponding to mode `f`.
The forcing required to do that, while not producing spurious QG potential vorticity,
should look like:

.. math::

   \mathbf{F}_{un} = (f_u, f_v) h(x,y, t) \phi_{f,n},

   \mathbf{F}_{hn} =  (f_v \partial_x h - f_u \partial_y h) \;
      f_0 H_n \phi_{f,n},

where (:math:`f_u`, :math:`f_v`) are constants, :math:`h(x,y)` describes the horizontal
spatial and temporal structure of the forcing, :math:`\phi_f` is the :math:`f` vertical mode number.
The spatiotemporal structure is chosen to match that of the expected wave, i.e. a sinusoidal-like
shape with wave number :math:`k` satisfying the following dispersion relationship

.. math::

   \omega^2 = f^2 + c_f^2 k^2.

In general: :math:`h=\sin[ k (y-y_c)-\omega t]`

A dirac in y and spatial uniformity in x may also work.


Forcing: vortices
-----------------

In order to force a turbulent eddy field, we apply a forcing that looks like:

.. math::

   \mathbf{F}_{un} = f_{uv} (-\partial_y h, \partial_x h ) \phi_{f,n} ,

   \mathbf{F}_{hn} =  f_h h(x,y,t) \frac{g H_n}{c_f^2} \phi_{f,n}.

where the spatial and temporal structure of the forcing is given by:

.. math::

   h(x,y,t) = s \Big [(x - x_c) \cos(\omega t) + (y-y_c) \sin(\omega t) \Big ],

where :math:`s` is a smoothly varying step function.
The QG potential vorticity will be generated at the following rate:

.. math::

   \D_t q_n = \Big \{ f_{uv} \nabla^2 h  - f_h h \frac{gf}{c_f^2} \Big \} \phi_{f,n}

We choose for now :math:`f_{uv}=0`.