1The inverse problem and what makes it hard

Let horizontal position be x = (λ, φ), longitude and latitude, and let d ≥ 0 be depth in metres, positive downward. The mandated output depths are

The task is to estimate T(x, t, d) for every d ∈ D from surface information alone. Write the causal encoder-decoder as

The history is strictly causal: no field dated after the issue time enters h_t. The encoder output h_t is the embedding-engine artifact the problem statement asks for, and is exported with its shape, stride, temporal support and weights.

Under squared error the optimal predictor is the conditional mean, and the error decomposes exactly:

This is an identity, not an estimate of either term. It is worth stating because it names the obstacle precisely: different subsurface states can share almost identical surface expressions. Two columns with the same SST, SSS, sea level and surface currents can differ by more than a degree at 100 m. The first term is what no surface-only method can remove; the second is the only one modelling effort can touch. We make no claim to have measured the split for this basin — a failed displacement regressor is not a measurement of an irreducible floor.

2Data and permitted information

Seven time-varying surface fields are harmonised onto one 0.25° daily grid of 101 × 241 cells, of which 11 736 are ocean.

The input tensor for one day carries 246 channels, composed exactly as:

• 234 temporal channels = 26 daily frames × 9 fields. Each frame holds the seven surface fields, a derived wind-curl channel, and the land/ocean mask, which is always last within a frame and never scaled.
• 12 static channels = 7 surface-climatology fields (one per input variable) + 2 normalised coordinate planes + 2 day-of-year channels (sin 2π·doy/365.25, cos 2π·doy/365.25) + 1 bathymetry plane.

The climatology is fit-only, 1993-2020. Any climatology estimated over a window that includes evaluation years leaks the answer, and this one does not. Enforcing it costs 1.43 % of skill, measured against a repaired control rather than assumed, and that cost is paid rather than avoided.

At inference no subsurface observation of any kind is read. Bathymetry is a static external prior and is disclosed as one. A surface field dated today but released days later cannot support a same-day operational claim, so the practical operating point is roughly today minus five days, set by each satellite product's own latency.

3The observation operator and its physics

3.1 From what a float measures to what the model predicts

A float reports adjusted pressure p, in-situ temperature t, and practical salinity S_P. None of those is metres, and none is the model's quantity. Using TEOS-10 throughout:

GSW takes pressure in dbar and position in degrees. Those unit boundaries are recorded explicitly because crossing one silently is the classic way this calculation goes wrong.

3.2 Two temperature conventions cannot be differenced

The model is fitted to in-situ ITS-90 Argo temperature. GLORYS12 ships potential temperature θ. Differencing one against the other is a units error wearing the same name, and it is not repaired by calling both arrays "temperature". The comparator is converted with the full chain:

Measured over our rows, t − θ reaches +0.108 °C at 1000 m (maximum 0.1379), which is comparable to the margins being argued about at depth. Section 4.1 of our protocol forbids borrowing the held-out Argo salinity to improve a comparator's conversion, so WOA23 monthly climatological salinity is used, present at all fifteen depths exactly with no vertical interpolation. The residual from using a climatological rather than the matched salinity field is bounded by redoing the conversion at S ± 0.5 psu: 0.0008 °C, two orders below the margins in question.

This conversion makes GLORYS worse (0.6941 → 0.7009), so it favours us and is held to a higher standard for that reason. It is performed because the physics requires it, and it would be required if it went the other way.

3.3 Interpolation support, and refusing to invent a label

Let I_i(d_k) be linear interpolation between the two valid native levels bracketing d_k. The observation validity indicator m^obs_ik is set to zero when the target has no bracket, fails QC, crosses an excessive gap, or lies below the seabed. No endpoint extrapolation is permitted, and the bracket limits are the INCOIS manual's own: <50 m spacing above 100 m depth, <100 m below.

This is not a technicality. One float with valid levels at 6.8 and 838.2 dbar was previously interpolated to produce a 125 m "observation" — a number with no measurement behind it. Sixteen such rows out of 278 carried 97 % of the squared error at 500 m. A method evaluated against invented labels is measuring its agreement with a straight line.

For a pairwise comparison against product c the mask is m_ik = m^obs_ik f^OV_ik f^c_ik, the product of observation validity and both products having a finite prediction there. The operator is held fixed between training variants and every benchmark candidate. Because each comparator has its own finite-row intersection, our own mean14 differs between contrasts — 0.6778 against INCOIS, 0.6885 against GLORYS12 — and that is a property of the row sets, not of two different models.

Observation, representation and model error all contribute to e. The common practice of subtracting an assumed observation-error variance to manufacture a "true model RMSE" requires zero covariance between the terms, which assimilating products violate by construction. We do not do it.

4The metric, and optimising the thing being judged

For the fourteen subsurface depths:

This is the mean of per-depth RMSEs, not the square root of a pooled MSE. They are different statistics and only the first is the quantity the problem statement describes. Pooled RMSE is dominated by whichever depths happen to be best sampled; thickness-weighted error answers a third question again. All three are kept distinct and only J₁₄ is reported as the headline. Bias is \(b_k = \sum_i m_{ik} e_{ik} / n_k\); correlation is Pearson over the same pairs, computed depth by depth, because a correlation pooled across depths is dominated by the mean vertical gradient and is a weak diagnostic.

With s₀ = 1 °C a fixed unit scale, the dimensionless sparse objective mirrors the metric rather than approximating it:

with ε = 10⁻⁴ a numerical stabiliser. Batches are sampled so deep levels appear regularly. Validation sufficient statistics are averaged over the full validation set before the square root is taken — averaging per-batch RMSE is a different statistic and a slightly optimistic one.

For dense GLORYS labels G, with wet masks q_jk and cell areas A_j ∝ cos φ_j:

Area weighting is correct for a grid objective and must not silently replace the declared float-sampled evaluation metric. The land mask is applied before any reduction. The full-data gradient of the reported metric for a valid row is

The 1/(n_kR_k) factor is the whole point: a pooled-MSE loss weights a row by nothing, whereas the reported metric weights it inversely by how well-sampled and how well-predicted its depth already is. Sparsely sampled deep levels therefore carry more gradient per row, which is what the metric asks for.

Objective conflict is monitored rather than assumed away, via the cosine between the two gradients on the shared encoder:

Repeated negative c indicates local disagreement between the two supervisions, not proof that the reanalysis is globally harmful.

5Architecture

The primary member is a 2-D encoder-decoder over the basin, treating the 246 channels as feature depth and the grid as spatial extent. It is spatial by design: a per-point mapper cannot see that a pixel sits on the warm flank of an eddy, and section 10 measures exactly what that costs.

Emitting all fifteen depths from one head lets the network share vertical structure across levels rather than solving fifteen unrelated regressions. The embedding h_t is the 56-channel activation immediately before the final 1×1 convolution — the representation a single linear readout turns into T(z).

That the embedding carries subsurface information beyond its own inputs is a measurement, not a design claim: a ridge readout on a held-out eastern longitude band scores R² +0.4749 from the 56 numbers against +0.4272 for a PCA control at equal width. An earlier version of this work quoted 0.4121 against 0.3216 from a docstring rather than a measurement; the conclusion survived, the margin halved.

6The algorithm

6.1 Training

Two properties of that loop are load-bearing. First, stop_metric is an explicit argument, so the P0 lineage — the problem statement's own supervision path — can be run with Argo touching neither the gradient, nor the argmin, nor the early-stop horizon. A run declaring --argo-weight 0 while still stopping on Argo is refused by a guard, because the exposure is otherwise invisible: it amounts to two epochs of selection signal, which is small but is not zero and must not be reported as zero. Second, the climatology used to build C is fit-only by construction, not by convention.

6.2 Inference

The curl line carries a real bug that was found and fixed. The helper that builds an inference window subtracted the climatology before assembling frames, so it computed the curl of the wind anomaly, while the dataset that trained the model used absolute wind. The difference is curl(u_clim): 0.2 % at the median but 34 % at the peak. Production reaches the network through the dataset path so the headline is unaffected, but two other entry points used the helper and were silently feeding the model a channel it was never trained on. A parity test between the two paths now runs in CI.

7Ensemble and blending

Per-depth nonnegative weights are fitted on a dedicated blend year, with shrinkage toward a prior w_0k:

A member trained on the blend year's labels is ineligible for this procedure. With Σ_k the covariance of member errors and b_k their bias vector,

which is why an individually weaker but decorrelated member can help while twenty correlated ones do not. The unconstrained minimum-variance solution Σ⁻¹1/(1^TΣ⁻¹1) is diagnostic only; the constrained regularised fit above is what ships.

8Dynamical features derived from the permitted inputs

Deriving features from the seven inputs preserves their information source; reading a new observed field at inference would change the task. Two families are used, and both are stated precisely because both are easy to get subtly wrong.

8.1 Wind forcing

Wind-velocity curl and wind-stress curl are different quantities with different units, and conflating them is common. With a documented bulk drag law,

Stress is in N/m². For eastward/northward components (A_λ, A_φ) the radial curl on a sphere is

With f = 2Ω sin φ, the large-scale Ekman pumping diagnostic, positive upward, is

Keeping f inside the derivative is not optional. Pulling it out discards its latitude dependence, which over a 25° meridional span is a large error. The domain starts at 5°N so f = 0 is avoided, but the low-latitude edge still needs one-sided or masked stencils, and invalid stencils are excluded rather than filled.

8.2 Sea level and currents

A current product derived from the same altimetry as η is not independent information, though its representation can still help learning. Vorticity \(\zeta = \partial_x v_s - \partial_y u_s\), divergence, and front magnitude \(\lvert\nabla \mathrm{SST}\rvert\) are computed on physical grid distances. These describe fronts, rotation and forcing; they do not uniquely determine the water column, which is equation (2) restated in physical terms.

Sea surface height contains steric, mass-loading and reference-depth contributions, so a temperature-only integral cannot be forced to equal observed SSH. No constraint of that kind is imposed.

8.3 Stratification, and why a monotonicity prior would be wrong here

In the Bay of Bengal, river discharge creates a fresh surface layer that can sit above warmer water, producing genuine temperature inversions beneath a shallow density mixed layer. A global rule forcing temperature to decrease with depth would be unphysical in exactly the region the problem statement cares most about. None is imposed. Where salinity is available the local buoyancy frequency is

and differentiating in-situ density with depth — thereby reading compressibility as stratification — is a known trap that is avoided. Without credible subsurface salinity labels any stability penalty is disabled: surface salinity alone cannot certify a density profile.

8.4 Vertical derivatives on an unequally spaced grid

The fifteen levels are not equally spaced — 5 m near the surface, 300 m at depth — so a uniform stencil is wrong by a factor of sixty across the column. With Δd_k = d_k+1 − d_k:

A linear profile must return exactly zero discrete curvature at interior levels; this is a unit test, not an assumption.

8.5 Why the thermocline is the hard band

A downward displacement ξ > 0 of a background profile has the first-order effect

Where |∂_dB| is large, a modest depth error becomes a large temperature error: a 10 m displacement across a 0.1 °C/m gradient is about 1 °C. That is the entire explanation for the shape of our error curve, for the shape of every comparator's error curve, and for why the gap to GLORYS2V4 sits almost entirely between 50 and 150 m. The approximation fails for large shifts or changing profile shape.

9Evaluation protocol

Five gates must pass before any number is reportable. They are adversarial by construction: each was written to make a specific way of fooling ourselves impossible, and each was verified to still reject on its failure modes after being made evidence-aware.

9.1 Simultaneous inference over a predeclared family

Fourteen depths × several comparators is a multiple-comparison problem. A per-depth 95 % interval repeated fourteen times is not a 95 % statement about the set. Define δ̂_j = R_j^c − R_j^OV over the predeclared family j, resample whole WMO floats with replacement retaining all their cycles and depths, and form

A win is claimed only where ℓ_j > 0. Floats are the resampling unit because profiles from one float are strongly dependent — treating rows as independent would inflate the sample by roughly the number of cycles per float and manufacture significance. This is an approximate simultaneous one-sided basic/pivotal construction, distinct from a percentile interval; for a single comparison it reduces to \(2\widehat{\delta} - Q_{.95}(\widehat{\delta}^{*})\). 5 000 resamples, seed 26066, one shared resample across the whole family.

9.2 Determining the comparator's analysis window without circularity

The INCOIS product ships no CF time bounds, so its timestamp convention was unknown and had been assumed. Tuning the offset against our own error would be circular — it would select whatever window flatters us. Instead the window was determined from the product's own Argo inputs: for each candidate offset, measure the misfit between the analysis field and the profiles collected in that window. The answer is then a property of the product.

Over offsets (−12 … +2) the minimum is at −8, i.e. a window of t−8 to t+2. Two controls: a permutation test pairing each field with another cycle's profiles gives p = 0.005, and the curve is resolved at 9.8× the standard error of the mean. The assumed centre was wrong. Residual exposure from the window choice is 0.0395 °C.

10Results

96 delayed-mode floats, 69 183 support-verified rows, 2023-24, identical rows per contrast.

The last two rows are from the external ladder on 10 091 pre-repair rows and are not directly comparable in absolute value with the rows above; the ranking is what they establish.

10.1 Decomposing the surface-only comparison

Running OSnet twice separates two claims that are usually conflated. Of its 0.3043 °C deficit under its published inputs, 0.2507 is information — two fields and one snapshot versus eight fields and a 26-day history — and the remaining 0.0536 is architecture. Handed identical information, the per-point mapper still loses at 13 of 14 depths with all thirteen supported, which is better-supported than our headline INCOIS result. That is the cleanest evidence in this work that the spatial treatment, not merely the input budget, is doing work.

10.2 Where the error lives

RMSE peaks at 1.24 °C at 100 m and falls to 0.27 °C at 1000 m. Via equation (19) this is the expected signature of displacement error acting on a steep background gradient. The same band is the whole of the gap to GLORYS2V4 (0.9572 against our 1.0824 at 75 m) and the gate on the hazard indicators of the previous page. One physical problem, three consequences.

11Limits and failure modes

11.1 Errors found by measuring what had been assumed

The following changed a number and are listed because a method that only reports its successful checks has not been checked.

• A comparison mixed potential and in-situ temperature. Correcting it helped us (6/14 → 8/14), which is why it is listed rather than absorbed.
• "Climatology leakage is essentially free" was withdrawn: the label cache was keyed without the climatology, so both arms saw identical labels. The broken run landed within 0.06 pp of the true answer by pure noise — the tell was one seed disagreeing in sign.
• A compliance report quoted an R² from a docstring rather than a measurement.
• A scorer had a network override but no blend override, so a locked weight file sat unread while production weights were scored.
• A required contract item could not fail, because the status line scanned only items that existed. Requirements are now predeclared.
• A monitoring command matched its own process, making a gate unfalsifiable from outside; and a memory gate was set eight times too high from an unattributed reading, idling a shared GPU for hours.

12Resources

12.1 Symbols

x, d, t	horizontal position; depth in metres positive downward; time
Xt, Mt	surface field stack and its validity mask at time t
ht	encoder output, the 56-channel embedding artifact
mik	joint validity: observation support × both products finite
eik, Rk	residual at profile i depth k; per-depth RMSE
J14	mean of the 14 subsurface per-depth RMSEs (≠ pooled RMSE)
ℓj	simultaneous one-sided lower bound for comparison j
θ, Θ, t	potential, Conservative and in-situ temperature
SP, SA	Practical and Absolute Salinity
ξ	vertical displacement of the background profile

12.2 Constants

ρ₀ = 1026 kg m⁻³ · c_p = 3985 J kg⁻¹ K⁻¹ · Ω = 7.2921 × 10⁻⁵ s⁻¹ · R = 6 371 km · s₀ = 1 °C · ε = 10⁻⁴ · bootstrap seed 26066, 5 000 resamples.

12.3 References

1. McDougall & Barker. Getting started with TEOS-10 and the Gibbs Seawater (GSW) Oceanographic Toolbox. SCOR/IAPSO WG127. Conversion APIs, unit conventions and the z_from_p, SA_from_SP, CT_from_pt, t_from_CT and Nsquared routines used in §3 and §8.3. teos-10.org
2. Pauthenet, Bachelot, Balem, Maze et al. Four-dimensional temperature, salinity and mixed-layer depth in the Gulf Stream, reconstructed from remote-sensing and in situ observations with neural networks. Ocean Science, 2022. The OSnet design reimplemented and retrained as the surface-only rival in §10.1. os.copernicus.org
3. Wang et al. Attention enhanced 3D-U-Net++ ocean temperature and salinity reconstruction in the northwestern Pacific based on transfer learning. Earth System Science Data, 2026. Precedent for observation pretraining followed by reanalysis fine-tuning. essd.copernicus.org
4. Chae, Donohue & Park. TS-Cast. Ocean Science, 2026. Surface-height ambiguity and observation-based reconstruction; the basis for §8.2's refusal to constrain a temperature integral to SSH. os.copernicus.org
5. Shay, Goni & Black. Effects of a warm oceanic feature on Hurricane Opal. Monthly Weather Review, 2000; and Mainelli et al., Weather and Forecasting, 2008. Tropical cyclone heat potential and its intensification threshold, used on the Hazards page.
6. Hobday et al. A hierarchical approach to defining marine heatwaves. Progress in Oceanography, 2016. The 90th-percentile, five-day definition.
7. Mechanisms of Barrier Layer Formation and Erosion from In Situ Observations in the Bay of Bengal. Journal of Physical Oceanography, 2019. Regional salinity stratification; why no monotonicity prior is imposed in §8.3. ametsoc.org
8. Kessler. The circulation of the eastern tropical Pacific: dynamics. NOAA PMEL. Ekman transport and pumping derivation used for equation (15). pmel.noaa.gov
9. Shalizi. Simulation for Inference I: The Bootstrap. Carnegie Mellon University, 2018. Basic/pivotal interval derivation; equation (20) applies that principle to a predeclared comparison family. stat.cmu.edu
10. Jean-Michel et al. The Copernicus Global 1/12° Oceanic and Sea Ice GLORYS12 Reanalysis. Frontiers in Earth Science, 2021. Dense training target and one comparator. frontiersin.org
11. Guinehut, Dhomps, Larnicol & Le Traon. High resolution 3D temperature and salinity fields derived from in situ and satellite observations (ARMOR3D). Ocean Science, 2012. The product that leads the ladder in §10, and which interpolates the Argo it is then scored against. os.copernicus.org
12. Argo. Argo float data and metadata from Global Data Assembly Centre. SEANOE. Delayed-mode profiles, QC ≤ 2, are the evaluation truth throughout. doi.org/10.17882/42182

12.4 Reproducibility

Every number on the Validation page is produced by a script that refuses to emit a result when its preconditions fail: the row cache must declare the post-repair baseline, the as-shipped arm must reproduce the published headline before any revision is accepted, per-depth verdict flips must reconcile against the aggregate, and the production bundle is checked file-by-file against recorded SHA-256 digests (10 files, 0 changed, 0 missing at time of writing). The suite is 295 tests, 0 failures, including the physical unit tests named in §8.4 and a round-trip test on the temperature conversion of §3.2.