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

Domain: earth_science
Published: 2024-01-16
Editorial outcome: accept

Abstract
Long-term changes in ocean heat content (OHC) represent a fundamental global warming indicator and are mostly caused by anthropogenic climate-altering gas emissions. OHC increases heavily threaten the marine environment, therefore, reconstructing OHC before the well-instrumented period (i.e., before the Argo floats deployment in the mid-2000s) is crucial to understanding the multi-decadal climate change in the ocean. Here, we shed light on ocean warming and its uncertainty for the 1961-2022 period through a large ensemble reanalysis system that spans the major sources of uncertainties. Results indicate a 62-year warming of 0.43 ± 0.08 Wm−2, and a statistically significant acceleration rate equal to 0.15 ± 0.04 Wm−2 dec−1, locally peaking at high latitudes. The 11.6% of the global ocean area reaches the maximum yearly OHC in 2022, almost doubling any previous year. At the regional scale, major OHC uncertainty is found in the Tropics; at the global scale, the uncertainty represents about 40% and 15% of the OHC variability, respectively before and after the mid-2000s. The uncertainty of regional trends is mostly affected by observation calibration (especially at high latitudes), and sea surface temperature data uncertainty (especially at low latitudes).

Introduction
Ocean warming is among the most notable threats to the marine environment, leading to e.g. sea-ice decline1,2, sea-level rise3, and an increase in frequency and amplitude of marine heatwaves4, which, together with other climate change effects (e.g., ocean acidification5), all endanger the marine ecosystem6. The ocean is characterized by non-uniform warming7, due, in turn, to the spatially and temporally varying heat uptake from the atmosphere and surface wind variability8, combined with multi-scale redistribution through advective processes and vertical mixing. Ocean heat content (OHC, that is, the total amount of heat stored in the ocean) is one of the most important indicators of climate change in the ocean. While inter-annual variations of OHC are shaped by both internal climate variability and external forcing9, its long-term increase is due to the human-induced increase of climate-altering gas concentrations10,11. Reliably quantifying ocean warming and its uncertainty has enormous importance for monitoring Earth’s energy budget and climate12,13.

Accounting for the main sources of uncertainty of OHC in the pre-Argo period is crucial to provide estimates with reliable uncertainty envelopes14. Before the Argo era (namely before the 2000s), uncertainty is dominated by the observations themselves, mostly the eXpendable BathyThermograph data15,16; however, horizontal mapping methods17,18 and vertical interpolation (i.e., the error length-scales19) are all non-negligible sources of errors20. Further to these, the accuracy of model-based reconstructions may also be affected by the vertical physics parametrization uncertainty and systematic model errors21, assumptions in the data assimilation systems, and the uncertainty in all other input datasets (notably, the atmospheric forcing and sea surface data22,23). Intrinsic (chaotic) ocean variability is also non-negligible for regional estimates24, making the overall picture quite complicated.

The analysis of ocean heat content is usually performed through (i) objective analyses17,25,26 which are statistical mapping of the observations, (ii) reanalyses that use an ocean general circulation model to project forward in time the ocean state, using observations to constrain the model trajectory through data assimilation27, (iii) a combination of them, for instance, using reanalyses to derive error statistics for unsampled regions28, (iv) use of proxy data29 or indirect data that infer OHC increase from satellite-derived steric sea level variations (i.e., the so-called geodetic estimates14), or (v) advanced data-driven techniques such as machine learning30. Reanalyses, while being very costly from a computational point of view, permit a multi-variate four-dimensional characterization of the ocean, useful in turn for process-oriented and cause–effect studies31; reanalyses are attractive because they can in principle ingest any measurement which is related to the ocean state (for instance, altimetry and gravimetry observations32,33), and for a much broader range of applications than objective analyses, notably to initialize long-range prediction systems34.

The present study aims to re-assess the OHC trends and their uncertainty as estimated by reanalyses, compare them with objective analyses, and draft a hierarchy of sources of uncertainty. It builds upon a large-ensemble reanalysis system (32 members) that includes, in the ensemble generation, the major sources of uncertainty. We also take advantage of an objective analysis system (i.e., a statistical mapping of the observations using the same analysis system as the reanalysis but without any numerical ocean model integration), which is used to test if any of the reanalysis signals are spuriously given by the interaction between the ocean model and the data assimilation35,36. The goal is thus to quantify the warming and assess and rank the sources of uncertainty in the OHC reconstruction.

Results
We summarize, in the next two sections, the main outcomes of the large ensemble reanalysis system (hereafter CIGAR, the Cnr Ismar Global historic Al Reanalysis, http://cigar.ismar.cnr.it) in terms of ocean warming distribution and uncertainties, respectively. For some uncertainty diagnostics, we divide the reanalysis period (1961–2022) into the observation-poor period 1961–2001 and the observation-rich period 2002–2022, considering that the Argo float deployment37 starts in the early 2000s and matures around the mid-2000s.

Ocean warming distribution
The total ocean warming from CIGAR in terms of energy per Earth’s area is equal to 0.43 ± 0.08 Wm−2 for the period 1961–2022 (Table 1), and it equals 0.41 ± 0.09 Wm−2 for the period 1961–2020 coinciding with the latest GCOS assessment38, which indicates 0.41 ± 0.10 Wm−2 for the same period. Note that we have followed the uncertainty definition of the GCOS assessment: assuming Gaussian distributions for the trend errors, the uncertainty is given as twice the ensemble standard deviation, which corresponds to the 95% confidence level38). Thus, the two estimates are identical in terms of trends, and very close in terms of uncertainty, confirming the two approaches are consistent and complementary, as expected. We found a significant global warming acceleration equal to 0.15 ± 0.04 Wm−2 dec−1 for the 1961–2022 period. During the recent period 2006–2018, the acceleration is equal to 0.20 ± 0.07 Wm−2 dec−1, in agreement within the error bars with the estimates 0.50 ± 0.47 Wm−2 dec−1 over mid-2005 to mid-201939 and 0.25 Wm−2 dec−1 over 2002–201913, although CIGAR shows much smaller uncertainty.

Timeseries for the global ocean heat content anomalies are shown in Fig. 1a, together with the GCOS20 assessment26. Additionally, we display estimates from ZANNA1940, which reconstructs OHC from sea surface data using Green’s functions, and from ARANN41, which uses autoregressive artificial neural networks. The figure shows that, although the full-period trend is similar, our reanalysis system shows enhanced interannual variability, and a pronounced non-linear OHC increase, with the first period (up to about 1997–2000) exhibiting a slower increase and a sharper increase afterward. The other timeseries show a steadier increase (see also Table 1 for quantitative diagnostics), similar also to CMIP simulations42. For the well-observed period 2007–2022, we also compare the reanalysis with other independent estimates, i.e., geodetic estimates43 and CERES-based timeseries28,44. This comparison, in terms of OHC yearly tendencies, shows a very good agreement between the datasets, except for some well-known differences in the CERES-based timeseries (e.g., all datasets but those based on CERES exhibit cooling in 201645). The high-frequency variability between CERES-based estimates and those from oceanic observations is known to be largely different, due to several weather and climate processes providing a different response in the TOA EEI (top-of-atmosphere Earth Energy Imbalance) compared to the ocean heat uptake12. The climate community is converging toward comparing these two complementary datasets at frequencies slower than 3 years28,43,46. However, an accurate understanding of the scale coherence between the two is still an open question, whose answer is complicated by the relatively short temporal record of the satellite data.

To gain confidence in the datasets, it is important to understand whether the enhanced interannual variations in the reanalysis system are somehow spurious; alternatively, the objective analysis methodology and the ensemble average operation performed over the quite diverse range of products of the GCOS assessment may flatten the ocean heat content tendency signals. To this end, we compare our ensemble mean realization with that of a corresponding objective analysis (same assimilation system but no ocean dynamical model, see “Methods” section for details), denoted OA in Fig. 1 and Table 1. This objective analysis has a trend, acceleration, and interannual variability very close to that of CIGAR. In one alternative objective analysis experiment, the background—namely, the prior estimate used in the statistical analysis—is taken from an external model simulation like CMIP (OA-BGSIM), mimicking to some extent OA methodologies relying on external model-based products as background47. This experiment shows a much more attenuated variability and a steady increase (see also Table 1), evident in the poorly sampled period. When the assimilation frequency and time window are extended from ten days to one month (experiment OA-MON), variability and OHC increase are damped during the most recent, well-observed, period; an attenuation of variability during the full reanalysis period occurs when horizontal correlation length scales in the data assimilation system are halved (experiment OA-SLS). In all these experiments, interannual variability, trend, and acceleration approach those of the GCOS20 assessment. This indicates that the combined choice of the background field, correlation length-scales, and temporal frequency in the objective analysis reconstructions affect the OHC interannual variability. Maps of interannual variability (Fig. 2a, d) reveal that areas of intense mesoscale activity (e.g., the western boundary currents, and the Antarctic Circumpolar Current) are those characterized by the largest interannual fluctuations. Here, objective analyses may struggle to reproduce such fluctuations due to their coarse temporal resolution and the lack of atmospheric forcing information. Likewise, observation-blind model simulations (as those used as background in OA-BGSIM), significantly under-estimate the inter-annual variability, especially in the Southern Ocean.

Timeseries for the three latitudinal bands (Northern and Southern Extra-Tropics, and Tropics) are shown in Fig. 1c. They indicate the largest interannual variability and trends in the extra-Tropics (particularly, in the Southern Ocean), while a steadier increase and attenuated interannual variability in the Tropics. In terms of trends, the Southern Extra-Tropics contribute more than other regions over the 1961–2022 period (0.63 Wm−2 against 0.35 and 0.37 Wm−2 for the Tropics and Northern Extra-Tropics, in units of the Earth’s surface), in agreement with previous studies48–50.

Maps of OHC’s significant trends and accelerations (1961–2022) are shown in Fig. 2, together with their uncertainty. Areas of large heat accumulation (>1 Wm−2) are visible in the Southern Ocean and the Arctic, while the Tropical band shows moderate warming, larger in the Atlantic Ocean than elsewhere. Although the in-situ observational sampling is limited in polar areas and the OHC uncertainty is therein large, CIGAR indicates the high latitudes as ocean warming hotspots, in agreement with many observation- and model-based studies51–54. Mid-latitudes show smaller warming than elsewhere. This agrees with previously estimated 1968–2019 trend maps55. The uncertainty generally follows the areas of the largest trends (high latitudes and Tropics), with notable uncertainty in the North Atlantic Ocean as well. The map of ocean heat content acceleration indicates that part of the Antarctic region (Weddel Sea) and the North Atlantic Ocean (Gulf of Mexico, western boundary currents, Labrador, Greenland, and Mediterranean Seas) suffer from ocean heat content acceleration exceeding 0.4 Wm−2 dec−1. The acceleration uncertainty is, similarly to the trend, the largest at high latitudes and near the Equator, but significantly smaller than the signal.

The large ensemble system allows us to quantify and understand recent changes and associate statistical significance with them. We show in Fig. 3 the areas that exhibit statistically significant OHC increase in 2022 compared to 2021 (last year’s OHC increase), which is of great importance for climate monitoring. Large portions of the mid-latitude Southern Hemisphere, Atlantic and Pacific Oceans, and Tropics, exhibit a significant increase in OHC, with patterns in agreement with a previous study56 (i.e., western Tropical Pacific Ocean, the southern part of the Indian Ocean, north-eastern Atlantic Ocean, etc.). The global OHC increase is driven by the accumulation of heat in many regions, rather than localized accumulation. To understand how global warming is affecting the OHC increase distribution, panel b of Fig. 3 shows (for each grid point of the large ensemble reanalysis system) the year owing the largest OHC. Large portions of all the main basins exhibit the latest two years as the warmest years, with some exceptions, for instance, located in the Eastern Tropical Pacific upwelling region, central North Atlantic Ocean, and Kuroshio extension. In terms of area-averaged values (Fig. 3c), on 11.6% of the global ocean, the year 2022 is the warmest year over the 1961–2022 record, almost doubling any previous year. These statistics confirm a robust acceleration of warming that concerns large areas of the global oceans. Other notable years (exceeding 5%) are 2021, 2016, and 2015.

Analysis of uncertainties
The ocean heat content ensemble standard deviation (spread) from CIGAR is shown in Fig. 4a, b for the global ocean and three latitudinal bands, together with the uncertainty of the GCOS20 global ocean heat content assessment. The top panel shows the absolute uncertainty, while the middle panel shows the percent uncertainty, normalized by the OHC interannual variability (temporal standard deviation of the detrended timeseries). Tropics appear as the most uncertain latitudinal band, while the Northern Extra Tropics are the least (with uncertainty values always <30%), mostly because of the maturity of its observing network before the Argo float deployment. Peaks in the Tropics correspond to strong El Nino events (e.g., 1998) and seem not much affected by the TAO/TRITON mooring array deployment in the 1990s. This is consistent with previous analyses57, which show that El Niño events lead to an increase and broadening of ensemble anomalies and covariances in reanalyses compared to neutral years, due to the enhanced thermocline variability. The Southern Ocean has very large uncertainty in the 1960s, which rapidly resembles that of the global ocean (35-40%) from the 1970s onward. Overall, the global ocean shows an uncertainty of about 40%, which drops to about 15% during the last decade (2013–2022). Compared to GCOS20, our reanalysis system has often a larger uncertainty but is steady at around 40% before the 2000s, while GCOS20’s fluctuates between 10% and 50%. The larger uncertainty in CIGAR is likely due to the many more sources of uncertainty in reanalyses than objective analyses (e.g., atmospheric forcing, model physics, etc.), which dominate the ensemble dispersion in periods with poor observational sampling. During the last decade, the normalized uncertainties of CIGAR and GCOS converged towards similar values of about 15%.

Temporal variations of warming uncertainty (ocean heat content trend), at the global scale, are evaluated by computing (Fig. 4c) the running ensemble standard deviation of trends, over both a 15 and a 30-year time window. The panel shows the running ensemble mean of the trend for comparison (in red, with the right-side axis). The trend uncertainty is high (more than 0.10 Wm−2) at the beginning of the timeseries and stabilizes between 0.05 and 0.10 Wm−2 until approximately 1995; then, it increases to more than 0.15 Wm−2 (2001) for the 15-year time interval and decreases afterward, up to values of <0.05 Wm−2. The increase starts earlier than the Argo float deployment; however, it is difficult to say how much of this is affected by the observational sampling, which only partially follows the 15-year running trend variability. Additional objective analysis (OA) experiments randomly withholding 50%, 75%, and 90% of Argo floats, or with no observation assimilation below 1000 m of depth, led to statistically insignificant differences in ocean warming and acceleration compared to the OA experiment with the full observing network; this suggests that our results are not significantly influenced by the change in observational sampling following the deployment of Argo floats. Figure 4c also reports diagnostics with 30-year time intervals, which show much steadier behavior than those with 15-year intervals.

Clustering the ensemble members by grouping their type of perturbations, i.e., their sources of uncertainty, allows us to understand the main sources of uncertainties, as explored in regional operational contexts58. In practice, this is done by considering the averaged spread for members with the same source of uncertainty (see the “Methods” section), in comparison with the total (full ensemble) spread. We show in Fig. 5 the outcomes of this assessment in terms of both global and regional diagnostics (Fig. 5a, b), and maps of significant prevailing sources of uncertainty (Fig. 5c, d; see the figure captions for the exact definitions).

The uncertainty of the global trend (top left panel) is a percentage estimate of the total ensemble standard deviation (which accounts for model physics and data assimilation uncertainty). In this case, the sources of uncertainty are not additive and saturate between 45% and 75%, meaning that they can all explain a large portion of the resulting ensemble dispersion. For the 1961–2022 global trend, the main source of uncertainty is the air–sea flux formulation and the observation dataset, provided that, in the global OHC trend estimates, the air–sea flux formulation directly drives the heat tendency because horizontal and vertical redistributions of heat do not affect the global mean signal. The sea surface temperature (SST) uncertainty exceeds 70% of the global trend uncertainty, while other sources show smaller values, with the initial conditions having the smallest percentage of 45%. The early period 1961–2001 shows similar results, with the atmospheric forcing uncertainty becoming the main source of uncertainty.

Regional trend uncertainty is defined as the percent area where a certain source of uncertainty is significantly prevailing. Thus, percent values are additive in this case. For the entire period, it is the largest for the observations, which account for about 22% of the uncertainty, followed by the SST uncertainty (13%) and initial conditions uncertainty (11%). Atmospheric forcing and the air–sea flux formulation together do not exceed 10% and represent more marginal sources of uncertainty than the others. In the observation poor period 1961–2001, the SST uncertainty, which plays a major role in the absence of other developed networks, is the largest source of uncertainty accounting for about 17% of the total regional uncertainty.

Maps of prevailing types of uncertainty indicate which uncertainty source is locally dominant—statistically significantly—in the trend uncertainty. Only areas where a statistically significant source occurs are shown. The uncertainty associated with the observation preprocessing (OBS) is found to prevail in the Southern Ocean and in the North Atlantic, initial conditions uncertainty dominates in the Indian Ocean, and SST uncertainty mostly manifests at low latitudes. The atmospheric boundary conditions uncertainty confirms more marginal and ascribed to small portions of the global oceans.

For the poorly sampled period 1961–2001, patterns are more complicated to interpret, with SST uncertainty extending its areas of dominance to mid-latitudes; the atmospheric boundary conditions become the prevailing source of uncertainty in large portions of the Southern Ocean.

Methods
Ocean model configuration
The ocean model is NEMO60, version 4.0.7, which includes the five-category sea-ice dynamic and thermodynamic model SI361. It is implemented at about 1° of horizontal resolution, with augmented meridional resolution in the Tropics (up to 1/3°), and 75 vertical depth levels with partial steps62. The model uses bulk formulas adapted from the AEROBULK package63, a 3-band RGB scheme for the penetrating component of the air–sea heat flux, with extinction coefficients that depend on a monthly climatology of chlorophyll64; for vertical mixing, the TKE scheme65 is used with re-tuned background eddy viscosity and diffusivity66, while the horizontal diffusion operator is composed by a Laplacian operator for temperature and salinity and a bi-Laplacian for momentum.

The system is forced at the surface by the ECMWF ERA5 reanalysis67. Hourly fields of near-surface temperature, humidity, wind, and pressure are used to calculate high-frequency turbulent fluxes, while daily means of radiation and precipitation fields are used, with an analytical diurnal modulation68 of downwelling solar radiation. At the sea surface, heat and freshwater fluxes are corrected through a relaxation scheme that nudges SST and sea surface salinity (SSS) to external analyses. The relaxation timescale is set to 1 month and 1 year for SST and SSS, respectively. SSS analyses are taken from the UKMO EN469, while SST analyses depend on the specific ensemble member (see the section “Sources of uncertainty and ensemble generation”).

The freshwater discharge onto the ocean from land and ice sheets is provided by the JMA JRA-55-do reanalysis70, as inter-annually varying daily mean freshwater inputs, re-gridded onto the NEMO grid through conservative remapping.

Assimilation scheme
An oceanic 3DVAR scheme is used for assimilation, which implements background-error covariances estimated from long-term climatological anomaly differences, re-tuned via a posterior calibration method71 within preliminary assimilation experiments. In particular, background-error vertical covariances are implemented through the application of multi-variate spatially varying EOFs, while a first-order recursive filter that uses spatially varying correlation length scales is adopted horizontally72. The set of assimilated observations includes all in-situ profiles (MBT, XBT, and CTD casts, moorings, floats and gliders, and animal-borne sensors), extracted from the UKMO EN4 dataset69. Observational errors were initially taken from previous estimates73 and calibrated employing the posterior method71. A variational quality control scheme74 is used for non-linear quality control of the assimilated observations. A 10-day assimilation time window and analysis increment application frequency are adopted.

In addition to the variational assimilation scheme, a large-scale model bias correction (LSMBC) scheme is applied75, where deep ocean waters (below 500 m of depth) are weakly relaxed towards external objective analyses, at a temporal and spatial scale of 10 years and 1000 km, respectively. Additional experiments with the withholding of the scheme, or in simulations using the LSMBC but without assimilating in-situ profiles, or comparing CIGAR with the objective analyses indicate that such a large-scale bias correction scheme has no significant impact on the resulting warming of CIGAR. Indeed, this is driven by the assimilation of in-situ profiles, in terms of long-term changes and interannual variability. The fact that the OA objective analysis experiment (see its description below), with no LSMBC, has a very close variability and trends to CIGAR, further confirms that LSMBC is negligible for ocean warming assessments.

For comparison with the reanalysis system, an objective analysis scheme (OA, hereafter) is built upon the same data assimilation scheme. The two analysis systems share the same data assimilation configuration, with an analysis frequency of 10 days, but OA does not include any dynamical ocean model (namely, the model time integration step), and simplifies to an observation statistical mapping, which uses the same variational data assimilation scheme as CIGAR. As background-error covariances are estimated from climatology anomalies, they are representative of the climate modes of variability rather than the background accuracy; thus, the same set is used for both reanalyses and objective analyses. In OA, background fields are taken from a 10-day climatology to which the 10-day (persistent) anomaly from the previous analysis cycle is added, mimicking the cumulative effect of data assimilation in reanalyses. The OA system is used to interpret the different behavior of the ocean heat content reconstructed from the ensemble reanalysis compared to that from objective analyses.

Sources of uncertainty and ensemble generation
The large-ensemble reanalysis system CIGAR accounts for the major sources of uncertainty in the system configuration. These are listed below, and correspond to the ensemble generation strategy:
- Observation bias correction (for MBT, XBT, and floats) is an important source of uncertainty for multi-decadal retrospective analyses15. Here, it is spanned using two different sets of observations (namely, two different algorithms for XBT corrections) which include, respectively, two different XBT corrections76,77. Both corrections are provided through the UKMO EN4 profile dataset.
- The SST dataset used in the surface relaxation scheme; here we consider two SST analysis realizations, given by the UKMO HadISST78 and the JMA COBE79 analyses. The differences between the two dataset exceed 1°C in areas of enhanced mesoscale activity (e.g., around the Antarctic Circumpolar Current), and show, in terms of globally averaged values, that COBE-SST is slightly colder than HadISST by about 0.15°C. After 2005, the two datasets tend to converge.
- Initial conditions at the beginning of the reanalyzed period are a possibly important source of uncertainty, especially for ocean heat content analyses;31 here, we use two sets of initial conditions, extracted from lagged simulation restarts in previous assimilation experiments. Specifically, we have taken the 1948 and 1968 initial conditions from previous pilot assimilation experiments as initial conditions valid in 1958; we consider the first three years (1958–1960) as an adjustment period. The pilot assimilation experiment is the last iteration of many iterative 1958–2021 reanalyses performed as spinup to stabilize the ocean model drift, with initial conditions taken from the latest restart of the previous iteration.
- Atmospheric forcing uncertainty is accounted for using two different ensemble members of the ERA5 reanalysis. We added the ensemble anomalies from members 1 and 4 of the ERA5 ensemble data assimilation (EDA) system to the deterministic ERA5 reanalysis fields. The two members have been chosen after clustering the atmospheric surface kinetic energy following a clustering method80.
- Air–sea flux formulation, whose uncertainty is modeled using two distinct bulk formulations, i.e., the NCEP/CORE formulas81 and that from ECMWF82. The two bulk formulas implement, among others, a different formulation of the Charnock coefficient and a different use of the sea surface temperature (bulk versus skin), respectively63,83; moreover, within the bulk formula implementation, the NCEP uses the absolute wind while the ECMWF the relative wind (namely, the surface currents are subtracted to the wind).

The large ensemble results from any possible combination of these sources of uncertainty, resulting in a 32-member reanalysis. On top of this configuration, we use a stochastic parameter perturbation (SPP) scheme84 to stochastically account for the model uncertainty (in particular, the uncertainty of the solar extinction coefficients, the TKE parameters, the surface nudging relaxation time, and the horizontal diffusivity and viscosity). Additionally, we perturb the observations before their ingestion in the variational data assimilation system, through a Gaussian random deviation, with the standard deviation equal to the representativeness error of the observations85. This allows us to consider stochastic perturbations for both the assimilation and forecast steps, in addition to the sources of uncertainty explicitly detailed above.

We estimated the sampling uncertainty of the ensemble mean of the global OHC trend (1961–2022) using previous definitions86 and found that the 32-member ensemble size reduces by 88% the sampling error compared to a 10-member ensemble. The sampling uncertainty may further decrease with a larger ensemble than 32; however, the size of 32 is chosen as a compromise between small sampling uncertainty and computational affordability.

This study is designed to account for the major sources of uncertainty in the reanalysis system. However, there may be other processes, that may induce additional uncertainties and are not sampled by our system, for instance: exchanges of heat between the ocean and the sea ice; small-scale energy exchanges that cannot be resolved by our reanalysis system due to its limited spatial resolution; systematic errors in the atmospheric forcing not spanned by the ERA5 members; observational sampling error not spanned by the use of one observation production system (EN4).

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The authors have adequately addressed my concerns and comments and I recommend the manuscript for publicafion.\nI have idenfified a couple minor edits that the authors might want to look into.\n-line 8: In my opinion it is the change in OHC that represents a climate change indicator. I recommend to rephrase to: \"Long-term changes in OHC represent ... and are mostly ... .\"\n-line 12: \"understanding\" instead of \"understand\".\n-line 36: threaten?\n-line 41: amount of\n-line 46: \"monitoring Earth's energy budget and climate change\"", "reviewer_2": "I find most of my comments have been addressed in the revised version of the manuscript. But I have sfill 2 reservafions about this manuscript\n1) One of my comments has been addressed by the authors but it is not reflected in the revised version of the manuscript. It is my comment about the befter results of recent reanalysis in simulafing OHC changes compared to previous versions of the reanalysis. So far, Palmer et al. 2017 was the reference in terms of assessment of the OHC in reanalysis. Here the authors results are significantly befter than in Palmer et al. 2007. The authors explained in the rebuftal why it was so. That is fine but these explanafions do not appear in the revised manuscript. I think they should appear in a way or another ( at least in a supplementary material) because they are important for the reader to understand what make the results befter here than in Palmer et al. 2017. It will also give a touch on the actual progress of reanalysis\n2) In the revised version of the manuscript the authors explain why their reanalysis-based esfimate of the OHC shows more variability that the objecfive analysis one. They argue this is due to a limitafions in objecfive analysis. Then why EBAF is also showing less interannual variability than the reanalysis? And CERES+ as well!. The author should answer this quesfion or at least discuss it to convince the readership that the interannual signal shown in reanalysis is not spurious\nI think that after correcfing these 2 points the manuscript can be accepted for publicafion\nI add hereafter a few minor comments\nL45-46: “Ocean Heat Content 45 (OHC), represenfing the total amount heat stored in the ocean, is one of the most important indicators for ocean warming.” This phrase is a tautology: ocean warming is evidently indicated by OHC changes. I suggest rephrasing\nL82: “The present study aims to re-assess the OHC trends…” you probably mean the” OHC trends as esfimated by reanalysis”\nL112: Habuba->Hakuba\nL127: I find the agreement is good with the geodefic solufion but less good with the GCOS solufion and not good with CERES. Can you comment on that? There is in general much more interannual variability in the OHU fime series from reanalysis in situ data and geodefic esfimate than in CERES. Why? Do you have hints on this? Note that the interannual variability in CERES+ is small as well. Why?\nL132: What about the comparison with CERES?\nL178: W.m-1 dec-1 -> W.m-2dec-1", "reviewer_3": "Most of my previous comments have been resolved to the safisfacfion. However, further revision is needed. See my comment below.\nOne of my previous comment is about the uncertainty in OHC changes that is associated with ocean observing system sampling error (e.g. changes in spafial coverage of global ocean observing system, in parficular before and after Argo). See my previous comment below\n\"One of the major uncertainfies associated with both ocean reanalysis and objecfive analysis when considering OHC changes, is the sampling errors in the heterogeneous ocean observing network, which has evolved dramafically over the last few decades. E.g., the Southern Ocean is only properly sampled after Argo deployment since mid-2000. It is possible that the accelerafion of a warming trend in global OHC after the 2000s is a result of befter sampling (both horizontally and verfically) of the global ocean. This signal is visible in Fig.1, but was not fully discussed in this manuscript. \"\nDespite effort from the authors to address the high uncertainfies OHC signals in the Southern Ocean, the above point (of changing observing network) is not really addressed in the revised manuscript. I am sfill a bit suspicious about this increased trend after 2000, related to my comment above. This is encouraging to see the new sub-panel in Figure 1-a that includes comparison with CERES data. However, inter-annual variabilifies from various products (ENS-REA and CERES related datasets) are really large. If you include also ZANNA19/ARANN/GCOS20 into this panel, how different are these datasets looks like?\nIn Fig.1-c panel, this increase of warming trend since 2000 is most pronounced in the Southern ET heat content anomaly, this is to me suggesfing that the accelerafion of warming trend may be a spafial sampling issue. A related quesfion is, If Argo data is excluded in your DA system, do you have the same increase of posifive trend in OHC anomalies after 2000 in global and in the Southern ET?", "reviewer_4": "REVIEWERS' COMMENTS\nReviewer #3 (Remarks to the Author):\nThanks to the authors who takes fime to resolve my concerns about uncertainfies in OHC trend related with observafion sampling errors. And I welcome the effort that the authors have made to address this issue, especially new content dedicated to discuss this topic (between L208-223). I recommend this manuscript for publicafion.\nOne minor comment\nWe need to be careful when discuss OA experiments by withholding a fracfion of Argo data because there is no Argo data before 2000. This approach (of withholding Argo) can be to sample uncertainfies in OA reconstrucfion of OHC trend in the Argo period, but can not tell much about observafion sampling error in the Pre-Argo period."}, "author_response": "Response to the Reviewers' comments. 2nd round of reviews\nManuscript NCOMMS-23-08971 by Storto and Yang\nDear Reviewers,\nFirst, we would like to thank you for the careful reading and the updated comments and suggestions. Below, we provide a point-by-point answer to your comments (reviewers in bold, our answer in normal font). Both authors agree with the changes. We have also corrected another couple of grammar issues within the manuscript. Additionally, in the revised version of the manuscript, we now refer to the updated name of our dataset as CIGAR (CIGAR: the Cnr Ismar Global historicAl Reanalysis), as data are released and advertised through a public website, and we have decided to update the name instead of the more generic and less appealing ENS-REA as before. Data and additional info are now released (see also http://cigar.ismar.cnr.it ).\nReviewer #1 (Remarks to the Author):\nThe authors have adequately addressed my concerns and comments and I recommend the manuscript for publication. I have identified a couple minor edits that the authors might want to look into.\n-line 8: In my opinion it is the change in OHC that represents a climate change indicator. I recommend to rephrase to: \"Long-term changes in OHC represent ... and are mostly ... .\"\nCorrected as suggested by the reviewer\n-line 12: \"understanding\" instead of \"understand\".\nCorrected as suggested by the reviewer\n-line 36: threaten?\nCorrected as suggested by the reviewer\n-line 41: amount of\nCorrected as suggested by the reviewer\n-line 46: \"monitoring Earth's energy budget and climate change\"\nCorrected as suggested by the reviewer\nReviewer #2 (Remarks to the Author):\nI find most of my comments have been addressed in the revised version of the manuscript. But I have still 2 reservations about this manuscript\n1) One of my comments has been addressed by the authors but it is not reflected in the revised version of the manuscript. It is my comment about the better results of recent reanalysis in simulating OHC changes compared to previous versions of the reanalysis. So far, Palmer et al. 2017 was the reference in terms of assessment of the OHC in reanalysis. Here the authors results are significantly better than in Palmer et al. 2007. The authors explained in the rebuttal why it was so. That is fine but these explanations do not appear in the revised manuscript. I think they should appear in a way or another ( at least in a supplementary material) because they are important for the reader to understand what make the results better here than in Palmer et al. 2017. It will also give a touch on the actual progress of reanalysis\nThanks a lot for this suggestion, we have been indeed hesitating whether to report part of the discussion of the previous response to the reviewers in the manuscript or not. The reviewer’s suggestion pushed us to add this discussion, finally. We have now included a section on this in the Methods (the new section 4.4), which summarizes all the previous arguments.\n2) In the revised version of the manuscript the authors explain why their reanalysis-based estimate of the OHC shows more variability that the objective analysis one. They argue this is due to a limitations in objective analysis. Then why EBAF is also showing less interannual variability than the reanalysis? And CERES+ as well!. The author should answer this question or at least discuss it to convince the readership that the interannual signal shown in reanalysis is not spurious\nWe thank the reviewer for raising this point; we believe, as demonstrated in previous works, that CERES-derived estimates contain high-frequency variability intrinsically different from that in the ocean heat uptake (and found smaller than the ocean heat uptake in many studies). For instance: i) Marti et al., 2022 (https://essd.copernicus.org/articles/14/229/2022/ ) removes signals lower than 3 years in the comparison with OHC tendencies; ii) Storto et al., 2022 (https://doi.org/10.1175/JCLI-D-21-0726.1) obtained coherence between the two datasets for the 3-year time-scale through cross-wavelet analysis; iii) Meyssignac et al., 2023, (https://www.nature.com/articles/s41558-023-01735-z) compares the EEI from the geodetic approach (inferred by altimetry minus gravimetry) to that of CERES after low-pass filtering the geodetic dataset, as this latter has more pronounced variability than CERES. This is because the signals at the top of the atmosphere respond to many atmospheric processes at very different scales (from tropical convection to decadal variability) which may be significantly different than the signals of the ocean heat uptake (see also the discussion in e.g. von Schuckmann et al., 2016, https://doi.org/10.1038/nclimate2876). As the effect of all these multi-scale processes on the different EEI estimates is known but far from being assessed in detail by the climate community, we prefer to discuss this issue in the revised version without being completely conclusive; indeed, it is still an open question which scales are fully comparable between the ocean heat uptake measured through the oceanic observing network and the EEI measured at the top of the atmosphere.\nWe have shown and proved that the interannual variability from reanalyses is not spurious compared to objective analyses, but it is beyond the scope of the present work to show how TOA or ocean data manifest different frequency variability.\nWe have added a discussion on this point in section 2.1: The high-frequency variability between CERES-based estimates and those from oceanic observations is known to be largely different, due to several weather and climate processes providing a different response in the TOA EEI compared to the ocean heat uptake (von Schuckmann et al., 2016). The climate community is converging toward comparing these two complementary datasets at frequencies slower than 3 years (Marti et al., 2022; Storto et al., 2022; Meyssignac et al., 2023). However, an accurate understanding of the scale coherence between the two is still an open question, whose answer is complicated by the relatively short temporal record of the data.\nI think that after correcting these 2 points the manuscript can be accepted for publication. I add hereafter a few minor comments. L45-46: “Ocean Heat Content 45 (OHC), representing the total amount heat stored in the ocean, is one of the most important indicators for ocean warming.” This phrase is a tautology: ocean warming is evidently indicated by OHC changes. I suggest rephrasing\nCorrected, thanks: “OHC… most important indicators of climate change in the ocean”\nL82: “The present study aims to re-assess the OHC trends…” you probably mean the” OHC trends as estimated by reanalysis”\nCorrected as suggested by the reviewer\nL112: Habuba->Hakuba\nCorrected\nL127: I find the agreement is good with the geodetic solution but less good with the GCOS solution and not good with CERES. Can you comment on that? There is in general much more interannual variability in the OHU time series from reanalysis in situ data and geodetic estimate than in CERES. Why? Do you have hints on this? Note that the interannual variability in CERES+ is small as well. Why?\nL132: What about the comparison with CERES?\nPlease see answer to the second major point above.\nL178: W.m-1 dec-1 -> W.m-2dec-1\nCorrected\nReviewer #3 (Remarks to the Author):\nMost of my previous comments have been resolved to the satisfaction. However, further revision is needed. See my comment below. One of my previous comment is about the uncertainty in OHC changes that is associated with ocean observing system sampling error (e.g. changes in spatial coverage of global ocean observing system, in particular before and after Argo). See my previous comment below\n\"One of the major uncertainties associated with both ocean reanalysis and objective analysis when considering OHC changes, is the sampling errors in the heterogeneous ocean observing network, which has evolved dramatically over the last few decades. E.g., the Southern Ocean is only properly sampled after Argo deployment since mid-2000. It is possible that the acceleration of a warming trend in global OHC after the 2000s is a result of better sampling (both horizontally and vertically) of the global ocean. This signal is visible in Fig.1, but was not fully discussed in this manuscript. \" Despite effort from the authors to address the high uncertainties OHC signals in the Southern Ocean, the above point (of changing observing network) is not really addressed in the revised manuscript. I am still a bit suspicious about this increased trend after 2000, related to my comment above. This is encouraging to see the new sub-panel in Figure 1-a that includes comparison with CERES data. However, inter-annual variabilities from various products (ENS-REA and CERES related datasets) are really large. If you include also ZANNA19/ARANN/GCOS20 into this panel, how different are these datasets looks like? In Fig.1-c panel, this increase of warming trend since 2000 is most pronounced in the Southern ET heat content anomaly, this is to me suggesting that the acceleration of warming trend may be a spatial sampling issue. A related question is, If Argo data is excluded in your DA system, do you have the same increase of positive trend in OHC anomalies after 2000 in global and in the Southern ET?\nWe thank the reviewer for clarifying this point, and we split the reply into two parts, one regarding the interannual variability, and one about the observational sampling\ni) Interannual variability\nRegarding the interannual variability, Figure 1a already includes GCOS20, and GCOS variability is very well aligned with that of ENS-REA (now called CIGAR). The only datasets showing damped variability are those related to CERES, and we argue that this is because CERES data measure a quantity intrinsically different from the OHC assessed by reanalyses or objective analyses, as it refers to the TOA EEI. For a detailed discussion, please see the answer to Reviewer 2, point 2). We have added a paragraph in section 2.1 to explain in detail the reasoning, corroborated by many studies in the recent literature, which suggest comparing the two datasets at a 3-year timescale, as timescales shorter than that respond to a different range of processes between OHC and TOA.\nii) Observational sampling\nRegarding the possible impact of the observational sampling on the estimates (global, or southern ocean), this is certainly an important point to discuss. We want to stress that in the second revision, we have already added (lines 101-105) a quantitative comparison of the OHC acceleration with independent estimates from Hakuba et al. (2021) (using TOA measurements), which already indicates the consistency of the two estimates. This in turn means that observational sampling is not critical in our estimates, and this is also indirectly confirmed by the bottom panel of Figure 4 in the manuscript (Figure 4c).\nAnother way to look at this is by withholding a large fraction of Argo data from the OA experiments (doing this in reanalysis experiments is too computationally demanding, therefore we use OA as a benchmark, provided that we have already demonstrated that the two have comparable long-term warming/acceleration signals).\nThe figure and table on the next page summarize the results at the global scale for the full period 1961-2022, with regional variations being negligible as well. There is no impact (statistically) of withholding a variable number of Argo floats: OA-50%, OA-75%, and OA-90% refer to the experiments where 50%, 75%, and 90%, respectively, of Argo floats have been randomly withheld, while OA-1000m when Argo float measurements below 1000m of depth have been rejected, to mimic the case where only upper ocean measurements are assimilated. Differences in terms of metrics of warming, acceleration, and interannual variability are very close to each other and indicate that the observational sampling has a negligible impact on the resulting warming and acceleration. We do not consider the case of 100% Argo rejection, which is not realistic considering that XBT and CTD data sampling decreased from the mid-1990s (see e.g. https://www.aoml.noaa.gov/phod/goos/xbtscience/data.php; Abraham et al., 2013, https://doi.org/10.1002/rog.20022; Meyssignac et al., 2019, https://doi.org/10.3389/fmars.2019.00432; Goes et al., 2020, https://doi.org/10.1175/JTECH-D-20-0027.1).\nWe have now mentioned this result while commenting on Figure 4c (lines 217-220 of the revised manuscript), but since the results of these new tests do not differ from OA, we think there is no need to include explicitly in the figures of the manuscript these new results.\nFigure R1. Global OHC anomaly for the experiments presented in the response above.\nExperiment Trend Acceleration Interannual Variability\nOA 0.42 0.12 0.18\nOA-50% 0.41 0.12 0.17\nOA-75% 0.40 0.11 0.17\nOA-90% 0.41 0.12 0.18\nOA-1000m 0.39 0.10 0.16\nTable R1. Reporting trend, acceleration, and interannual variability (as defined in the manuscript) for the additional OA experiments presented above. Note, values are like in Table 1 of the manuscript, except the period here is the full 1961-2022 period (therefore numbers may vary slightly)."}]

Decision letter