1.1. Primitive Equations
At resolutions larger than 1 km (more marginally above 100 m), The ocean is a fluid that can
be described as a good approximation by the primitive equations (PE). The PE equations are
simplifications from the Navier-Stokes equations made from scale considerations, along
with a nonlinear equation of state, which couples the two active tracers (temperature and salinity):
Hydrostatic hypothesis: the vertical momentum equation is reduced to a balance between
the vertical pressure gradient and the buoyancy force (non-hydrostatic processes such as
convection must be parametrized)
Boussinesq hypothesis: density variations are neglected except in their contribution
to the buoyancy force
Incompressibility hypothesis (stemming from the former): the three-dimensional divergence
of the velocity vector is assumed to be zero.
Spherical earth approximation: the geopotential surfaces are assumed to be spheres so
that gravity (local vertical) is parallel to the earth’s radius
Thin-shell approximation: the ocean depth is neglected compared to the earth’s radius
Turbulent closure hypothesis: the turbulent fluxes (which represent the effect of
small-scale processes on the large scale) are expressed in terms of large-scale features
By default (#undef NBQ), CROCO solves the primitive equations. However, it also
has the ability to relax the first 3 hypotheses (#define NBQ). When SOLVE3D is
not activated, CROCO can be used as a classical shallow water model.
1.1.1. Equations in Cartesian coordinates
\[ \begin{align}\begin{aligned}\frac{\partial u}{\partial t}
+ \vec{\bf \nabla} . \left ( \vec{\textbf v} u \right ) - f v =
- \frac{\partial \phi}{\partial x} + \mathcal{F}_u + \mathcal{D}_u\\\frac{\partial v}{\partial t}
+ \vec{\bf \nabla} . \left ( \vec{\textbf v} v \right ) + f u =
- \frac{\partial \phi}{\partial y} + \mathcal{F}_v + \mathcal{D}_v\\Turbulent closure schemes are applied to parametrized subgrid-scale vertical fluxes.\end{aligned}\end{align} \]
The time evolution of a scalar concentration field, \(C(x,y,z,t)\) (e.g., salinity,
temperature, or nutrients), is governed by the advective-diffusive equation :
\[\frac{\partial C}{\partial t}
+ \vec{\bf \nabla} . \left ( \vec{\textbf v} C \right ) =
\mathcal{F}_C + \mathcal{D}_C\]
\[\rho = \rho(T,S,P)\]
In the Boussinesq approximation, density variations are neglected in the momentum
equations except in their contribution to the buoyancy force in the vertical momentum
equation. Under the hydrostatic approximation, it is further assumed that the vertical
pressure gradient balances the buoyancy force :
\[\frac{\partial \phi}{\partial z} = - \frac{\rho g}{\rho_0}\]
\[\vec{\bf \nabla} . \vec{\textbf v} =
\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y}
+ \frac{\partial w}{\partial z} = 0\]
The variables used are :
\(\mathcal{D}_u, \mathcal{D}_v, \mathcal{D}_C\) : diffusive terms
\(\mathcal{F}_u, \mathcal{F}_v, \mathcal{F}_C\) : forcing terms
\(f(x,y)\) : Traditional Coriolis parameter \(2 \Omega sin \phi\)
\(g\) : acceleration of gravity
\(\phi(x,y,z,t)\) : dynamic pressure \(\phi=P/\rho_0\), with P the total pressure
\(\rho_0+\rho(x,y,z,t)\) : total in situ density
\(u,v,w\) : the (x,y,z) components of vector velocity \(\vec{\textbf v}\)
1.1.2. Equations in terrain-following coordinates
We first introduce a generalized stretched vertical coordinate system (\(s\)),
which sets the variable bottom flat at \(z=-h(x,y)\). \(s\) spans the range
from -1 (bottom) to 0 (surface), and the transformation rules are:
\[\begin{split} \begin{aligned}
{\left(\frac{\partial}{\partial x}\right)_{z}=\left(\frac{\partial}{\partial x}\right)_{s}-\left(\frac{1}{H_{z}}\right)\left(\frac{\partial z}{\partial x}\right)_{s} \frac{\partial}{\partial s}} \\
{\left(\frac{\partial}{\partial y}\right)_{z}=\left(\frac{\partial}{\partial y}\right)_{s}-\left(\frac{1}{H_{z}}\right)\left(\frac{\partial z}{\partial y}\right)_{s} \frac{\partial}{\partial s}} \\
{\frac{\partial}{\partial z}=\left(\frac{\partial s}{\partial z}\right) \frac{\partial}{\partial s}=\frac{1}{H_{z}} \frac{\partial}{\partial s}} \\
\text{where } {H_{z} \equiv \frac{\partial z}{\partial s}}
\end{aligned}\end{split}\]
The vertical velocity in \(s\) coordinate is:
\[\begin{split}\begin{aligned} \Omega(x, y, s, t) &=\frac{1}{H_{z}}\left[w-(1+s) \frac{\partial \zeta}{\partial t}-u \frac{\partial z}{\partial x}-v \frac{\partial z}{\partial y}\right] \\
w &=\frac{\partial z}{\partial t}+u \frac{\partial z}{\partial x}+v \frac{\partial z}{\partial y}+\Omega H_{z} \end{aligned}\end{split}\]
\(\Omega=0\) at both surface and bottom.
Next, the requirement for a laterally variable grid resolution can also be met for
suitably smooth domains by introducing an appropriate orthogonal coordinate
transformation in the horizontal. Let the new coordinates be \(\xi(x,y)\)
and \(\eta(x,y)\) where the relationship of horizontal arc length to the differential
distance is given by:
\[ \begin{align}\begin{aligned}\qquad(d s)_{\xi}=\left(\frac{1}{m}\right) d \xi\\\qquad(d s)_{\eta}=\left(\frac{1}{n}\right) d \eta\end{aligned}\end{align} \]
Here \(m(\xi,\eta)\) and \(n(\xi,\eta)\) are the scale factors which relate the
differential distances \((\Delta \xi, \Delta \eta)\) to the actual (physical) arc lengths.
\[\begin{split} \begin{aligned} \frac{\partial}{\partial t}\left(\frac{H_{z} u}{m n}\right)+\frac{\partial}{\partial \xi}\left(\frac{H_{z} u^{2}}{n}\right)+\frac{\partial}{\partial \eta}\left(\frac{H_{z} u v}{m}\right)+\frac{\partial}{\partial s}\left(\frac{H_{z} u \Omega}{m n}\right) \\-\left\{\left(\frac{f}{m n}\right)+v \frac{\partial}{\partial \xi}\left(\frac{1}{n}\right)-u \frac{\partial}{\partial \eta}\left(\frac{1}{m}\right)\right\} H_{z} v=\\-\left(\frac{H_{z}}{n}\right)\left(\frac{\partial \phi}{\partial \xi}+\frac{g \rho}{\rho_{o}} \frac{\partial z}{\partial \xi}+g \frac{\partial \zeta}{\partial \xi}\right)+\frac{H_{z}}{m n}\left(\mathcal{F}_{u}+\mathcal{D}_{u}\right) \\ \frac{\partial}{\partial t}\left(\frac{H_{z} v}{m n}\right)+\frac{\partial}{\partial \xi}\left(\frac{H_{z} u v}{n}\right)+\frac{\partial}{\partial \eta}\left(\frac{H_{z} v^{2}}{m}\right)+\frac{\partial}{\partial s}\left(\frac{H_{z} v \Omega}{m n}\right) \\+\left\{\left(\frac{f}{m n}\right)+v \frac{\partial}{\partial \xi}\left(\frac{1}{n}\right)-u \frac{\partial}{\partial \eta}\left(\frac{1}{m}\right)\right\} H_{z} u=\\-\left(\frac{H_{z}}{m}\right)\left(\frac{\partial \phi}{\partial \eta}+\frac{g \rho}{\rho_{o}} \frac{\partial z}{\partial \eta}+g \frac{\partial \zeta}{\partial \eta}\right)+\frac{H_{z}}{m n}\left(\mathcal{F}_{v}+\mathcal{D}_{v}\right) \\
\frac{\partial}{\partial t}\left(\frac{H_{z} T}{m n}\right)+\frac{\partial}{\partial \xi}\left(\frac{H_{z} u T}{n}\right)+\frac{\partial}{\partial \eta}\left(\frac{H_{z} v T}{m}\right)+\frac{\partial}{\partial s}\left(\frac{H_{z} \Omega T}{m n}\right) &=\frac{H_{z}}{m n}\left(\mathcal{F}_{T}+\mathcal{D}_{T}\right) \\ \frac{\partial}{\partial t}\left(\frac{H_{z} S}{m n}\right)+\frac{\partial}{\partial \xi}\left(\frac{H_{z} u S}{n}\right)+\frac{\partial}{\partial \eta}\left(\frac{H_{z} v S}{m}\right)+\frac{\partial}{\partial s}\left(\frac{H_{z} \Omega S}{m n}\right) &=\frac{H_{z}}{m n}\left(\mathcal{F}_{S}+\mathcal{D}_{S}\right) \\ \rho &=\rho(T, S, P) \\ \frac{\partial \phi}{\partial s} &=-\left(\frac{g H_{z} \rho}{\rho_{o}}\right) \\ \frac{\partial}{\partial t}\left(\frac{H_{z}}{m n}\right)+\frac{\partial}{\partial \xi}\left(\frac{H_{z} u}{n}\right)+\frac{\partial}{\partial \eta}\left(\frac{H_{z} v}{m}\right)+\frac{\partial}{\partial s}\left(\frac{H_{z} \Omega}{m n}\right) &=0 \end{aligned}\end{split}\]
