============ Wave Forcing ============ Breaking waves push water shoreward. The wave field's momentum flux, the radiation stress [LonguetHiggins1964]_, varies in space, and its divergence is a force on the water column that sets up the water level where waves break and draws it down just seaward of the breaking zone. In a hurricane this wave setup is a large share of the surge in the surf zone [Dietrich2011]_. Cocoa takes that force from one of two sources, chosen by ``forcing.wave.source``: - **A file** (`External Wave Forcing`_) carrying the force itself, from any wave model (an ADCIRC+SWAN ``rads.64`` converts directly), or the radiation stress, whose divergence Cocoa takes. - **Cocoa's own spectral wave model** (`Internal Wave Model`_), run on the circulation mesh and coupled to it every few minutes. Either way the force enters the equations exactly as the wind stress does. Every key is listed in the configuration reference (:ref:`wave-forcing-config`). .. contents:: On This Page :local: :depth: 2 Theory ------ For a spectrum of action density :math:`N(\sigma, \theta)` with group and phase speeds :math:`c_g` and :math:`c`, and :math:`n = c_g / c`, the radiation stress divided by the water density (units m\ :sup:`3` s\ :sup:`-2`) is .. math:: :label: wave-radiation-stress \frac{S_{xx}}{\rho_0} &= g \iint \left(n \cos^2\theta + n - \tfrac12\right) \sigma N \, d\sigma\, d\theta,\\ \frac{S_{xy}}{\rho_0} &= g \iint n \sin\theta\cos\theta \, \sigma N \, d\sigma\, d\theta,\\ \frac{S_{yy}}{\rho_0} &= g \iint \left(n \sin^2\theta + n - \tfrac12\right) \sigma N \, d\sigma\, d\theta. The force per unit density, in m\ :sup:`2` s\ :sup:`-2`, the units of the kinematic wind stress :math:`\boldsymbol{\tau}/\rho_0`, is .. math:: :label: wave-force \mathbf{F} = -\nabla \cdot \frac{\mathbf{S}}{\rho_0} = -\left(\frac{\partial S_{xx}}{\partial x} + \frac{\partial S_{xy}}{\partial y},\; \frac{\partial S_{xy}}{\partial x} + \frac{\partial S_{yy}}{\partial y}\right) / \rho_0. The two are one surface stress: they enter the GWCE at time level :math:`n` (averaged over an element under the wind limiter) and the momentum equations as the Crank-Nicolson pair .. math:: :label: wave-force-momentum \frac{\Delta t}{2}\left(\frac{\mathbf{F}^n}{H^n} + \frac{\mathbf{F}^{n+1}}{H^{n+1}}\right). ADCIRC does the same by adding the force into its wind stress array. Applying the Force ------------------ These conventions hold for both sources. **One surface stress.** The wave force is summed with the wind stress where the equations read it, so it inherits the wind stress's wet/dry treatment: the wind limiter near drying nodes and the depth floor in the momentum division. A run without meteorological forcing still applies the wave force. **Divergence of a stress.** Where Cocoa differentiates a radiation stress (a stress file, or the internal model's), the element-mean gradient .. math:: :label: wave-stress-gradient \frac{\partial S}{\partial x}\bigg|_e = \frac{1}{2A}\sum_i S_i\,b_i, \qquad \frac{\partial S}{\partial y}\bigg|_e = \frac{1}{2A}\sum_i S_i\,a_i uses the same Green coefficients :math:`b_i`, :math:`a_i` (:eq:`quad-green`) and momentum scale factors as every other element term, which is exact for the linear stress on a triangle and the bilinear stress on a quadrilateral. A node then takes the area-weighted mean of its elements' values. The force is zero at a node on the mesh boundary and at a node touching any element that is not wet (not active, or with a dry vertex), as in SWAN: a wave model sets :math:`S = 0` at dry nodes, so a gradient there is an artifact of the waterline, not wave forcing. **Ramp.** ``forcing.wave.ramp`` is the wave force's own ramp (ADCIRC's ``RampWRad``), the same hyperbolic-tangent form as the other ramps. It is on by default over one day and does not inherit ``forcing.ramp``. **Cap.** ``forcing.wave.cap`` bounds the magnitude of the ramped force the equations see, keeping its direction (ADCIRC's ``WaveStressGrad_Cap``). The output keeps the uncapped value. **Precision.** The force is stored in single precision at two time levels, the precision in which wave models hand over their stresses; all arithmetic is double. **Checkpoints.** The force's current level is checkpointed (the next step stages it as time level :math:`n`), so a restart reproduces the uninterrupted run. A restart from a checkpoint written without wave forcing evaluates the force at the resume time. External Wave Forcing --------------------- .. code-block:: yaml forcing: wave: enabled: true source: file filename: "wave_force.nc" quantity: force # force (default) or stress ramp: enabled: false # a converted rads.64 is already ramped cap: 0.05 # optional, m2 s-2 A file needs ``source: file``, since the internal model is the default. ``quantity`` says what the file carries: the force, or the radiation stress, whose divergence Cocoa takes each step. File Format ^^^^^^^^^^^ A NetCDF file with dimensions ``time`` and ``node`` (the mesh's node count, in mesh order): .. list-table:: :header-rows: 1 :widths: 25 20 55 * - Variable - Units - Content * - ``time`` - CF, e.g. ``seconds since 2011-08-25 12:00:00`` - Strictly ascending snapshot times covering the run. * - ``force_x``, ``force_y`` - ``m2 s-2`` - Force per unit density, geographic east/north (``quantity: force``). * - ``sxx``, ``sxy``, ``syy`` - ``m3 s-2`` - Radiation stress per unit density, mesh frame (``quantity: stress``). Values must be finite; dry nodes carry zero. A ``mesh_id`` global attribute, when present, is checked against the run's mesh. The file is validated at startup, and every failure names the file and what to fix. **Time.** Snapshots are interpolated linearly in time, ADCIRC's ``NRS = 1`` reading, and the file must cover the whole run. A stress is interpolated before its divergence is taken, so every step uses its own wet mask. **Frames.** A force file is geographic east/north and is rotated into the solver frame on a rotated-pole mesh. A stress file is in the mesh frame and is refused on a rotated-pole mesh, where the tensor would need rotating before it is differentiated. Converting ADCIRC Output ^^^^^^^^^^^^^^^^^^^^^^^^ ``utils/adcirc_rads_to_cocoa.py`` converts an ADCIRC+SWAN ``rads.64`` (ASCII or NetCDF) into a force file: dry-node values become zero and the units attribute is corrected. .. code-block:: bash python utils/adcirc_rads_to_cocoa.py rads.64.nc wave_force.nc --prepend-zero python utils/adcirc_rads_to_cocoa.py rads.64 wave_force.nc \ --start "2011-08-25 12:00:00" --prepend-zero ``rads.64`` holds ADCIRC's ramped force, so run Cocoa with ``forcing.wave.ramp.enabled: false``. Its first record is one output interval after the cold start; ``--prepend-zero`` adds a zero record at the cold start, where a ramped ADCIRC run's force is zero, so the file covers a Cocoa run beginning at the same time. Write ``rads.64`` at the coupling interval (``NSPOOLGW``) for the interpolation to follow the forcing. Deviations from ADCIRC ^^^^^^^^^^^^^^^^^^^^^^ - **Waterline and boundary nodes.** Cocoa follows SWAN's ``SwanComputeForce`` and ADCIRC's ``MARCELSWAN`` build option: the force is zero at boundary nodes and at nodes touching a dry element. ADCIRC's default build averages each node's force over every element around it, dry ones included, which puts an artificial shoreward gradient into waterline nodes. Interior nodes agree with ADCIRC exactly. - **Nodal gradient.** Cocoa keeps ADCIRC's area-weighted mean of element gradients rather than SWAN's Green's-theorem gradient over the polygon joining the centroids of the surrounding cells. - **Dryness.** An element is dry when Cocoa's wet/dry scheme marks it inactive or any of its vertices dry, not when a depth falls below SWAN's ``DEPMIN``. - **MARCELSWAN partial sum.** At a node touching a dry element, ADCIRC's ``MARCELSWAN`` branch zeroes the output array instead of its accumulator (``couple2swan.F``), so the node keeps the undivided sum of the elements visited before the dry one. Cocoa implements the intended zero. - **Units attribute.** ADCIRC's NetCDF ``rads.64`` labels the force ``m-2 s-2``; the values are m\ :sup:`2` s\ :sup:`-2`, which is what the converter writes. Internal Wave Model ------------------- With ``source: model`` Cocoa runs its own spectral wave model on the same mesh. Its equations, numerics and source terms, and where each came from, are in :doc:`../theory/spectral_waves`; this section covers running it. .. code-block:: yaml forcing: wave: enabled: true source: model # the default Every key has a default, and the defaults are the SWAN+ADCIRC hurricane configuration, so this is a complete setup. The waves take their wind from ``forcing.meteorological``; without it they have no wind to grow from. The configuration reference shows every key at its default (:ref:`wave-forcing-config`). Coupling ^^^^^^^^ At every coupling time :math:`t_k` the model takes the circulation as it stands (water depth, depth-averaged current, the wet mask), advances the spectrum by one coupling interval (``model.coupling_interval``, 10 minutes by default), and forms the radiation stress :eq:`wave-radiation-stress` as SWAN's sum over the bins, .. math:: :label: wave-stress-quadrature \iint (\cdot)\,\sigma N\,d\sigma\,d\theta \;\to\; \sum_f \sum_d (\cdot)\,\sigma_f N_{fd}\,\Delta\sigma_f\,\Delta\theta. Its divergence (`Applying the Force`_) is the force :math:`\mathbf F_k`. The next force is not known until :math:`t_{k+1}`, so between couplings the circulation sees ADCIRC+SWAN's schedule, a linear ramp from :math:`\mathbf F_k` toward the extrapolated :math:`2\mathbf F_k - \mathbf F_{k-1}`: .. math:: :label: wave-force-schedule \mathbf F(t) = \mathbf F_k + \frac{t - t_k}{t_{k+1} - t_k} \left(\mathbf F_k - \mathbf F_{k-1}\right), \qquad t_k \le t < t_{k+1} (:math:`\mathbf F_0` is held over the first interval). The couplings run every interval from the start, or, on a restart whose checkpoint carries no wave state, from the restart time. A checkpoint carries the spectrum, the force bracket, the depth of the last coupling and the relaxed sink rate, so a restart continues the schedule exactly. The spectrum lives in the solver frame, so the force needs no rotation on a rotated-pole mesh; its output copy is rotated back to geographic like every other vector. The progress table shows the largest Hs and each coupling's iterations (:doc:`../getting_started/quick_start`). A step that reaches ``stop.max_iterations`` is logged and the run continues; an Hs that is not finite, or that reaches ``stop.max_height`` (30 m by default), stops the run as a diverged circulation does. Choosing the Spectral Grid ^^^^^^^^^^^^^^^^^^^^^^^^^^ ``directions`` sets the direction bins, :math:`360^\circ/n` wide. ``frequencies`` sets how many frequencies are spaced geometrically from ``lowest_frequency`` to ``highest_frequency``, both included. This is SWAN's ``CGRID ... CIRCLE mdc flow fhigh msc``, whose ``msc`` counts the intervals between frequencies, one fewer than ``frequencies``: the default grid is SWAN's ``CIRCLE 36 0.031384 1.420416 40`` (:ref:`spectral-waves-grid`). - **Frequency spacing.** Keep the ratio between neighbors near 1.1, .. math:: :label: wave-frequency-ratio r = \left(\frac{f_{max}}{f_{min}}\right)^{1/(N_f - 1)}, with :math:`N_f` the ``frequencies`` count. The quadruplet approximation interpolates across bins and is tuned for that spacing, as SWAN recommends. To widen the range, add frequencies rather than stretching the spacing. - **Lowest frequency.** Below the longest swell the domain can generate: 0.031 Hz is a 32 s period. - **Highest frequency.** At least an octave above the peak of the lightest wind sea that matters: a young sea under 4 m/s peaks near 0.6-0.7 Hz. The spectrum is extended above the grid as :math:`f^{-4}`, but only the bins on the grid grow, and a peak close to the top grows too slowly. - **Directions.** 36 (10 degrees) is the SWAN+ADCIRC standard. Operational unstructured hurricane runs use 24 to 36; fewer directions spread swell into rays (the garden-sprinkler effect) and resolve refraction more coarsely. Use an even count: an odd one makes the sweep's snapshot hold every node instead of the block-boundary rows (:ref:`spectral-waves-memory`), and the minimum is 3. Cost and memory scale with the number of bins, directions times frequencies. The model holds two single-precision spectra (the iterate and the step's start, 4 bytes per node per bin each) and a snapshot of the sweep blocks' boundary rows. The full CPRA Louisiana mesh (1.57M nodes, 36 x 41 bins) needs 9.3 GB per spectrum: the iterate stays on a 16 GB V100 and ``sweep.host_spectra: auto`` puts the starting spectrum in pinned host memory, costing about a quarter more time per iteration and leaving the results unchanged. The log reports where each spectrum was placed and how much device memory remained. Choosing the Source Terms ^^^^^^^^^^^^^^^^^^^^^^^^^ ``sources.physics`` picks the wind input and whitecapping (:ref:`spectral-waves-sources`); breaking, friction and the quadruplets are the same under both. - ``default``: SWAN's ``GEN3 KOMEN AGROW`` with ``WCAP KOMEN``, the configuration of the SWAN+ADCIRC hurricane studies. Its whitecapping coefficients are exposed (``sources.whitecapping``). - ``st6``: SWAN's ST6 package as Day and Dietrich recommend it for ADCIRC+SWAN storm runs (:ref:`spectral-waves-st6`). Its coefficients are fixed, and setting ``sources.whitecapping`` with it is refused. The wind the waves feel is relative to the current, :math:`\mathbf U_{10} - s\,\mathbf U_c` :eq:`wave-relative-wind`, with :math:`s` set by ``sources.wind.current_coupling``; the circulation's own wind stress uses the absolute wind (:ref:`spectral-waves-wind`). The default, 1, is SWAN's form; the atmosphere partly adjusts to a moving surface, so a smaller value may be more realistic. :math:`u_*` comes from the wave model's own drag law at that relative speed: ``sources.wind.drag_law``, chosen from the laws the circulation offers and capped by ``sources.wind.drag_limit``. ``swell`` reads the spread of the spectrum at the start of each step, where SWAN re-forms it every iteration. The defaults, Garratt capped at 0.0025 with SWAN's constant :math:`1.2875\times10^{-3}` at and below 7.5 m/s (``sources.wind.low_wind_drag``), are how ADCIRC+SWAN couples when the circulation keeps its own defaults. The two caps are independent here, whereas ADCIRC+SWAN hands SWAN the circulation's cap: a run that changes ``physics.wind_drag_limit`` and wants the coupled behavior must set ``sources.wind.drag_limit`` to the same value. To run as standalone SWAN does, set a ``drag_limit`` no law reaches and either ``wu`` for the default package, whose law already holds :math:`1.2875\times10^{-3}` at and below 7.5 m/s so ``low_wind_drag`` changes nothing, or ``hwang`` with ``low_wind_drag: false`` for ST6. The waves take the meteorological wind with its ``wind_scale`` divided back out, as ADCIRC hands SWAN an ``NWS = 12`` wind without fort.22's multiplier, times ``sources.wind.multiplier`` (ADCIRC's ``waveWindMultiplier``). That needs every gridded source to share one ``wind_scale``, and a vortex to be composed only with unscaled sources; the run refuses to start otherwise. Without meteorological forcing there is no wind input. ``sources.friction.roughness: manning`` derives the friction roughness from the mesh's Manning's :math:`n` at every coupling, as ADCIRC does for SWAN's ``READINP ADCFRIC``. The value ADCIRC passes is a roughness length :eq:`wave-manning-kn`, which SWAN uses as :math:`k_N` unchanged (:ref:`spectral-waves-shared`). A roughness needs ``sources.friction.enabled: true``; one given with friction off is refused. Limitations ^^^^^^^^^^^ - **No incoming boundary spectra.** Mesh boundaries let waves out and none in, so the model sees only waves generated inside the domain; swell arriving from outside it is absent. - **Triads and wave-induced currents** beyond the radiation-stress force are not modeled; the Stokes drift is omitted from continuity, as in SWAN+ADCIRC. - **Waves start from rest** at the first coupling (or at a restart whose checkpoint has no wave state), so allow a spin-up before the waves matter. Output ------ ``output.variables.wave_force: true`` writes ``wave_force_x`` and ``wave_force_y`` (m\ :sup:`2` s\ :sup:`-2`): the ramped force, uncapped and geographic, as ADCIRC writes ``rads.64``. Off by default, like the other forcing diagnostics. The internal model also writes its integral parameters at the last coupling: ``wave_hs``, ``wave_dir`` and ``wave_tm01`` by default, ``wave_tmm10`` and ``wave_tps`` when set to ``true`` (:doc:`output_files`). They follow SWAN's definitions with CF metadata (:ref:`spectral-waves-output`). The peak file carries ``wave_hs_max`` and, for each period and the direction written, its value at the coupling of that maximum (``wave_tm01_max``, ...), so the periods describe the sea state of the largest waves rather than maxima of their own (:ref:`output-peak-variables`). Two conventions matter when comparing with ADCIRC+SWAN: the periods are absolute, as SWAN's ``TM01`` and ``TMM10`` (the frequency a fixed point sees on a current), and ``wave_dir`` is where the waves come from, clockwise from north, where SWAN's default ``DIR`` is where they travel, counterclockwise from east. References ---------- .. [LonguetHiggins1964] Longuet-Higgins, M. S., and R. W. Stewart (1964). Radiation stresses in water waves; a physical discussion, with applications. *Deep Sea Research*, 11(4), 529-562. .. [Dietrich2011] Dietrich, J. C., et al. (2011). Modeling hurricane waves and storm surge using integrally-coupled, scalable computations. *Coastal Engineering*, 58(1), 45-65.