\chapter{Evolution equations of incompressible 2D fluids} 
%
\section{Classical non-rotating case} \label{sec_2Dcase}
%
The dimensional evolution equations of incompressible
homogeneous 2D-fluids on the plane (see e.g.\ \cite{batchelor1967}
or \cite{canutoetal1988}) in vorticity velocity form
\begin{equation} \label{eq_2Dvortvel}
  \zeta_{t} + \mathbf{u} \cdot \nabla \zeta = F + D,
\end{equation}
with $\mathbf{u} = \left( u,v \right)$ the vector of velocity fields 
in $x$ and $y$-direction, $F$ a forcing and $D$ a dissipation term. 
From homogeneity and incompressibility of the fluid we get
\begin{equation} \label{eq_conti}
  \mathbf{\nabla \cdot} \rho \mathbf{u} 
    = 
  \rho \mathbf{\nabla \cdot u} = 0,
\end{equation} 
with $\rho$ the fluid density. Using this property we can introduce 
a volume (mass) stream function measuring the volume (mass) flux  
across an arbitrary line from the point $(x_{0},y_{0})$ to a point $(x,y)$ 
via the path integral
\begin{equation} \label{eq_psizetadef}
  \psi(x,y) - \psi(x_{0},y_{0}) = - \int_{(x_{0},y_{0})}^{(x,y)}
  \left[
   \left(
    \begin{array}{c}
     u  
     \\ 
     v
    \end{array}
   \right)
    \mathbf{\cdot}
   \left(
    \begin{array}{c}
     -dy  
      \\ 
      dx 
    \end{array}
   \right)
  \right].
\end{equation}
The minus sign in front of the integral just changes the direction
of positive massflux across the line and is chosen to make 
the classical stream function of 2D fluids compatible to 
the stream function (see e.g.\ \cite{danilovandgurarie2000}) 
typically used for 2D rotating fluids in geosciences. 
In terms of the stream function $\psi$ the integrated mass flux 
$\mathcal{M}$ across a line joining the points $(x_{0},y_{0})$ and 
$(x,y)$ is given by
\begin{equation} \label{eq_calM}
  \mathcal{M}(x_{0},y_{0}|x,y)
   = - \rho H_{0} 
       \left[\psi(x,y) -  \psi(x_{0},y_{0}) \right],
\end{equation}
with $H_{0}$ the depth of the fluid. At the same time
the stream function $\psi$ is connected to the velocity 
field $(u,v)$ and the (relative) 
vorticity $\zeta = v_{x} - u_{y}$ via
\begin{equation} \label{eq_psiuv}
  (u,v) = (- \psi_{y},\psi_{x}) \ \mbox{and} \ \ \zeta = \Delta \psi.
\end{equation}
Using the stream function $\psi$, the evolution equation (\ref{eq_2Dvortvel}
can be also written in the form
\begin{equation} \label{eq_2Dvortstream}
  \zeta_{t} + J(\psi,\zeta) = F + D.
\end{equation}
Since the velocity field is divergence free (equation \ref{eq_conti})
we can further derive the equivalent equation in flux form
\begin{equation} \label{eq_2Dflux}
  \zeta_{t} + \nabla \cdot \left(\mathbf{u} \zeta \right) = F + D.
\end{equation}
Starting from the flux form one can after repeated application of the
continuity equation (\ref{eq_conti}) finally derive an equation where
the Jacobian only depends on the velocity two fields $(u,v)$
\begin{equation} \label{eq_2Duv}
  \zeta_{t} 
   + \partial_{x} \partial_{y} \left( v^{2} - u^{2} \right) 
   + \left(\partial^{2}_{x} - \partial^{2}_{y} \right) uv
   = 
  F + D.
\end{equation}
This form allows in the pseudo-spectral method (see \ref{ssec_evolfourier}) 
to reduce the number of Fourier transforms per time-step from 
$3$, i.e.\ ($u,v,\zeta$) to $2$, i.e.\ ($u,v$).  

%
\section{Quasi-two-dimensional rotating case} \label{sec_quasi2Dcase}
%
Starting from the shallow water equation on the $\beta$-plane one can
derive (see e.g. \cite{danilovandgurarie2000}) an equation
describing a rotating barotropic quasi-two-dimensional fluid which is a
generalization of the 2D-equation (\ref{eq_2Dvortvel}). For small Rossby 
numbers $\mathrm{Ro} = U/Lf$, with $U$ a typical horizontal fluid velocity
at the length scale $L$ considered and $f$ the local coriolis paramter 
we get the potential vorticity (PV) $q$ in quasi-geostrophic (QG) approximation 
\begin{equation} \label{eq_qdef}
  q = \left( \nabla^{2}- \alpha^{2} \right) \psi + f,
\end{equation}
where we have the modification parameter $\alpha = 1/L^{2}_{\mathrm{R}}$ 
with the Rossby-Obukhov radius of deformation 
$L_{\mathrm{R}} = \sqrt{g H_{0}}/f$ and the stream function 
$\psi = gh/f$. Here $g$ is the gravitational acceleration and $h(x,y)$ 
the deviation of the mean fluid depth $H_{0}$ of the original 
shallow water layer. As in the $2$-dimensional case the stream function $\psi$ 
(remind the different definition) is related to the velocity fields $(u,v)$ 
as defined in equation (\ref{eq_psiuv}). In an unforced non-dissipative fluid
the QG PV is materially conserved
\begin{equation} 
  q_{t} + J(\psi,q) = 0. 
\end{equation}  
Using the linear approximation of the coriolis parameter $f = f_{0} + \beta y$
and introducing again forcing and dissipation we can write the evolution
equation in the form
\begin{equation} \label{eq_quasi2Dbaro}
  q_{t} + J(\psi,q) + \beta \psi_{x} = F + D,
\end{equation}
with the vorticity $q$ given by
\begin{equation} \label{eq_vortquasi2Dbaro}
  q = \left(\nabla^2 -\alpha^{2} \right) \psi 
    = \zeta - \alpha^{2} \psi.
\end{equation}
Further one can introduce a variable mean fluid depth $H(x,y)$, 
which in the simple case of a linear slope in $y$-direction leads 
to a topographic $\beta$-effect (see e.g.\ \cite{vanheist1994}). 

In the form (\ref{eq_quasi2Dbaro}) and (\ref{eq_vortquasi2Dbaro}) 
one can simulate incompressible 2D fluids and rotating quasi-2D fluids 
with the same set of equations using different parameters. 
In this more general frame the simplest case of a non-rotating
2D incompressible fluid is characterized by a vanishing 
ambient vorticity gradient, i.e.\ $\beta = 0$, and the limit of
an infinite Rossby radius $L_{\mathrm{R}} \longrightarrow \infty$
or a vanishing modification parameter $\alpha \longrightarrow 0$.

One only have to keep in mind that the stream functions different in 
the two cases. In the non-rotating case $\psi$ is defined by 
equation (\ref{eq_calM}). In the rotating case we get 
\begin{equation} \label{eq_hquasi2Dbaropsi}
  h(x,y) = \frac{f}{g} \ \psi(x,y),    
\end{equation} 
so that $\psi$ is proportional to pressure deviations, which is not 
the case in the non-rotating 2D case where the relation is more
complex, see e.g.\ \cite{johnstonandliu2004}.

Using the property of the Jacobian $J(f,f) = 0$ for all fields
$f(x,y)$ on the fluid domain equation (\ref{eq_quasi2Dbaro})
is equivalent to
\begin{equation} \label{eq_quasi2Dbarozeta}
  q_{t} + J(\psi,\zeta) + \beta \psi_{x} = F + D,
\end{equation}
where the vorticity $q$ is still defined by equation 
(\ref{eq_vortquasi2Dbaro}). From this from it directly follows that
that the form of the 2D Jacobian in equations 
(\ref{eq_2Dflux}) and (\ref{eq_2Duv}) can be also applied in the 
quasi-two-dimensional rotating case.

\section{Non-adiabatic terms}

\subsection{Laplacian based Viscosity and friction}

Internal viscosity and external friction of the fluid are described by 
the dissipation term $D$ on the left hand side of the equation 
(\ref{eq_quasi2Dbaro}). A classical way - the default reference in our
 model - to describe this term is to use a linear operator wich is a
superposition of powers of the Laplacian. 
The default dissipation in the fluid simulator is defined by
\begin{equation} \label{eq_Laplace_dissip}
  D \ q = - \left[\sigma \left(-1 \right)^{p_{\sigma}} \Delta^{p_{\sigma}}
                     +
                  \lambda \left(-1 \right)^{p_{\lambda}} \Delta^{p_{\lambda}}
            \right] q.
\end{equation}
Introducing dissipation time and length-scales we can write the 
dissipation operator also in the form 
\begin{equation} \label{eq_Laplace_dissip_02}
  D \ q = - \left[ 
              \frac{\left(-1 \right)^{p_{\sigma}}}{t_{\sigma}} 
              \left( 
               \frac{L_{\sigma}}{2 \pi}
              \right)^{2 p_{\sigma}}
              \Delta^{p_{\sigma}}
               +
              \frac{\left(-1 \right)^{p_{\lambda}}}{t_{\lambda}} 
              \left( 
                \frac{L_{\lambda}}{2 \pi}
              \right)^{2 p_{\lambda}}
              \Delta^{p_{\lambda}}
            \right] q.
\end{equation}
Here $L_{\sigma}$ and $L_{\lambda}$ are the small and 
large-scale cut-off length scales. The corresponding small and 
large-scale "damping" time scales are given by $t_{\sigma}$ 
and $t_{\lambda}$. The powers $p_{\sigma}$ of small-scale 
viscosity are in the range of $p_{\sigma} \in [1,2,3, \ \dots \ ]$ 
and the powers $p_{\lambda}$ of large-scale friction are in the 
range $p_{\lambda} \in [0,-1,-2,-3, \ \dots \ ]$. 
For $p_{\sigma} = 1$ and $p_{\lambda} = 0$ we speak of viscosity
and linear drag (friction), for $p_{\sigma} > 1$ and $p_{\lambda} < 0$
of hyperviscosity and hypofriction (see also \cite{danilovandgurarie2001}).
From equation (\ref{eq_Laplace_dissip}) it follows that the coefficients
$\sigma$ and $\lambda$ can be written by
\begin{equation} \label{eq_siglam}
  \sigma  = \left(\frac{L_{\sigma}}{2 \pi} \right)^{2 p_{\sigma}}
            \ \frac{1}{t_{\sigma}} 
          = \left(\frac{1}{k_{\sigma}}\right)^{2 p_{\sigma}}
            \ \frac{1}{t_{\sigma}} 
  \ \ \mbox{and} \ \
  \lambda = \left(\frac{L_{\lambda}}{2 \pi} \right)^{2 p_{\lambda}}
            \ \frac{1}{t_{\lambda}}
          = \left(\frac{1}{k_{\lambda}} \right)^{2 p_{\lambda}}
            \ \frac{1}{t_{\lambda}}.
\end{equation}
We have chosen the additional factor of $2 \pi$ since finally the
fluid domain is rescaled to multiples of $2 \pi$. Using this scaling
the coefficients $\sigma$ and $\lambda$ are characterized respectively 
by the damping the time scales $t_{\sigma}$ and $t_{\lambda}$ as well
as the cut-off wave numbers $k_{\sigma}$ and $k_{\lambda}$.

Above dissipation operator is a special case of dissipation operators 
which are polynomials with positive and negative powers of the Laplacian
\begin{equation} \label{eq_Laplace_dissip_poly}
  D q = \sum_{n = 1}^{n_{max}} D_{n} q,
  \ \mbox{with} \ \
  D_{n} q 
   = 
 -\left[
   \sigma_{n-1} \left(-1 \right)^{n-1} \ \Delta^{n-1}
   \ +
   \lambda_{n} \left(-1 \right)^{-n} \ \Delta^{-n}
  \right] \ q.
\end{equation}
In Fourier space (see subsection \ref{ssec_evolfourier}) the dissipation 
operators of this class reduce to the multiplication with polynomials 
in positive and negative powers of the wave numbers. More general 
dissipation operators can be constructed directly in Fourier space 
(see subsection \ref{ssec_Dschemes}).

\subsection{Forcing}
The forcing term $F$ in equations (\ref{eq_2Dvortvel}) and 
(\ref{eq_quasi2Dbaro}) describes forcings due to either 
external processes as a wind-stress or a moving plate or
non-resolved internal processes as, e.g.\ baroclinic instability 
or diabatic heating. Both types of forcings can be described 
in physical or spectral space. For a constant external forcing, 
e.g.\ a wind stress or drag of a moving plate $(\tau^{u},\tau^{v})$ 
given in $[N/m^{2}]$ and acting in $x$ and $y$ direction at the surface 
of the fluid the forcing term is given by 
\begin{equation} \label{eq_Fstressdrag}
 F(x,y) = \frac{1}{\rho H_{0}} \left( \tau_{x}^{v} - \tau_{y}^{u} \right),
\end{equation}
where in the rotating quasi-2D case the height deviations $h$ of the
fluid are neglected. The forcing can also be defined in spectral
space, see subsection \ref{ssec_Fschemes} below.
%---
\section{Geometry and boundary conditions}
%
\subsection{Doubly periodic boundary condition}
%
The evolution equations have to be completed by boundary conditions.  
The default geometry of the fluid domain is a square with an edge of length
$L_{x} = L_{y} = L$ and the default boundary conditions are doubly periodic,  
i.e.\ $(f(x,y) = f(x+L,y+L)$ for all functions $f$ on the fluid domain.
%
\subsection{Channel boundary condition}
%
In $x$-direction (zonal direction) we have periodic boudary conditions, 
i.e.\ $f(x,y) = f(x+L,y)$. In $y$-direction (meridional direction) 
at $y=0$ and $y=L$ we introduce walls with no-slip boundary conditions, 
i.e.\ the meridional velocity at the walls is zero $v(x,0) = v(x,L) =  0$.

%
\subsection{Box boundary condition}
%
In $x$-direction (zonal direction) we introduce walls with no-slip boundary
conditions, i.e.\ the zonal velocity at the walls is zero 
$u(0,y) = u(L,y) = 0$. In $y$-direction (meridional direction) 
at $y=0$ and $y=L$ we introduce walls with no-slip boundary conditions, 
i.e.\ the meridional velocity at the walls is zero $v(x,0) = v(x,L) =  0$.