Momentum Equations
Cocoa solves the non-conservative form of the 2D depth-integrated momentum
equations in spherical coordinates.
Governing Equations
The momentum equations in spherical coordinates on an equidistant cylindrical
projection (\(\lambda, \phi\) for longitude and latitude):
\[\begin{split}\begin{split}
\frac{\partial U}{\partial t}
&+ \frac{U}{R\cos\phi}\frac{\partial U}{\partial \lambda}
+ \frac{V}{R}\frac{\partial U}{\partial \phi}
- f_{\text{eff}} V \\
&= -\frac{g}{R\cos\phi}\frac{\partial \zeta}{\partial \lambda}
+ \frac{\tau_{sx} - \tau_{bx}}{\rho_0 H} + M_x
\end{split}\end{split}\]
\[\begin{split}\begin{split}
\frac{\partial V}{\partial t}
&+ \frac{U}{R\cos\phi}\frac{\partial V}{\partial \lambda}
+ \frac{V}{R}\frac{\partial V}{\partial \phi}
+ f_{\text{eff}} U \\
&= -\frac{g}{R}\frac{\partial \zeta}{\partial \phi}
+ \frac{\tau_{sy} - \tau_{by}}{\rho_0 H} + M_y
\end{split}\end{split}\]
where \(M_x\) and \(M_y\) represent lateral stress terms, and the
effective Coriolis parameter includes a spherical metric correction:
\[f_{\text{eff}} = 2\Omega\sin\phi + \frac{U \tan\phi}{R}\]
The second term accounts for the convergence of meridians and is evaluated
using the velocity from the previous timestep.
Coordinate Systems
Cocoa supports multiple cylindrical map projections for spherical coordinates.
The default is the equidistant cylindrical (CPP) projection, corresponding to
ADCIRC’s ICS=21 setting. See Coordinate Systems
Key points for momentum assembly:
The mesh precomputes shape function derivatives in projected coordinates
During assembly, x-derivatives are multiplied by the momentum scale factor
\(S_x^{mom}\) (stored in elem_scale_momentum_x
The y-derivatives require no correction for the equidistant cylindrical
projection since latitude spacing is uniform
The effective Coriolis parameter includes a spherical metric correction
using \(\tan\phi/R\) (see Equidistant Cylindrical (Default)
Time Discretization
Bottom friction and Coriolis use Crank-Nicolson weighting (\(\theta = 0.5\) ).
Writing the x-momentum equation:
\[\begin{split}\begin{split}
\frac{U^{n+1} - U^n}{\Delta t}
- f_{\text{eff}}\left(\frac{V^{n+1} + V^n}{2}\right)
&= -\tau_{xx}\left(\frac{U^{n+1} + U^n}{2}\right) \\
&\quad - \tau_{xy}\left(\frac{V^{n+1} + V^n}{2}\right) + F_x
\end{split}\end{split}\]
Rearranging:
\[\begin{split}\begin{split}
\left(1 + \frac{\Delta t}{2}\tau_{xx}\right) U^{n+1}
&+ \frac{\Delta t}{2}\left(\tau_{xy} - f_{\text{eff}}\right) V^{n+1} \\
&= \left(1 - \frac{\Delta t}{2}\tau_{xx}\right) U^n
- \frac{\Delta t}{2}\left(\tau_{xy} - f_{\text{eff}}\right) V^n
+ \Delta t \cdot F_x
\end{split}\end{split}\]
Per-Node 2x2 System
digraph momentum {
rankdir=LR;
node [shape=box, style="filled,rounded", fontname="Helvetica", fontsize=10];
edge [color="#555555"];
bgcolor="transparent";
// RHS contributions (inputs)
subgraph cluster_rhs {
label="RHS Contributions (element loop)";
style="filled,rounded";
fillcolor="#F5F5F5";
color="#BDBDBD";
fontname="Helvetica";
fontsize=10;
pres [label="Pressure\nGradient", fillcolor="#E3F2FD", color="#1E88E5"];
adv [label="Advection", fillcolor="#E3F2FD", color="#1E88E5"];
lat [label="Lateral\nStress", fillcolor="#E3F2FD", color="#1E88E5"];
tipg [label="Tide Potential\nGradient", fillcolor="#E3F2FD", color="#1E88E5"];
wind [label="Wind Stress\n(nodal)", fillcolor="#E3F2FD", color="#1E88E5"];
}
// LHS contributions
subgraph cluster_lhs {
label="LHS Matrix (per node)";
style="filled,rounded";
fillcolor="#F5F5F5";
color="#BDBDBD";
fontname="Helvetica";
fontsize=10;
fric [label="Bottom Friction\n(TKM diagonal)", fillcolor="#FFF3E0", color="#FB8C00"];
cor [label="Coriolis\n(off-diagonal)", fillcolor="#FFF3E0", color="#FB8C00"];
}
// Assembly
assemble [label="Assemble 2x2\nSystem per Node", fillcolor="#E8F5E9", color="#43A047"];
// Solve
solve [label="Cramer's Rule\nDirect Solve", fillcolor="#FFEBEE", color="#E53935"];
// Output
vel [label="U, V\n(velocity)", fillcolor="#F3E5F5", color="#8E24AA"];
flux [label="Qx, Qy\n(flux = H * U)", fillcolor="#F3E5F5", color="#8E24AA"];
pres -> assemble;
adv -> assemble;
lat -> assemble;
tipg -> assemble;
wind -> assemble;
fric -> assemble;
cor -> assemble;
assemble -> solve -> vel;
solve -> flux;
}
Fig. 21 Momentum solver: per-node 2x2 system assembly and solve
The momentum solve at each node:
(9) \[\begin{split}\begin{bmatrix}
1 + \frac{\Delta t}{2}\tau_{xx} & \frac{\Delta t}{2}(\tau_{xy} - f_{\text{eff}}) \\[6pt]
\frac{\Delta t}{2}(\tau_{xy} + f_{\text{eff}}) & 1 + \frac{\Delta t}{2}\tau_{yy}
\end{bmatrix}
\begin{bmatrix} U^{n+1} \\ V^{n+1} \end{bmatrix}
=
\begin{bmatrix} R_x \\ R_y \end{bmatrix}\end{split}\]
The right-hand side:
\[\begin{split}\begin{split}
R_x = \left(1 - \frac{\Delta t}{2}\tau_{xx}\right) U^n
&- \frac{\Delta t}{2}(\tau_{xy} - f_{\text{eff}}) V^n \\
&+ \Delta t \cdot F_x
\end{split}\end{split}\]
\[\begin{split}\begin{split}
R_y = \left(1 - \frac{\Delta t}{2}\tau_{yy}\right) V^n
&- \frac{\Delta t}{2}(\tau_{xy} + f_{\text{eff}}) U^n \\
&+ \Delta t \cdot F_y
\end{split}\end{split}\]
Coriolis appears only in the off-diagonal entries of Equation (9) .
The system is solved via Cramer’s rule, keeping all four matrix entries in GPU
registers.
Bottom Friction (TKM Tensor)
The TKM tensor stores the linearized friction coefficient. For Manning’s n
(see (4) ):
(10) \[C_d = \frac{g n^2}{H^{1/3}}\]
(11) \[\tau = \min\left(\frac{C_d |\mathbf{U}|}{H}, \frac{2}{\Delta t}\right)\]
The upper bound \(2/\Delta t\) prevents the implicit friction term from
overwhelming the solution in very shallow water.
For isotropic friction, \(\tau_{xx} = \tau_{yy} = \tau\) and
\(\tau_{xy} = 0\) . When internal tide friction is enabled, its
contribution is added to the TKM tensor after Manning’s friction is computed
(see Internal Tide Friction \(\tau_{xy} \neq 0\) and the full
\(2 \times 2\) system in Equation (9) is required.
These values are stored directly in PhysicsData
RHS Forcing Terms (Continuous)
The forcing terms \(F_x\) and \(F_y\) include:
Pressure Gradient (time-centered):
\[\begin{split}\begin{aligned}
F_x^{pres} &= -\frac{g}{2} \alpha_L
\left(\frac{\partial \zeta^{n-1}}{\partial x}
+ \frac{\partial \zeta^{n+1}}{\partial x}\right) \\[4pt]
F_y^{pres} &= -\frac{g}{2} \alpha_L
\left(\frac{\partial \zeta^{n-1}}{\partial y}
+ \frac{\partial \zeta^{n+1}}{\partial y}\right)
\end{aligned}\end{split}\]
where \(\alpha_L\) is the slope limiter from wet/dry processing
(Equation (24) ) that reduces pressure gradients in shallow water.
Tide Potential Gradient:
\[\begin{split}\begin{aligned}
F_x^{tip} &= \frac{g}{2} \alpha_L
\left(\frac{\partial \eta_{tip}^{n-1}}{\partial x}
+ \frac{\partial \eta_{tip}^{n+1}}{\partial x}\right) \\[4pt]
F_y^{tip} &= \frac{g}{2} \alpha_L
\left(\frac{\partial \eta_{tip}^{n-1}}{\partial y}
+ \frac{\partial \eta_{tip}^{n+1}}{\partial y}\right)
\end{aligned}\end{split}\]
The tide potential acts as a body force opposing the free surface gradient.
Advection (non-conservative form, time-weighted):
\[F_x^{adv} = -\left(\bar{U} \frac{\partial U}{\partial x}
+ \bar{V} \frac{\partial U}{\partial y}\right), \quad
F_y^{adv} = -\left(\bar{U} \frac{\partial V}{\partial x}
+ \bar{V} \frac{\partial V}{\partial y}\right)\]
Cocoa uses 2-level time centering with weights between \(n-1\) and \(n\) .
Wind Stress (nodal application):
\[F_x^{wind} = \frac{f_w \, \tau_{sx}}{\rho_0 H}, \quad
F_y^{wind} = \frac{f_w \, \tau_{sy}}{\rho_0 H}\]
where \(f_w\) is the depth-dependent wind stress limiter that smoothly
ramps wind forcing to zero near the drying threshold. Wind stress is computed
from 10-meter wind velocity using the Garratt drag law and applied per node
(not per element). See Meteorological Forcing
Lateral Stress Divergence:
\[F_x^{ls} = \frac{1}{H}\left(\frac{\partial \sigma_{xx}}{\partial x}
+ \frac{\partial \sigma_{xy}}{\partial y}\right), \quad
F_y^{ls} = \frac{1}{H}\left(\frac{\partial \sigma_{yx}}{\partial x}
+ \frac{\partial \sigma_{yy}}{\partial y}\right)\]
The stress tensor uses the Smagorinsky eddy viscosity model.
Discretized RHS Assembly
The forcing terms are discretized using linear triangular finite elements.
For an element with nodes \((1, 2, 3)\) , area \(A_e\) , and shape function
derivatives stored as \(FDX_i = 2 A_e \frac{\partial N_i}{\partial x}\) .
Pressure Gradient (Discretized):
\[\begin{split}\begin{aligned}
F_x^{pres,e} &= -\frac{g \Delta t}{4} \alpha_L
\sum_{i=1}^{3} (\zeta_i^{n-1} + \zeta_i^{n+1}) \cdot FDX_i \\[4pt]
F_y^{pres,e} &= -\frac{g \Delta t}{4} \alpha_L
\sum_{i=1}^{3} (\zeta_i^{n-1} + \zeta_i^{n+1}) \cdot FDY_i
\end{aligned}\end{split}\]
Note: The \(\Delta t\) factor is included because the RHS represents
\(\Delta t \cdot F\) , not \(F\) alone. The factor of \(1/4\)
combined with \(FDX = 2A_e \partial N/\partial x\) (and similarly
\(FDY = 2A_e \partial N/\partial y\) ) yields the continuous
\(g/2\) coefficient for the time-centered average.
Advection (Discretized):
The velocity gradients are computed as:
\[\frac{\partial U}{\partial x} = \frac{1}{2A_e}
\sum_{i=1}^{3} U_i \cdot FDX_i, \quad
\frac{\partial V}{\partial x} = \frac{1}{2A_e}
\sum_{i=1}^{3} V_i \cdot FDX_i\]
and similarly for the y-derivatives using \(FDY_i\) .
The advection contribution uses element-averaged velocities:
\[\begin{split}\begin{aligned}
F_x^{adv,e} &= -\Delta t \left(\bar{U} \frac{\partial U}{\partial x}
+ \bar{V} \frac{\partial U}{\partial y}\right) \\[4pt]
F_y^{adv,e} &= -\Delta t \left(\bar{U} \frac{\partial V}{\partial x}
+ \bar{V} \frac{\partial V}{\partial y}\right)
\end{aligned}\end{split}\]
where \(\bar{U} = \frac{1}{3} \, \sum_{i=1}^{3} U_i\) . Cocoa applies time-weighting
between levels \(n-1\) and \(n\) with coefficients matching ADCIRC’s
2-level centered scheme.
Lateral Stress (Discretized):
The Smagorinsky eddy viscosity is computed per-element using the deviatoric strain rate
(matching ADCIRC’s formulation):
\[\nu_e = C_s \cdot A_e \cdot |S|, \quad
|S| = \sqrt{
\left(\frac{\partial U}{\partial x}
- \frac{\partial V}{\partial y}\right)^{\!2}
+ \left(\frac{\partial U}{\partial y}
+ \frac{\partial V}{\partial x}\right)^{\!2}
}\]
The stress tensor uses per-node depth \(H_i\) . Cocoa uses ADCIRC’s
symmetric IBPSV formulation (CME_LS_IBPSV
\[\sigma_{xx,i} = \nu_e H_i \frac{\partial U}{\partial x}, \quad
\sigma_{yy,i} = \nu_e H_i \frac{\partial V}{\partial y}\]
\[\sigma_{xy,i} = \sigma_{yx,i}
= \tfrac{1}{2} \nu_e H_i
\left(\frac{\partial U}{\partial y} + \frac{\partial V}{\partial x}\right)\]
The lateral stress divergence contribution (symmetric IBPSV formulation):
\[F_x^{ls,i} = -\frac{\Delta t}{H_i} \cdot 1.5 \cdot
(\sigma_{xx,i} \cdot FDX_i + \sigma_{xy,i} \cdot FDY_i)\]
\[F_y^{ls,i} = -\frac{\Delta t}{H_i} \cdot 1.5 \cdot
(\sigma_{yx,i} \cdot FDX_i + \sigma_{yy,i} \cdot FDY_i)\]
Element-to-Node Assembly:
Element contributions are accumulated to nodes using atomic operations (shown for
x-momentum; y-momentum follows the same pattern with \(F_y\) terms):
\[\begin{split}\begin{aligned}
RHS_{x,i} &\mathrel{+}= F_x^{pres,e} + F_x^{adv,e} + F_x^{ls,i} \\
RHS_{y,i} &\mathrel{+}= F_y^{pres,e} + F_y^{adv,e} + F_y^{ls,i}
\end{aligned}\end{split}\]
After assembly, the accumulated RHS is normalized by the node’s control volume:
\[RHS_{x,i} \leftarrow \frac{RHS_{x,i}}{A_{total,i}}, \quad
RHS_{y,i} \leftarrow \frac{RHS_{y,i}}{A_{total,i}}\]
where \(A_{total,i} = \sum_{e \ni i} A_e/3\) is the sum of one-third of each
element area surrounding node \(i\)
Land Boundary Conditions
At land boundaries, Cocoa enforces a no-normal-flow condition:
\[n_x U + n_y V = 0\]
This is implemented by transforming the \(2 \times 2\) system to tangential-normal coordinates:
Transform RHS: \((RHS_t, RHS_n) = R \cdot (RHS_x, RHS_y)\)
Transform matrix to rotated coordinates
Set \(RHS_n = 0\) (no normal flow)
Solve transformed system
Transform back to \((U, V)\)
The transformation uses:
Flux Computation
After solving for velocity, the flux is computed:
\[Q_x = H \cdot U, \quad Q_y = H \cdot V\]
where \(H = h + \zeta\) is the total water depth.
A minimum depth threshold is applied to avoid division issues in very shallow water.
Implementation Files
MomentumSolver.hpp
MomentumSolveKernels.hpp \(2 \times 2\) solve kernel
MomentumRhsKernels.hpp
BottomFriction.hpp
LateralStress.hpp