Mathematical description of the grid_swath_ssha.py pipeline: the expanded swath grid, the observation operators, the optional background / coherent internal-tide decomposition, and the optimal interpolation of the residual. Equations are written in plain LaTeX-style notation.

1. Coordinate frame and notation

All fitting and mapping happens in a local azimuthal-equidistant (aeqd) plane centered on the grid centroid: T: (lon, lat) -> (x, y) in km, shared by the grid nodes and every observation, so distances are consistent across all stages. Times t are measured in days relative to the analysis time t0 (the SWOT segment mid time, or --target-time; on the regional cube, each frame time t_k). Each observation carries its native measurement time t_i — not a granule-center, pass-center, or analysis-frame stamp.

An observation i is (x_i, y_i, t_i, eta_i, sigma_i): position, time, SSHA value (m), and observational error (m). Weights are w_i = 1 / sigma_i^2 throughout.

2. Expanded swath grid

One grid row per along-track line of the SWOT segment. The nadir centerline c_j = (lon_j, lat_j) and the per-line along-track azimuth alpha_j (circular mean of the geodesic azimuths to the neighboring lines) define the cross-track axis: grid nodes are placed geodesically (WGS84) at

g_{j,k} = fwd(c_j, alpha_j + 90 deg * sign(chi_k), |chi_k|),

where chi_k = k * dx is the signed cross-track coordinate (dx = 2 km, negative = left of flight, SWOT convention). By default chi spans +/- (ceil(edge/dx)*dx + expansion), i.e. the native swath edge (~68 km) widened by --expansion-km (100 km) on both sides. --target-width-km W instead places nodes at chi = -W/2 ... +W/2 (300 km total → ±150 km, 151 columns at 2 km posting).

3. Observations and their errors

SWOT KaRIn (Expert): eta = ssha_karin_2 + height_cor_xover + internal_tide_hret. The HRET internal tide is added back because the L2 product removes it while the RADS sla formula never does — both sources then carry the same internal-tide signal. Quality masks: ssha_karin_2_qual == 0, crossover quality <= 1, ocean depth < 0.

Nadir altimeters (RADS): eta = ssha (all standard corrections applied; barotropic ocean tide FES14 removed; internal tide retained). Each along-track sample keeps the NetCDF time value. With no RADS thinning (regional-cube default; also the 300 km pass mapper) that timestamp is unchanged.

Native timestamps. SWOT Expert L2 time is along-track (num_lines): one clock per line, broadcast across that line's cross-track pixels. A KaRIn line is acquired nearly simultaneously; Expert has no per-pixel clock. The OI and the harmonic tide use that per-line t_i (Δt = t_i − t0 in the covariance; ω t_i in A cos(kr − ω t) + B sin(...)). The 2nd-order background polynomial discards time by design. Map times (t0, or the cube's 6-hourly t_k) never overwrite observation times.

Thinning (both sources): non-overlapping blocks — N x N pixels for KaRIn (--swot-median-window), N consecutive along-track points per pass for nadir (--rads-thin-window). Each block becomes one observation

eta_block = median(eta_members)

and, by default (--swot-block-error floor),

sigma_block = max(sigma_floor, std(eta_members)),

with sigma_floor = --swot-obs-error (KaRIn) or --rads-thin-error-floor (nadir). --swot-block-error std uses std(eta_members) itself (sample std, ddof=1), skipping blocks with fewer than two valid pixels; uncertainty is max(std, 1e-4 m) so a numerically zero std cannot blow up the OI. The production 300 km pass mapper uses 5 x 5 SWOT blocks in std mode, nadir uncertainty 0.02 m with no RADS thinning, a 4-day OI time half-width, and a ±5.5-day observation search that can include other SWOT Expert passes. A SWOT block is placed at the center pixel when that pixel is valid, otherwise at the spherical-mean position, and tagged with the center line's time (or the median of valid times in the block if the center is invalid). A RADS block is placed at the spherical-mean position and median time. The scatter std is a pixel-based observational error: a block straddling a front or noisy region is trusted less.

Instrument bias corrections (optional, applied first): per-(mission, cycle, pass) constant median(eta) removed from nadir data (--fit-rads-pass-bias); per-pass, per-side polynomial in the scaled cross-track coordinate removed from KaRIn (--fit-swot-cross-swath-poly).

4. Large-scale background (--fit-internal-tides, stage 1)

A 2nd-order 2D polynomial, time discarded, fit to all observations in the wide window |t_i| <= T_IT (--time-half-width-days-internal-tides, 10 days):

B(x, y) = c0 + c1 x~ + c2 y~ + c3 x~^2 + c4 x~ y~ + c5 y~^2,
x~ = x / 100 km,  y~ = y / 100 km,

with coefficients from the weighted least-squares problem

min_c  sum_i w_i ( eta_i - B(x_i, y_i) )^2 .

The 100 km scaling keeps the normal equations well conditioned. Subtracting B leaves mesoscale + submesoscale + internal tides.

5. Coherent internal tide (stage 2)

Model

S point sources (--tide-n-sources, 6) at fixed, equally spaced azimuths theta_s = theta_0 + s * 360/S deg around the region center, each at a fitted distance D_s:

source position   p_s = D_s (cos theta_s, sin theta_s),
range             r_s(x) = | x - p_s |,
per-source field  eta_s(x, t) = a_s cos(k_s r_s - omega t) + b_s sin(k_s r_s - omega t),

with per-source wavenumber k_s = 2 pi / lambda_s (each source has its own wavelength) and the fixed constituent frequency omega = 2 pi / T (M2: T = 12.4206012 h; S2/N2/K1/O1 selectable). Equivalently

eta_s = A_s cos(k_s r_s - omega t - phi_s),   A_s = sqrt(a_s^2 + b_s^2),
phi_s = atan2(b_s, a_s).

The circular wavefronts encode the geometry the model is built for: at the region center the wavefront curvature radius equals D_s, so a nearby source produces visibly curved tidal beams across the region while a distant source degenerates to a plane wave (direction -theta_s, and D_s becomes unidentifiable — only the wave direction, wavelength, and complex amplitude matter, which is why D_s is bounded above).

Fitting (variable projection)

Free parameters: (D_s, lambda_s, a_s, b_s) for each source — 4S = 24 for the default 6. The amplitudes enter linearly: for fixed distances and wavelengths, build the design matrix

G(D, lambda)  with columns  cos(k_s r_s - omega t), sin(k_s r_s - omega t),

and solve the ridge-regularized weighted normal equations

c^ = ( G^T W G + eps I )^{-1} G^T W y,     eps = 1e-3 * mean(diag(G^T W G)),

where y is the background-removed SSHA and W = diag(w_i). The nonlinear search then runs only over the 2S = 12 parameters (D_s, lambda_s):

  1. Greedy initialization (matching pursuit): source by source, scan a coarse grid (8 log-spaced distances x 12 wavelengths within bounds) against the residual of the other sources, keep the best-scoring pair; two sweeps.
  2. Joint refinement: bounded trust-region least squares (scipy.optimize.least_squares, method trf) on the packed (D_1..D_S, lambda_1..lambda_S) with the inner linear solve embedded in the residual (so amplitudes are always optimal for the trial geometry).

Bounds: lambda_s in --tide-wavelength-bounds (default 80–250 km; mode-1 M2 is ~120–160 km in mid-latitudes) and D_s in --tide-distance-bounds (default 1.5 x region radius to 3000 km, keeping sources outside the data region — the model has no near-field physics).

Diagnostics and identifiability

  • at_bounds per source: a fitted D_s or lambda_s within 1% of a bound signals a degenerate or absorbing fit. In particular lambda pinned at the upper bound means the least-oscillatory basis available is soaking up non-tidal (mesoscale) signal — tighten the wavelength bounds.
  • Temporal observability: the weighted phasor concentration R = | sum_i w_i e^{i omega t_i} | / sum_i w_i. If R > 0.95 the sampling times cluster at one tidal phase (single snapshot, or repeat-orbit aliasing) and amplitude/phase are weakly constrained; a warning is emitted.
  • Variance explained is reported for both stages; note that with flexible bounds part of the "tide" variance is mesoscale leakage.

6. Optimal interpolation of the residual (stage 3)

The OI maps eta_res,i = eta_i - B(x_i, y_i) - eta_IT(x_i, t_i) for observations inside the OI window |t_i| <= dt (--delta-t-days). Without --fit-internal-tides it maps eta_i directly.

Covariance model. Space-time separable, isotropic:

C(r, tau) = sigma_sig^2 * rho(r; L) * rho_t(tau; L_t),

with the shape selectable per axis (--covariance-model, --time-covariance-model):

gaussian:     rho(r) = exp( - r^2 / (2 L^2) )
exponential:  rho(r) = exp( - r / L )
markov:       rho(r) = (1 + r/L) exp( - r/L )        (2nd-order autoregressive)

The scale L is set from the half-width flags with a shape-independent convention: rho(h) = 1/2 at the half-width h, i.e.

gaussian:     L = h / sqrt(2 ln 2)
exponential:  L = h / ln 2
markov:       L = h / 1.67835        (solves (1+x) e^{-x} = 1/2)

sigma_sig^2 is the variance of the observations about their median.

Local ordinary kriging. For each grid node, the --max-neighbors nearest observations (in L-scaled space-time distance) form the system

[ K + R    1 ] [ w  ]   [ c ]
[ 1^T      0 ] [ mu ] = [ 1 ],

where K_ij = C(r_ij, tau_ij) is the obs-obs covariance, R = diag(sigma_i^2) the observational error variances, c_i = C(r_i0, tau_i0) the obs-node covariances, and the Lagrange multiplier mu enforces the unbiasedness constraint sum w_i = 1 (unknown constant mean). Then

prediction        eta^(x0) = w^T y,
error variance    sigma_OK^2(x0) = sigma_sig^2 - c^T w - mu     (>= 0),

written to ssh_anomaly_oi and oi_error_variance. Because the obs-obs covariance is used, clustered redundant observations are down-weighted — ten co-located points carry roughly one point's worth of information (verified against a brute-force solve). Keep the neighborhood's data span at 2–3x the spatial half-width (thin the data accordingly) or the node-to-node switching of neighbor sets shows up as edges in the map.

7. Total reconstruction

eta_total(x) = B(x) + eta_IT(x, t = t0) + eta_OI-residual(x),

written to ssh_anomaly_total — the field directly comparable to ssha_swot. The components ssh_background and ssh_internal_tide are stored separately, and the fitted per-source parameters (azimuth, distance, wavelength, amplitude, phase), the variance explained by each stage, and the phasor diagnostic are in the NetCDF attributes.

8. Windows

Two half-windows govern the data: the internal-tide/background fits use +/- T_IT (10 days — long enough to see most of an M2 alias cycle under ~daily sampling), while the OI uses +/- dt (4–5 days) about t0. When the tide stage is off, only the OI window is loaded.

9. Regional lon/lat time cube

map_region_ssha.py (odsl-multi-altimeter-region) uses the same loaders, native t_i, background polynomial, point-source tide, and OI as above, on a regular lon/lat grid at many analysis times t_k (default 6-hourly). The background and tide are fit once on the padded window; each frame then maps the residual at t_k with Δt = t_i − t_k and stores

eta_total(x, t_k) = B(x) + eta_IT(x, t_k) + eta_OI-residual(x, t_k).

The 6-hourly axis is the map time. Observation clocks are unchanged.