Baroclinic Coupling

Cocoa can couple a depth-averaged (2DDI) hydrodynamic solve to an external 3D ocean climatology through three terms: a horizontal baroclinic pressure gradient (BPG) with its free-surface correction, an internal-tide wave-drag friction enhancement, and a dispersion friction correction. All three read a single preprocessor-generated nodal forcing file (Configuration) and share a common activation ramp.

The method follows the 2DDI approach of Pringle et al. (2019) [Pringle2019] and its global application [Pringle2021], as refined in ADCIRC by C. Blakely, whose global tidal dissipation analysis [Blakely2022] underpins the internal-tide drag’s role in the global model. The internal-tide and dispersion terms are carried over from ADCIRC’s nodalattr.F and internaltide.F90.

Global depth-averaged velocity after 30 days of tidal plus baroclinic forcing — Fig. 6 Depth-averaged velocity magnitude on a 3.7M-node global shell, forced by tides plus a single frozen (diagnostic-mode) RTOFS density snapshot, at day 30 of spin-up. The western boundary currents and the Antarctic Circumpolar Current are resolved.

Overview 

A depth-averaged model has no vertical structure, so it cannot represent the pressure gradient a real 3D density field produces, the drag a stratified ocean exerts on the barotropic tide through internal-wave generation, or the extra bottom friction implied by unresolved vertical shear. Cocoa closes these gaps with fields computed offline from a 3D ocean climatology (HYCOM GOFS 3.1) by the cocoa_ocean3d Python preprocessor and interpolated onto mesh nodes as a CF NetCDF time series:

The BPG and its free-surface correction (Baroclinic pressure gradient) enter GWCE and momentum as an extra pressure-gradient forcing, in the same way as the barotropic term.
Internal-tide wave drag (Internal-tide wave drag) adds a symmetric 2x2 friction-tensor contribution, active only where internal tides can propagate freely.
Dispersion as friction (Dispersion as friction) adds a further friction-tensor contribution representing depth-averaging error.

All three are gated independently under forcing.baroclinic (see Configuration) and share the file’s BPG, sigma-t, buoyancy-frequency, and mixed-layer products; enabling internal-tide drag or dispersion without the parent enabled: true is rejected at startup.

Theory 

Baroclinic pressure gradient 

The 2DDI BPG at a mesh node is a depth-averaged, horizontal derivative of the baroclinic pressure anomaly

\[\mathrm{BPG}_x = \frac{g}{\rho_0 H}\int_{-H}^{0} \int_{-H}^{z}\frac{\partial \rho}{\partial x}\,dz'\,dz,\]

computed by the preprocessor from a TEOS-10 (GSW) in-situ density profile via one-sided horizontal differences and a trapezoidal vertical integral (cocoa_ocean3d.bpg). \(\mathrm{BPG}_y\) is defined identically in the meridional direction. Both are staggered, with BPGX at longitude midpoints and BPGY at latitude midpoints, and the vertical integral is truncated at the shallower of any two neighboring wet columns, never extrapolated past the shallower column’s seabed. The preprocessor clamps the horizontal density-anomaly bracket to \(\pm 0.1\) kg/m ³ before scaling by \(g/\rho_0\), bounding \(|\mathrm{BPG}| \le g/\rho_0 \times 0.1 \approx 9.8\times10^{-4}\) m/s ².

The BPG enters GWCE and momentum with the same sign convention as the barotropic pressure gradient (Physics::BaroclinicGradient::compute_momentum_contribution / compute_gwce_contribution): the momentum contribution is area-weighted and Crank-Nicolson time-centered, each node of element \(e\) receiving its own nodal value,

\[\Delta(\text{momentum RHS})_i = -\frac{\Delta t}{2}\, A_e\, \bigl(\mathrm{BPG}_i^{\,n} + \mathrm{BPG}_i^{\,n+1}\bigr),\]

so that dividing by the node’s wet area recovers \(\mathrm{BPG}_i\) exactly, as ADCIRC’s per-node DBCPDX1A/2A/3A do (a quadrilateral weights each corner by its lumped area weight in place of \(A_e\)). The GWCE contribution is the negative element average of the nodal products, \(-\overline{H_i\,\mathrm{BPG}_i}\). Time-centering the BPG matches Cocoa’s other forcing terms (atmospheric pressure, tide potential, wind stress).

Free-surface correction 

The BPG above is computed from a fixed reference density profile and does not know the model’s own free surface. Pringle et al. [Pringle2019] correct for this by adding the elemental gradient of the surface potential-density anomaly (sigma-t) times the live free-surface elevation,

\[\Delta\mathrm{BPG}_{\mathrm{FS}} = \frac{g}{\rho_0}\, \nabla\bigl(\sigma_t\,\eta\bigr),\]

evaluated with the element gradient coefficients the pressure-gradient term already uses (the element tables’ spatial_grad_a/spatial_grad_b, scaled by the momentum map factors). The term is recomputed every step from the current eta and the time-interpolated sigma_t and added into the effective bpg_x/bpg_y fields before either consumer reads them, so the momentum and GWCE contributions both include it. A node whose surface sigma-t is the missing-value sentinel (see Nodal forcing file) zeroes the whole element’s correction.

Internal-tide wave drag 

A barotropic tide flowing over rough topography in a stratified ocean radiates internal waves, extracting momentum from the barotropic flow. Cocoa represents this as an additional symmetric 2x2 contribution to the friction tensor (Momentum Equations), driven by a user-supplied topographic nodal attribute and the file’s runtime bottom/mean buoyancy frequencies (\(N_b\), \(N_m\)). Three forms are supported, selected by the attribute’s column count (physics.internal_tide_friction: mesh, an existing NodalAttributeReader variable named internal_tide_friction):

Nycander (2005) [Nycander2005], 5 columns (“nonlocal”): columns are the coefficient \(C\) and two topographic-slope pairs \((H_x, H_y)\), \((J_x, J_y)\),

\[\mathrm{Fit} = C\,\frac{N_b}{H}\sqrt{1 - f^2/\omega^2}, \qquad \gamma^2 = \frac{H_x^2 + H_y^2}{(\omega^2-f^2)/(N_b^2-\omega^2)}, \qquad \mathrm{Fit} \leftarrow \frac{\mathrm{Fit}}{\max(1,\gamma^2)},\]

with tensor add \((2\,\mathrm{Fit}\,H_xJ_x,\ 2\,\mathrm{Fit}\,H_yJ_y,\ \mathrm{Fit}(H_xJ_y+H_yJ_x))\) for \((xx, yy, xy)\).
Lyard et al. (2006) [Lyard2006], 3 columns (“local”), and its Zaron & Egbert (2006) [ZaronEgbert2006] 2-column scalar form share one drag coefficient,

\[\mathrm{Fit} = C\,\frac{\sqrt{(N_b^2-\omega^2)(N_m^2-\omega^2)}}{\omega} \Big/ \max\Bigl(1,\ \frac{\omega^2-f^2}{N_b^2-\omega^2}\Bigr),\]

isotropic (\(C_{xx}=C_{yy}=\mathrm{Fit}\cdot\text{slope}^2\), \(C_{xy}=0\)) for the scalar form, or the full tensor (\(C_{xx}=\mathrm{Fit}H_x^2\), etc.) for the 3-column form.

Both forms act only where internal tides can propagate freely: \(\omega < N_b\) (supercritical stratification) and \(\omega > |f|\) (above the local inertial frequency), \(\omega\) being the M2 tidal frequency (Constants::omega_m2(), \(1.4052\times10^{-4}\) rad/s) throughout. Outside those bounds the tensor add is exactly zero. A 1/3-column attribute without baroclinic coupling (forcing.baroclinic.internal_tide.enabled: false) is the static form: the raw attribute values are added to the friction tensor directly, with no \(N_b\)/\(N_m\)/de-tiding dependence. A 2- or 5-column attribute always requires the coupling (there is no static Zaron & Egbert or Nycander equivalent) and is rejected at startup otherwise; with the coupling enabled, a 1-column attribute is rejected in turn, since the scalar dynamic form needs the slope column.

Two conventions differ from the published formulation so that an ADCIRC fort.13 converted by cocoa-mesh-tools keeps working with the values it already has.

First, the coefficient: [Pringle2019] (Eq. 5) writes the drag tensor with a \(4\pi\) in the denominator and reports a tuned \(C_{it} = 4.0\). Neither ADCIRC nor Cocoa carries the \(4\pi\); both fold it into the attribute value, so the number in column 1 is \(C_{it}/4\pi\) (about 0.318 for the paper’s 4.0), and taking \(C_{it}\) straight from the paper overstates the drag by a factor of about 12.6.

Second, the shelf cutoff: [Pringle2019] sets \(\mathcal{C} = 0\) where \(h < 100\) m to prevent excessive dissipation on the continental shelves, while Cocoa applies no such depth gate (the only cutoffs are the propagation bounds above). internal_tide_friction is a user-supplied mesh attribute that cocoa_ocean3d does not generate, so zeroing it on the shelf is the responsibility of whatever produced the attribute. An ADCIRC fort.13 prepared for the published method already carries that zeroing; an attribute built from raw topographic slopes does not, and will dissipate energy across the shelf.

De-tiding 

The internal-tide and dispersion formulas above act on the non-tidal (subtidal) velocity, since the raw depth-averaged velocity is dominated by the barotropic tide itself. Cocoa extracts the non-tidal component with a 25-hour, 12-minute-sampled lagged boxcar mean (equivalent to ADCIRC’s UNTIDE, including its partial mean during the first 25 hours of spin-up), then forms a capped de-tiding ratio

\[\mathrm{adju} = \min\!\Bigl(\frac{D_{\mathrm{tidal}}}{D_{\mathrm{total}}},\ 100\Bigr), \qquad D = C_{xx}u^2 + 2C_{xy}uv + C_{yy}v^2\ \text{(evaluated at raw and de-tided velocity)},\]

that scales the internal-tide tensor contribution: a steady (non-tidal) flow has \(D_{\mathrm{tidal}}=0\), so the drag reduces to bare friction, while a pure oscillation has the de-tided velocity go to zero, so \(\mathrm{adju}\to1\). The ring state and the effective BPG pair are persisted, so a mid-window restart reproduces an uninterrupted run exactly (see Checkpoint format).

The ring buffer lives in device memory: float ring samples and double running accumulators take about 1.6 GB and 25.6 MB at 1.6M nodes.

Dispersion as friction 

Depth averaging discards a vertical shear dispersion term. Cocoa adds a further diagonal friction-tensor contribution, active only in deep water (dispersion.deep_water_depth, default 1000 m, tested against the local bathymetric depth rather than total depth so the gate does not chatter with the tide):

\[\mathrm{add}_{xx} = \mathrm{ramp}\cdot c_d\cdot \mathrm{adj}_u\cdot \frac{\lVert(\bar u,\bar v)\rVert}{H}, \qquad \mathrm{add}_{yy} = \mathrm{ramp}\cdot c_d\cdot \mathrm{adj}_v\cdot \frac{\lVert(\bar u,\bar v)\rVert}{H},\]

with \(\mathrm{adj}_u,\mathrm{adj}_v\) the same capped de-tiding ratio form as above, under a symmetric \([-100,100]\) cap since this ratio, unlike the internal-tide ratio, is a raw mean-over-instantaneous quotient that can be negative during a tidal reversal. Both ratios are applied under a single joint \(|u|,|v|>10^{-3}\) guard: a large \(u\) with a near-zero \(v\) still gets \(\mathrm{adj}_u=1\), not the ratio. The off-diagonal entry is never touched by this term. \(c_d\) comes from one of two modes (dispersion.mode):

mld (mixed-layer-depth power law): \(c_d = a\,(x^b - 1)\), where \(x\) is the file’s dimensionless mixed_layer_ratio, floored at \(4\times10^{-3}\) before the power. a/b are required config keys with no default; appropriate values are case-specific.
cdisp (direct coefficient): \(c_d\) is the file’s dispersion_cd value directly, clamped to \(\le 1.0\) at read time.

mixed_layer_ratio is the mixed-layer depth divided by the local wet depth and is at most 1 by construction (see Nodal forcing file), so with positive a and b the power law gives \(c_d \le 0\) and the term reduces friction. The positive-semidefinite clamp below keeps the combined tensor physical when this reduction, together with a tidal-reversal sign flip in \(\mathrm{adj}_u\)/\(\mathrm{adj}_v\), would otherwise drive a diagonal entry negative.

Positive-semidefinite clamp 

Internal-tide drag and dispersion are both added on top of the base Manning-friction tensor; either can, in principle, push the combined symmetric 2x2 tensor non-positive-semidefinite, an unphysical friction that would inject energy. Cocoa applies a closed-form eigenvalue floor (physics/FrictionTensorPsd.hpp) after every additive contribution:

\[\lambda_{1,2} = \frac{\mathrm{tr}}{2} \pm \sqrt{\left(\frac{xx-yy}{2}\right)^2 + xy^2},\]

with any negative eigenvalue floored at zero and the tensor reconstructed from the unchanged eigenvectors. The clamp is an exact no-op on the pure Manning tensor (\(xy=0\) makes both eigenvalues equal the non-negative diagonal entries), so a run with no internal tide and no dispersion reproduces Cocoa’s pre-baroclinic references exactly.

Preprocessor pipeline 

The cocoa_ocean3d Python package (part of Cocoa Tools (Python), in utils/cocoa_tools) produces the nodal forcing file in three stages tied together by pipeline.SnapshotRegridder:

Compute (compute.compute_snapshot): TEOS-10 (via the gsw package) in-situ density from practical salinity and in-situ temperature, using gsw.p_from_z for the depth-to-pressure conversion. BPGX/BPGY, surface sigma-t, \(N_b\)/\(N_m\) (bottom and depth-mean buoyancy frequency via gsw.Nsquared), and the mixed-layer ratio (a de Boyer Montegut et al. 2004 [DeBoyerMontegut2004] fixed 0.03 kg/m³ sigma-t-anomaly criterion) are all computed here.
Regrid (regrid.compute_weights/apply_weights): bilinear by default, falling back to a box average when the mesh’s local neighboring-element-centroid footprint spans more source grid cells than a configurable threshold. A fill guard (any contributing source cell non-finite) NaNs the regridded node value.
Write (writer.NodalForcingWriter): the CF nodal file (see Nodal forcing file), one call per snapshot. Appending to an existing file validates and extends its time axis instead of overwriting it.

The cocoa-ocean3d CLI wraps download/compute/regrid/all subcommands against HYCOM GOFS 3.1. regrid and all take the Cocoa mesh with --mesh and stamp its mesh id on the file. all --with-velocity also fetches the velocity fields and writes dispersion_cd, which dispersion.mode: cdisp requires; without it the file carries only the mld products.

Configuration 

Key	Default	Description
`forcing.baroclinic.enabled`	`false`	Master switch. There is no scalar-constant fallback: the nodal file is either present and read, or the feature is off entirely.
`forcing.baroclinic.filename`	(required if enabled)	Nodal forcing NetCDF path (see Nodal forcing file).
`forcing.baroclinic.ramp`	enabled, `1d`	One-sided cold-start spinup ramp (normalized tanh, 0 at `simulation.start_time` rising to 1 over `duration`), shared by the BPG/free-surface correction, the dynamic internal-tide drag, and dispersion, so `ramp_factor == 0` reproduces the pre-baroclinic friction tensor exactly. Independent of the global `forcing.ramp` (disabling one does not disable the other). A static `internal_tide_friction` attribute is not a baroclinic contribution (it is legal with no baroclinic file at all), so it is left unramped.
`forcing.baroclinic.internal_tide.enabled`	`false`	Activates the dynamic (\(N_b\)/\(N_m\)-driven) internal-tide drag. Requires `physics.internal_tide_friction: mesh` naming a 2/3/ 5-column `internal_tide_friction` nodal attribute. Gated separately from the BPG, but still requires the parent `enabled: true` (the shared file supplies \(N_b\)/\(N_m\)).
`forcing.baroclinic.dispersion.enabled`	`false`	Activates dispersion-as-friction. Requires `internal_tide.enabled: true` (dispersion rides on the same de-tided `adju`/`adjv` machinery internal-tide drag needs).
`forcing.baroclinic.dispersion.mode`	`mld`	`mld` or `cdisp`; see Dispersion as friction.
`forcing.baroclinic.dispersion.deep_water_depth`	`1000.0` [m]	Deep-water gate on local bathymetric depth.
`forcing.baroclinic.dispersion.mld_coefficients.a`/`b`	(required if mode: mld)	Power-law coefficients applied to the dimensionless `mixed_layer_ratio`; no default exists, and appropriate values are case-specific. Any positive pair produces \(c_d \le 0\) (see Dispersion as friction).

forcing:
  baroclinic:
    enabled: true
    filename: "baroclinic_forcing.nc"
    ramp:
      enabled: true
      duration: 12h
    internal_tide:
      enabled: true
    dispersion:
      enabled: true
      mode: mld
      deep_water_depth: 1000.0
      mld_coefficients:
        a: 0.01
        b: 0.5
  tide:
    # ... tidal forcing, unaffected by baroclinic coupling ...

physics:
  internal_tide_friction: mesh   # required when internal_tide.enabled

Nodal forcing file 

CF-1.8 NetCDF, dims (time, node), written by the preprocessor. Cocoa interpolates linearly in time between the two snapshots that bracket the model time.

Variable	Units	Required when
`time`	CF numeric (seconds since an epoch the file itself declares)	always
`bpg_x`, `bpg_y`	m s^-2	always
`sigma_t_surface`	kg m^-3	always
`n_bottom`, `n_mean`	s^-1	`internal_tide.enabled`
`mixed_layer_ratio`	1 (dimensionless)	`dispersion.enabled` and `mode: mld`
`dispersion_cd`	1 (dimensionless)	`dispersion.enabled` and `mode: cdisp`

One file covers every product; the conditional variables above let a run enable any subset without forking the file format. A frozen file (exactly one time value) applies for the whole simulation window (diagnostic mode) and is exempt from the reader’s time-coverage validation.

Every variable declares _FillValue = NaN, but only sigma_t_surface may carry it in a well-formed file. The preprocessor zero-fills BPG/N/dispersion products at dry columns (so a coastal bilinear blend fades toward zero forcing), while surface sigma-t is NaN at any node whose regridding stencil touches a dry/land source cell. NaN in any other variable is a fail-loud file error naming the variable and the first offending node; the free-surface-correction kernel is the one consumer of the sigma_t_surface NaN sentinel (see Free-surface correction).

Global attributes: Conventions (“CF-1.8”), title/institution/ source/history (CF provenance), mesh_num_nodes (validated against the reader’s own count, fail loud on mismatch), mesh_id (the mesh’s cm1_... fingerprint, which cocoa-ocean3d always stamps; a file written for a different mesh is rejected, and a file without the attribute is accepted on the node count alone), source_dataset, utility_version, generation_command.

The file stores every product in fp64. The model stores the nodal fields as float and does all arithmetic on them in double, including the free-surface-correction accumulation. Float storage loses nothing the source data holds: an OGCM analysis differentiated onto the mesh carries four to five significant digits.

Regression checks 

Zero forcing reproduces the pre-baroclinic baseline exactly: an all-zero nodal forcing file with baroclinic enabled matches the same configuration with baroclinic disabled (Integration_WNAT_baroclinic_zero_forcing), and so do a static internal-tide attribute of zeros and a disabled dispersion path. The composed integration test (Integration_WNAT_baroclinic_composed) checks the other direction: with the BPG, free-surface correction, internal-tide drag, and dispersion all active, the run must differ measurably from its baroclinic-disabled counterpart, and it is compared against a committed reference at the 1e-6 cross-run tolerance of the other WNAT implicit references.

Checkpoint format 

The internal-tide de-tiding ring buffer’s state (ring samples as float with two velocity components, the double running accumulator, the ring write position, whether the buffer has completed a full cycle, and the last-sample simulation time) is persisted in Cocoa’s checkpoint format (Checkpoint / Restart) when internal_tide.enabled is set. The effective BPG pair and its previous level are persisted whenever baroclinic forcing is on, as f8 regardless of the in-model float storage; GWCE consumes the value composed during the previous step and momentum centers against the step before that, so both levels carry across a step boundary. These fields arrived in checkpoint format version 4, and a checkpoint from a run without the corresponding forcing carries none of them. A mid-window restart reproduces an uninterrupted run exactly.

Full citations are collected in References.