Wetting and Drying

Coastal flooding simulations must track a moving shoreline, with elements passing between wet and dry states [Medeiros2012] [Carrier1958]. Cocoa uses an element-based wetting/drying algorithm that conserves mass and remains numerically stable.

Algorithm Overview

Fig. 48 Wet/dry algorithm phases executed after each GWCE solve

Cocoa uses a 9-phase wet/dry algorithm executed after each GWCE solve:

Initialize cycle: copy committed node status to working status
Nodal drying (D1): mark nodes dry if depth falls below threshold
Nodal wetting (W1): mark dry nodes wet if neighbors can supply water
Barrier wetting: force-wet dry neighbors of barrier (NIBNODECODE) nodes
Elemental drying (DE1): deactivate elements acting as shallow control sections (see Phase Details)
Active element count: count wet elements per node
Landlocked drying (D2): dry nodes with no active elements
Status reconciliation: commit working status to final status
Slope limiter (ALPHA): compute pressure gradient limiting factors

Under MPI, elements have no ghosts, so a rank sees only its own elements’ verdicts on a node it shares with a neighbor. Every element-to-node result is therefore combined over all ranks before any rank reads it (WetDry.cpp), including at cold start:

the W1 wetting verdict and the post-wetting friction boost, by maximum (a node wets, and its friction rises, if any rank holding one of its elements says so);
the barrier node marks, by maximum, and the per-node barrier counts, by sum;
the active-element counts and wet patch areas of phase 6, by sum.

The nodal phases (D1, D2, reconciliation) read only values every copy of a node holds alike, so every copy decides alike and no exchange is needed.

Node Status

Each node has a wet/dry status stored in WetDryData:

enum class NodeStatus : int {
  Dry = 0,
  Wet = 1
};

Two status arrays are maintained:

m_node_wet_status: Committed status used in computations
m_node_wet_working: Working status during algorithm phases

Element Status

Elements have an active/inactive status:

enum class ElementStatus : int {
  Inactive = 0,
  Active = 1
};

Inactive elements are excluded from GWCE and momentum assembly.

Depth Thresholds

The algorithm uses several depth thresholds, which can be configured via the physics section of the configuration file (see Configuration):

Parameter	Default	Config Key	Description
\(H_0\)	0.1 m	`h0`	Minimum depth for wet classification
\(H_{off}\)	0.12 m	(derived)	Wetting threshold (\(1.2 \times H_0\))
\(H_{abs}\)	0.08 m	(derived)	Absolute floor (\(0.8 \times H_0\))
\(V_{min}\)	0.01 m/s	`velmin`	Minimum velocity for wetting

Number line showing wet/dry depth thresholds and hysteresis — Fig. 49 Wet/dry depth thresholds. Nodes dry when \(H \le H_0\) and re-wet when \(H > H_{off}\), creating hysteresis that prevents oscillation at the shoreline.

Phase Details

Nodal Drying (D1)

A node becomes dry when its total depth falls below the threshold:

\[H_i = h_i + \zeta_i \leq H_0 \implies \text{Dry}\]

A node that dries with \(H_i < H_{abs}\) has its elevation raised to \(H_{abs} - h_i\). The same pass records the lowest elevation of every node that was wet on entry, before any of them dries, so the stability check can bound elevation from below even where a runaway drop is hidden by drying.

Nodal Wetting (W1)

A dry node becomes wet when:

It is the one dry corner of an element whose other two corners are wet, and both of those wet corners are at least \(H_{off}\) deep
Water can flow from wet to dry (potential gradient favorable)
The estimated velocity of the flow onto it exceeds \(V_{min}\).

The estimate balances the free-surface slope \(S = \Delta\zeta / \Delta x\) from the wet source node (the wet neighbor with the higher surface, total depth \(H\), Manning drag coefficient \(C_d\) floored at cf_lower_limit) toward the dry node against a resistance chosen by physics.wetting_velocity. On a global mesh the projected-x part of \(\Delta x\) is reduced to its shortest equivalent across the longitude branch cut (Longitude Branch Cut). The three forms are ADCIRC’s wetDryControl settings (physics/WettingFlow.hpp):

Form	Velocity \(V\)	TKNF	ADCIRC equivalent
`linear` (default)	\(\dfrac{g S}{\max(C_d V_{min} / H,\ 10^{-4})}\)	1	default (`directvelWD = F`)
`quadratic`	\(\sqrt{g H S / C_d}\)	1	`directvelWD = T`, `useHF = F`
`inertial`	\(\sqrt{g H S}\)	100	`directvelWD = T`, `useHF = T`

linear holds the friction at its value for a flow of \(V_{min}\), so a thin, rough source still wets readily; quadratic solves the full quadratic friction balance; inertial drops friction and takes the free shallow-water speed, the most reluctant of the three on a gentle slope. linear is ADCIRC’s default; inertial is the form Global STOFS runs.

When a node wets, the friction at its two wet neighbors is raised to damp the spurious velocity that follows:

\[\tau = \min\left(\frac{TKNF \cdot C_d V}{H}, \frac{2}{\Delta t}\right)\]

The 100-fold boost of inertial compensates for a velocity that carries no friction of its own.

Elemental Drying (DE1)

DE1 deactivates elements that act as shallow control sections. Ordering the element’s nodal elevations \(\zeta_A \geq \zeta_B > \zeta_C\), the element is deactivated when either of the two higher-elevation nodes is shallower than the wetting threshold:

\[H_A < H_{off} \quad \text{or} \quad H_B < H_{off} \implies \text{Inactive}\]

Elements touching a node with a barrier contribution (NIBCNT > 0) are exempt, so weir overflow can keep its receiving elements active.

Active Element Count

For each node, count elements that are:

Active (not deactivated by DE1), AND
Have all nodes wet (three on a triangle, four on a quadrilateral)

This determines if a node is “landlocked.”

Landlocked Drying (D2)

A node with zero active elements is marked dry, even if its depth suggests wet. This prevents isolated wet nodes from persisting.

Quadrilateral Rules (W1, DE1)

The rules (numeric/cg/wetdry/QuadNodeSelection.hpp) are stated once, and both are the triangle rule [Luettich1999] [Dietrich2004] read over the four half-triangles of the quadrilateral’s two diagonal splits. The element states no diagonal, so deferring to one would make the answer depend on a choice the mesh never made; the four halves are also the four ways to omit one corner, which is how the code enumerates them.

W1 (nodal wetting). Every half-triangle with exactly one dry corner is a candidate triple, and each runs the identical triangle wetting test (detail::apply_wetting_candidate) on its (dry, wet, wet) vertices. The pair a triple names is the pair that half contains: with the dry corner’s diagonal omitted it is the two edge neighbors, in the order (d-1, d+1); with one edge neighbor omitted it is the surviving edge neighbor and the diagonal corner. Three wet corners give three triples, all for the one dry corner d, naming (d-1, d+1), (d+1, d+2) and (d-1, d+2); the diagonal two-wet pattern gives one triple per dry corner, and so does the adjacent two-wet pattern; every other pattern gives none. Barrier friction boosting applies to the same two wet corners a triple names, exactly as on a triangle. The wetting itself is a compare-and-swap on the dry node’s status, so a corner named by several triples is tested several times.

DE1 (elemental drying). The barrier gate runs first, exactly as on triangles. Then each half-triangle is judged by the triangle rule itself (its lowest corner by elevation is a unique minimum and one of the other two is below the wetting threshold), and the quadrilateral dries when all four halves would. A half that still holds water keeps the whole element active, so the element never landlocks a corner that the mesh’s own split would have kept active. W1 admits what any half admits and DE1 keeps while any half keeps: one reading of the same four halves, in the same direction.

The wet flag, the slope limiter and the active-element counts read all four corners: each is the three-corner rule generalized to \(N = 4\) with no other change.

Where this differs from a mesh split along a fixed diagonal. For W1 the set of dry corners admitted differs in five of the sixteen wet/dry patterns: the four adjacent two-wet patterns, and the diagonal two-wet pattern whose wet corners are not the split’s own diagonal. In each of those the quadrilateral admits both dry corners where the fixed split admits one or none. The other eleven patterns agree. For DE1 the answers differ wherever the two halves of the fixed split disagree with each other: the split can deactivate one half and keep the other, and a quadrilateral has no half to deactivate alone, so it stays active. Both differences follow from an element that states no diagonal; neither discretization is wrong. The unit tests (test_wetdry_quads.cpp) assert agreement on the patterns where the rules predict agreement and divergence on the patterns where they predict that, with the candidate triples of all sixteen patterns pinned by an explicit table.

The sloping-beach validation case measures what the two rules produce together against ADCIRC on the same node set (see the front study): on a beach whose shoreline moves 1035 m per half cycle, a quadrilateral mesh tracks ADCIRC’s front to 0.01% in speed and 49 m in position, against 0.06% and 75 m for the triangle mesh.

Slope Limiter (ALPHA)

The slope limiter \(\alpha_L\) reduces pressure gradients in shallow water to prevent spurious flows. The limiter is controlled by the wet_dry_slope_limiter configuration parameter (ADCIRC’s SLIM), with a default of \(10^9\), which never limits in practice.

When limiting is enabled (e.g., wet_dry_slope_limiter: 0.0004), the scaling factor is the ratio of SLIM to the element’s water-surface slope. Where the bathymetry has a defined gradient (\(|\nabla h| > 10^{-5}\)), the surface slope is measured along the downslope bathymetric direction; otherwise it is the full surface gradient:

(31)\[\begin{split}\alpha_L = \begin{cases} \dfrac{\text{SLIM}\,|\nabla h|}{-\nabla\zeta \cdot \nabla h} & |\nabla h| > 10^{-5} \\[8pt] \dfrac{\text{SLIM}}{|\nabla \zeta|} & \text{otherwise} \end{cases}\end{split}\]

and \(\alpha_L = 1\) whenever this value falls outside \((0, 1)\). The limiter is only evaluated for active elements whose nodes are all wet and with at least one node of total depth at or below 0.5 m; every other element keeps \(\alpha_L = 1\). This factor multiplies the free-surface gradient terms in both GWCE and momentum (in the GWCE also the atmospheric pressure and tide potential gradients):

\(\alpha_L = 1\): Full pressure gradient (deep water or limiting disabled)
\(\alpha_L < 1\): Reduced pressure gradient (shallow water with limiting)

Data Structures

WetDryData

The WetDryData class stores all wet/dry state (abridged; the per-element members hold one view per element type):

class WetDryData {
  NodeStatusView m_node_wet_status;      // Committed status
  NodeStatusView m_node_wet_working;     // Working during algorithm
  OrdinalView m_active_element_count;    // Per-node count
  DoubleView m_active_element_area;      // Per-node wet patch area

  PerElementType<ElementStatusField> m_active;    // Active/Inactive, 2 levels
  PerElementType<FloatView> m_slope_limiter;      // ALPHA values (float storage)
  PerElementType<FlagView> m_wet_flag;            // Active and all corners wet
  PerElementType<OrdinalView> m_wet_indices;      // Compacted wet-element list

  ElevationScalarView m_min_wet_elevation;        // Lowest wet elevation (D1)
};

Element Iteration

Every element phase is an element-parallel kernel that scatters to its corners with atomics; no phase walks a node’s element list. (DeviceMesh keeps node-to-element CSR connectivity on the host only, for sizing the GWCE matrix rows at construction.)

Implementation Files

All in numeric/cg/wetdry/ unless noted:

WetDry.hpp: Algorithm orchestrator with static compute() method
WetDryKernels.hpp: Individual phase kernels. The nodal phases take no element type; each element phase is one function templated on the element-type tag, with a definition per type in WetDryTriangleKernels.cpp and WetDryQuadrilateralKernels.cpp (the corner rules of the two are different physics, not the same body at two vertex counts)
NodeSelection.hpp, QuadNodeSelection.hpp: Triangle and quadrilateral corner rules for W1 and DE1, and the shared wetting test
physics/WettingFlow.hpp: The three wetting-velocity forms
data/WetDryData.hpp: Data container class