Locally stationary Argo ocean heat content estimates: Modeling, validation and uncertainty quantification

Source: arXiv:2606.31957 · Published 2026-06-30 · By Thea Sukianto, Mikael Kuusela, Donata Giglio, Anirban Mondal, Pulong Ma, Douglas W. Nychka

TL;DR

This paper addresses the challenging problem of producing statistically rigorous and computationally tractable global maps of Ocean Heat Content (OHC) from Argo float data, along with principled uncertainty quantification. OHC is a critical climate metric reflecting Earth's energy imbalance and ocean warming, but Argo data are complex in structure with spatio-temporal sparsity and heterogeneity in vertical sampling. The authors propose a new end-to-end framework modeling vertically integrated Argo temperature profiles as a locally stationary Gaussian process with anisotropic space-time covariance. This local Gaussian process approach is combined with a climatological mean field including seasonal and long-term trends, leveraging data-driven decorrelation scales. A key novel contribution is a computationally efficient conditional simulation ensemble methodology that generates spatially and temporally correlated uncertainty ensembles, fully propagating uncertainty into integral OHC estimates over space and time. They validate modeling choices with a novel paired cross-validation method and demonstrate that models including climatological time trends and time in the covariance outperform simpler baselines without these features. They present new Argo-based OHC anomaly maps with uncertainty for 2004–2022, and show the framework’s modular open-source software enables reproducibility and extensibility. The resulting maps and uncertainty estimates improve accuracy and enable downstream climatological analyses previously impossible due to inadequate uncertainty quantification.

Key findings

• Including a climatological time trend in the mean field nearly doubles the estimated OHC trend relative to models omitting it (consistent with Cheng and Zhu 2015).
• Modeling time dependence in the covariance structure significantly improves cross-validated monthly median absolute prediction errors over models using only spatial covariance.
• Estimates use 1,731,330 quality-controlled Argo profiles from 2004–2022, divided into two vertical sections yielding 1,540,593 profiles for 15-975 dbar and 1,128,932 for 975-1850 dbar.
• Local Gaussian process regression with data-driven anisotropic decorrelation scales produces more accurate spatial interpolations than fixed-scale global covariance models.
• The proposed local conditional simulation algorithm efficiently generates ensembles that reproduce observed spatio-temporal covariance and lag-1 autocorrelation structures in Argo data.
• Paired cross-validation shows that uncertainty quantification from the conditional simulations reliably matches empirical prediction errors across space and time.
• The modular open-source code reproduces the entire pipeline from raw Argo GDAC data to final OHC maps and uncertainties on a 1°×1°×monthly grid.
• Uncertainty quantification for the full OHC vertical integral uses a conservative upper bound due to no vertical covariance modeled between sections.

Methodology — deep read

1. Threat Model & Assumptions: The problem is to estimate global and regional Ocean Heat Content (OHC) from noisy, irregularly spaced Argo float observations sampling the upper 2000 m ocean. The adversary analogy is not relevant here; rather, the statistical challenge is to estimate a spatio-temporally correlated field with heterogenous vertical sampling and produce accurate uncertainty quantification.

2. Data: The input data consist of 1,731,330 Argo profiles from 2004–2022, quality-controlled using flags and filtering to remove duplicate profiles within 15 minutes at same spatial locations. The profiles vary in vertical sampling density due to early Argos telemetry limitations vs recent Iridium floats. Profiles with insufficient depth coverage are excluded. The final analysis divides depth into two vertical sections (15-975 dbar and 975-1850 dbar) to maximize profile inclusion.

3. Framework Overview: Vertical integration of temperature profiles over the two vertical sections is performed with shape-preserving PCHIP interpolation and trapezoidal integration over depth. A local polynomial regression estimates the climatological mean field μ(x,y,t) including seasonal harmonics and quadratic time trends within moving spatial windows.

4. Local Gaussian Process Anomaly Mapping: The residual anomaly a(x,y,t) = OHC - μ is modeled as a zero-mean locally stationary Gaussian process with anisotropic exponential covariance including space (longitude, latitude) and time decorrelation scales. The parameters (variance, nugget variance, decorrelation lengths) are estimated by maximum likelihood within moving spatio-temporal windows. Temporal seasonal dependence of parameters is small, so median parameters over months are used for robustness.

5. Prediction: For each 1°×1° grid cell and month, the local GP conditional mean is computed using data inside the window with plug-in covariance parameters. The final OHC map at that grid cell and time is the sum of the mean field and anomaly prediction.

6. Uncertainty Quantification: To capture spatial and temporal dependence in prediction uncertainty, the authors develop a local conditional simulation algorithm. Gaussian white noise is simulated over the full spatio-temporal grid. For each grid point, local predictive covariance calculated from the local GP is decomposed (eigendecomposition), square root applied to the local white noise subset, and the simulated anomaly at the grid point is obtained by convolution with the kernel defined from the covariance. This is repeated for all grid points to produce correlated conditional realization ensembles.

7. Ensemble outcomes capture full predictive covariance structure implicitly, enabling uncertainty estimates for the spatially integrated OHC and any linear or nonlinear functionals by transforming the ensemble members.

8. Cross-validation: Paired cross-validation compares observed anomaly pairs with simulated realizations to validate uncertainty quantification and shows good agreement in predictive variances and persistence statistics.

9. Software Implementation: The entire pipeline from raw Argo data download, quality control, vertical integration, local polynomial regression, local GP fitting, prediction, and conditional simulation ensemble creation is reproducible with modular open-source code released on GitHub.

Example end-to-end: For a given month in the North Atlantic, Argo profiles are vertically integrated using PCHIP, mean field estimated from local polynomial regression over a spatial window, residual anomalies computed and modeled as local GPs with window-specific covariance parameters estimated by MLE, Gaussian conditional mean predictions computed at each 1° grid cell, then local conditional simulations generate an ensemble of anomaly maps reflecting uncertainty and spatial-temporal covariance, which are combined with mean estimates to yield UQ-enabled OHC anomaly maps.

Technical innovations

• Extension of locally stationary Gaussian process regression to modeling vertically integrated Argo temperature profiles with space-time anisotropic covariance.
• Development of a novel local conditional simulation algorithm that efficiently generates ensembles capturing spatially and temporally correlated uncertainty from local GP models.
• Introduction of a paired cross-validation methodology to jointly validate prediction accuracy and the quality of uncertainty quantification for OHC fields.
• Integration of a climatological mean field model with seasonal harmonics and quadratic time trends within a moving local polynomial regression framework.

Datasets

• Argo GDAC profiles (2004–2022) — 1,731,330 profiles quality-controlled — publicly available at Argo global data assembly centres

Baselines vs proposed

• Model without climatological time trend: OHC trend estimate roughly half that obtained when including time trend.
• Kuusela and Stein (2018) model without time in covariance: higher cross-validated median absolute prediction error compared to proposed local space-time anisotropic model.
• Predictive uncertainties validated by paired cross-validation agreed well with empirical prediction errors, indicating reliable uncertainty quantification.

Figures from the paper

Figures are reproduced from the source paper for academic discussion. Original copyright: the paper authors. See arXiv:2606.31957.

Fig 1: Illustration of the local conditional simulation algorithm for a region in the North Atlantic (left panel)

Fig 2: 4-step software pipeline for mapping Argo data.

Fig 3: Ocean heat content time series and trends for different modeling choices.

Fig 4: Cross-validated monthly median prediction errors (midocean) for different modeling choices.

Fig 5: OHC anomalies for different modeling choices. The insets show the estimated lag-1 autocorrelations.

Fig 6: Cross-validated monthly median absolute prediction errors (upper ocean) for different modeling choices.

Fig 7 (page 21).

Fig 8 (page 21).

Limitations

• Vertical covariance between depth sections (15–975 dbar and 975–1850 dbar) is not modeled; total uncertainty uses a conservative upper bound.
• Conditional simulation ensembles rely on Gaussian assumptions; non-Gaussian features or nonlinearities are not explicitly addressed.
• Seasonal variation of covariance parameters is ignored by using median values for robustness, possibly missing finer temporal variability.
• Early Argos profiles are sparser and more unevenly sampled in depth, potentially reducing integration accuracy compared to more recent data.
• Validation limited to statistical cross-validation within Argo-sampled time period; no explicit out-of-sample or adversarial robustness testing.

Open questions / follow-ons

• How to incorporate vertical covariance between ocean depth layers to reduce uncertainty conservatism in full-depth OHC estimates.
• Extension to non-Gaussian modeling frameworks or incorporation of nonlinear oceanographic dynamics in anomaly modeling.
• Applicability and adaptation of the local conditional simulation approach to other oceanographic variables with different spatial-temporal correlation structures.
• Robust prediction and uncertainty quantification under changing Argo array configurations or future data gaps.

Why it matters for bot defense

This work is focused on oceanographic data mapping and spatio-temporal uncertainty quantification, not on security or bot defense mechanisms such as CAPTCHA or bot detection. However, the advanced statistical modeling and conditional simulation techniques presented demonstrate rigorous approaches to estimating spatially and temporally correlated fields with uncertainty from sparse, irregular observations. Bot-defense engineers might draw analogies in handling large-scale spatio-temporal behavioral data, incorporating local nonstationarity, and quantifying uncertainty in probabilistic models. The paired cross-validation method for validating both predictive accuracy and uncertainty calibration could inspire analogous validation techniques in bot detection systems. Overall, while not directly applicable, the detailed modeling insights and computational methods for scalable Gaussian process approximations have potential indirect relevance for large-scale, uncertainty-aware security analytics.

Cite

@article{arxiv2606_31957,
  title={ Locally stationary Argo ocean heat content estimates: Modeling, validation and uncertainty quantification },
  author={ Thea Sukianto and Mikael Kuusela and Donata Giglio and Anirban Mondal and Pulong Ma and Douglas W. Nychka },
  journal={arXiv preprint arXiv:2606.31957},
  year={ 2026 },
  url={https://arxiv.org/abs/2606.31957}
}

Read the full paper

• arXiv:2606.31957 (abstract)
• arXiv:2606.31957 (PDF)