Governing Equations

Cocoa solves the 2D depth-integrated shallow water equations [Vreugdenhil1994] [Toro2001].

Derivation

The shallow water equations assume:

1. Horizontal length scales are much larger than vertical scales

2. Vertical accelerations are negligible (hydrostatic pressure)

3. Density is constant (barotropic flow)

Conservation Laws

Mass Conservation

The continuity equation expresses conservation of mass:

(1)\[\frac{\partial \zeta}{\partial t} + \frac{\partial (UH)}{\partial x} + \frac{\partial (VH)}{\partial y} = 0\]

where:

• \(\zeta\) is the water surface elevation above the geoid

• \(H = h + \zeta\) is the total water depth

• \(h\) is the bathymetric depth (positive downward)

• \(U, V\) are the depth-averaged velocity components

Fig. 16 Water column variables: bathymetry \(h\) (positive downward from geoid), surface elevation \(\zeta\), and total depth \(H = h + \zeta\).

Momentum

The momentum equations express Newton’s second law:

x-momentum:

(2)\[\begin{split}\begin{split} \frac{\partial U}{\partial t} + U\frac{\partial U}{\partial x} + V\frac{\partial U}{\partial y} - fV &= -g\frac{\partial \zeta}{\partial x} - \frac{1}{\rho_0}\frac{\partial p_s}{\partial x} \\ &\quad + \frac{\tau_{sx}}{\rho_0 H} - \frac{\tau_{bx}}{\rho_0 H} + \nu_h \nabla^2 U \end{split}\end{split}\]

y-momentum:

(3)\[\begin{split}\begin{split} \frac{\partial V}{\partial t} + U\frac{\partial V}{\partial x} + V\frac{\partial V}{\partial y} + fU &= -g\frac{\partial \zeta}{\partial y} - \frac{1}{\rho_0}\frac{\partial p_s}{\partial y} \\ &\quad + \frac{\tau_{sy}}{\rho_0 H} - \frac{\tau_{by}}{\rho_0 H} + \nu_h \nabla^2 V \end{split}\end{split}\]

Term Descriptions

x-momentum Term	y-momentum Term	Physical Meaning
\(\partial U/\partial t\)	\(\partial V/\partial t\)	Local acceleration
\(U \frac{\partial U}{\partial x} + V \frac{\partial U}{\partial y}\)	\(U \frac{\partial V}{\partial x} + V \frac{\partial V}{\partial y}\)	Advective acceleration
\(-fV\)	\(+fU\)	Coriolis acceleration
\(-g \frac{\partial \zeta}{\partial x}\)	\(-g \frac{\partial \zeta}{\partial y}\)	Barotropic pressure gradient
\(-\frac{1}{\rho_0}\frac{\partial p_s}{\partial x}\)	\(-\frac{1}{\rho_0}\frac{\partial p_s}{\partial y}\)	Atmospheric pressure gradient (see Meteorological Forcing)
\(\tau_{sx}/(\rho_0 H)\)	\(\tau_{sy}/(\rho_0 H)\)	Wind stress (see Meteorological Forcing)
\(-\tau_{bx}/(\rho_0 H)\)	\(-\tau_{by}/(\rho_0 H)\)	Bottom friction
\(\nu_h \nabla^2 U\)	\(\nu_h \nabla^2 V\)	Lateral diffusion

Bottom Friction

Cocoa uses a quadratic bottom friction formulation with Manning’s n:

\[\tau_{bx} = \rho_0 C_f U \sqrt{U^2 + V^2}\]

\[\tau_{by} = \rho_0 C_f V \sqrt{U^2 + V^2}\]

where the friction coefficient \(C_f\) is derived from Manning’s n:

(4)\[C_f = \frac{g n^2}{H^{1/3}}\]

Here \(n\) is the Manning’s roughness coefficient and \(H\) is the total water depth. This formulation accounts for the depth-dependence of bottom friction, with shallower water experiencing stronger frictional effects.

Fig. 17 Manning’s friction coefficient \(C_f\) as a function of depth for several roughness values \(n\). The coefficient increases sharply in shallow water and is floored at \(C_{f,\min} = 0.001\).

Internal Tide Friction

In regions where barotropic tidal energy is converted to internal (baroclinic) tides through interaction with topographic features, the energy loss from the barotropic tide can be parameterized as an additional linear friction term. This internal tide friction is added to the bottom friction TKM tensor after the Manning’s friction has been computed:

\[\tau_{xx} \leftarrow \tau_{xx} + \tau_{IT,xx}, \quad \tau_{yy} \leftarrow \tau_{yy} + \tau_{IT,yy}, \quad \tau_{xy} \leftarrow \tau_{xy} + \tau_{IT,xy}\]

The internal tide friction coefficient is typically derived from tidal energy dissipation estimates (e.g., from satellite altimetry or numerical models) and is specified as a spatially-varying nodal attribute in the mesh file. Two storage modes are supported:

• Scalar (1 value per node): isotropic friction added to the diagonal components (\(\tau_{IT,xx} = \tau_{IT,yy} = \tau_{IT}\), \(\tau_{IT,xy} = 0\))

• Tensor (3 values per node): anisotropic friction with independent \(\tau_{IT,xx}\), \(\tau_{IT,yy}\), and \(\tau_{IT,xy}\) components

For global meshes using coordinate rotation, the tensor form of internal tide friction must be rotated from the geographic frame to the computational frame. This is done once during initialization using a rank-1 tensor decomposition: the tensor is decomposed into a vector \(\mathbf{H} = (H_x, H_y)\), rotated using the per-node velocity transformation matrix, and reconstructed. This matches ADCIRC’s Apply2DInternalWaveDrag rotation for the IFSPROTS=1 case.

See Configuration for configuration details.

Coriolis Parameter

The Coriolis parameter is computed as:

\[f = 2\Omega \sin(\phi)\]

where \(\Omega = 7.292 \times 10^{-5}\) rad/s is Earth’s angular velocity and \(\phi\) is the latitude.