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

The default projection is equidistant cylindrical (CPP), ADCIRC’s ICS=21 setting; Coordinate Systems covers the others and their scale factors. In the 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

Fig. 40 Momentum solver: per-node 2x2 system assembly and solve

The momentum solve at each node:

(14)\[\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 (14). 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 (9)):

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

(16)\[\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. \(C_d\) is floored at cf_lower_limit (default 0.001), \(H\) at min_depth_friction, and the speed \(|\mathbf{U}|\) is formed in double precision from the time-level \(n\) velocity (stored as float) and floored at min_velocity_friction, following ADCIRC’s Apply2DBottomFriction.

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 in the governing equations). With tensor internal tide friction, \(\tau_{xy} \neq 0\) and the full \(2 \times 2\) system in Equation (14) is required.

These values are stored directly in PhysicsData and used as shown in the matrix above.

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}}{\partial x} + \frac{\partial \zeta^{n+1}}{\partial x}\right) \\[4pt] F_y^{pres} &= -\frac{g}{2} \alpha_L \left(\frac{\partial \zeta^{n}}{\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 (31)) that reduces pressure gradients in shallow water.

Tide Potential Gradient:

\[\begin{split}\begin{aligned} F_x^{tip} &= \frac{g}{2} \left(\frac{\partial \eta_{tip}^{n}}{\partial x} + \frac{\partial \eta_{tip}^{n+1}}{\partial x}\right) \\[4pt] F_y^{tip} &= \frac{g}{2} \left(\frac{\partial \eta_{tip}^{n}}{\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)\]

The advecting velocity and its gradients are evaluated explicitly at time level \(n\), using ADCIRC’s half-scaled gradient convention.

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 for the drag law formulation, wind limiter definition, and discretized form.

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 [Smagorinsky1963].

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} + \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} + \zeta_i^{n+1}) \cdot FDY_i \end{aligned}\end{split}\]

The \(\Delta t\) factor appears 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\), with all velocities taken at time level \(n\).

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) [Luettich2004], so the off-diagonal shear components are equal:

\[\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 nodal stresses are formed in two passes, following ADCIRC. The element pass scatters each wet element’s \(\nu_e\) times its velocity gradients to its corners, weighted by area; the node pass divides each node’s sum by its wet patch area \(A_{total,i}\) (below) and multiplies by \(H_i\). Where ADCIRC divides by the static total area, Cocoa divides by the wet patch so that dry neighbors do not dilute the stress. Under MPI the element sums are added over all ranks before the division, so a node on a partition cut sees every element around it.

ADCIRC’s default CME_LS_IBPV is asymmetric, \(\sigma_{xy} = \nu H\,\partial U/\partial y\) and \(\sigma_{yx} = \nu H\,\partial V/\partial x\). The two forms agree wherever the flow is nearly unidirectional and part where there is cross-stream shear, so a cross-implementation comparison has to match the form anywhere there is shear to resolve: a bend, and the entry corner of a flux boundary, both measured under Ideal Channel Validation.

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\), summed over the active elements around node \(i\) whose corners are all wet, is the area of the node’s wet patch (active_element_area, built in the wet/dry pass and summed over ranks). Each element term above carries one factor of \(A_e\), so the division yields the area-weighted average of \(\Delta t\,F\) over the patch. Wind stress is added after this normalization (it is a nodal body force).

On quadrilaterals. The pressure gradient, advection, and lateral-stress RHS terms above generalize with \(N = 4\), the quadrilateral Green coefficients in place of \(FDX_i\)/\(FDY_i\), and the nodal mean \(\tfrac{1}{4}\sum_{i=1}^{4}\) in place of \(\bar{U} = \tfrac{1}{3}\sum_{i=1}^{3} U_i\) and similar averages. momentum_rhs_element_kernel<Quadrilateral>, the tag-templated launcher the triangles also go through, computes them and accumulates into the same \(RHS_{x,i}\), \(RHS_{y,i}\) atomics, launched immediately after the triangle instantiation. The factor 1.5 in \(F^{ls}_{x,i}\), \(F^{ls}_{y,i}\) above is the same constant for both element types (the lumping convention). The control-volume normalization is where the two types differ: \(A_{total,i}\) is \(\sum_{e \ni i} 3\int_{\Omega^e}\phi_i\,d\Omega\), which is \(A_e\) per triangle but three times the lumped mass row, \(3\,ML_i\) (about \(\tfrac{3}{4}A_e\), unequal between corners), per quadrilateral corner, so on a hybrid mesh it is not a plain area sum at any node both types touch. The quadrilateral’s element-uniform terms are rescaled per corner from its area to that weight, so the division still recovers an average.

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:

1. Transform RHS: \((RHS_t, RHS_n) = R \cdot (RHS_x, RHS_y)\)

2. Transform matrix to rotated coordinates

3. Set \(RHS_n = 0\) (no normal flow)

4. Solve transformed system

5. Transform back to \((U, V)\)

The transformation uses the normalized physical land normal, not the raw projected one:

• \(\cos\theta = \text{SFMX} \cdot n_x / |\mathbf{n}_{phys}|\) (CSII in ADCIRC notation)

• \(\sin\theta = \text{SFMY} \cdot n_y / |\mathbf{n}_{phys}|\) (SIII in ADCIRC notation)

with \((n_x, n_y)\) the projected chord normal and \(|\mathbf{n}_{phys}|\) the norm of \((\text{SFMX}\,n_x, \text{SFMY}\,n_y)\). The momentum solve zeroes \(u \cdot n\) using physical velocity components, so it needs the physical normal: on an equidistant map \(\text{SFMY} = 1\) and this reduces to ADCIRC’s stored \((\text{SFMX}\,n_x, n_y)\), but on Mercator \(\text{SFMX} = \text{SFMY}\) so the physical normal is the raw projected one, and on an equal-area map neither identity holds. ADCIRC’s stored form rotates a 45-degree wall by 0.65 degrees at 29 N under Mercator.

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.

The flux is not floored; a dry node is written as zero velocity and zero flux instead. Where the momentum path divides by depth (the Manning friction coefficient, (15) and (16), and the wind stress term), \(H\) is floored at physics.min_depth_friction (default 0.1 m) first.

Implementation Files

• MomentumSolver.hpp: High-level solver class

• MomentumSolveKernels.hpp: Per-node \(2 \times 2\) solve kernel

• MomentumRhsKernels.{hpp,cpp}: RHS assembly kernels, one tag-templated element launcher instantiated per element type into the same nodal accumulators (Quadrilateral Elements)

• BottomFriction.hpp: TKM tensor computation

• LateralStress.hpp: Smagorinsky turbulence model