Acceleration of the ocean warming from 1961 to 2022 unveiled by large-ensemble reanalyses

Ocean warming is a major threat to the marine environment, contributing to sea-ice decline, sea-level rise, and more frequent and intense marine heatwaves. Ocean heat content (OHC), the total heat stored in the ocean, is a key indicator of climate change. While interannual OHC variations reflect internal variability and external forcing, its long-term increase is primarily driven by human-induced greenhouse gas concentrations. Accurately quantifying OHC change and its uncertainty is essential for monitoring Earth’s energy budget. Before the Argo era, OHC uncertainty is dominated by observational issues (e.g., XBT biases) and by reconstruction choices such as horizontal mapping and vertical interpolation; additional uncertainties arise from model physics, data assimilation assumptions, atmospheric forcing, sea-surface datasets, and intrinsic ocean variability. Analyses of OHC have used objective analyses (statistical mapping), ocean reanalyses with data assimilation, hybrid approaches, indirect geodetic estimates from steric sea level, and data-driven techniques (e.g., machine learning). This study aims to reassess OHC trends and uncertainty using a large-ensemble reanalysis (CIGAR), compare with objective analyses, and rank sources of uncertainty over 1961–2022, with particular attention to the pre-Argo period.

The study situates itself within multiple streams of prior work: (i) GCOS assessments of global OHC providing community benchmarks; (ii) objective analyses that statistically map historical in situ data; (iii) ocean reanalyses that assimilate multi-platform observations into dynamical models; (iv) geodetic/indirect estimates inferring OHC from satellite altimetry and gravimetry via steric sea level; and (v) advanced data-driven methods (e.g., ARANN neural networks, Green’s function reconstructions such as Zanna et al. 2019). Prior comparisons indicate varying interannual variability and trend estimates among products. Satellite TOA energy imbalance products (e.g., CERES EBAF) and optimized CERES series have been used for consistency checks with ocean-based OHC tendencies, though differences exist at high frequencies. Literature also highlights non-uniform ocean warming with high-latitude hotspots and significant Southern Ocean contributions, as well as lingering observational biases (e.g., XBT) and methodological uncertainties affecting pre-Argo reconstructions.

Reanalysis system: The large-ensemble reanalysis CIGAR (CNR ISMAR Global historical Reanalysis) comprises 32 members designed to span major uncertainty sources. The ocean model is NEMO v4.0.7 with a five-category sea-ice model, ~1° horizontal resolution with enhanced meridional resolution in the Tropics (to ~1/3°), and 75 vertical levels with partial steps. Vertical mixing uses a TKE scheme; lateral diffusion is Laplacian for tracers and bi-Laplacian for momentum. Air–sea fluxes use bulk formulations from either NCEP/CORE or ECMWF; a 3-band RGB scheme represents penetrating solar radiation with climatological chlorophyll-dependent extinction coefficients. The system is forced by ECMWF ERA5: hourly near-surface meteorology for turbulent fluxes, daily radiation and precipitation (with diurnal modulation). Surface relaxation nudges SST (1-month timescale) and SSS (1-year) to external analyses (EN4 for SSS; SST dataset depends on ensemble member). Freshwater inputs from JRA55-do are applied. Data assimilation: A 3DVAR scheme assimilates in situ profiles (MBT, XBT, CTD, moorings, floats, gliders, animal-borne sensors) from EN4 with calibrated observation errors and variational quality control. Background-error covariances employ multivariate, spatially varying EOFs vertically and a first-order recursive filter with spatially varying horizontal length scales. The assimilation window and update frequency are 10 days. A large-scale model bias correction (LSMBC) weakly relaxes deep waters (>500 m) toward external objective analyses over 10-year/1000-km scales; tests indicate negligible impact on long-term warming metrics. Objective analyses for comparison: An OA system uses the same assimilation configuration and 10-day cycle but no model integration, mapping observations onto a background built from a 10-day climatology plus the persistent anomaly from the previous analysis. Sensitivity OA experiments include OA-BGSIM (background from free model simulations), OA-MON (1-month window/frequency), and OA-SLS (halved horizontal correlation length scales). Ensemble generation and uncertainty sampling: Five main uncertainty axes are combined to form 32 members: (1) Observation bias correction using two XBT-corrected EN4 realizations (different XBT correction schemes); (2) SST analysis for surface relaxation: HadISST versus JMA COBE (COBE slightly colder globally by ~0.15 °C; convergence after 2005); (3) Initial conditions (two sets from lagged restarts, nominally 1948 and 1968 origins, with 1958–1960 as adjustment); (4) Atmospheric forcing variants by adding ERA5 EDA ensemble anomalies (members 1 and 4) to deterministic ERA5; (5) Air–sea flux bulk formulations: NCEP/CORE vs ECMWF (differing Charnock, SST usage—bulk vs skin—and absolute vs relative wind). Additional stochastic parameter perturbations (solar extinction coefficients, TKE parameters, SST nudging timescale, horizontal diffusivity/viscosity) and stochastic observation perturbations (Gaussian with representativeness error) are applied to represent model and assimilation uncertainties. Ensemble sampling error analysis indicates a 32-member size reduces sampling error by 88% relative to a 10-member ensemble. Diagnostics: Global and regional OHC time series, trends (linear slope), accelerations (quadratic term), and interannual variability (standard deviation of detrended annual means) are computed over 1961–2022 and subperiods. Uncertainty is quantified as twice the ensemble standard deviation (95% confidence). Running trend uncertainties use 15- and 30-year windows. Regional uncertainty attribution identifies statistically prevailing sources (OBS, SST, INI, atmospheric forcing ATF, bulk flux BLK) by clustering members by perturbation type and comparing spreads. Additional OA sensitivity tests with random Argo withholding and no assimilation below 1000 m assess observing system impacts.

- Global OHC increase (units per Earth’s surface area): 0.43 ± 0.08 W m−2 over 1961–2022. For 1961–2020: 0.41 ± 0.09 W m−2, virtually identical to GCOS (0.41 ± 0.10 W m−2). - Statistically significant acceleration: 0.15 ± 0.04 W m−2 dec−1 (1961–2022). During 2006–2018: 0.20 ± 0.07 W m−2 dec−1, consistent within uncertainties with other estimates (e.g., ~0.50 ± 0.47 W m−2 dec−1 for mid-2005–mid-2019; ~0.25 W m−2 dec−1 for 2002–2019). - Interannual variability: CIGAR shows enhanced variability and a pronounced nonlinear increase (slower pre-1997/2000, sharper thereafter), relative to GCOS, ZANNA19, and ARANN. OA configured like CIGAR reproduces similar variability, trend, and acceleration; OA-BGSIM (external model background), OA-MON (monthly window), and OA-SLS (shorter length scales) damp variability and yield steadier increases, approaching GCOS characteristics. - Regional contributions (1961–2022): Southern Extra-Tropics contribute most to warming (0.63 W m−2), versus Tropics 0.35 W m−2 and Northern Extra-Tropics 0.37 W m−2. High latitudes (Southern Ocean, Arctic) are warming hotspots; mid-latitudes show weaker trends. Acceleration hotspots (>0.4 W m−2 dec−1) occur in parts of the Weddell Sea and North Atlantic (Gulf of Mexico, western boundary currents, Labrador, Greenland, Mediterranean Seas). - Recent record-setting years: 2022 exhibits significant OHC increases over vast regions relative to 2021; 11.6% of global ocean area records its maximum annual OHC in 2022 (nearly double any previous year). Other notable years (>5% area) include 2021, 2016, and 2015. - Uncertainty levels and distribution: Global normalized OHC uncertainty ~40% before the 2000s, decreasing to ~15% in 2013–2022. Tropics are most uncertain; Northern Extra-Tropics least (<30%). Running trend uncertainty (global) exceeds 0.10 W m−2 early in the record, stabilizes 0.05–0.10 W m−2 until ~1995, peaks (>0.15 W m−2 around 2001 for 15-year windows), then declines to <0.05 W m−2 in recent years. OA tests withholding large fractions of Argo or excluding data below 1000 m yield no significant changes in warming/acceleration. - Sources of uncertainty: For the global trend, non-additive contributions from all sources saturate at 45–75% of the total spread; SST uncertainty explains >70%, with air–sea flux formulation and observation dataset also important; initial conditions contribute least (~45%). Regionally (additive metric), observation preprocessing dominates (~22% of ocean area), followed by SST (~13%) and initial conditions (~11%); atmospheric forcing and bulk flux formulation together contribute <10%. In 1961–2001, SST uncertainty is the leading regional source (~17%). Spatially, OBS dominates in the Southern Ocean and North Atlantic, INI in the Indian Ocean, SST at low latitudes; ATF and BLK are marginal, though ATF gains importance over parts of the Southern Ocean in the early period. - Consistency with independent estimates: For 2007–2022 OHC tendencies, strong agreement is found among CIGAR, geodetic estimates, and optimized CERES-based series, acknowledging known high-frequency discrepancies with CERES EBAF.

The study addresses the central question of how much and how rapidly the ocean has warmed since 1961, and with what confidence, by constructing a 32-member ensemble reanalysis that explicitly spans key uncertainties. The derived global warming and acceleration closely match community assessments (GCOS), supporting the robustness of reanalysis-based estimates. The enhanced interannual variability seen in CIGAR is shown not to be an artifact of model–assimilation interactions; instead, objective analysis design choices (backgrounds, error length-scales, assimilation frequency) can suppress variability and push reconstructions toward steadier increases, explaining differences with some objective analyses and multi-product averages. Regionally, results clarify that high latitudes and the Southern Ocean are dominant contributors to multi-decadal warming and acceleration, while the Tropics display steadier increases but higher relative uncertainties. The finding that over 11% of the global ocean reached record OHC in 2022 provides compelling evidence of ongoing acceleration and spatially widespread heat accumulation. Uncertainty attribution reveals that while multiple factors contribute to global trend uncertainty, regional uncertainties are most sensitive to observation preprocessing (bias correction, particularly at high latitudes) and SST analysis accuracy (especially at low latitudes). These insights directly inform the design of future reanalyses and guide where improvements in observations and surface datasets can most effectively reduce uncertainty.

This work delivers a comprehensive, uncertainty-aware reassessment of ocean warming from 1961 to 2022 using a 32-member large-ensemble reanalysis (CIGAR). Key contributions include: (i) precise estimates of global OHC increase (0.43 ± 0.08 W m−2) and significant acceleration (0.15 ± 0.04 W m−2 dec−1), consistent with independent assessments; (ii) identification of high-latitude hotspots and stronger Southern Ocean contributions; (iii) evidence that 2022 was exceptionally warm over extensive ocean areas; and (iv) a quantitative ranking of uncertainty sources showing regional dominance of observation preprocessing and SST analysis uncertainty, with all sources contributing to global trend uncertainty. The multi-perturbation ensemble approach proves reliable and is well-suited for extending reconstructions further back in time and for initializing long-range predictions. Future research should (a) expand ensemble sizes and resolution, (b) incorporate additional uncertainty sources (e.g., sea-ice–ocean heat exchanges, alternative observation products), (c) refine surface flux and SST analyses, and (d) further integrate satellite-derived energy budget constraints at coherent temporal scales.

- Some processes and uncertainties are not fully sampled: ocean–sea-ice heat exchanges; unresolved small-scale energy pathways due to model resolution limits; systematic errors in atmospheric forcing not represented by the selected ERA5 ensemble members; observational sampling errors beyond those captured by a single observation production system (EN4). - Reliance on historical bias corrections (e.g., XBT) and SST analyses introduces structural uncertainties that vary regionally and temporally, especially pre-Argo. - Although 32 members substantially reduce sampling error compared to smaller ensembles, larger ensembles could further constrain uncertainty. - Observational sparsity in high-latitude regions and early decades leads to higher uncertainties and potentially limits the robustness of regional trend attribution there. - Objective analysis comparisons show sensitivity to methodological choices (backgrounds, correlation length-scales, assimilation windows), which may influence apparent interannual variability in non-reanalysis products.