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 (Wave Forcing).

Theory 

For a spectrum of action density \(N(\sigma, \theta)\) with group and phase speeds \(c_g\) and \(c\), and \(n = c_g / c\), the radiation stress divided by the water density (units m³ s^-2) is

(1)\[\begin{split}\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.\end{split}\]

The force per unit density, in m² s^-2, the units of the kinematic wind stress \(\boldsymbol{\tau}/\rho_0\), is

(2)\[\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 \(n\) (averaged over an element under the wind limiter) and the momentum equations as the Crank-Nicolson pair

(3)\[\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

(4)\[\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 \(b_i\), \(a_i\) ((40)) 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 \(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 \(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 

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

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.

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² s^-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 Spectral Wave Model; this section covers running it.

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 (Wave Forcing).

Coupling 

At every coupling time \(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 (1) as SWAN’s sum over the bins,

(5)\[\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 \(\mathbf F_k\). The next force is not known until \(t_{k+1}\), so between couplings the circulation sees ADCIRC+SWAN’s schedule, a linear ramp from \(\mathbf F_k\) toward the extrapolated \(2\mathbf F_k - \mathbf F_{k-1}\):

(6)\[\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}\]

(\(\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 (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, \(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 (The Spectral Grid).

Frequency spacing. Keep the ratio between neighbors near 1.1,

(7)\[r = \left(\frac{f_{max}}{f_{min}}\right)^{1/(N_f - 1)},\]

with \(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 \(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 (Precision and 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 (Source Terms); 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 (ST6). Its coefficients are fixed, and setting sources.whitecapping with it is refused.

The wind the waves feel is relative to the current, \(\mathbf U_{10} - s\,\mathbf U_c\) (64), with \(s\) set by sources.wind.current_coupling; the circulation’s own wind stress uses the absolute wind (Wind). The default, 1, is SWAN’s form; the atmosphere partly adjusts to a moving surface, so a smaller value may be more realistic. \(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 \(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 \(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 \(n\) at every coupling, as ADCIRC does for SWAN’s READINP ADCFRIC. The value ADCIRC passes is a roughness length (59), which SWAN uses as \(k_N\) unchanged (Shared Terms). 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² s^-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 (Output Files). They follow SWAN’s definitions with CF metadata (Output Parameters). 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 (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.