===================
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
------------------

.. graphviz::
   :align: center
   :caption: Wet/dry algorithm phases executed after each GWCE solve

   digraph wetdry {
     rankdir=TB;
     node [shape=box, style="filled,rounded", fontname="Helvetica", fontsize=10];
     edge [color="#555555"];
     bgcolor="transparent";

     init [label="1. Initialize Cycle\nCopy committed -> working status", fillcolor="#E0E0E0", color="#616161"];
     d1 [label="2. Nodal Drying (D1)\nH < H0 -> Dry", fillcolor="#FFEBEE", color="#E53935"];
     w1 [label="3. Nodal Wetting (W1)\nWet neighbor + flow potential\n-> Wet", fillcolor="#E3F2FD", color="#1E88E5"];
     barrier [label="4. Barrier Wetting\nForce-wet barrier neighbors", fillcolor="#E3F2FD", color="#1E88E5"];
     de1 [label="5. Elemental Drying (DE1)\nShallow control section -> Inactive", fillcolor="#FFEBEE", color="#E53935"];
     count [label="6. Active Element Count\nPer-node wet element tally", fillcolor="#FFF3E0", color="#FB8C00"];
     d2 [label="7. Landlocked Drying (D2)\nZero active elements -> Dry", fillcolor="#FFEBEE", color="#E53935"];
     commit [label="8. Status Reconciliation\nWorking -> committed", fillcolor="#E8F5E9", color="#43A047"];
     alpha [label="9. Slope Limiter\nCompute pressure gradient\nlimiting factors", fillcolor="#F3E5F5", color="#8E24AA"];

     init -> d1 -> w1 -> barrier -> de1 -> count -> d2 -> commit -> alpha;
   }

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``:

.. code-block:: cpp

   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:

.. code-block:: cpp

   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
:doc:`../getting_started/configuration`):

.. list-table::
   :header-rows: 1

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

.. figure:: ../_static/images/diagrams/depth_thresholds.svg
   :alt: Number line showing wet/dry depth thresholds and hysteresis
   :width: 550px
   :align: center

   Wet/dry depth thresholds. Nodes dry when :math:`H \le H_0` and re-wet when
   :math:`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:

.. math::

   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:

.. math::

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

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

.. math::

   \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 :math:`\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:

.. math::

   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.

.. _quad-wetdry-rules:

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
:math:`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
:ref:`the front study <channel-validation-wetdry>`): 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 :math:`\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 :math:`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:

.. math::
   :label: slope-limiter

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

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

Data Structures
---------------

WetDryData
^^^^^^^^^^

The ``WetDryData`` class stores all wet/dry state:

.. code-block:: cpp

   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:

.. code-block:: cpp

   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