Wetting and Drying

Coastal flooding simulations require robust handling of moving shorelines where elements transition between wet and dry states [Medeiros2012] [Carrier1958]. Cocoa implements an element-based wetting/drying algorithm that maintains mass conservation and numerical stability.

Algorithm Overview

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

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

1. Initialize Cycle: Copy committed node status to working status

2. Nodal Drying D1: Mark nodes dry if depth falls below threshold

3. Nodal Wetting W1: Mark dry nodes wet if neighbors can supply water

4. Barrier Wetting: Force-wet dry neighbors of barrier (NIBNODECODE) nodes

5. Elemental Drying DE1: Deactivate elements acting as shallow control sections (see Phase Details)

6. Active Element Count: Count wet elements per node

7. Landlocked Drying D2: Dry nodes with no active elements

8. Status Reconciliation: Commit working status to final status

9. Slope Limiter (ALPHA): Compute pressure gradient limiting factors

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

Fig. 48 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}\]

Nodal Wetting (W1)

A dry node becomes wet when:

1. It has at least one wet neighbor

2. Water can flow from wet to dry (potential gradient favorable)

3. Estimated velocity exceeds minimum:

\[V = \sqrt{g H \frac{\Delta \zeta}{\Delta x}} > V_{min}\]

When a node wets, the TKM friction tensor is updated with TKNF scaling:

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

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 just 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 simply 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.

W2 reactivation, 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. This is the price of an element that states no diagonal, and it is not a defect of either discretization – a fixed-diagonal mesh answers a question the quadrilateral was never asked. 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 effectively disables limiting.

When limiting is enabled (e.g., wet_dry_slope_limiter: 0.0004), the scaling factor is computed based on the ratio of bathymetric gradient to water surface gradient:

(31)\[\alpha_L = \frac{\text{SLIM} \cdot |\nabla h|}{|\nabla \zeta|}\]

The limiter is only evaluated for active elements whose nodes are all wet and that contain shallow water (total depth at or below the shallow-water threshold). This factor multiplies the pressure gradient terms in both GWCE and momentum:

• \(\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:

class WetDryData {
  NodeStatusView m_node_wet_status;      // Committed status
  NodeStatusView m_node_wet_working;     // Working during algorithm

  ElementStatusField m_element_active;   // TemporalField<ElementStatus, 2, 1>
  OrdinalView m_active_element_count;    // Per-node count

  FloatView m_slope_limiter;             // ALPHA values (float storage)
};

Node-to-Element Connectivity

For efficient iteration, DeviceMesh stores CSR (Compressed Sparse Row) format connectivity:

OrdinalView m_node_to_elem_offsets;  // [num_nodes + 1]
OrdinalView m_node_to_elem_list;     // [total_connections]

Implementation Files

• 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)

• WetDryData.hpp: Data container class