======================= Meteorological Forcing ======================= Cocoa applies atmospheric forcing through two mechanisms: wind stress at the water surface and atmospheric pressure gradients. Both enter the governing equations as body forces and are ramped up during spinup to suppress initial transients. .. contents:: On This Page :local: :depth: 2 Wind Stress ----------- Surface Wind Drag ^^^^^^^^^^^^^^^^^ Wind stress is computed from 10-meter wind velocity using the Garratt (1977) bulk drag law [Garratt1977]_. The drag coefficient is a linear function of wind speed, capped at a configurable maximum: .. math:: :label: garratt-cd C_d = \min\bigl((0.75 + 0.067 \, U_{10}) \times 10^{-3},\; C_{d,\max}\bigr) where :math:`U_{10} = \sqrt{u_{10}^2 + v_{10}^2}` is the 10-meter wind speed magnitude and :math:`C_{d,\max}` is the drag coefficient cap (default 0.0025). The wind stress vector is then: .. math:: :label: wind-stress \frac{\boldsymbol{\tau}_s}{\rho_0} = \frac{\rho_a}{\rho_0} \, C_d \, U_{10} \, \begin{pmatrix} u_{10} \\ v_{10} \end{pmatrix} with :math:`\rho_a = 1.293` kg/m\ :sup:`3` (air density at 0 |deg|\ C, ADCIRC convention) and :math:`\rho_0 = 1000` kg/m\ :sup:`3` (reference seawater density). The result is a kinematic stress with units m\ :sup:`2`/s\ :sup:`2`. .. note:: The cap :math:`C_{d,\max}` prevents unrealistically large stress at extreme wind speeds, where observational evidence suggests the drag coefficient levels off or decreases. The default value of 0.0025 follows ADCIRC convention. .. figure:: ../_static/images/diagrams/garratt_drag.svg :alt: Garratt drag coefficient as a function of wind speed :width: 500px :align: center Garratt (1977) drag coefficient showing linear growth capped at :math:`C_{d,\max} = 0.0025`. The transition occurs near 26 m/s. .. _wind-limiter: Wind Stress Limiter ^^^^^^^^^^^^^^^^^^^ To prevent abrupt wind stress changes at the wet/dry boundary, Cocoa applies a depth-dependent limiter that smoothly ramps wind stress to zero as water depth approaches the drying threshold :math:`H_0`. This is always active and improves stability during wetting/drying cycles with strong winds. The limiter computes a scaling factor :math:`f_w \in [0, 1]`: .. math:: :label: wind-limiter f_w(H) = \begin{cases} 1 & H > 4 H_0 \\[4pt] \dfrac{1}{2}\tanh\!\left(\dfrac{8\,(H - 2.5\,H_0)}{4\,H_0}\right) + \dfrac{1}{2} & H_0 < H \le 4 H_0 \end{cases} Since both GWCE and momentum kernels iterate only over wet nodes/elements (where :math:`H > H_0` by construction), the :math:`H \le H_0` case is never reached. The early return for :math:`H > 4 H_0` avoids the ``tanh`` evaluation for the vast majority of the domain, where all threads in a GPU warp take the same branch with no divergence. This formulation follows ADCIRC's ``windLimiter`` subroutine (DMW 2022). .. figure:: ../_static/images/diagrams/wind_limiter.svg :alt: Wind stress limiter as a function of water depth :width: 500px :align: center Wind stress limiter :math:`f_w` smoothly transitions from 0 near the drying threshold :math:`H_0` to 1 at :math:`4H_0`, suppressing wind forcing in nearly-dry cells. Application in the Momentum Equations ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ In the momentum equations (:eq:`x-momentum`, :eq:`y-momentum`), wind stress appears as a surface forcing divided by the water column depth: .. math:: :label: wind-momentum F_x^{wind} = \frac{f_w \, \tau_{sx}/\rho_0}{H}, \quad F_y^{wind} = \frac{f_w \, \tau_{sy}/\rho_0}{H} where :math:`\tau_s/\rho_0` is the kinematic wind stress (:eq:`wind-stress`), :math:`H` is the total water depth, and :math:`f_w` is the wind stress limiter (:eq:`wind-limiter`). **Discretized form (nodal):** Wind stress is applied per node after the element-based RHS assembly, with Crank-Nicolson centering between the two time levels: .. math:: RHS_{x,i} \mathrel{+}= \frac{\Delta t}{2} \, f_{w,i} \left( \frac{R_{met}^{n} \, \tau_{sx,i}^{n}/\rho_0}{H_i^{n}} + \frac{R_{met}^{n+1} \, \tau_{sx,i}^{n+1}/\rho_0}{H_i^{n+1}} \right) and identically for :math:`RHS_{y,i}`, where :math:`R_{met}` is the meteorological ramp factor (see :ref:`meteorological-ramp`; the stress fields are stored pre-ramped) and :math:`f_{w,i}` is evaluated at the :math:`n+1` depth. Application in the GWCE ^^^^^^^^^^^^^^^^^^^^^^^^ In the GWCE (:eq:`gwce`), wind stress enters the :math:`\mathbf{J}` vector as a depth-integrated surface forcing, scaled per node by the wind limiter before element averaging: .. math:: :label: wind-gwce J_x^{wind} = \frac{R_{met}}{3} \sum_{i=1}^{3} f_{w,i} \, \frac{\tau_{sx,i}}{\rho_0}, \quad J_y^{wind} = \frac{R_{met}}{3} \sum_{i=1}^{3} f_{w,i} \, \frac{\tau_{sy,i}}{\rho_0} where :math:`f_{w,i}` is the wind limiter evaluated at node :math:`i` using that node's total depth :math:`H_i = h_i + \zeta_i`. Unlike the momentum equation, wind stress in the GWCE is not divided by depth. The GWCE operates on depth-integrated quantities, and wind stress is a surface force per unit area. Atmospheric Pressure Gradient ----------------------------- The atmospheric pressure gradient force accelerates water from high to low atmospheric pressure. This is the primary mechanism by which storm surge is generated, with the "inverted barometer" effect producing roughly 1 cm of surge per 1 hPa of pressure deficit. Application in the GWCE ^^^^^^^^^^^^^^^^^^^^^^^^ The atmospheric pressure gradient enters the GWCE through the :math:`\mathbf{J}` vector: .. math:: :label: atm-pressure-gwce J_x^{atm} = -\frac{\alpha_L \, \bar{H}}{\rho_0 \cdot 2A_e} \sum_{i=1}^{3} p_{s,i} \cdot FDX_i, \quad J_y^{atm} = -\frac{\alpha_L \, \bar{H}}{\rho_0 \cdot 2A_e} \sum_{i=1}^{3} p_{s,i} \cdot FDY_i where: - :math:`p_{s,i}` is the atmospheric pressure at node :math:`i` [Pa] - :math:`\alpha_L` is the wet/dry slope limiter - :math:`\bar{H}` is the element-averaged total depth - :math:`FDX_i, FDY_i` are the shape function derivative products The slope limiter :math:`\alpha_L` is applied consistently with the tide potential gradient (see :doc:`../theory/gwce`), reducing pressure gradients in shallow water to suppress spurious oscillations. **Pressure ramping:** During spinup, the atmospheric pressure is blended from the background reference pressure toward the actual forcing data: .. math:: :label: pressure-ramp p_{s}^{ramped} = p_{bg} + R_{met} \cdot (p_s - p_{bg}) where :math:`p_{bg} = 101300` Pa is the background atmospheric pressure and :math:`R_{met}` is the meteorological ramp factor, so the pressure forcing starts from hydrostatic equilibrium and grows smoothly. Application in the Momentum Equations ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ In the momentum equations, the atmospheric pressure gradient appears alongside the barotropic pressure gradient: .. math:: F_x^{atm} = -\frac{1}{\rho_0}\frac{\partial p_s}{\partial x}, \quad F_y^{atm} = -\frac{1}{\rho_0}\frac{\partial p_s}{\partial y} This term is already present in the momentum equation formulation (see :eq:`x-momentum`). For storm surge simulations, the atmospheric pressure gradient typically dominates over wind stress in deep water, while wind stress becomes more important in shallow coastal areas. .. _meteorological-ramp: Meteorological Ramp ------------------- To prevent impulse responses at initialization, meteorological forcing is multiplied by a spin-up ramp :math:`R_{met}(t)` that rises from 0 to 1 (:doc:`../user_guide/ramps`). It starts at the earlier of the cold start and the first meteorological snapshot. By default, the meteorological ramp uses the same configuration as the tidal ramp (``forcing.ramp``). An optional ``forcing.meteorological.ramp`` section allows a different ramp duration for meteorological forcing, which is useful when starting from a pre-spun-up tidal state that needs only a short met ramp. Internally, ``ForcingManager`` maintains separate ramp objects for tidal and meteorological forcing. .. _meteo-temporal-interp: Temporal Interpolation ---------------------- Meteorological forcing data is provided at discrete time snapshots. Between snapshots, Cocoa linearly interpolates the forcing fields: .. math:: \phi(t) = (1 - \alpha) \, \phi_{k} + \alpha \, \phi_{k+1} where :math:`\alpha = (t - t_k) / (t_{k+1} - t_k)` is the interpolation weight and :math:`\phi` represents any forcing field (pressure, wind_u, wind_v). The interpolation is performed at each mesh node after spatial interpolation from the meteorological grid. Cocoa maintains a ring buffer of pre-read snapshots to amortize I/O costs across timesteps. Spatial Interpolation --------------------- Meteorological data is typically provided on a regular rectangular grid at coarser resolution than the computational mesh. Cocoa uses bilinear interpolation to transfer forcing data from the meteorological grid to mesh nodes. For a mesh node located at :math:`(x, y)` within a grid cell with corners at :math:`(x_0, y_0)` and :math:`(x_1, y_1)`: .. math:: :label: bilinear \phi(x, y) = w_{00} \, \phi_{00} + w_{10} \, \phi_{10} + w_{01} \, \phi_{01} + w_{11} \, \phi_{11} where the bilinear weights :math:`w_{ij}` are products of the normalized distances: .. math:: s = \frac{x - x_0}{x_1 - x_0}, \quad t = \frac{y - y_0}{y_1 - y_0} .. math:: w_{00} = (1-s)(1-t), \quad w_{10} = s(1-t), \quad w_{01} = (1-s)t, \quad w_{11} = st The interpolation weights are precomputed once when the meteorological data is initialized and cached for reuse at each timestep. For regular (rectangular) grids, this is an index computation. For irregular grids (OWI NetCDF with moving vortex tracking), the weights are recomputed when the grid moves. Multi-Domain Composition ^^^^^^^^^^^^^^^^^^^^^^^^ Cocoa can compose several meteorological sources into one forcing field. Each source is a ``forcing.meteorological.domains`` entry of any format, or the flat single-source form. During interpolation, Cocoa applies a priority rule: the innermost valid domain takes precedence over outer domains, where "innermost" ranks first by configured source order (later entries win at overlap) and then, within a source that is itself multi-domain, by that source's own domain index. Mesh nodes that fall outside all domains receive background values (101300 Pa pressure, zero wind). OWI's internal nesting of several grids within one file (a coarse outer grid plus fine inner grids) is the same last-wins rule applied one level deeper, so composition does not depend on format. This approach allows, for example, a coarse basin-scale field to cover the full model domain while a fine hurricane-following nest (from the same source or a different one, e.g. a GRIB regional model nested inside a coarser OWI field) provides high-resolution forcing in the region of interest. .. _meteo-cadence-interp: Cadence Interpolation onto the Union Time Axis ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Composed sources need not share a common time step. Cocoa builds a single union time axis as the sorted, deduplicated set of every domain's own sample times, and reads each domain against it. Where a union time coincides with one of a domain's own samples, that domain's data is used verbatim, without interpolation. Where a union time falls strictly between two of a domain's own samples (a coarser-cadence source relative to the axis), the domain's fields (and, for a moving grid, its coordinates) are linearly interpolated between its own bracketing snapshots: .. math:: \phi_d(t) = (1 - \beta) \, \phi_{d,k} + \beta \, \phi_{d,k+1}, \qquad \beta = \frac{t - t_{d,k}}{t_{d,k+1} - t_{d,k}} where :math:`\phi_d` is domain :math:`d`'s field, :math:`t_{d,k}` and :math:`t_{d,k+1}` are its own bracketing sample times, and :math:`t` is the union-axis time being synthesized. This source-side (grid-space) interpolation is distinct from the node-side temporal blend of :ref:`meteo-temporal-interp`: it lets a coarse-cadence source (e.g. a 3-hourly field) coexist with a fine-cadence one (e.g. hourly) on the same union axis without forcing every composed source to share a time step. .. _meteo-activation-window: Activation Window Fade ^^^^^^^^^^^^^^^^^^^^^^^ A source active only over part of the run (a moving nest, or any source with a configured activation window) can be given an activation fade, a ramp block whose weight :math:`w(t)` rises from 0 at the window's open time and falls back to 0 at its close time on the union axis (:ref:`ramp-activation-fades`). Where the source covers a node, the composed field cross-fades toward it, :math:`\text{node} = (1 - w) \cdot \text{underlying} + w \cdot \text{source}`, in place of the hard override of the priority rule above. A zero ``duration`` (the default) disables the fade, recovering hard window edges. The same weight gates the Eulerian per-domain combine, the Lagrangian moving-grid overlay, and the GAHM vortex overlay's per-storm temporal factor. The moving-grid and vortex overlays evaluate it every time step; the Eulerian combine evaluates it at each union snapshot time while composing that slot, so between snapshots the weight is interpolated linearly along with the fields. Wind Reduction -------------- Cocoa optionally reduces wind velocity over land before applying the drag law. Two reductions apply in sequence. Both are applied per meteorological source, folded into that source's wind at composition time (before it is blended into the combined field) and gated by the source's ``wind_reduction`` flag, so co-located sources can be treated differently: for example, a parametric vortex reduced while an already-surface GFS field is left unreduced. A node covered by a single source gets the same stress a post-blend reduction would give; only feather rings between a reduced and an unreduced source differ. Directional Roughness Reduction ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ The directional roughness reduction modifies wind components based on the ratio of upwind land surface roughness to marine roughness. This follows ADCIRC's ``ApplyDirectionalWindReduction`` formulation. Marine roughness is derived from the Garratt drag coefficient via Charnock's relation: .. math:: z_{0m} = \frac{0.018}{g} \, C_d \, U_{10}^2 where :math:`C_d` is the drag coefficient computed from the pre-reduction wind speed. Land roughness :math:`z_{0l}` is obtained by linear interpolation from the two nearest directional bins. The mesh stores 12 roughness lengths per node (one per 30-degree bin, starting from East and proceeding counter-clockwise). For overland flooding, when the total water depth :math:`H = h + \zeta` exceeds :math:`2 h_0` (twice the minimum depth threshold), the effective land roughness is reduced: .. math:: z_{0l}^{eff} = z_{0l} - \frac{H}{30} This accounts for the reduced influence of land surface features when they are submerged by floodwater. When :math:`z_{0l} > z_{0m}`, the wind reduction factor is: .. math:: f_r = \min\left(1,\; \left(\frac{z_{0l}}{z_{0m}}\right)^{0.0706} \frac{\ln(10 / z_{0l})}{\ln(10 / z_{0m})} \right) When :math:`z_{0l} \le z_{0m}` (land is smoother than open water), no reduction is applied (:math:`f_r = 1`). The reduced wind components are: .. math:: u_{10}' = f_r \, u_{10}, \quad v_{10}' = f_r \, v_{10} Canopy Coefficient ^^^^^^^^^^^^^^^^^^ The canopy coefficient :math:`C_c \in \{0, 1\}` is a binary multiplicative factor representing wind attenuation under forest canopy: .. math:: u_{10}'' = C_c \, u_{10}', \quad v_{10}'' = C_c \, v_{10}' A value of :math:`C_c = 1` (default) means no canopy effect; :math:`C_c = 0` suppresses the wind (and hence the stress) at that node. The coefficient is folded into the source wind here rather than scaling the post-drag stress. For a binary mask this is exact (:math:`0` and :math:`1` are unchanged by the drag law's squaring, so :math:`C_c = 0` gives zero stress and :math:`C_c = 1` is a no-op either way), which is why fractional values are rejected at startup. After both reductions the corrected source winds are blended, and the combined field is passed to the Garratt drag law (:eq:`garratt-cd`, :eq:`wind-stress`) once, after the meteorological ramp. Physical Constants ------------------ .. list-table:: :header-rows: 1 :widths: 35 20 45 * - Constant - Value - Description * - :math:`\rho_0` (seawater density) - 1000 kg/m\ :sup:`3` - Reference density (ADCIRC convention) * - :math:`\rho_a` (air density) - 1.293 kg/m\ :sup:`3` - Dry air density at 0 |deg|\ C (ADCIRC convention) * - :math:`p_{bg}` (background pressure) - 101300 Pa - Ambient pressure away from any storm (1013 mb, the ADCIRC/ATCF convention). The parametric vortex resolves its own far-field :math:`p_n` per track snapshot from the ATCF ``RADP`` column and falls back to this value only when that column is blank; see :ref:`parametric-vortex`. * - :math:`g` (gravity) - 9.81 m/s\ :sup:`2` - Gravitational acceleration * - :math:`C_{d,\max}` (drag cap) - 0.0025 - Default maximum drag coefficient Implementation Files -------------------- - ``forcing/meteorological/DragLaw.hpp``: ``wind_drag()`` Garratt drag coefficient and ``wind_ice_drag()`` ice attenuation (device-callable) - ``forcing/meteorological/WindStress.hpp``: ``wind_stress()`` free function template for kinematic wind stress (device-callable) - ``forcing/meteorological/WindReduction.hpp``: Directional wind reduction kernel (device-callable) - ``physics/WindLimiter.hpp``: Depth-based wind stress limiter for drying cells (device-callable) - ``numeric/cg/terms/WindStressTerm.hpp``: ``GwceWindStress`` wind stress GWCE contribution - ``numeric/cg/terms/AtmPressureGradient.hpp``: ``GwceAtmPressureGradient`` atmospheric pressure gradient GWCE contribution - ``numeric/cg/momentum/MomentumRhsKernels.cpp``: Wind stress in momentum RHS - ``forcing/meteorological/MeteoForcingProvider.hpp``: Ring buffer, temporal interpolation, and overlay orchestration - ``forcing/meteorological/MeteoSnapshotSource.hpp``: Rank-0 reader thread; union time axis, per-domain windows, cadence interpolation - ``cocoa_meteo/CompositeMeteoReader.hpp``: Composes multiple configured sources into one flat, priority-ordered domain space - ``cocoa_meteo/interpolation/MultiDomainInterpolator.hpp``: Host spatial combine, priority rule, winner mask, activation weights - ``forcing/WindowedRamp.hpp``: Activation/deactivation window fade (shares the ramp function with ``Ramp``) - ``forcing/ForcingManager.hpp``: Ramp computation, wind reduction, and forcing orchestration