Scientific Background

This section provides a comprehensive explanation of the scientific principles, mathematical foundations, and computational algorithms underlying the ADCIRC subgrid method.

Theoretical Foundation

Historical Development

Subgrid corrections in hydrodynamic modeling evolved from the need to represent high-resolution topographic features on computationally efficient coarse meshes:

Early Development: Defina (2000) introduced subgrid corrections for partially wet computational cells and advection enhancements in tidal models
Finite Volume Extensions: Casulli (2009) improved mass balance and wetting/drying in unstructured finite volume models
Bottom Friction Improvements: Volp et al. (2013) incorporated spatial variability in bottom roughness and depth
Formalization: Kennedy et al. (2019) formalized closure approximations and added water surface gradient corrections
Storm Surge Integration: First comprehensive implementation in hurricane storm surge models (ADCIRC) achieved 10-50× computational speedup while improving accuracy

Scale Separation Problem

The fundamental challenge in coastal hydrodynamic modeling is the scale separation between:

Model scale: Tens to hundreds of meters (mesh resolution)
Subgrid scale: 1 meter or smaller (DEM resolution)
Feature scale: Critical topographic features (channels, barriers, dunes)

Computational Trade-offs:

High-resolution uniform meshes: Accurate but computationally prohibitive
Coarse meshes: Fast but lose critical flow pathways and barriers
Subgrid solution: Use coarse computational mesh with fine-scale topographic corrections

Subgrid-Scale Modeling Principles

Subgrid modeling addresses scale separation by:

Acknowledging sub-mesh heterogeneity: Each computational cell contains unresolved topographic complexity
Parameterizing sub-mesh effects: Correction terms capture aggregate impacts of unresolved features
Preserving computational efficiency: Avoids prohibitive computational costs of uniform high resolution
Improving accuracy-efficiency trade-offs: Achieves better accuracy on coarse meshes than conventional fine meshes

Performance Characteristics:

Typical coarsening factor: 3-20× (up to 50× in some applications)
Computational savings: 1-2 orders of magnitude reduction in runtime
Accuracy improvement: Subgrid coarse models often outperform conventional fine models

The Shallow Water Equations

The standard 2D shallow water equations form the basis for ADCIRC:

Continuity Equation:

\[\frac{\partial \zeta}{\partial t} + \frac{\partial (uH)}{\partial x} + \frac{\partial (vH)}{\partial y} = 0\]

Momentum Equations:

\[\frac{\partial u}{\partial t} + u\frac{\partial u}{\partial x} + v\frac{\partial u}{\partial y} - fv = -g\frac{\partial \zeta}{\partial x} - \frac{\tau_{bx}}{\rho H}\]

\[\frac{\partial v}{\partial t} + u\frac{\partial v}{\partial x} + v\frac{\partial v}{\partial y} + fu = -g\frac{\partial \zeta}{\partial y} - \frac{\tau_{by}}{\rho H}\]

Where:

\(\zeta\) = water surface elevation
\(u, v\) = depth-averaged velocities
\(H = h + \zeta\) = total water depth
\(h\) = bathymetric depth
\(f\) = Coriolis parameter
\(g\) = gravitational acceleration
\(\tau_b\) = bottom stress
\(\rho\) = water density

Subgrid-Modified Equations

The subgrid method uses averaged versions of the shallow water equations with closure terms.

Averaged Continuity Equation:

\[\phi \frac{\partial \langle \zeta \rangle_W}{\partial t} + \frac{\partial \langle UH \rangle_G}{\partial x} + \frac{\partial \langle VH \rangle_G}{\partial y} = 0\]

Averaged Conservative Momentum Equations:

x-direction:

\[\frac{\partial \langle UH \rangle_G}{\partial t} + \frac{\partial C_{UU} \langle U \rangle \langle UH \rangle_G}{\partial x} + \frac{\partial C_{VU} \langle V \rangle \langle UH \rangle_G}{\partial y}\]

\[+ gC_\zeta \langle H \rangle_G \frac{\partial \langle \zeta \rangle_W}{\partial x} = -\frac{C_{M,f} \langle U \rangle \langle UH \rangle_G}{\langle H \rangle_W} - f \langle VH \rangle_G\]

\[- g \langle H \rangle_G \frac{\partial P_A}{\partial x} + \frac{\tau_{sx}}{\rho_0 \phi} + \frac{\partial}{\partial x}\left(E_h \frac{\partial \langle UH \rangle_G}{\partial x}\right) + \frac{\partial}{\partial y}\left(E_h \frac{\partial \langle UH \rangle_G}{\partial y}\right)\]

y-direction:

\[\frac{\partial \langle VH \rangle_G}{\partial t} + \frac{\partial C_{UV} \langle U \rangle \langle VH \rangle_G}{\partial x} + \frac{\partial C_{VV} \langle V \rangle \langle VH \rangle_G}{\partial y}\]

\[+ gC_\zeta \langle H \rangle_G \frac{\partial \langle \zeta \rangle_W}{\partial y} = -\frac{C_{M,f} \langle V \rangle \langle VH \rangle_G}{\langle H \rangle_W} + f \langle UH \rangle_G\]

\[- g \langle H \rangle_G \frac{\partial P_A}{\partial y} + \frac{\tau_{sy}}{\rho_0 \phi} + \frac{\partial}{\partial x}\left(E_h \frac{\partial \langle VH \rangle_G}{\partial x}\right) + \frac{\partial}{\partial y}\left(E_h \frac{\partial \langle VH \rangle_G}{\partial y}\right)\]

Where:

\(\langle \cdot \rangle\) denotes averaged variables
Subscripts \(W\) and \(G\) indicate wet-averaging and grid-averaging respectively
\(U, V\) are depth-averaged velocities in x and y directions
\(H\) is total water depth
\(\zeta\) is water surface elevation
\(\phi\) is the wet area fraction (0 ≤ φ ≤ 1)
\(P_A\) is atmospheric pressure
\(\tau_{sx}, \tau_{sy}\) are surface wind stresses
\(\rho_0\) is reference water density
\(f\) is Coriolis parameter
\(E_h\) is horizontal eddy viscosity
\(g\) is gravitational acceleration

Closure Coefficients:

The subgrid method introduces five key closure coefficients:

\(C_\zeta\) - Surface gradient correction (typically set to 1.0)
\(C_{UU}, C_{VU}\) - Advection correction terms (typically 1.06-1.10)
\(\phi\) - Wet area fraction (0 ≤ φ ≤ 1)
\(C_{M,f}\) - Bottom friction coefficient correction

Closure Approximation Levels:

Level 0 Closure (Basic Implementation):

All advection corrections set to unity: \(C_{UU} = C_{VU} = C_{UV} = C_{VV} = 1\)
Bottom friction: \(C_{M,f} = gn^2/H^{1/3}\) (conventional Manning’s formula)
Surface gradient: \(C_\zeta = 1\)

Level 1 Closure (Enhanced Implementation):

Bottom friction correction:

\[C_{M,f} = \langle H \rangle_W R_v^2 \quad \text{where} \quad R_v = \frac{\langle H \rangle_W}{\langle H^{3/2}C_f^{-1/2} \rangle_W}\]

Advection corrections:

\[C_{UU} = \frac{\langle U^2H \rangle_G}{\langle U \rangle \langle UH \rangle_G}\]

Level 1 corrections provide significant accuracy improvements, particularly in shallow regions with complex topography and varying friction.

The Phi (φ) Factor

Physical Interpretation

The φ factor represents the wet area fraction of a computational element or vertex at a given water surface elevation. This concept encapsulates several important physical processes:

Partial Wetting: As water levels rise, different portions of a mesh element become submerged at different rates
Complex Topography: Sub-mesh topographic variability controls the wetting pattern
Flow Obstruction: Dry areas within an element impede flow, effectively reducing the wet cross-sectional area
Discrete Values: φ ranges from 0 (fully dry) to 1 (fully wet) in discrete increments

Vertex-Averaged Approach:

Current implementations use vertex-averaged areas rather than element-averaged areas. This approach:

Aligns more closely with conventional ADCIRC’s numerical scheme
Reduces lookup table size by approximately factor of 4
Improves computational performance and memory usage

Mathematical Definition

For a computational element with area \(A_{total}\), the φ factor is:

\[\phi(\zeta) = \frac{A_{wet}(\zeta)}{A_{total}}\]

Where \(A_{wet}(\zeta)\) is the area that is wet (submerged) when the water surface elevation is \(\zeta\).

Key properties of φ:

Bounds: \(0 \leq \phi \leq 1\)
Monotonic: \(\phi\) increases monotonically with increasing \(\zeta\)
Boundary Conditions:
- \(\phi \to 0\) as \(\zeta \to -\infty\) (completely dry)
- \(\phi \to 1\) as \(\zeta \to +\infty\) (completely wet)

Computational Algorithm

The φ factor is computed using high-resolution digital elevation models (DEMs):

Step 1: Element Elevation Sampling

For each computational mesh element, extract all DEM pixels that fall within the element boundaries:

\[E = \{z_i : (x_i, y_i) \in \text{element}, i = 1, 2, ..., n\}\]

Step 2: Water Level Discretization

Create a series of water surface elevations spanning the range of interest:

\[\zeta_j = \zeta_{min} + j \cdot \Delta\zeta, \quad j = 0, 1, ..., N_{levels}\]

Step 3: Wet Fraction Calculation

For each water level \(\zeta_j\), compute the fraction of pixels that are submerged:

\[\phi_j = \frac{1}{n} \sum_{i=1}^{n} H(\zeta_j - z_i)\]

Where \(H(\cdot)\) is the Heaviside step function:

\[\begin{split}H(x) = \begin{cases} 1 & \text{if } x > 0 \\ 0 & \text{if } x \leq 0 \end{cases}\end{split}\]

Bottom Friction Modifications

Manning’s Roughness Integration

Subgrid corrections also account for sub-mesh variability in bottom friction. The effective Manning’s coefficient is computed as an area-weighted average:

\[n_{eff} = \frac{\sum_{i=1}^{n} n_i A_i}{\sum_{i=1}^{n} A_i}\]

Where \(n_i\) is the Manning’s coefficient for land cover class \(i\), and \(A_i\) is the area covered by that class within the element.

Friction Coefficient Calculation

The friction coefficient \(C_f\) is related to Manning’s coefficient:

\[C_f = \frac{g n_{eff}^2}{H^{1/3}}\]

In the subgrid formulation, this becomes:

\[C_f = \frac{g n_{eff}^2}{(\phi H)^{1/3}}\]

This modification accounts for the reduced effective depth due to partial wetting.

Water Level Distribution Methods

Linear Distribution

The simplest approach distributes water levels evenly across a specified range:

\[\zeta_j = \zeta_{min} + \frac{j}{N-1}(\zeta_{max} - \zeta_{min})\]

Advantages: - Simple implementation - Uniform coverage of water level range

Disadvantages: - May not capture critical transitions - Inefficient for topographically diverse areas

Histogram-Based Distribution

A more sophisticated approach bases water level distribution on the elevation histogram within each element:

Step 1: Elevation Histogram

Compute the histogram of elevations within each element:

\[H(z_k) = \text{number of DEM pixels with elevation } z_k\]

Step 2: Cumulative Distribution

Calculate the cumulative distribution function:

\[F(z) = \frac{\sum_{z_i \leq z} H(z_i)}{\sum_{i} H(z_i)}\]

Step 3: Quantile-Based Levels

Distribute water levels based on elevation quantiles:

\[\zeta_j = F^{-1}\left(\frac{j}{N-1}\right)\]

Advantages: - Adaptive to local topography - Focuses resolution on critical elevation ranges - More efficient representation

Lookup Table Optimization

Table Size Reduction:

Modern subgrid implementations optimize lookup table sizes for computational efficiency:

Original Table Sizes:

\[N_\zeta \times N_{VE} \times N_E \quad \text{(element-averaged)}\]

\[N_\zeta \times N_V \quad \text{(vertex-averaged)}\]

Optimized Table Size:

\[N_\phi \times N_V\]

Where:

\(N_\zeta\) = number of water surface elevations (e.g., 401 for -20m to +20m at 0.1m increments)
\(N_{VE}\) = number of sub-elements per vertex
\(N_E\) = number of elements in mesh
\(N_V\) = number of vertices in mesh
\(N_\phi\) = number of φ values (typically 11)

Key Improvements:

Elimination of element-averaged tables - reduces memory requirements
Discrete φ values - \(N_\phi = 11\) instead of \(N_\zeta = 401\)
Memory reduction factor - approximately 160× for high-resolution meshes

Implementation Details

Wetting and Drying Criteria

Updated Approach:

Modern implementations use total water depth criteria instead of minimum wet area fraction:

\[\langle H \rangle_G > \langle H \rangle_{G,min}\]

Where \(\langle H \rangle_{G,min} = 0.1\) m (typical threshold).

Benefits:

Improved numerical stability
Avoids extremely small water depths that cause computational instabilities
Uses robust conventional ADCIRC wetting/drying scheme
For triangular elements: all 3 vertices must meet minimum criteria

Numerical Algorithms

Raster Processing Pipeline:

Mesh Element Identification: For each mesh node, identify the surrounding element(s)
DEM Intersection: Extract DEM pixels within element boundaries using spatial indexing
Land Cover Assignment: Map land cover classes to Manning’s coefficients
Statistical Computation: Calculate elevation statistics and φ relationships
Lookup Table Generation: Create optimized φ-based lookup tables

Memory Management:

The algorithm processes large raster datasets efficiently using:

Windowed Reading: Process raster data in memory-efficient chunks
Optimized Lookup Tables: Reduced memory footprint through φ-based indexing
Spatial Indexing: Use R-tree or similar structures for fast spatial queries
Caching: Cache frequently accessed data to minimize I/O

Parallel Processing:

Subgrid calculations are inherently parallel:

Vertex Independence: Each mesh vertex can be processed independently
Thread Safety: Algorithms designed for multi-threaded execution
Load Balancing: Distribute computational load based on vertex complexity

Quality Control and Validation

Physical Consistency Checks

Monotonicity Verification:

Ensure that φ increases monotonically with water level:

\[\phi(\zeta_{j+1}) \geq \phi(\zeta_j) \quad \forall j\]

Boundary Condition Verification:

Check that φ approaches appropriate limits:

\[\lim_{\zeta \to \zeta_{min}} \phi(\zeta) = 0\]

\[\lim_{\zeta \to \zeta_{max}} \phi(\zeta) = 1\]

Conservation Properties

Volume Conservation:

The subgrid method should preserve volume conservation in the discrete system:

\[\frac{d}{dt} \int_{\Omega} \zeta \, dA + \oint_{\partial\Omega} \phi uH \cdot \mathbf{n} \, dl = 0\]

Mass Conservation:

Verify that the modified equations maintain mass conservation properties of the original system.

Error Analysis

Sources of Error

Discretization Error: Finite resolution of DEM and computational mesh
Averaging Error: Spatial averaging within computational elements
Interpolation Error: Temporal interpolation of φ values during simulation
Data Quality Error: Uncertainties in input DEM and land cover data

Error Quantification

Convergence Analysis:

Study convergence behavior as: - DEM resolution increases - Number of φ levels increases - Computational mesh resolution changes

Validation Against High-Resolution Results:

Compare subgrid results with high-resolution simulations:

\[\text{Error} = \frac{|\zeta_{subgrid} - \zeta_{high-res}|}{\zeta_{high-res}}\]

Sensitivity Analysis:

Evaluate sensitivity to: - Manning’s coefficient variations - DEM uncertainty - φ level distribution choices

Performance Characteristics

Computational Complexity

Preprocessing Cost:

Time Complexity: O(N_elements × N_pixels × N_levels)
Space Complexity: O(N_elements × N_levels)

Where: - N_elements = number of mesh elements - N_pixels = average DEM pixels per element - N_levels = number of φ levels

Runtime Cost in ADCIRC:

Additional Memory: ~20-30% increase for subgrid arrays
Computational Overhead: ~10-20% increase in runtime
I/O Overhead: Initial reading of subgrid file

Scaling Properties

Mesh Resolution Scaling:

Subgrid benefit increases with decreasing mesh resolution:

\[\text{Benefit} \propto \left(\frac{\text{DEM resolution}}{\text{Mesh resolution}}\right)^2\]

Domain Size Scaling:

Processing time scales approximately linearly with domain size for fixed resolution ratios.

Applications and Limitations

Optimal Application Domains

Subgrid corrections are most effective for:

Coastal Storm Surge: Complex nearshore bathymetry and topography
Real-time Forecasting: Requiring rapid computation with maximum accuracy
Urban Flooding: Built environment with significant sub-mesh features
Wetland Modeling: Areas with complex channel and vegetation patterns
Large-Scale Simulations: Where uniform high resolution is computationally prohibitive
Hurricane Storm Surge: Combining operational time constraints with accuracy requirements

Validated Performance Ranges:

Coarsening factors: 3-20× (up to 50× in some applications)
Minimum nearshore resolution: 60-400m (depending on coastal complexity)
Domain scales: From regional bays (O(10 km)) to ocean basins (O(1000+ km))
Computational speedup: 10-50× compared to equivalent-accuracy fine meshes

Resolution Guidelines and Limitations

Critical Resolution Thresholds:

Barrier Island Aliasing: Primary limitation occurs when mesh resolution becomes coarser than critical flow-blocking features
Recommended minimum: Resolve primary barrier islands and major flow channels in underlying DEM
Practical limits: Beyond 400-1000m nearshore resolution, important coastal features become aliased
Connectivity issues: Over-connectivity in protected areas when barriers are under-resolved

Physical Process Limitations:

Unidirectional Approximation: Assumes φ depends only on local water level
Flow Direction Independence: Does not account for flow direction effects on wetting patterns
Static Topography: Does not account for morphodynamic changes
Simplified Wetting/Drying: May not capture complex hysteresis effects
Sub-element Assumptions: Requires all three sub-elements to be wet for proper channel flow representation

Future Research Directions

Potential improvements include:

Dynamic φ Factors: Account for flow-dependent wetting patterns
Momentum Subgrid Terms: Develop correction terms for advection and viscosity
Morphodynamic Coupling: Integrate with sediment transport and bed evolution
Machine Learning Enhancement: Use ML techniques for improved parameterizations

This scientific foundation provides the theoretical understanding necessary for effective application of subgrid methods in ADCIRC modeling.