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

  • The momentum balance in zonal x and meridional y directions, written in terms of grid-scale (resolved) and subgrid-scale velocity components:

\[ \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\end{aligned}\end{align} \]

Turbulent closure schemes are applied to parametrized subgrid-scale vertical fluxes.

  • 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\]
  • The equation of state is given by :

\[\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}\]
  • The final equation expresses the continuity equation. For an incompressible fluid (Boussinesq approximation):

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