Recent acceleration in global ocean heat accumulation by mode and intermediate waters

Domain: earth_science
Published: 2023-10-28
Editorial outcome: accept

Abstract
The ocean absorbs >90% of anthropogenic heat in the Earth system, moderating global atmospheric warming. However, it remains unclear how this heat uptake is distributed by basin and across water masses. Here we analyze historical and recent observations to show that ocean heat uptake has accelerated dramatically since the 1990s, nearly doubling during 2010–2020 relative to 1990–2000. Of the total ocean heat uptake over the Argo era 2005–2020, about 89% can be found in global mode and intermediate water layers, spanning both hemispheres and both subtropical and subpolar mode waters. Due to anthropogenic warming, there are significant changes in the volume of these water-mass layers as they warm and freshen. After factoring out volumetric changes, the combined warming of these layers accounts for ~76% of global ocean warming. We further decompose these water-mass layers into regional water masses over the subtropical Pacific and Atlantic Oceans and in the Southern Ocean. This shows that regional mode and intermediate waters are responsible for a disproportionate fraction of total heat uptake compared to their volume, with important implications for understanding ongoing ocean warming, sea-level rise, and climate impacts.

Introduction
The ocean directly impacts the Earth’s climate by absorbing and redistributing large amounts of heat, freshwater, and carbon, and by exchanging these properties with the atmosphere1. About 91% of the excess heat trapped by greenhouse gases1 and 31% of human emissions of carbon dioxide2 are stored in the ocean, shielding humans from even more rapid changes in climate. However, warmer oceans result in sea-level rise, ice-shelf melt, intensified storms, tropical cyclones, and marine heatwaves, as well as more severe marine species and ecosystem damage1,3,4. These effects depend on the pattern of ocean warming; it is thus critical to quantify the dynamics and distribution of ocean warming to better understand its consequences and predict its implications5.

The observed distribution of ocean warming is not uniform. About 90% of total ocean warming is found in the upper 2000 m, with over two-thirds concentrated in the upper 700 m since the 1950s5–8, and an increase of warming rates at both intermediate depths of 700–2000 m5,7,9–12, and in the deeper ocean below 2000 m5,7,10,11. The Southern Ocean south of 30°S has been estimated to account for 35–43% of global ocean warming from 1970 to 201713, and an even greater proportion in recent years12,13, while Northern Hemisphere ocean warming appears to be concentrated in the Atlantic Ocean6,14,15. Due to the accumulated excess heat in ocean basins, an acceleration of total ocean warming has become more evident from recent observational-based studies5,6,16,17. While much past work has focused on the distribution of ocean warming as a function of depth and basin15,18–20, relatively little analysis has been undertaken of the distribution as a function of water-mass layers and within specific water masses. This is the focus of the present study.

Ocean warming is often measured by tracking changes in ocean heat content (OHC). A useful coordinate to analyze OHC changes is a density- or temperature-based framework because this is a more natural coordinate to explore the thermodynamics of ocean warming21–23. Deepening of isopycnals is associated with ocean heat uptake, with past work suggesting the maximum deepening can be found in the mode water density range over the Southern Hemisphere extratropics and in the North Atlantic, implying a volumetric increase of mode waters and associated subduction and lateral spreading of heat from high-latitude well-ventilated regions14,18–20,24. A temperature-based framework further suggests that about half of the surface heat uptake during 1970–2014 is confined to about a quarter of the ocean’s surface area in the subpolar regions, which is in turn capable of exchanging heat with the coldest 90% of the global ocean volume25.

The above studies have provided insights and identified where we lack observational-based constraints in understanding OHC variability. However, most past studies of OHC changes are depth- or basin-integrated, obscuring the OHC evolution in mode and intermediate waters and the related processes at play. Mode and intermediate waters are distinguished by their properties of density, salinity and low stratification by wintertime ventilation, and their importance in storing and redistributing heat, carbon and oxygen26–28, yet their heat content changes associated with global ocean warming are yet to be quantified. Furthermore, relatively little attention has been paid to their temperature and volumetric evolution and the associated role of subduction, spreading, and mixing in accumulating oceanic heat, despite observations highlighting a significant upper ocean volumetric increase of mode waters18 and showing striking local warming signals at depths of 700–2000 m in the North Atlantic and Southern Oceans19,29.

In this study, we examine the global and regional variability of OHC in the upper 2000 m of the ocean using products of subsurface temperature and OHC that have been constructed from historical measurements and Argo observations. We explore heat accumulation within global mode and intermediate water layers and their role in increasing global OHC. We further define a set of regional mode and intermediate waters using specific geographic and density constraints within those two layers (“Methods”). Using a classification of water masses allows us to quantify observed changes in the heat content of individual mode and intermediate waters and to calculate their contribution to the global OHC increase. We separately consider OHC increases due to volume and temperature changes to better understand the processes driving ocean heat uptake and redistribution by these water masses.

Results
Accelerating global ocean warming
Observations of the interior ocean thermal structure relied on measurements from research vessels prior to the advent of expendable bathythermographs (XBTs), which have been widely used since the 1970s. During the 1990s, unprecedented coverage of the global ocean using high-quality hydrographic sections was achieved as part of the World Ocean Circulation Experiment30 (WOCE). Subsequently, the Argo float array began to dominate ocean observations since around 2005, measuring the upper 2000 m of the water column. These observing systems have evolved to provide a valuable means to improve OHC estimation, leading to major advancements in global climate studies6. The community has also made considerable progress in developing gridded products of subsurface temperature for the global ocean by merging all available observing systems8,10,31,32. Due to instrumental biases in bathythermograph measurements33 and limited data coverage prior to Argo (before 2005), as well as geographic and depth limitations of conventional Argo floats, various correction schemes and objective analysis methods have been applied to mitigate these biases and increase our confidence in historical OHC estimates.

To provide a more robust analysis of the global and regional OHC changes in this study, we present the ensemble average of OHC estimates derived from various ocean data products (Fig. 1; “Methods”). The uncertainty of this estimate is taken to be ±2 standard deviations (±2σ), to indicate the degree of agreement across the products.

Although some differences in the global OHC estimates are evident among the different datasets, the ensemble mean of OHC estimates for the upper 2000 m shows that ocean warming is not only continuing but accelerating5,16,17 (Fig. 1). The average linear trend of global OHC over the pre-Argo period 1955–2004 is 3.9 ± 0.7 × 1021 J yr−1 (Fig. 1a; Table 1). For the Argo era 2005–2020, we estimate this trend to have increased almost threefold to 10.2 ± 1.8 × 1021 J yr−1. Further comparison with previous estimates of OHC trends (Table 1) and the latest assessment report of the Intergovernmental Panel on Climate Change (IPCC) Working Group I7 also reveals significant increases in the ocean warming rate across the measurement record5,8–12,17,18,20. This increasing ocean warming is not just limited to the upper layer of the ocean but also occurs at intermediate depths (700–2000 m) according to a further breakdown of OHC trends by depth over the upper 2000 m5,7.

When evaluating the ocean heat uptake for each decade (“Methods”), analysis of the past three decades reveals that the ocean heat uptake during 2010–2020 has increased more than 25% relative to 2000–2010 and has nearly doubled relative to the 1990’s WOCE era, as seen in Fig. 1b, where we highlight the decadal ocean heat uptake since the 1960s. Note that there has been both increased ocean sampling and a shift of the observational network from a ship-based system to the Argo network since the initiation of the global Argo array (2001–2003)34. This may impact the estimated increase in global ocean warming over the past three decades (Fig. 1). However, the rate of global mean sea-level rise has also been increasing since 1993 based on an independent estimate from satellite altimeter data1,35, providing confidence in our results given that half of the global sea surface height increase is due to thermal expansion of the ocean since altimeter measurements began1. Significant ocean warming and accelerating OHC changes are also consistent with the increase in net radiative energy absorbed by Earth detected in satellite observations5,36, something that is likely to continue throughout the 21st century6,16 in the absence of substantial greenhouse gas emissions reductions.

The increased ocean warming is non-uniformly distributed across ocean basins. Overall, in each ocean basin, an increase in OHC is observed (values indicated in Fig. 2a, b), with stronger warming in the mid-latitude Atlantic Ocean and the Southern Ocean compared with other basins6,12. Total warming in the Southern Ocean is estimated to account for ~31% of the global upper 2000-m OHC increase from 1980–2000 to 2000–2010 (Fig. 2a), and almost half of the global OHC increase from 2000–2010 to 2010–2020 (values indicated in parentheses of Fig. 2b). Hence the Southern Ocean has seen the largest increase in heat storage over the past two decades, holding almost the same excess anthropogenic heat as the Atlantic, Pacific, and Indian Oceans north of 30°S combined (Fig. 2d). The most striking warming in the Southern Ocean is concentrated on the northern flank of the Antarctic Circumpolar Current, the location of deep mixed layers and subduction hotspots for Subantarctic Mode Water and Antarctic Intermediate Water, as well as the location of subtropical mode waters formation further equatorward (Fig. 3). The well-ventilated regions near western boundary current extensions in the North Atlantic and North Pacific also reveal large warming over the past two decades. These hotspots of ocean warming are likely linked to enhanced uptake, subduction, and lateral spreading of heat associated with mode and intermediate waters that warrant further investigation.

Even though all products show robust ocean warming rates averaged in the upper 2000 m over the past few decades, reaching record warming over 2010–20205,17 (Table 1), when decomposed into water-mass trends these products exhibit better agreement after ~2005 once the Argo network was widely deployed. Thus, for the remainder of this study, we focus on the Argo era 2005–2020, a period common to all data products analyzed. We also limit the water-mass decomposition to be between 65°S and 65°N, because of sparse observations in the polar regions even after the Argo array was established.

Heat accumulation in mode and intermediate water layers
The global OHC change of the upper 2000 m is next decomposed into the tropical water layer, the mode water layer, and the intermediate water layer to study the distribution of ocean warming by global water-mass layers (Fig. 4). The most remarkable OHC change over the Argo era is the increase within the mode water layer (Fig. 4c, f). The mode water layer occupies only 20% of the upper 2000 m of the ocean, yet plays a dominant role in total ocean heat content increase (Fig. 4a, Supplementary Table 1). The intermediate water layer occupies 42% of the upper 2000 m, but exhibits a decreasing heat content trend (Fig. 4d). These opposing signals are mainly due to the mode water layer accumulating heat content via increasing its volume, while the intermediate water layer loses heat content due to a volume decrease (Fig. 5c, f). After factoring out volumetric effects, both the global mode and intermediate water layers show a near-continuous and monotonic warming across the Argo era (Fig. 5, Supplementary Fig. 2), revealing a robust footprint of global warming penetrating into the ocean interior. Taken together, the total heat accumulated in the mode and intermediate water layers is 89% of the net global ocean warming during the Argo era, despite their total volume just occupying 62% of the upper 2000 m. The net warming of these two layers, separate from volumetric effects, accounts for ~76% of global ocean warming, while ~12% is due to their combined volume change (Fig. 5a, Supplementary Table 2).

Heat content within the tropical water layer reveals large variations associated with El Niño-Southern Oscillation (ENSO) variability (Figs. 4b, f and 5b). In particular, the heat content of the upper 100 m increases significantly during El Niño years (2009/10, 2014/15, 2015/16, and 2018/19), offset partially by opposing variability below 100 m, confirming previous studies12,37. This ENSO-driven variability in OHC is related to net air-sea heat exchanges and the redistribution of heat laterally across different latitudes and ocean regions, and vertically between surface and subsurface layers37,38. In this way, tropical warming exhibits strong interannual variability, mainly due to El Niño, whereas warming in global mode and intermediate water layers is less variable and characterized by a near-continuous accumulation of heat. Because of the tremendous heat accumulation within the mode and intermediate water layers, we next focus on the OHC increase in a few key water masses found within these two layers.

Regional intensification in ocean warming
To further explore how the upper 2000-m ocean warming is distributed by basin during the Argo era 2005–2020, we next compute heat content trends over the upper 2000 m and find the increased OHC to be concentrated in the tropical Pacific, subtropical oceans and the Southern Ocean (Fig. 6a). Note that trends within this relatively short period can indicate both the forcing signal due to anthropogenic warming as well as signals superimposed due to internal climate variability. For example, the OHC change in the tropical Pacific shows an El Niño-like spatial structure, with warming in the east and cooling in the west (Fig. 6a). This is due to the prevalence of El Niño events during the latter half of the analyzed time period, relative to the start of that period, which was characterized by three strong La Niña events (2007/08, 2010/11 and 2011/12). This is linked to a transition of the Interdecadal Pacific Oscillation (IPO) from its negative to positive phase in around 201439 (note that our analysis period does not include the recent triple La Niña spanning 2020–2023). Warming in the extratropical latitudes is mainly confined to the subtropical Atlantic and Pacific Oceans in both hemispheres and the Southern Ocean. From the zonally integrated OHC increase, we find that this extratropical ocean warming can be observed vertically from the surface down to at least 300–1500 m depth, covering the mode and intermediate water layers (Fig. 7).

An independent measure of OHC changes can be found in altimeter measurements of sea level if salinity adjustments remain small and the ocean mass component of sea-level change is known and accounted for. This is because the ocean expands when it is warmed, thus directly affecting sea level. We find that the locations of strong increase in OHC and thermal expansion (Fig. 6a, c) correspond closely to regions where the sea level has also increased40,41 (Fig. 6b). These locations include the eastern tropical Pacific, the subtropical oceans, and the Southern Ocean. The effects of changed ocean salinity are small (Fig. 6d), mainly limited to the North Atlantic where Greenland ice melt and adjustments in the Atlantic Meridional Overturning Circulation have had an impact42,43. The mass component of sea-level change, due to land ice melt (including glaciers, and the Greenland and Antarctic ice sheets) and land water storage change, has also significantly increased across the altimeter era, but displays a more uniform pattern of sea-level rise compared to the thermal expansion component40,41.

Increased heat uptake by mode and intermediate waters
Anthropogenic warming has stratified most regions of the world’s oceans44, thereby, reducing ocean ventilation and slowing vertical exchanges of heat and carbon. However, it has been estimated that the total ocean from the surface to the seafloor has been warming at an accelerating rate11 (Fig. 1), and there is an indication that mode waters, on the basin average, have increased in their volume and depth14,18–20,24,45. Here we show that the accelerated OHC increase is concentrated in the subtropical and Southern Ocean (Figs. 2 and 6), from the surface to 300–1500 m depth (Fig. 7), spanning the mode and intermediate layers (Figs. 4, 5 and 7). These regions are characterized by deep wintertime mixed layers and vigorous mode water formation46 (Fig. 3), which act to imprint a signature of air-sea interaction and climate conditions in the ocean interior. Therefore, a component of this accelerated OHC increase must be due to increased heat uptake by mode and intermediate waters via wintertime ventilation, in combination with deepening of the warmed water by mixing. In this way warming from the surface can propagate into the ocean interior and spread laterally and vertically via circulation and mixing with surrounding water masses. This motivates us to study the heat content changes within individual mode and intermediate waters to better understand their role in warming the ocean. We do this by first defining 13 regional water masses using density and geographic limits within the mode and intermediate layers and over the subtropical and subpolar oceans respectively (“Methods”; Figs. 3 and 7, Supplementary Table 3), and then estimating their OHC changes over the Argo era 2005–2020 (Fig. 8).

The total heat content of a water mass can increase because it warms, or because the total volume of that water mass increases, or both (“Methods”). By decomposing the OHC increase over the Argo era 2005–2020, we find that the actual warming (within a fixed volume) of all regionally defined mode and intermediate waters can account for 48% of the global ocean warming (Figs. 8e and 9g, Supplementary Tables 3 and 4). The volume of each of these water masses also changes markedly but non-uniformly across this period, with their combined volume change modestly offsetting (−4%) the net mode and intermediate water warming calculated over a fixed volume. The combined OHC increase of the sum of all these regionally defined mode and intermediate waters, due to both warming (fixed isopycnal depth and volume) and volume changes, amounts to a total contribution of 44% to the global OHC increase (using the values shown in Fig. 8e and Supplementary Table 4). The total volume of the regionally defined water masses within the subtropical and subpolar oceans occupies just 24% of the upper 2000 m averaged over 2005–2020. Therefore, this disproportionately large heat uptake per volume of mode and intermediate waters reveals their outsized role in absorbing and storing heat for the world ocean under global warming. We next separately consider the OHC increased due to both volume and temperature changes.

Warming and thickening of water masses
The most striking ocean warming and thickening of the mode and intermediate water layers is over the subtropical Atlantic and Pacific Oceans and in the Southern Ocean, dominating the global OHC increase (Figs. 5, 7, 8 and 9). The Subtropical Mode Water layer has shown the largest thickening of the upper 2000-m ocean over the Argo era (Fig. 5, Supplementary Fig. 2). This thickening and warming of the layer are concentrated mainly over the subtropical Atlantic and Pacific Oceans, with a weak increase over the subtropical South Indian Ocean (Figs. 8b and 9a, d), which likely corresponds to ventilation regions of regional Subtropical Mode Waters26,47 (Fig. 3). The Subpolar Mode Water layer and the intermediate water layer both show broad warming (volume effects removed) over the Southern Ocean (Fig. 9e, f). In the Southern Ocean the 0–1500-m depth range is occupied by Subantarctic Mode Water and Antarctic Intermediate Water26,48,49. Warming from the sea surface to 800 m has been measured along a single transect in the Southern Ocean across mode and intermediate water layers29; this can also be seen in our analysis of the full zonally integrated OHC trends (Fig. 7). In particular, our analysis reveals that Subantarctic Mode Water has warmed and deepened substantially within and north of its ventilation region in the Southern Ocean, while Antarctic Intermediate Water has warmed and thinned due to the deepening of lighter isopycnals (Figs. 8c, d and 9e, f, Supplementary Table 4). Total warming (volume effects removed) of Subantarctic Mode Water and Antarctic Intermediate Water accounts for ~36% of global ocean warming over the Argo era (Fig. 9g).

Among the mode and intermediate waters defined in this study, we find that OHC trends of mode waters are dominated by their volumetric changes. Instead warming (fixed volume) exhibits a larger signal in the heat content trend for intermediate waters, including both Antarctic and North Pacific Intermediate Water (Figs. 3 and 9, Supplementary Fig. 2 and Table 4). This is likely because while the warming of water masses propagates downward, the density-defined boundary of water masses shifts both vertically and laterally as a response to ocean warming (Supplementary Fig. 3). In particular, the subtropical oceans and the Southern Ocean reveal a significant deepening and poleward shift of mode water isopycnals. This indicates that the Subtropical/Subantarctic Mode Water has been expanding and shifting toward locations that once were classified as Subpolar Mode Water and/or intermediate water. Note that although the mode water has shifted, this does not necessarily indicate a net warming of the water mass. Instead, the intermediate water layer could potentially be lightened and warmed due to mixing with the warmed and deepened mode waters over the subtropical oceans and the Southern Ocean19 (Fig. 9). In reality, a combination of these two effects has likely played out. In addition to density-based constraints, alternative criteria like stratification or salinity within specific locations are often used to define mode and intermediate waters26. These criteria also allow the examination of changes in the density and temperature of water masses. Within this framework, the poleward and downward migration of isopycnals over the subtropical oceans and the Southern Ocean can serve as an indication of the warming and lightening of mode and intermediate waters (Supplementary Fig. 3).

There are several possible mechanisms driving the increased ocean heat uptake, particularly over the ventilation regions of mode and intermediate waters. Anthropogenic warming provides the excess heat to warm the ocean in regions where sea surface temperature is cooler than the overlying surface air temperature—for example, an acceleration and poleward shift of the westerly winds50 has kept surface waters relatively cool in the Southern Ocean, and thus favorable for ongoing ocean heat uptake51. The circulation-related changes in heat transport also play a profound role in heat accumulation over the subtropical Atlantic Ocean and the Southern Ocean15,52. In addition, western boundary currents accumulate and feed heat toward the formation regions of Subtropical Mode Waters. Deep mixed layers at the mode water formation sites pump heat into the ocean interior via subduction53–55 (Fig. 3). In the Southern Ocean, strengthening westerly winds and Ekman pumping steepen isopycnal surfaces and increase the Subantarctic Mode Water volume via subduction56,57, while also pushing cold surface Antarctic water northward58, enabling increased heat uptake and subduction into mode and intermediate layers51, as noted above. Enhanced mesoscale eddy activity over the subtropical oceans and the Southern Ocean59 might also contribute to increased subduction and vertical transport of heat for Subtropical and Subantarctic Mode Waters and Antarctic Intermediate Water54,60,61.

While the globally integrated OHC of the upper 2000 m and in the mode and intermediate water layers both show a robust and near-monotonic increase (Figs. 1 and 5), variations in the formation mechanisms, properties, and heat content of localized mode and intermediate waters can be influenced by interannual to decadal variability (e.g., ENSO and the Indian Ocean Dipole62 on interannual time-scales, and the IPO on decadal time-scales). However, our analyses show that the main footprint of internal climate variability occurs over the tropical Pacific Ocean, whereas the rest of the Pacific as well as the Indian, Atlantic, and Southern Oceans all tend to show a consistent and robust warming, in particular over the subduction and well-ventilated regions near western boundary current extensions and north of the ACC (compare Figs. 2a, b and 6a). The near-monotonic warming in mode and intermediate water layers across the Argo era (Figs. 4 and 5a), suggests a robust ocean heat gain forced by anthropogenic climate change. In contrast, OHC changes within the tropical water layer and in the tropical Pacific Ocean reveal variations likely associated with ENSO and/or IPO variability.

Methods
Observational products
A total of 5 state-of-the-art gridded ocean datasets and OHC products are used in this study to estimate global ocean warming rates and three of those products are used to estimate OHC changes within the global water-mass layers as well as in regionally defined mode and intermediate waters. These datasets can be grouped into two different types: (1) gridded ocean datasets based on Argo only63 and (2) four products merged from Argo and other ocean measurements8,10,31,32,64. The estimate of global and basin-wide OHC trends among these selected observational-based products was previously shown to exhibit good agreement in the upper 2000-m ocean over the Argo era45 (Fig. 1). Several factors can affect the estimate of ocean warming rate and should be factored in, particularly over the pre-Argo period (Fig. 1, Supplementary Fig. 1). Such factors include uncertainties due to data scarcity and data quality prior to the late 1960s10, mapping artifacts from horizontal and vertical interpolation of the sparse historical data65, the increased spatial sampling coverage over time (particularly in the tropical and Southern Hemisphere oceans), and the shift of the observational network from a ship-based system to the Argo network during 2001–200334.

We use monthly means of in-situ temperature and practical salinity from a gridded Scripps Institution of Oceanography (SIO) Argo product63 (hereafter referred to as “SIO RG-Argo”), which is provided on 56 standard pressure levels (0–2000 dbar) and a 1° regular grid. This gridded Argo product was constructed from Argo data only. We only use the SIO RG-Argo data between 65°S–65°N because of low observational data coverage in high-latitude regions, particularly ice-covered areas and in marginal seas. As such, the resultant global OHC trends from SIO RG-Argo in Fig. 1 will be underestimated if these excluded polar ocean regions are also warming5. For the group of gridded datasets merged from both Argo and other historical observations, we include in-situ temperature and salinity observations from the Institute of Atmospheric Physics (IAP) data10,64 (monthly averaged from 1955 to 2020 with 1° horizontal resolution; referred to as “IAP data”) and four ensemble members of EN4.2.2 data from the Met Office Hadley Center32 (monthly averaged from 1955 to 2020 with 1° horizontal resolution). The EN4.2.2 data used in this study includes all four ensemble members available from the latest version of the EN4 product, with four different sets of XBT and mechanical bathythermograph (MBT) bias corrections applied (see ref. 32 for more details). When presenting results using the EN4.2.2 data, we choose to show the ensemble average from those four EN4.2.2 ensembles (referred to as the “EN4.2.2 ensemble mean”), to mitigate the uncertainty arising from different approaches applied in correcting biases in XBT and MBT data (Supplementary Fig. 1). Due to the large inconsistency within the four EN4.2.2 ensembles over the 1990s, and due to a large shift of global OHC around 2000–2004 in the EN4.2.2 ensemble mean, we have only included the EN4 data over the Argo period 2005–2020 in this study. The Gibbs-Sea Water (GSW) Oceanographic Toolbox from the TEOS-10 software66,67 is applied to convert the SIO RG-Argo data and EN4.2.2 data to Conservative Temperature68 and Absolute Salinity69, as well as the IAP in-situ temperature to Conservative Temperature, and achieves static stability on time scales of months by applying vertical stabilization software70. The thermal expansion and haline contraction coefficients, the surface-referenced potential temperature, θ, and potential density, ρ, are then calculated for SIO RG-Argo, IAP data and EN4.2.2 data, using functions from the GSW Oceanographic Toolbox. The neutral density71, γn, is calculated using functions from PreTEOS-10 software.

The 5-year running mean of global OHC during 1955–2020 from ref. 8 and the annual average of OHC during 2005–2020, updated by the NOAA/NESDIS/NCEI Ocean Climate Laboratory using temperature and salinity fields from the World Ocean Database (both are referred to as “Levitus et al. (2012)”), as well as the annual average of OHC during 1955–2020 from ref. 10 (referred to as “Cheng et al. (2017)”) and ref. 31 (referred to as “Ishii et al. (2017)”) are applied to reconcile the OHC estimate from SIO RG-Argo and EN4.2.2 ensemble mean in Figs. 1 and 2. These three OHC fields are all provided in a single depth layer as they were vertically integrated over the upper 2000 m. As for the estimate of OHC changes within the mode and intermediate water layers and in the key water masses in Figs. 3–9 and Supplementary Tables 1–4, and in the expansion and haline contraction in Fig. 6, in which the temperature and salinity fields in the longitude, latitude, and vertical directions are needed, we show the ensemble average of OHC estimates from SIO RG-Argo, IAP data, and EN4.2.2 ensemble mean. The uncertainty range (±2σ) is derived using the ensemble standard deviation (σ) of the corresponding estimates to indicate the degree of agreement across the products.

The basin mask used to distinguish ocean basins of the Southern Ocean (south of 30°S), the Atlantic, Pacific, and Indian Oceans in Figs. 2 and 7 is obtained from ref. 8. We also consider changes in sea-level anomaly in Fig. 6, using the E.U. Copernicus Marine Service Information (CMEMS) gridded multi-mission altimeter sea-level anomaly data with 1/4 × 1/4° spatial resolution, covering the Argo era (2005–2020).

OHC computation and decomposition
This study focuses on the warming of the upper 2000 m of the ocean in part because this is the volume of the ocean in closest contact with the atmosphere to experience global warming. This 2000-m constraint is also due to the depth limit of the Argo profiling floats, with only sparse data coverage deeper than 2000 m. Waters deeper than 2000 m are generally not ventilated rapidly enough to have experienced the last 50 years of global warming, but note that both the abyssal overturning slowdown72,73 and oceanic internal waves74 can result in bottom ocean warming over these time scales, as seen in observations75.

From the observational-based ocean datasets, we compute the monthly and annual mean OHC (J m−2) across the measurement period 1955–2020 for each isopycnal layer as follows:

OHC = Cp ∫ z2(x,y,z,t) z1(x,y,z,t) ρ(x,y,z,t) θ(x,y,z,t) dz, (1)

where Cp = 3992 J kg−1 K−1 is the specific heat capacity of seawater, θ is the potential temperature with respect to 0°C, and z1 and z2 (m) are the depth of the lower and upper isopycnal surfaces for integration, respectively. The OHC is partitioned into Δγn = 0.1 kg m−3 density bins, spanning between 20 < γn < 28 kg m−3. The OHC trend located at layers denser than γn > 27.75 kg m−3 is not discussed because these layers are not fully sampled within the 2000 m-depth Argo profiles, which also precludes us from producing analyses for deeper depths.

Since the OHC represents the vertical integration of temperature for each isopycnal layer, we separate the temporal change in OHC into two contribution terms24:

∂(OHC)/∂t = Cp ∫ z2(x,y,z) z1(x,y,z) ρ(x,y,z) ∂θ(x,y,z,t)/∂t dz + Cp ∫ z2(x,y,z,t) z1(x,y,z,t) ∂ρ(x,y,z,t)/∂t θ'(x,y,z) dz

where ρ' and θ' are the climatologically averaged potential density and temperature over the Argo era 2005–2020 respectively, and the last two terms represent changes in OHC associated with changes in thickness and temperature over the climatologically fixed temperature and density field, respectively. The terms z'1 and z'2 indicate the fixed location and depth/thickness of a density layer over the Argo era, based on the climatologically averaged density field. Thus the warming component of OHC changes is estimated based on fixed (Argo era averaged) geographic locations and depths of isopycnals, with a constant volume being used to integrate the warming effects of each water-mass layer and regionally defined water mass over the Argo period.

Note that the calculation of OHC is relative to a fixed reference temperature, making the estimate of trends in OHC and its thickness dependent on the somewhat arbitrary choice of reference temperature when the volume of a water mass changes76. In contrast, the temperature component of OHC change is not dependent on the choice of reference temperature. In the calculations presented here, we take the reference temperature to be 0°C for the global ocean, following similar choices made in past work20,24. We also repeated our calculations of OHC trends using −2°C (the approximate freezing point of seawater) as the fixed reference temperature and found the OHC trends of water masses to be overall robust to this choice, although the relative changes are more sensitive at cooler temperatures. Nonetheless, our findings are overall robust to the choice of reference temperature, for reasonable choices of this parameter over the global ocean. Another simplification is that the effects of salinity changes on the thickness change of a water-mass layer, and thus the heat content change of the layer, have not been separately considered in Eq. (2). Salinity changes can sometimes be large enough to impact the ocean stratification and density structure, and thus influence the ocean circulation and layer thickness77, impacting the OHC within a given density layer. These salinity-driven impacts on ocean sequestration of heat require further investigation but are beyond the scope of this study.

Water-mass framework
To study heat content changes within the global water-mass layers and in regional mode and intermediate waters, we first divide the upper 2000 m of the ocean into the tropical water layer (TW, γn < 25.0 kg m−3), mode water layer (MW, 25.0 ≤ γn < 27.05 kg m−3), and intermediate water layer (IW, 27.05 ≤ γn < 27.75 kg m−3), as sketched in Fig. 4e and given in Supplementary Tables 1 and 2. The mode water layer is further divided into the Subtropical Mode Water layer (STMW, 25.0 ≤ γn < 26.45 kg m−3) and the Subpolar Mode Water layer (SPMW, 26.45 ≤ γn < 27.05 kg m−3), to separately consider heat content changes in specific Subtropical Mode Waters and Subantarctic Mode Water within these two layers.

The regional water masses of individual mode and intermediate waters are defined within the mode and intermediate water layers, using specific longitude, latitude, and neutral density constraints in Supplementary Table 3 (see also in Figs. 3 and 7). For instance, in the Subtropical Mode Water layer, we separately define regional Subtropical Mode Waters at the subtropical South Atlantic26,47,78,79 (SA-STMW), North Atlantic26,47,80–82 (NA-STMW, also known as the North Atlantic’s Eighteen Degree Water), western South Pacific26,47,83,84 (SP-WSTMW), eastern South Pacific26,47,85 (SP-ESTMW), western North Pacific26,47,86–88 (NP-STMW), eastern North Pacific26,47,86,88 (NP-ESTMW), and Indian Ocean26,47,89,90 (IO-STMW), by further using a set of more specific density and geographic limit (Figs. 3, 7, and 8b, e). The Central Mode Water in the North Pacific26,47,86,88 (NP-CMW) is also defined within the Subtropical Mode Water layer. In the Subpolar Mode Water layer, we separately define the Subantarctic Mode Water in the Southern Ocean south of 30°S26,49, and the Madeira Mode Water in the subtropical North Atlantic26,47,82 (NA-MMW) just beneath the NA-STMW (Figs. 7 and 8c, e). In the intermediate water layer, we define the Antarctic Intermediate Water in the South Atlantic26,27,61 due to sustained warming of the intermediate water layer in the Atlantic basin8 and in the Indian and Pacific sectors of the Southern Ocean (Figs. 7 and 8d, e). We also define the Subpolar Mode Water in the subpolar North Atlantic26,47,82,91 (NA-SPMW) and the North Pacific Intermediate Water26,92,93 (NP-IW), which occupies both the Subpolar Mode Water layer and the intermediate water layer. The heat content changes of all these regionally defined modes and intermediate waters are provided in Figs. 8e and 9g, and Supplementary Tables 3 and 4.

Note that the estimate of OHC changes within the mode and intermediate layers and in regionally defined water masses in Figs. 3–9 and Supplementary Tables 1–4 are limited to the Argo era 2005–2020 when the subsurface thermal structure of the ocean was widely observed by the Argo network from 2005. This estimation is based on the ensemble average of SIO RG-Argo, IAP data, and EN4.2.2 ensemble mean. For simplicity, this study uses a set of fixed longitude, latitude, and density or depth criteria to define these water masses and estimate how their OHC has changed over the Argo era, even though the oceanographic features that define them (e.g., temperature, density, salinity and stratification) have also evolved as the upper ocean has warmed and stratified over the Argo era. However, we do nonetheless separate out the volumetric effects from water-mass transformation in attributing total OHC change to each water mass. By decomposing the OHC change into its volumetric and warming components (Eq. 2), the warming component of OHC change within the water-mass layers and in regionally defined water masses (seen in Figs. 5 and 9, and Supplementary Fig. 2 and Tables 2 and 4) is estimated based on a set of fixed longitudes, latitudes, and isopycnal depths (referred to as “fixed volume” or “volume effects removed” in this study) averaged over the Argo era, while the volumetric component is based on a set of density ranges whose depths and volume can evolve over time. This decomposition allows us to distinguish OHC changes of a water mass due to warming effects separate from isopycnal layer thickness (or volume) changes.

OHC trend and heat uptake
Linear trends of OHC are calculated within each selected time period for spatially integrated annual mean OHC timeseries, to indicate the increase or decrease in OHC across each selected period. For the trend maps, the fields are first vertically integrated over depth or isopycnal layers, and then the linear trends over the selected period are computed for each grid point. The ocean heat uptake in Fig. 1b is calculated as the linear trend of annual mean OHC for every 11-year period across the measurement era (e.g., 1960–1970, 1970–1980, 1980–1990, 1990–2000, 2000–2010, 2010–2020) multiplied by the length of period. For linear trends within a period of less than 15 years, decadal variability can be part of the signal as well as long-term trends due to anthropogenic warming17. Another factor is start- and end-point sensitivity in estimating linear trends and capturing nonlinearity in rates of OHC change17. Overall, the monotonic increase in ocean heat uptake has remained clear since the 1980s.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Li et al. invesfigated the ocean heat content changes in the global mode and intermediate water\nlayers, found that 94% of the total OHC increase are stored in these layers. They further invesfigated\nthe contribufion from thickness change and warming effect, suggesfing a ~40% of warming\ncontribufion in these layers. The novelty of this study is a comprehensive analysis of OHC change in\ndensity layers associated with water masses, which is a valuable perspecfive and a fresh insight.\nThus, I do think the paper is worthy a publicafion in NC. I have some comments listed below, mainly\nsupport to increase the robustness of the quanfificafion and strengthen the conclusion. Some of the\ncomments are about terminology and presentafion.\nComments:\n1) The use of data. This paper mainly relies on Levitus data for long-term changes and Roemmich\ndata for Argo-period analysis. Using single dataset is potenfial risky because they are suffered from\nbiases, an example is Cheng&Zhu, 2014 GRL study idenfified a shift in Levitus yearly fime series\nassociated with data sampling (2000~2005 during the observafion system transfer). It is increasingly\nclear that sampling is a serious issue. Thus, I would suggest to at least use ensemble mean of several\nbefter datasets (Levitus pentadal esfimates, Ishii et al. 2017 and Cheng et al. 2017 data). This could\nprovide you more reliable quanfificafion. Roemmich dataset might also not be ideal because of not\nglobal coverage (ITF regions, polar regions are not included) and other issues (for instance Trenberth\net al. 2016 indicates a quite larger noise level).\n2) Figure 1, 3 and other places, it is not proper to label 2010-2020 data as Argo because of two\nreasons: first, Levitus (together with Ishii/Cheng data in my first comment) also used Argo data;\nsecond, there are several Argo gridded products, so not only RG09 data is “Argo data”. The benefit of\nLevitus/Ishii/Cheng data is they have real global coverage by merging Argo with other plafforms. This\ncould be a thing that impact your quanfificafion. Note that some previous analyses indicate the\ndifference of ocean area could contribute to ~10% of trend difference in OHC.\n3) Figure 2. It is a bit misleading to add shading area around the fime series because it looks like error\nbars. I suggest to have some rectangle shading in the panel to denote the fime periods and add\nactual error bars if possible.\n4) Uncertainty assessment should be given for any quanfificafions in this study. One way is to use the\nspread among different products.\n5) “Accelerafing global ocean warming” secfion. Please find a recent study of Cheng et al. 2022 for\nquanfifying the rate change of OHC, this paper could strongly support this part as they do a more\nthorough analysis to deal with the short-term fluctuafions.\n6) “Levitus period” is defined in this study, I don’t support this definifion because Levitus does not\ndefine any period, it starts from 1955 because of befter data coverage.\n7) Lines 111-117. ENSO: this is not only because of the wind-driven variability, but also air-sea net\nheat exchanges (see Cheng et al. 2019 ENSO related changes of OHC paper). The thermocline change\nwould not cause net OHC change when averaging over the tropics, the only mechanism can lead to\nnet OHC change is heat exchange at interfaces (e.g. air-sea, tropical vs. middle-lafitudes, upper vs\ndeep layers).\n8) Lines 120-onward. I’m wondering why 2010-2020 period is the focus, the short period means lots\nof internal variability mixed into anthropogenic signals as you discussed within 133-147.\n9) Figure 3 and 4 can be combined to save space: add zonal mean plot to the right side of the spafial\nmaps.\n10) Figure 3 vs. Figure 5. Similar to Fig.5, Figure 3 is also a trend map (difference between 2000-2010\nand 1980-2000). Comparing the two trend maps indicate quite a significant regional inconsistency,\nindicafing to me that ENSO/IPO could be important at different fime periods. Do these variabilifies\nimpact your conclusion (i.e. mode/intermediate waters contribufion to total OHC)?\n11) Table 1. Von Schuckmann et al. 2020 ESSD had a figure nicely summarizing previous efforts and\nalso show the increasing trend of OHC based on a more complete list of literature. I would also\nremove this table to save space (or move to the SI).", "reviewer_2": "Review of “Recent accelerafion in global ocean heat accumulafion by mode\nand intermediate waters” by Li, England, Groeskamp.\nGeneral Comments:\nI have enjoyed reading this paper. It has a wide appeal to oceanographers and to a wider\ncommunicafions audience. This appeal is because the manuscript documents in a more thorough\nway than previously the changes in a key water-mass type that is found widely across the global\noceans. The global perspecfive of this paper means that it details a coherent set of changes in both\nnorthern and southern hemispheres. It also shows, quite remarkably in my mind, how this key water\nmass has actually dominated the total heat uptake in the global oceans, a point that I had known\nintuifively but had not quite appreciated the scale of that uptake. Finally the heat uptake increase in\nthe most recent period (in both hemispheres, e.g. figure 4 of the manuscript) is also extra-ordinary.\nThe result that these mode waters have an uptake of about 94% of the total ocean heat uptake in\nthese periods from the 0-2000m layer and they occupy just 16% of the global ocean volume.\nIt really means that from the context of ocean heat uptake, these mode (and intermediate)\nwater are the crifical water mass.\nWhile not stated, it also implies that the other 84% of the global ocean volume has only received 6%\nof the net warming. In some ways its also good news, because these mode waters are far from the\nPoles and therefore “buffering” to some extent the very negafive changes that could occur by the\nmelt of the Polar Ice Cap (Antarcfica and Greenland), and in terms of the oceanic carbon cycle\nprobably helpful in sequestering anthropogenic carbon in these same regions. Also not stated, but\nrelevant to how we might sample the oceans to monitor it and to “predict” how it might change,\nthus tracking the earth’s energy imbalance.\nSo this paper is very interesfing for its relevance in the global climate system. It uses relafively\nstandard techniques and approaches to assess the ocean heat content changes for this key (but\nvaried water mass) including the contribufions from density (or volume change) and heat storage.\nMy belief is that this is a very relevant and fimely paper for nature communicafions.\nThat said, the paper is not without gaps, some minor and some more major (but solvable). This\npaper is really a statement of the changing state of the upper 2000 metre of the ocean. Its not a full\ncensus of the full depth of the ocean. This point should be made more clearly. Unlike modern\nsimulafions of the climate system, or re-analysis of the climate system, this paper basically relies on a\nsingle gridded data product (RG09) and as support the Satellite Alfimetry product (DUACS and\ndistributed by EU CMEMS) and for earlier periods the Levitus product from WODC (reference 5).\nThere are actually at least 6 different gridded products and they do vary between each other when\nthese OHC calculafions are done. IPCC would look to provide an ensemble view of the observafions\nrather than relying on a single product. The second gap is that not a single table, or a single graph has\nany sense of an error bar on them. What is the uncertainty in these 10 year trends for each of the\nten year periods provided in this manuscript? Are these trends real at the different layer averages, and\ndifferent spafial scales, grid point to global averages? Is the heat uptake as extreme as it seems\nfrom the abstract and results. Some would say that no errors bars means the authors don’t know the\nuncertainty and so it is very large (rather than very small) and therefore none of the results are valid.\nMy familiarity with these data sets tells me that the results would be replicated across the suite of\nthe products. These two points are part of the solufion though. By using the 6 or so available gridded\nproducts one can immediately get the ensemble esfimate of the changes in these key water masses,\nand also get from each data set 6 esfimates of every quanfity presented in the paper and hence the\nspread and a completely jusfifiable error bar on every result presented. Perhaps the most important\nresults that would be befter for confidence intervals would be Figure 1, Figure 2 panels a,b,c,d,\nFigure 4, Figure 7 a, Figure 8 g. Table 1, Table 2. The mapped fields/panels in the figures could include\nthe places where change is greater than some confidence interval, but actually I think that is not so\nimportant.\nOtherwise the figures are quite clear.\nWhile the message of this paper is clear it does leave a number of quesfions unresolved. It implies\nthat these changes are a simple accelerafing trend and are anthropogenically forced. The fitle says\nthe changes are accelerafing. But there is no evidence in this paper that it is an accelerafion that\ncan't be explained by say internal variability and climate change (say). The conclusions (line 248 to\n254) could be easily rewriften to address whether the increase in uptake relafive 1990-2000 is\nanthropogenically induced or internal variability and climate change. I totally agree its a key water\nmass in the climate system and the results from this paper are compelling. BUT it could be\nimplausibly high and in 2020-2030 period or the 2030-2040 period we could be wrifing a similar\npaper but on the decelerafion of mode and intermediate water uptake. I think these are the types of\nquesfions this paper should be addressing in lines 248 to 254 (as it sums up the conclusions), and\nperhaps as a locafion on the planet that is very important for monitoring the state of the climate.\nSpecific Comments:\nI like the decomposifion of Ocean Heat Content in equafion 2. It follows earlier work (that is referred\nto). And in general this equafion is spot on, but there have been some papers talking very specifically\nabout the role of ocean heat content in sub-tropical gyres where the salinity changes within the\nwater column are large enough to induce a thickness change and thus on density surface cause\ncooling signal see on the equatorward side of the gyres in the Indian Ocean and presumably in the\npacific ocean. So there should be a note that the thickness change calculafion has not taken into\naccount the effects of salinity changes on thickness change on OHC in the methods secfion.\nLevitus Era or Levitus period (Figure 1 capfion). I dont like this terminology. Syd Levitus should have\nmore fame but the period 1955-2010 is not a period of fime that we should name after Levitus. We\ncertainly could use the term Levitus Atlas. I would prefer simply “pre Argo”, or simply 1955-2010.\nThe IPCC is wary of the ocean thermal data before the 1970’s typically because these earlier data are\nunreliable (see the discussions about xbt and fall rates), and the trends 1960 to 1970 and 1970 to\n1980 have clear variability in term. We could also eliminate the “Roemmich period” as well (Figure 2\ncapfion)", "reviewer_3": "his manuscript analyzes ocean heat content using yearly maps of ocean heat content from 0-2000 m\nfrom the World Ocean Atlas 2018 and an updated monthly maps of temperature (as a funcfion of\npressure, as well as density) from Roemmich and Gilson (2009). It focuses on contribufions of mode\nand intermediate water warming to ocean warming as a whole. While this work is interesfing, and\nshould eventually be suitable for publicafion, there are substanfial issues that preclude publicafion of\nthe manuscript in its present form. Some of the major issues are expanded upon below, and then a\nfew specific comments are made.\n1. First, neither the Roemmich and Gilson (2009) nor the World Ocean Atlas maps are immune to the\neffects of changes in oceanographic sampling (e.g., Lyman and Johnson, 2014, J. Climate). The\nRoemmich and Gilson (2009) monthly maps relax back to the long-term seasonal climatolgy in\nlocafions and fimes where data are non-existent or sparse. This means that the first few years of\nmaps are somewhat suspect, as are the global integrals from it in those years as shown in Figures 1\nand 2. Similarly, the World Ocean Atlas anomalies in the earlier years are based on very, very sparse\ndata (especially prior to the late 1960s when the XBT came into widespread use), and analyses of\nthese maps should be treated with caufion. This means the first decade in figure 1 is very suspect,\nand really the increasing spafial (and depth of sampling) coverage as fime goes on plays a role in the\nshape of the curve. This should at least be acknowledged in the claims of dramafic accelerafion. This\nis not to say that there is not accelerafion, just that some of it may be owing to increases in sampling\n(both in terms of ocean regions and depths) as the decades go on. It also means that the pafterns in\nthe first few years in the panels of Figure 2 are suspect.\n2. Second, mode and intermediate waters are defined by oceanographic features that are expected\nto change density over fime. For instance, the temperature, salinity, density, and depth of the salinity\nminimum that defines the core of the Antarcfic Intermediate water has evolved over the decades as\nthe upper ocean has warmed and upper ocean salinity pafterns have intensified over the decades\n(e.g., Schmidtko and Johnson, 2012, J. Climate). Similarly, mode waters are defined as a layer of low\nstrafificafion. The water properfies of those mode waters (including their temperature, salinity,\ndepth, or density) will change over fime as upper ocean water properfies change over fime, so their\nisopycnal boundaries might also change with fime These facts need to be acknowledged, and the\ninterpretafion of results (especially in Figure 2f, but also in the maps in Figures 7 and 8) should keep\naccount for this possibility.\n3. Third, as the authors are aware, but do not emphasize enough in the manuscript, analyses of\nocean heat content in isopycnal layers can be dominated by volume changes between those layers,\nrather than temperature/salinity changes in those layers. Since the density boundaries defining these\nmode waters may change with fime, one is left asking what these volume changes actually mean?\nFor instance, if the density of the actual mode water shifted across an isopyncal boundary (say it\nbecame substanfially less dense in a warming world), what would a volume change for that layer\nmean? Including the volumetric change without changing the isopyncal boundary changes\nwould overemphasize the heat content change.\nFor the reasons discussed in the preceding two paragraphs it might actually be simpler, and perhaps\nbefter, to define a set of fixed depths (they could be defined by isopycnal depths at a specific fime, or\nfor a decadal average) to analyse different mode waters, rather than allowing those depths to follow\nisopycnals as they migrate with fime.\n5. It is not good pracfice to name fime periods after people. Please find other names than “The\nRoemmich period” or “The Levitus period”.\n6. As presented, the units of J yr-1 in Figures 3, 4, 5a, 5c-d, 6, 7, and 8 are all resolufion dependent,\nand need to include area units for the maps (e.g., J m-2 yr-1) and a length unit for the line plots in\nFigure 4 (J m-1 yr-1). Also, the units of Figure 2f are also resolufion dependent as presented, and\nneed to be in J kg-1 m3.\n7. Lines 78-89 and elsewhere: There are a number of different groups esfimafing ocean warming as is\ndocumented in von Shuckman et al. (2020, Earth Syst. Sci. Data), where the accelerafion noted here\nis also discussed in some detail. In addifion Loeb et al. (2021, Geophys. Res. Left.) also document the\naccelerafion of ocean warming 2005-2019, and show that it agrees with increases in net-top-of the\natmosphere energy fluxes observed by satellite.\n8. Lines 104-108. Does this perhaps signal actually a change in the density of the mode\nwaters (becoming lighter with fime in a warming ocean), rather than a change in volume of intermediate\nversus mode waters?\n9. Line 125. Here it would be appropriate to note the work of Roemmich, who noted previously on\nthat recent warming has been concentrated in the southern Ocean.\n10. Lines 170-186. This explicafion (and others elsewhere which include volume changes) really\noveremphasizes the importance of warming in the mode and intermediate waters. Volume changes\ncan result in huge “heat” changes, but since heat changes really should be evaluated in a constant\nvolume, they don’t mean much and need to be explained very carefully."}, "author_response": "RESPONSE TO REVIEWER COMMENTS:\nBelow reviewer comments have been copy-pasted in black font, and our\nresponses are indicated in blue font. Line numbers in our responses refer to the\nrevised manuscript (except where indicated).\nReviewer #1 (Remarks to the Author):\nLi et al. investigated the ocean heat content changes in the global mode and\nintermediate water layers, found that 94% of the total OHC increase are stored in these\nlayers. They further investigated the contribution from thickness change and warming\neffect, suggesting a ~40% of warming contribution in these layers. The novelty of this\nstudy is a comprehensive analysis of OHC change in density layers associated with\nwater masses, which is a valuable perspective and a fresh insight. Thus, I do think the\npaper is worthy a publication in NC. I have some comments listed below, mainly support\nto increase the robustness of the quantification and strengthen the conclusion. Some of\nthe comments are about terminology and presentation.\nResponse: We thank Reviewer #1 for the positive assessment and detailed comments\nin this round of review. We believe these comments have significantly helped us to\nimprove the paper.\nReviewer comments:\n1) The use of data. This paper mainly relies on Levitus data for long-term changes and\nRoemmich data for Argo-period analysis. Using single dataset is potential risky because\nthey are suffered from biases, an example is Cheng&Zhu, 2014 GRL study identified a\nshift in Levitus yearly time series associated with data sampling (2000~2005 during the\nobservation system transfer). It is increasingly clear that sampling is a serious issue.\nThus, I would suggest to at least use ensemble mean of several better datasets (Levitus\npentadal estimates, Ishii et al. 2017 and Cheng et al. 2017 data). This could provide you\nmore reliable quantification. Roemmich dataset might also not be ideal because of not\nResponse to Reviewers\nglobal coverage (ITF regions, polar regions are not included) and other issues (for\ninstance Trenberth et al. 2016 indicates a quite larger noise level).\nResponse: Thanks for pointing this out and providing the insightful suggestions. This\ncomment picks up on points also raised by Reviewer 2. We agree that this is a possible\nsource of uncertainty by relying only on Levitus data for the pre-Argo period analysis\n(1955–2004) and Roemmich Argo data for the Argo-period analysis (2005–2020). To\naddress this issue and improve the robustness of our results, we have now included\nobservational-based ocean heat content (OHC) products from Levitus et al. (2012),\nCheng et al. (2017), and Ishii et al. (2017), and also included the estimate of global and\nregional OHC using temperature and salinity observations from SIO RG-Argo, IAP data\n(Cheng and Zhu 2016; Cheng et al. 2017), and the ensemble mean of Met Office\nHadley Centre EN4.2.2 data (Good et al., 2013), into the manuscript. We have now\nadded a detailed description of this point in Methods of the revised manuscript (lines\n356–404) and updated all results and display diagrams in the manuscript accordingly.\nWe have also included a statement of the uncertainty that may come from using\nSIO RG-Argo, which does not have a global coverage, on lines 367–371: “We only use\nthe SIO RG-Argo data between 65°S-65°N because of the low observational data\ncoverage in high latitude regions, particularly ice-covered areas and in marginal seas.\nAs such, the resultant global OHC trends from SIO RG-Argo in Fig. 1 will be\nunderestimated if these polar ocean regions are warming5.” Also, a clarification of the\nshift in OHC estimates that is due to the change in locations of the observational\nnetwork has been added to lines 128–131: “For example, increased spatial sampling\ncoverage, particularly over the tropical and Southern Hemisphere oceans due to the\nshift in location of the observational network from a ship-based system to the Argo\nnetwork, can affect the estimate of ocean warming during 2001–200334 (Fig. 1,\nSupplementary Fig. 1).”.\n2) Figure 1, 3 and other places, it is not proper to label 2010-2020 data as Argo\nbecause of two reasons: first, Levitus (together with Ishii/Cheng data in my first\nResponse to Reviewers\ncomment) also used Argo data; second, there are several Argo gridded products, so not\nonly RG09 data is “Argo data”. The benefit of Levitus/Ishii/Cheng data is they have real\nglobal coverage by merging Argo with other platforms. This could be a thing that impact\nyour quantification. Note that some previous analyses indicate the difference of ocean\narea could contribute to ~10% of trend difference in OHC.\nResponse: Thanks for the suggestions. We have now referred to the gridded Argo\nproduct from Roemmich and Gilson (2009) as Scripps Institution of Oceanography\n(SIO) RG-Argo, and have additionally included the EN4.2.2 ensemble mean from the\nMet Office Hadley Center (Good et al., 2013), the IAP data (Cheng and Zhu 2016;\nCheng et al. 2017), and the Ishii et al. (2017) data into our analyses. These datasets are\nconstructed by merging both Argo and other historical observations and are all global in\ncoverage. As such, the inclusion of these datasets provides a more robust estimate of\nglobal OHC changes. We have clarified this point on lines 367–371 “This gridded Argo\nproduct was constructed from Argo data only. We only use the SIO RG-Argo data\nbetween 65°S-65°N because of the low observational data coverage in high latitude\nregions, particularly ice-covered areas and in marginal seas. As such, the resultant\nglobal OHC trends from SIO RG-Argo in Fig. 1 will be underestimated if these polar\nocean regions are warming5.”\nIt may also be noted that in Figs. 2–10 of the manuscript, where we show OHC\nchanges in global mode and intermediate layers and within localized mode and\nintermediate waters, we focus on the region between 65°S-65°N due to the low\nsampling coverage at higher latitudes and for better consistency of the data’s spatial\ncoverage, as now clarified in the caption of Figs. 2, 3 and Supplementary Table 1, and\non line 135–138: “Therefore, this study mainly focuses on the Argo era 2005–2020, an\noverlapping period common to all products analyzed, and on the ocean domain\nbetween 65°S and 65°N, because of the sparse observations in the polar regions even\nafter the Argo array was established.”.\n3) Figure 2. It is a bit misleading to add shading area around the time series because it\nResponse to Reviewers\nlooks like error bars. I suggest to have some rectangle shading in the panel to denote\nthe time periods and add actual error bars if possible.\nResponse: We thank the reviewer for pointing this out. We have now included the OHC\nestimate based on temperature and salinity observations from SIO RG-Argo, the IAP\ndata (Cheng and Zhu 2016; Cheng et al. 2017), and EN4.2.2 ensemble mean from the\nMet Office Hadley Center (Good et al., 2013), into the Argo period analyses, and have\nupdated Figs. 2–5, 9, 10, Supplementary Fig. 1, and Supplementary Tables 1–4,\nregarding this point (see more details in Methods, lines 356–404). Now both the OHC\nestimate from each of these three products and their ensemble mean are provided in\nFig. 2. The uncertainty range (±2 standard deviations) from these three products is\nprovided in Figs. 3a-d and Supplementary Tables 1, 2.\n4) Uncertainty assessment should be given for any quantifications in this study. One\nway is to use the spread among different products.\nResponse: Thanks for the suggestion. We have now included observational-based\nOHC products from Levitus et al. (2012), Cheng et al. (2017), and Ishii et al. (2017), and\nalso included the estimate of global and regional OHC using temperature and salinity\nobservations from SIO RG-Argo, IAP data (Cheng and Zhu 2016; Cheng et al. 2017),\nand the ensemble mean of Met Office Hadley Centre EN4.2.2 data (Good et al., 2013),\ninto the manuscript. We have now updated all results and figures in the manuscript\nbased on the inclusion of these new datasets. Now both the ensemble mean and the\nuncertainty range (±2σ) from all products are provided in the manuscript. We have now\nincluded a statement of this point on lines 101–104: “To provide a more reliable analysis\nof the global and regional OHC changes in this study, we present the ensemble average\nof OHC estimates derived from various ocean data products (Fig. 1; Methods). The\nuncertainty range is given by ±2 standard deviations (±2σ) to indicate the degree of\nagreement across the products.”, in figure and table captions, and on lines 397–404.\n5) “Accelerating global ocean warming” section. Please find a recent study of Cheng et\nResponse to Reviewers\nal. 2022 for quantifying the rate change of OHC, this paper could strongly support this\npart as they do a more thorough analysis to deal with the short-term fluctuations.\nResponse: Thanks for the suggestions. We have now added a citation to Cheng et al.\n2022 (Journal of Climate) in the manuscript (e.g., lines 50, 107) and also added a\nclarification of this point in this section (lines 114–117): “Further comparison with\nprevious estimates of OHC trends (Table 1) and the latest assessment report of the\nIntergovernmental Panel on Climate Change (IPCC) Working Group I7 also reveals\nsignificant increases in the ocean warming rate across the measurement record5,8–\n11,16,17,19.”, and on lines 126–127: “Another factor is start- and end-point sensitivity in\nestimating linear trends and capturing nonlinearity in rates of OHC change16.”\n6) “Levitus period” is defined in this study, I don’t support this definition because Levitus\ndoes not define any period, it starts from 1955 because of better data coverage.\nResponse: Agreed – thanks, good suggestion. We have now referred to the period\n1955–2004 as the ‘pre-Argo period’ and the period 2005–2020 as the ‘Argo era’ in the\nmanuscript. The terminologies of ‘Levitus period’ and ‘Roemmich period’ have all been\nremoved from the manuscript.\n7) Lines 111-117. ENSO: this is not only because of the wind-driven variability, but also\nair-sea net heat exchanges (see Cheng et al. 2019 ENSO related changes of OHC\npaper). The thermocline change would not cause net OHC change when averaging over\nthe tropics, the only mechanism can lead to net OHC change is heat exchange at\ninterfaces (e.g. air-sea, tropical vs. middle-latitudes, upper vs deep layers).\nResponse: Agreed - thanks for raising this. We have now rephrased the description of\nENSO related changes of OHC in this section (lines 162–164 of revised manuscript):\n“This ENSO-driven variability in OHC is related to net air-sea heat exchanges and the\nredistribution of heat laterally across different latitudes and ocean regions, and vertically\nbetween surface and subsurface layers35,36.”, and added a citation to Cheng et al.\n(2019) in the manuscript (lines 162, 164).\nResponse to Reviewers\n8) Lines 120-onward. I’m wondering why 2010-2020 period is the focus, the short period\nmeans lots of internal variability mixed into anthropogenic signals as you discussed\nwithin 133-147.\nResponse: Agreed, part of this came about because we initially focused on the decade-\nby-decade accumulation of heat in the ocean, but the downside is a signal of internal\nvariability, as the reviewer notes. To reduce the impacts from internal variability and\ntake advantage of available observations, we have now changed our analysis period to\nfocus on the Argo era 2005–2020 (an overlapping period of all datasets used), when\nArgo data became mature in spatial coverage from 2005. We have now updated Figs.\n2–4, 6–10 and Supplementary Tables 1–4 accordingly. We have now added a\nstatement in regard to this point on lines 131–138: “Even though all products show\nrobust ocean warming rates for the upper 2000 m over the past few decades, reaching\nrecord warming over 2010–20205,16 (Table 1), these products reach better agreement\nafter ∼2005 when the Argo network was widely deployed. To minimize the impact of\ninternal climate variability and focus on longer-term trends, we consider the longest\nreliable OHC record available. Therefore, this study mainly focuses on the Argo era\n2005–2020, an overlapping period common to all products analyzed, and on the ocean\ndomain between 65°S and 65°N, because of the sparse observations in the polar\nregions even after the Argo array was established.”.\n9) Figure 3 and 4 can be combined to save space: add zonal mean plot to the right side\nof the spatial maps.\nResponse: Agreed - good suggestion, thanks. We have now combined the previous\nFigs. 3 and 4 as the new Fig. 5. Due to the inclusion of additional observational-based\ndatasets, the ensemble average of decadal mean OHC change is now shown in the\nspatial maps (Figs. 5a-b), and both the ensemble mean and the uncertainty range (±2σ)\nare provided in the zonal mean plot (Figs. 5c-d). This is a nice improvement to these\nanalyses.\nResponse to Reviewers\n10) Figure 3 vs. Figure 5. Similar to Fig.5, Figure 3 is also a trend map (difference\nbetween 2000-2010 and 1980-2000). Comparing the two trend maps indicate quite a\nsignificant regional inconsistency, indicating to me that ENSO/IPO could be important at\ndifferent time periods. Do these variabilities impact your conclusion (i.e.\nmode/intermediate waters contribution to total OHC)?\nResponse: Thanks for raising this. We agree that there are regional inconsistencies in\nOHC changes between Fig. 3 (now Fig. 5) and Fig. 5 (now Fig. 6). This occurs mainly\nover the tropical Pacific Ocean, whereas the rest of the Pacific Ocean and the Indian,\nAtlantic, and Southern Oceans tend to show a consistent and robust warming, in\nparticular over subduction/ventilation regions near the western boundary current\nextensions and the ACC (Figs. 5a-b, 6a, 8; Rathore et al. 2020; Zika et al. 2020). This\nrobust long-term warming over the subtropical Atlantic and Pacific Oceans and in the\nSouthern Ocean, as well as the near-monotonic warming in the mode and intermediate\nwater layers across the Argo era (Figs. 2, 3a; line 150–153), suggest a robust ocean\nheat gain by mode and intermediate waters. In contrast, OHC changes within the\ntropical water layer and in the tropical Pacific Ocean reveal large variations associated\nwith ENSO and/or IPO variability (Figs. 2b, 3b, 5a-b, 6a; lines 158–169 and 187–194).\nWe further agree that the formation of mode and intermediate waters can be\ninfluenced by interannual to decadal variability of the climate system. For example,\ninterannual variations in the formation of South Pacific SAMW have been linked to\ntropical teleconnections from the Indo-Pacific (Li and England 2020). In order to reduce\nthe footprint from climate modes to our conclusion, we have now changed our main\nwater-mass analyses to focus on a longer period over the Argo era 2005–2020 (an\noverlapping period of all the datasets analyzed). We have now added a clarification of\nthis point on lines 131–138 and lines 294–307: “While the globally integrated OHC of\nthe upper 2000 m and in the mode and intermediate water layers both show a robust\nand near-monotonic increase (Figs. 1–3), variations in the formation mechanisms,\nproperties, and heat content of localized mode and intermediate waters can be\ninfluenced by interannual to decadal variability (e.g. ENSO and the Indian Ocean\nResponse to Reviewers\nDipole60 on interannual time-scales, and the IPO on decadal time-scales). However, our\nanalyses show that the main footprint of internal climate variability occurs over the\ntropical Pacific Ocean, whereas the rest of the Pacific as well as the Indian, Atlantic,\nand Southern Oceans all tend to show a consistent and robust warming, in particular\nover the subduction and well-ventilated regions near western boundary current\nextensions and north of the ACC (compare Figs. 5a-b and 6a). The near-monotonic\nwarming in mode and intermediate water layers across the Argo era (Figs. 2, 3a),\nsuggests a robust ocean heat gain forced by anthropogenic change. In contrast, OHC\nchanges within the tropical water layer and in the tropical Pacific Ocean reveal\nvariations likely associated with ENSO and/or IPO variability.”. Text has also been\nadded on this point at lines 328–331: “Although our results are robust within a set of\nobservational-based datasets, important questions remain, such as separating out the\nanthropogenic signal of OHC change from internal variability, particularly in mode and\nintermediate waters where decadal variability is poorly constrained by sparse\nmeasurements prior to the Argo era.”\nWe have also added text on the point of continuing establishing and improving\nsubtropical, polar, and deep ocean observing systems, which is of fundamental\nimportance for better monitoring the state of the ocean and understanding the\nmechanisms that are driving ongoing ocean warming, at lines 325–342 and 343–347,\nas per the suggestion from Reviewer 2.\n11) Table 1. Von Schuckmann et al. 2020 ESSD had a figure nicely summarizing\nprevious efforts and also show the increasing trend of OHC based on a more complete\nlist of literature. I would also remove this table to save space (or move to the SI).\nResponse: Thanks for the suggestion. We have now updated Table 1 and decided to\nkeep it in the main manuscript for readers to have ready access to this table showing\naccelerating ocean warming built on past studies. So we feel it’s worth keeping the\nTable. We have nonetheless now included a citation to von Shuckman et al. (2020) and\nthe latest assessment report of the IPCC Working Group I in the manuscript text, lines\n114–123: “Further comparison with previous estimates of OHC trends (Table 1) and the\nResponse to Reviewers\nlatest assessment report of the Intergovernmental Panel on Climate Change (IPCC)\nWorking Group I7 also reveals significant increases in the ocean warming rate across\nthe measurement record5,8–11,16,17,19. This increasing ocean warming rate is not just\nlimited to the upper layer of the ocean, but also occurs at intermediate depths (700–\n2000 m) according to a further breakdown of OHC trends by depth over the upper 2000\nm5,7. Significant ocean warming and accelerating OHC changes are consistent with the\nincrease in net radiative energy absorbed by Earth detected in satellite observations5,33,\nsomething that is likely to continue throughout the 21st Century6,15 in the absence of\nsubstantial greenhouse gas emissions reductions.”.\nResponse to Reviewers\nReviewer #2 (Remarks to the Author):\nReview of “Recent acceleration in global ocean heat accumulation by mode\nand intermediate waters” by Li, England, Groeskamp.\nGeneral Comments:\nI have enjoyed reading this paper. It has a wide appeal to oceanographers and to a\nwider communications audience. This appeal is because the manuscript documents in a\nmore thorough way than previously the changes in a key water-mass type that is found\nwidely across the global oceans. The global perspective of this paper means that it\ndetails a coherent set of changes in both northern and southern hemispheres. It also\nshows, quite remarkably in my mind, how this key water mass has actually dominated\nthe total heat uptake in the global oceans, a point that I had known intuitively but had\nnot quite appreciated the scale of that uptake. Finally the heat uptake increase in the\nmost recent period (in both hemispheres, e.g. figure 4 of the manuscript) is also extra-\nordinary. The result that these mode waters have an uptake of about 94% of the total\nocean heat uptake in these periods from the 0-2000m layer and they occupy just 16% of\nthe global ocean volume. It really means that from the context of ocean heat uptake,\nthese mode (and intermediate) water are the critical water mass.\nWhile not stated, it also implies that the other 84% of the global ocean volume has only\nreceived 6% of the net warming. In some ways its also good news, because these\nmode waters are far from the Poles and therefore “buffering” to some extent the very\nnegative changes that could occur by the melt of the Polar Ice Cap (Antarctica and\nGreenland), and in terms of the oceanic carbon cycle probably helpful in sequestering\nanthropogenic carbon in these same regions. Also not stated, but relevant to how we\nmight sample the oceans to monitor it and to “predict” how it might change, thus\ntracking the earth’s energy imbalance.\nSo this paper is very interesting for its relevance in the global climate system. It uses\nrelatively standard techniques and approaches to assess the ocean heat content\nchanges for this key (but varied water mass) including the contributions from density (or\nResponse to Reviewers\nvolume change) and heat storage.\nMy belief is that this is a very relevant and timely paper for nature communications.\nResponse: We thank Reviewer #2 for this favourable assessment of our paper, and for\ntheir insightful comments in this round of review. We believe these comments have\nhelped us to improve the paper. Thank you.\nThat said, the paper is not without gaps, some minor and some more major (but\nsolvable). This paper is really a statement of the changing state of the upper 2000 metre\nof the ocean. Its not a full census of the full depth of the ocean. This point should be\nmade more clearly.\nResponse: We thank the reviewer for pointing this out. We have now added text to\nclarify that the focus of this study is the changing state of the upper 2000 m of the ocean\nin Methods (lines 412–418): “This study focuses on warming of the upper 2000 m of the\nocean in part because this is the volume of the ocean in closest contact with the\natmosphere, to experience global warming. This 2000-m constraint is also due to the\ndepth limit of the Argo profiling floats, with only sparse data coverage below 2000 m.\nWaters below 2000 m, apart from recently overturned Antarctic Bottom Water and\nLower North Atlantic Deep Water, are generally not ventilated rapidly enough to have\nexperienced the last 50 years of global warming.”.\nUnlike modern simulations of the climate system, or re-analysis of the climate system,\nthis paper basically relies on a single gridded data product (RG09) and as support the\nSatellite Altimetry product (DUACS and distributed by EU CMEMS) and for earlier\nperiods the Levitus product from WODC (reference 5). There are actually at least 6\ndifferent gridded products and they do vary between each other when these OHC\ncalculations are done. IPCC would look to provide an ensemble view of the\nobservations rather than relying on a single product. The second gap is that not a single\ntable, or a single graph has any sense of an error bar on them. What is the uncertainty\nResponse to Reviewers\nin these 10 year trends for each of the ten year periods provided in this manuscript? Are\nthese trends real at the different layer averages, and different spatial scales, grid point\nto global averages? Is the heat uptake as extreme as it seems from the abstract and results.\nSome would say that no errors bars means the authors don’t know the\nuncertainty and so it is very large (rather than very small) and therefore none of the\nresults are valid. My familiarity with these data sets tells me that the results would be\nreplicated across the suite of the products. These two points are part of the solution\nthough. By using the 6 or so available gridded products one can immediately get the\nensemble estimate of the changes in these key water masses, and also get from each\ndata set 6 estimates of every quantity presented in the paper and hence the\nthe spread and a completely justifiable error bar on every result presented. Perhaps the most important\nresults that would be better for confidence intervals would be Figure 1, Figure 2 panels\na,b,c,d, Figure 4, Figure 7 a, Figure 8 g. Table 1, Table 2. The mapped fields/panels in\nthe figures could include the places where change is greater than some confidence\ninterval, but actually I think that is not so\nimportant.\nOtherwise the figures are quite clear.\nResponse: Thanks for pointing this out and providing these useful suggestions. This\ncomment picks up on points also raised by Reviewer 1. We agree that this is a possible\nsource of uncertainty by relying only on Levitus data for the pre-Argo period analysis\n(1955–2004) and Roemmich Argo data for the Argo-period analysis (2005–2020). To\naddress this issue and improve the robustness of our results, we have now included\nobservational-based ocean heat content (OHC) products from Levitus et al. (2012),\nCheng et al. (2017), and Ishii et al. (2017), and also included the estimate of global and\nregional OHC using temperature and salinity observations from SIO RG-Argo, IAP data\n(Cheng and Zhu 2016; Cheng et al. 2017), and the ensemble mean of Met Office\nHadley Centre EN4.2.2 data (Good et al., 2013), into the manuscript. We have now\nadded a detailed description of this point in Methods of the revised manuscript (lines\n356–404) and updated all results and display diagrams in the manuscript accordingly.\nResponse to Reviewers\nIn addition, the ensemble mean estimate of global OHC changes and OHC\nchanges in key water masses from these products is now provided throughout the\nmanuscript, instead of the previous approach, and the uncertainty range (±2σ) is\nprovided in Figs. 1, 3a-d, 5c-d, 9a, 10g, and S1 and Tables 1, S1–4.\nWhile the message of this paper is clear it does leave a number of questions\nunresolved. It implies that these changes are a simple accelerating trend and are\nanthropogenically forced. The title says the changes are accelerating. But there is no\nevidence in this paper that it is an acceleration that can't be explained by say internal\nvariability and climate change (say). The conclusions (line 248 to 254) could be easily\nrewritten to address whether the increase in uptake relative 1990-2000 is\nanthropogenically induced or internal variability and climate change. I totally agree its a\nkey water mass in the climate system and the results from this paper are compelling.\nBUT it could be implausibly high and in 2020-2030 period or the 2030-2040 period we\ncould be writing a similar paper but on the deceleration of mode and intermediate water\nuptake. I think these are the types of questions this paper should be addressing in lines\n248 to 254 (as it sums up the conclusions), and perhaps as a location on the planet that\nis very important for monitoring the state of the climate.\nResponse: We agree that OHC trends within this relatively short time period can\nindicate a combination of both an anthropogenic forced warming signal as well as a\nsignal due to internal variability (lines 187–188). For example, the OHC changes within\nthe tropical water layer and in the tropical Pacific Ocean reveal large variations\nassociated with ENSO and/or IPO variability (Figs. 2b, 3b, 5a-b, 6a; lines 158–167 and\n187–194), while the mode and intermediate water layers exhibit less variability and\nmore steady warming across the Argo era (see Figs. 2, 3a; line 150–153), especially for\nwarming of the intermediate water layer (Figs. 3d, e). Overall, the subtropical Pacific,\nIndian, and Atlantic Oceans, and the Southern Ocean, all show a consistent and robust\nwarming, in particular over subduction/ventilation regions near the western boundary\ncurrent extensions and the ACC (Figs. 5a-b, 6a, 8; Rathore et al. 2020; Zika et al.\n2020). This robust long-term warming over the subtropical Atlantic and Pacific Oceans\nResponse to Reviewers\nand in the Southern Ocean, as well as the near-monotonic warming in the mode and\nintermediate water layers across the Argo era (Figs. 2, 3a; line 150–153), suggest a\nrobust ocean heat gain by mode and intermediate waters.\nWe further agree that the formation of mode and intermediate waters can be\ninfluenced by interannual to decadal variability of the climate system. For example,\ninterannual variations in the formation of South Pacific SAMW have been linked to\ntropical teleconnections from the Indo-Pacific (Li and England 2020). In order to reduce\nthe footprint from climate modes to our conclusion, we have now changed our main\nwater-mass analyses to focus on a longer period over the Argo era 2005–2020 (an\noverlapping period of all the datasets analyzed). We have now added a clarification of\nthis point on lines 131–138 and lines 294–307: “While the globally integrated OHC of\nthe upper 2000 m and in the mode and intermediate water layers both show a robust\nand near-monotonic increase (Figs. 1–3), variations in the formation mechanisms,\nproperties, and heat content of localized mode and intermediate waters can be\ninfluenced by interannual to decadal variability (e.g. ENSO and the Indian Ocean\nDipole60 on interannual time-scales, and the IPO on decadal time-scales). However, our\nanalyses show that the main footprint of internal climate variability occurs over the\ntropical Pacific Ocean, whereas the rest of the Pacific as well as the Indian, Atlantic,\nand Southern Oceans all tend to show a consistent and robust warming, in particular\nover the subduction and well-ventilated regions near western boundary current\nextensions and north of the ACC (compare Figs. 5a-b and 6a). The near-monotonic\nwarming in mode and intermediate water layers across the Argo era (Figs. 2, 3a),\nsuggests a robust ocean heat gain forced by anthropogenic change. In contrast, OHC\nchanges within the tropical water layer and in the tropical Pacific Ocean reveal\nvariations likely associated with ENSO and/or IPO variability.”.\nText on this point has also been added at lines 343–347: “The world’s oceans hold\nby far the largest reservoir of excess anthropogenic heat in the Earth system1. Knowing\nexactly where and how much heat is stored in the ocean is thus fundamental to\nunderstanding how climate and sea-level rise will change in the future. This knowledge\nis also vital to designing the optimal observing networks for monitoring the state of the\nResponse to Reviewers\nocean and climate system.”, and lines 325–342: “Our analysis of heat accumulation\nwithin these key water masses focused on the Argo era 2005–2020, when the\nsubsurface thermal structure of the ocean is well-constrained by the Argo network.\nAlthough our results are robust within a set of observational-based datasets, important\nquestions remain, such as separating out the anthropogenic signal of OHC change from\ninternal variability, particularly in mode and intermediate waters where decadal\nvariability is poorly constrained by sparse measurements prior to the Argo era.\nIn addition, due to sparse data coverage over the ice-covered oceans, in marginal seas,\nin the deep and abyssal oceans, and in the years prior to the Argo network, not all regions\nof the ocean or time-scales are well-constrained by a suitably dense array of\nmeasurements. It is thus essential to continue establishing and improving coastal, polar,\nand deep ocean observing systems to better monitor changes of the ocean and track\nthe global energy balance more reliably. Nonetheless, our work reveals accelerating\nwarming in the upper 2000 m of the ocean over the past several decades and highlights\nthe increasing heat uptake by mode and intermediate waters, as well as the drivers\nregulating the formation of, and heat uptake in, these water masses. Regional\nobservations of subsurface temperature and salinity in the subtropical and subpolar\noceans are thus of fundamental importance to better understand the mechanisms that\nare driving ongoing warming and predict the state of climate in the future.”.\nSpecific Comments:\nI like the decomposition of Ocean Heat Content in equation 2. It follows earlier work\n(that is referred to). And in general this equation is spot on, but there have been some\npapers talking very specifically about the role of ocean heat content in sub-tropical\ngyres where the salinity changes within the water column are large enough to induce a\nthickness change and thus on density surface cause cooling signal see on the\nequatorward side of the gyres in the Indian Ocean and presumably in the pacific ocean.\nSo there should be a note that the thickness change calculation has not taken into\naccount the effects of salinity changes on thickness change on OHC in the methods\nsection.\nResponse to Reviewers\nResponse: Agreed – thanks for the suggestions. We have now added a statement of\nthis point on lines 449–455: “Another simplification is that the effects of salinity changes\non the thickness change of a water-mass layer, and thus the heat content change of the\nlayer, have not been separately considered in Eq. (2). Salinity changes can sometimes\nbe large enough to impact the ocean stratification and density structure, and thus\ninfluence the ocean circulation and layer thickness68, impacting the OHC within a given\ndensity layer. These salinity-driven impacts on ocean sequestration of heat require\nfurther investigation, but are beyond the scope of this study.”.\nLevitus Era or Levitus period (Figure 1 caption). I dont like this terminology. Syd Levitus\nshould have more fame but the period 1955-2010 is not a period of time that we should\nname after Levitus. We certainly could use the term Levitus Atlas. I would prefer simply\n“pre Argo”, or simply 1955-2010. The IPCC is wary of the ocean thermal data before the\n1970’s typically because these earlier data are unreliable (see the discussions about xbt\nand fall rates), and the trends 1960 to 1970 and 1970 to 1980 have clear variability in\nterm. We could also eliminate the “Roemmich period” as well (Figure 2 caption)\nResponse: Agreed – thanks for pointing this out. We have now updated the\nterminology for the period 1955-2004 as “pre-Argo period” and the period 2005-2020 as\n“Argo era” in the manuscript. The terminologies of “Levitus period” and “Roemmich\nperiod‘’ in the previous version of the manuscript have now all been removed.\nResponse to Reviewers\nReviewer #3 (Remarks to the Author):\nThis manuscript analyzes ocean heat content using yearly maps of ocean heat content\nfrom 0-2000 m from the World Ocean Atlas 2018 and an updated monthly maps of\ntemperature (as a function of pressure, as well as density) from Roemmich and Gilson\n(2009). It focuses on contributions of mode and intermediate water warming to ocean\nwarming as a whole. While this work is interesting, and should eventually be suitable for\npublication, there are substantial issues that preclude publication of the manuscript in its\npresent form. Some of the major issues are expanded upon below, and then a few\nspecific comments are made.\nResponse: We thank Reviewer #3 for the insightful and detailed comments that helped\nus to significantly improve the paper.\n1. First, neither the Roemmich and Gilson (2009) nor the World Ocean Atlas maps are\nimmune to the effects of changes in oceanographic sampling (e.g., Lyman and\nJohnson, 2014, J. Climate). The Roemmich and Gilson (2009) monthly maps relax back\nto the long-term seasonal climatolgy in locations and times where data are non-existent\nor sparse. This means that the first few years of maps are somewhat suspect, as are\nthe global integrals from it in those years as shown in Figures 1 and 2. Similarly, the\nWorld Ocean Atlas anomalies in the earlier years are based on very, very sparse data\n(especially prior to the late 1960s when the XBT came into widespread use), and\nanalyses of these maps should be treated with caution. This means the first decade in\nfigure 1 is very suspect, and really the increasing spatial (and depth of sampling)\ncoverage as time goes on plays a role in the shape of the curve. This should at least be\nacknowledged in the claims of dramatic acceleration. This is not to say that there is not\nacceleration, just that some of it may be owing to increases in sampling (both in terms\nof ocean regions and depths) as the decades go on. It also means that the patterns in\nthe first few years in the panels of Figure 2 are suspect.\nResponse: Thanks for pointing this out. We have now added a clarification of this point\non lines 95–107: “The community has also made considerable progress in developing\nResponse to Reviewers\ngridded products of subsurface temperature for the global ocean by merging all\nobserving systems8,10,30,31. Due to instrumental biases in bathythermograph\nmeasurements32 and insufficient data coverage prior to Argo (before 2005), as well as\ngeographic and depth limitations of conventional Argo floats, various correction\nschemes and reanalysis methods have been applied to mitigate theses biases and\nincrease our confidence in the OHC estimate. To provide a more reliable analysis of the\nglobal and regional OHC changes in this study, we present the ensemble average of\nOHC estimates derived from various ocean data products (Fig. 1; Methods). The\nuncertainty range is given by ±2 standard deviations (±2σ) to indicate the degree of\nagreement across the products. Although some discrepancies of the global OHC\nestimates are evident among the different datasets, the ensemble mean of OHC\nestimates for the upper 2000 m shows that ocean warming is not only continuing, but\naccelerating5,8,16 (Fig. 1).”, and lines 124–134: “While this is a striking acceleration of\nglobal ocean warming that warrants further investigation, uncertainties due to data\nscarcity and data quality prior to the late 1960s10 should be factored in when interpreting\nthese results. Another factor is start- and end-point sensitivity in estimating linear trends\nand capturing nonlinearity in rates of OHC change16. For example, increased spatial\nsampling coverage, particularly over the tropical and Southern Hemisphere oceans due\nto the shift in location of the observational network from a ship-based system to the\nArgo network, can affect the estimate of ocean warming during 2001–200334 (Fig. 1,\nSupplementary Fig. 1). Even though all products show robust ocean warming rates for\nthe upper 2000 m over the past few decades, reaching record warming over 2010–\n20205,16 (Table 1), these products reach better agreement after ∼2005 when the Argo\nnetwork was widely deployed.”\nTo reduce the uncertainty from only using SIO RG-Argo and Levitus data, and\nimprove the robustness of our results, we have now additionally included observational-\nbased OHC products from Cheng et al. (2017) and Ishii et al. (2017), and also included\nthe estimate of global and regional OHC using temperature and salinity observations\nfrom IAP data (Cheng and Zhu 2016; Cheng et al. 2017), and the ensemble mean of\nMet Office Hadley Centre EN4.2.2 data (Good et al., 2013), into the manuscript. These\ndatasets are constructed from Argo and other historical observations, and are based on\nResponse to Reviewers\nvarious correction approaches by the community to correct biases in XBT and MBT\ndata. We have added a detailed description of this point in Methods (lines 356–404) and\nupdated all results and figures in the revised manuscript accordingly. Now the ensemble\nmean estimate of global OHC changes and OHC changes in key water masses from\nthese products is provided throughout the manuscript, and the uncertainty range (±2σ)\nis provided in Figs. 1, 3a-d, 5c-d, 9a, 10g, S1 and Tables 1, S1–4.\n2. Second, mode and intermediate waters are defined by oceanographic features that\nare expected to change density over time. For instance, the temperature, salinity,\ndensity, and depth of the salinity minimum that defines the core of the Antarctic\nIntermediate water has evolved over the decades as the upper ocean has warmed and\nupper ocean salinity patterns have intensified over the decades (e.g., Schmidtko and\nJohnson, 2012, J. Climate). Similarly, mode waters are defined as a layer of low\nstratification. The water properties of those mode waters (including their temperature,\nsalinity, depth, or density) will change over time as upper ocean water properties\nchange over time, so their isopycnal boundaries might also change with time These\nfacts need to be acknowledged, and the interpretation of results (especially in Figure 2f,\nbut also in the maps in Figures 7 and 8) should keep account for this possibility.\nResponse: Thanks for the suggestions. We agree that the properties of mode and\nintermediate waters (including their temperature, salinity, density, PV and stratification)\nhave changed as the upper ocean has warmed and stratified over the past decades. So\nwe have now included clarifications of this point on lines 319–325: “Mode and\nintermediate waters are traditionally defined by their oceanographic features such as\ntemperature, density, salinity, and stratification. These features have evolved as the\nupper ocean has warmed and stratified over the last few decades, and as precipitation\nminus evaporation fields have changed. Here for simplicity, we use a set of fixed\nlongitude, latitude, and density or depth criteria to define these water masses and\nestimate how their OHC has changed over the Argo era, although we do separate out\nthe effects of volumetric changes from water-mass transformation in attributing total\nOHC change to each water mass.”, and in Methods (lines 491–497): “For simplicity this\nResponse to Reviewers\nstudy uses a set of fixed longitude, latitude, and density or depth criteria to define these\nwater masses and estimate how their OHC has changed over the Argo era, even\nthough the oceanographic features that define them (e.g. temperature, density, salinity\nand stratification) have also evolved as the upper ocean has warmed and stratified over\nthe Argo era. However, we do nonetheless separate out the volumetric effects from\nwater-mass transformation in attributing total OHC change to each water mass.”.\n3. Third, as the authors are aware, but do not emphasize enough in the manuscript,\nanalyses of ocean heat content in isopycnal layers can be dominated by volume changes between those layers,\nrather than temperature/salinity changes in those layers. Since the density boundaries defining these\nmode waters may change with time, one is left asking what these volume changes actually mean?\nFor instance, if the density of the actual mode water shifted across an isopyncal boundary (say it\nbecame substantially less dense in a warming world), what would a volume change for that layer\nmean? Including the volumetric change without changing the isopyncal boundary changes\nwould overemphasize the heat content change.\nResponse: Thanks for pointing this out. We agree that OHC changes of a mode water\ndensity layer can often be dominated by its volumetric changes. Volume expansion of\nmode waters suggests more oceanic heat and carbon has been injected into the ocean\ninterior, while warming is of fundamental importance as the subducted mode water can\nwarm the surrounding water masses in the ocean interior via mixing (lines 261–264). By\ndecomposing the OHC change into its volumetric and warming components (Eq. 2), we\nhave decomposed OHC changes due to volumetric changes of the density layer (e.g.,\nblue dashed lines in Fig. 3) separately from net warming of the layer (calculated using\nfixed volume based on the Argo era average density field; e.g., red dashed lines in Fig.\n3; Figs. 4, 10, and Tables S2, S4). We have now added further clarification of this point\non lines 433–438 and lines 497–505 of the revised manuscript: “By decomposing the\nOHC change into its volumetric and warming components (Eq. 2), the warming\ncomponent of OHC change within the water-mass layers and in regionally defined water\nmasses (seen in Figs. 3, 4, 10, and Supplementary Tables 2, 4) is estimated based on a\nResponse to Reviewers\nset of fixed longitudes, latitudes, and isopycnal depths (referred to as “fixed volume” or\n“volume effects removed” in this study) averaged over the Argo era, while the volumetric\ncomponent is based on a set of density ranges whose depths and volume can evolve\nover time. This decomposition allows us to distinguish OHC changes of a water mass\ndue to warming effects separate from isopycnal layer thickness (or volume) changes.”.\nThen we have rephrased statements in the manuscript to highlight the role of net\nwarming (fixed volume) in the OHC change and to avoid overemphasizing the heat\ncontent change in mode and intermediate waters, such as on lines 150–157 and 226–\n240.\n4. For the reasons discussed in the preceding two paragraphs it might actually be\nsimpler, and perhaps better, to define a set of fixed depths (they could be defined by\nisopycnal depths at a specific time, or for a decadal average) to analyse different mode\nwaters, rather than allowing those depths to follow isopycnals as they migrate with time.\nResponse: Thanks for the suggestions. Actually by decomposing the OHC change into\nits volumetric and warming components (Eq. 2), the warming component of OHC\nchange within the water-mass layers and in localized water masses (Figs. 3–4, 10, and\nTables S2, S4) can be estimated based on a set of fixed longitudes, latitudes, and\nisopycnal depths averaged over the Argo era, while the volumetric component is based\non a set of density ranges whose depth and volume can migrate with time. This\ndecomposition of the OHC change into its warming component at fixed isopycnal depths\nand volumetric component at fixed temperatures provides a more complete analysis of\nOHC changes in water masses. We have now added a statement of this point at lines\n433–438: “The terms z̅ and z̅ indicate the fixed location and depth/thickness of a\n1 2\ndensity layer over the Argo era, based on the climatologically averaged density field.\nThus the warming component of OHC changes is estimated based on fixed (Argo era\naveraged) geographic locations and depths of isopycnals, with a constant volume being\nused to integrate the warming effects of each water-mass layer and regionally defined\nwater mass over the Argo period.”, and on lines 497–505: “By decomposing the\nOHC change into its volumetric and warming components (Eq. 2), the warming component of\nResponse to Reviewers\nOHC change within the water-mass layers and in regionally defined water masses\n(seen in Figs. 3, 4, 10, and Supplementary Tables 2, 4) is estimated based on a set of\nfixed longitudes, latitudes, and isopycnal depths (referred to as “fixed volume” or\n“volume effects removed” in this study) averaged over the Argo era, while the volumetric\ncomponent is based on a set of density ranges whose depths and volume can evolve\nover time. This decomposition allows us to distinguish OHC changes of a water mass\ndue to warming effects separate from isopycnal layer thickness (or volume) changes.”.\nWe have also estimated OHC trends per decade, for the upper 2000 m of the\nocean and for the mode and intermediate water layers respectively, based on a set of\ndecadal averaged density fields from 1960 to 2020, as per the reviewer suggestion.\nThis initial test indicates that heat accumulation by the global mode and intermediate\nwater layers can account for 81% and 77% of the total heat accumulation in the upper\n2000 m of the ocean over 2000–2010 and 2010–2020, respectively. These values are\nthus close to the estimate that ~76% of ocean heat accumulation is accounted for by the\nnet warming (fixed volume) of mode and intermediate water layers over the Argo era\n(see added text on lines 150–157).\n5. It is not good practice to name time periods after people. Please find other names\nthan “The Roemmich period” or “The Levitus period”.\nResponse: Agreed – thanks for pointing this out. We have now referred to the period\n1955–2004 as “pre-Argo period” and the period 2005–2020 as “Argo era” in the\nmanuscript. The terminologies of “Levitus period” and “Roemmich period‘’ in the\nprevious version of the manuscript have now all been removed.\n6. As presented, the units of J yr-1 in Figures 3, 4, 5a, 5c-d, 6, 7, and 8 are all resolution\ndependent, and need to include area units for the maps (e.g., J m-2 yr-1) and a length\nunit for the line plots in Figure 4 (J m-1 yr-1). Also, the units of Figure 2f are also\nresolution dependent as presented, and need to be in J kg-1 m3.\nResponse to Reviewers\nResponse: Thanks for pointing this out. We have now corrected the units in Figs. 2–10\nto be resolution independent units, such as [𝐽 𝑚−2] for Figs. 5a-b, [𝐽 per degree latitude]\nfor Figs. 5c-d, and [𝐽 𝑚−2 𝑦𝑟−1] for Figs. 6a and 7, 9–10. We have also added a\nclarification to indicate that the unit of Figs. 2f and 3e-f is in [𝐽 per 0.1 𝑘𝑔 𝑚−3 density\nbin] in the caption of Figs. 2 and 3.\n7. Lines 78-89 and elsewhere: There are a number of different groups estimating ocean\nwarming as is documented in von Shuckman et al. (2020, Earth Syst. Sci. Data), where\nthe acceleration noted here is also discussed in some detail. In addition Loeb et al.\n(2021, Geophys. Res. Lett.) also document the acceleration of ocean warming 2005-\n2019, and show that it agrees with increases in net-top-of the atmosphere energy fluxes\nobserved by satellite.\nResponse: Thanks for the suggestions. We have now added a citation to von\nShuckman et al. (2020), the latest assessment report of the IPCC Working Group I,\nand Loeb et al. (2021) in the manuscript text, and also included clarifications of this point on\nlines 48–50: “Due to the accumulated excess heat in ocean basins, an acceleration of\ntotal ocean warming has become more evident from recent observational-based\nstudies5,6,15,16.”, lines 105–107: “Although some discrepancies of the global OHC\nestimates are evident among the different datasets, the ensemble mean of OHC\nestimates for the upper 2000 m shows that ocean warming is not only continuing, but\naccelerating5,8,16 (Fig. 1).”, and lines 114–123: “Further comparison with previous\nestimates of OHC trends (Table 1) and the latest assessment report of the\nIntergovernmental Panel on Climate Change (IPCC) Working Group I7 also reveals\nsignificant increases in the ocean warming rate across the measurement record5,8–\n11,16,17,19. This increasing ocean warming rate is not just limited to the upper layer of the\nocean, but also occurs at intermediate depths (700–2000 m) according to a further\nbreakdown of OHC trends by depth over the upper 2000 m5,7. Significant ocean\nwarming and accelerating OHC changes are consistent with the increase in net\nradiative energy absorbed by Earth detected in satellite observations5,33,\nsomething that is likely to continue throughout the 21st Century6,15 in the absence of\nsubstantial greenhouse gas emissions reductions.”.\nResponse to Reviewers\n8. Lines 104-108. Does this perhaps signal actually a change in the density of the mode\nwaters (becoming lighter with time in a warming ocean), rather than a change in volume\nof intermediate versus mode waters?\nResponse: Yes, warming of the mode waters can contribute to its volume increase due\nto changes in the density field. Volume of the intermediate water layer has been\ndecreasing at a rate of –0.4×1014 m3 year−1, and the mode water layer volume has\nbeen increasing at a rate of 1.9×1014 m3 year−1 over the Argo era 2005–2020 (Figs.\n10a-c). By decomposing the temporal OHC change into volumetric and warming (fixed\nvolume) components (Eq. 2), we find that the increased OHC of the mode water layer is\ndue to both net warming (fixed volume) and a volume increase of the layer, while\ndecreased OHC in the intermediate water layer is mainly due to its volume loss (Figs.\n3c-f, Table S2). Both the mode and intermediate water layers have shown a near-\ncontinuous warming (at fixed depths) that can account for about 76% of the total ocean\nwarming measured across the Argo era (Figs. 3a, c-e, Table S2). We have now\nincluded a plot of this point into the main manuscript as Fig. 3, and have rephrased\nthese sentences on lines 148–157: “These opposing signals are mainly due to the mode\nwater layer accumulating heat content via increasing its volume, while the intermediate\nwater layer loses heat content due to volume decreases (Figs. 3c-f). After factoring out\nthese volumetric effects, both the global mode and intermediate water layers have\nshown a near-continuous and monotonic warming across the Argo era (Figs. 3 and 4),\nrevealing a robust footprint of global warming penetrating into the ocean interior. Taken\ntogether, the total heat accumulated in the mode and intermediate water layers is 89%\nof the net global ocean warming during the Argo era, despite their total volume just\noccupying 62% of the upper 2000 m. The net warming of these two layers, separate\nfrom volume effects, accounts for ~76% of this global ocean warming, while ~12% is\ndue to their combined volume change (Fig. 3a, Supplementary Table 2).”.\nWe have also included a statement to clarify the uncertainty of using a fixed density\nrange across the Argo era at lines 319–325 and in Methods (lines 491–505).\nResponse to Reviewers\n9. Line 125. Here it would be appropriate to note the work of Roemmich, who noted\npreviously on that recent warming has been concentrated in the southern Ocean.\nResponse: Yes. Done – thanks.\n10. Lines 170-186. This explication (and others elsewhere which include volume\nchanges) really overemphasizes the importance of warming in the mode and\nintermediate waters. Volume changes can result in huge “heat” changes, but since heat\nchanges really should be evaluated in a constant volume, they don’t mean much and\nneed to be explained very carefully.\nResponse: Thanks for pointing this out - agreed. To separately address the role of\nwarming (fixed volume) and volumetric components in total OHC changes, we have\nnow rephrased the statements in this section (now on lines 226–238): “The total heat\ncontent of a water mass can increase both because it warms, or because the total\nvolume of that water mass increases, or both (Methods). By decomposing the OHC\nincrease over the Argo era 2005–2020, we find that the actual warming (within a fixed\nvolume) of all regionally defined mode and intermediate waters can account for 48% of\nthe global ocean warming (Fig. 9a, 10g; Supplementary Tables 3, 4). The volume of\nthese water masses also changes markedly but non-uniformly across this period, with\ntheir net volume change overall explaining a further −4% of the global OHC increase.\nThe combined OHC increase of the sum of all these regionally defined mode and\nintermediate waters, due to both warming (fixed isopycnal depth and volume) and\nvolume changes, amounts to a total contribution of 44% to the global OHC increase\n(using the values shown in Fig. 9a and Supplementary Table 4). The total volume of\nthose regionally defined water masses within the subtropical and subpolar oceans\noccupies just 24% of the upper 2000 m averaged over 2005–2020.”. Text on this point\nhas also been rephrased at lines 148–157, which can be seen in our response to\ncomment #8 above.\nResponse to Reviewers\nWe have also now added clarifications such as “fixed volume” and “volume effects\nremoved” when we refer to the contribution from warming of water masses on lines 156,\n228, 234, 252, 270, and 314, and removed the sentence on lines 170–172 of the\noriginal manuscript, to avoid overemphasizing the role of warming in mode and\nintermediate waters. The warming component of OHC change (seen in Figs. 3–4, 10,\nand Tables S2, S4) is estimated based on a set of fixed longitudes, latitudes, and\nisopycnal depths averaged over the Argo era (Eq. 2), and therefore, a set of constant\nvolumes have been used to integrate the warming effects of each water mass. We have\nadded a statement on this point at lines 433–438 and 497–505.\n__________________________________________________\nREVIEWER COMMENTS\nReviewer #1 (Remarks to the Author):\nIt is nice that the authors have addressed most of my previous comments. I have some follow-on\nsuggesfions, which are not substanfial, mainly help the authors to correct some minor errors and\nimprove the presentafion of the manuscript, so I’d classified my remaining comments as “minor\nrevisions”. I hope my suggesfions (somefimes very detailed and specific) could help the authors to\nimprove the paper and make it more powerful and impacfful.\n1) Line 87-138: this secfion can be shortened. Many of the texts within line 124-135 can be moved\ninto the Methods secfion when jusfifying your data/method.\n2) Line 90: “since the late 1970s”, should be “since the late 1960s”\n3) Line 99: “reanalysis”, should change to “objecfive analysis”\n4) Line 105: “discrepancies” befter for “differences”\n5) Line 93-95: ref.5 did not do this, you should cite ref.6 here.\n6) Line 107: here ref.8 does not provide any analysis for rate change, you should replace ref.8 by\nref.15 here, for the very first claim of warming accelerafion.\n7) Fig.1, a very important figure of this study, can be improved. 1), the terminology “ocean heat\naccumulafion” should be changed or clearly explained (please check Cheng et al. 2022 NREE paper\nfor terminology related to ocean warming, this term of “ocean heat accumulafion” is rarely used and\nnot clear). 2) what are the blue rectangles? Why the central points of the blue bars are located at the\ntop of the rectangles? These are not clear. Please make it a self- explanatory figure, which is very\nhelpful for readers to quickly get your points. 3) in panel-a, the start and end points of the three\nperiods are not apparent (one needs to read the figure capfion to know this), probably you can add\nlight (maybe different level of grey) background shadings to befter illustrate the three periods.\n8) Fig. 1b, Cheng et al. 2022 J. Climate indicates that linear trend is only valid for ~15 years segments,\nthus please consider change your segments accordingly to make the plot more meaningful.\n9) I strongly object using EN4, as explicitly stated regarding the EN4 data paper Good et al. 2013: “It\nis important to note that the analyses will relax to climatology in the absence of any observafions.\nCare must therefore be taken if using the analyses for applicafions such as idenfifying trends in\ntemperature or salinity, because a trend may be unrealisfic if analyzing periods when there were no\nobservafions.”. It is not recommended to use EN4 for trend analysis as the current study.\nFurthermore, if you look at Supplementary Fig.1, the big shift around 2000 is not physical, and the\nchanges before 2000 are too liftle to be true.\n10) Fig. 4: can the red bars be more transparent so one can see the dark blue bars behind?\n11) Line 172-183 and Fig.5, I actually don’t find these discussions needed in this paper. The\ndifference of the trends for the two periods are unclear, as you correctly stated in line 187-188 “Note\nthat trends within this relafively short period can indicate both the forcing signal due to\nanthropogenic warming as well as signals superimposed due to internal variability”, the “regional\nintensificafion” of trend for 2010-2020 compared with 2000-2010 period is not very meaningful and\nout of the scope of this study (sfick to the OHC changes in different water masses would be\ninteresfing enough, and sfick to 2005-2020 would be sufficient).\n12) Fig. 8, probably you can merge this figure with Fig.6, to make each plot more informafive.\n13) Fig. 8, are you able to befter separate SP-ESTMW and SP-WSTMW, NP-CMW and NP-STMW in the\nplot? Where are the edges of these water masses?\n14) Fig. 6, please consider moving panel-d into the SI, because it is not directly linked to the main\nflow.\n15) Fig. 9 capfion, you don’t need to repeat the full name of all water masses.\n16) Line 243: not all water masses are warming and thickening, so please consider revise the fitle.\n17) Line 308: this line probably starts a new secfion of “Summary and Discussion”?\n18) Line 337-340: this has repeated what you have said before. This paragraph needs to be\nshortened.\n19) Line 347-351: this should be shown combined of the quanfifively statements within 308-318.\n20) Lin 351-353: in previous paragraph, you already stated something similar about “befter measure\nof heat”. Please improve the organizafion of the final two paragraphs.\n21) Line 356-357: “a total of 8”: different members of EN4 should be regarded as 1 dataset.\n22) Line 363: “relafively minor”, please consider to change it to a more objecfive presentafion,\nmaybe “befter consistency over XX than YY”.\n23) Line 366-369: this choice means the global OHC esfimate will be biased because of the smaller\narea.\n24) Line 382-384: IAP salinity is “absolute salinity”, so you should not need to do conversion for this\ndataset.\n25) Line 419-421: just to check: is this following TEOS-10 standard?\n26) Line 507-512: how did you consider the reducfion in degree of freedom when calculafing least-\nsquare linear trend?\n27) Line 540-542, IAP data link is: hftp://www.ocean.iap.ac.cn/\n28) Table 1. Because you aim to show the “accelerafing global ocean warming”, probably you can re-\narrange the order of rows, starfing from longest-period trend and end up with most-recent-period\ntrend, in this way, one can nicely observe the increase in trends from up to boftom.\n29) Figure 10: Please consider to make it similar to the style of Fig. 9. At least, you can make bar plots\nat top and then spafial maps below, so one can befter compare Fig.9 and Fig.10. To further improve\nthe illustrafion, Fig.9 and 10 can be combined into a powerful/informative plot (deserve to do so for\nthis high-level journal).\n30) Fig. 3, the x-axis is not consistent between the first four panels and the last two, please make it\nconsistent, so one can befter compare upper four panels with the lowest two.\n31) Fig. 4: Actually I think you can combine Fig.4 into Fig.3, by make Fig.4 small and put it to the right\nside the Fig.3f. Actually in Fig.4, doing the stafisfics for every 0.1 kg m^-3 is not necessary, you can do\nit for 0.5 (or 0.25 at least), to simplify the plot without losing the key informafion you want to convey.\nReviewer #2 (Remarks to the Author):\nGeneral Comments.\nAs I said in the first review “I have enjoyed reading this paper. It has a wide appeal to oceanographers\nand to a wider communicafions audience. This appeal is because the manuscript documents in a\nmore thorough way than previously the changes in a key water-mass type that is found widely across\nthe global oceans. The global perspecfive of this paper means that it details a coherent set of\nchanges in both northern and southern hemispheres. “ and that “My belief is that this is a very\nrelevant and fimely paper for nature communicafions.\n“\nThe reviewers have responded to my comments very thoroughly. They have adopted an ensemble\napproach through out all of the figures, provided esfimates of error bars in figures and tables.\nFrom my perspecfive they have thoroughly addressed my original comments.\nThe only quibble I now have is the use of the word accelerafion in heat content. The\naccelerafion that is presented is often based on differing periods of fime, and trends calculated over shorter and\nlong periods will always be different. Thus this asserfion of an accelerafion can really only be made\nwhen the trend periods are presented for the same length of fime. The text should be altered to be\nclearer on this point.\nReviewer #3 (Remarks to the Author):\nThis manuscript analyzes ocean heat content using eight (slightly) different g fime-series of global\nocean heat content and the mean of three (slightly) different gridded maps for ocean temperature\nand salinity fields. It focuses on contribufions of mode and intermediate water warming to ocean\nwarming as a whole. The revision is improved from the original, but the one of main issues with the\nmanuscript brought up in the previous review remain largely unaddressed, and another is not\nadequately addressed in the revision.\n1. The potenfial impact of increases in ocean sampling on the accelerafion of ocean warming\nbetween the 1990s and 2010-2020 remains under-acknowledged. For instance, the claim of doubling\nrate over those fime periods is made in the abstract without any menfion of possible mapping\narfifacts. After a full paragraph (L105-121) on the accelerafion, only one oblique sentence (l124-126)\nbriefly menfions this issue. Linking to accelerafion of sea level rise might provide a bit more context.\n2. Second, the manuscript sfill mis-interprets the volume changes in isopycnal classes For instance,\nthe long and elaborate descripfion of possible reasons for changes in L261-293 glosses over the\nsimplest explanafion, that warmer surface condifions propagate down, and density horizons shift for\nthe various water masses. What was once a subtropical mode water density class has shifted towards\nwhat is now a subpolar mode water density class and what was once a subpolar mode water density\nclass shifts to an intermediate water density class. This interpretafion is very clear in the figures, and\nprefty clear in the previous literature on the subject. For instance, look as the quadrapoles at the\nboundaries between those water masses in Figure 2f and you can see the all the water masses\nshifting to lighter density classes. The same effect can be seen even more clearly in Figure 3f. Figures\n9d and 9e also show these changes really clearly, with the SMPW supposedly shifting southward and\neastward in the Indian and Southwest Pacific Ocean, and intermediate water supposedly\ndisappearing in the in the Southeast Pacific Ocean. It seems much more likely (and supported by\nother analyses in the literature) that the density of intermediate water has lessened over fime, as has\nthe density of SPMW, and STMW, as the surface condifions have warmed (and in some cases\nfreshened). This is reflected in the change of volumes in fixed density layers, as the water masses\nthemselves shift to lighter layers.\n3. This a more minor point, but on L218-220 and elsewhere the manuscript seem to grapple with\nhow heat builds up in the mode and intermediate waters. These waters warm mostly because\nsurface waters warm, even in winter, and thus the mode and intermediate waters warm when\nwinterfime mixing deepens the reach of that warmer layer.\nResponse to Reviewers\nRESPONSE TO REVIEWER COMMENTS:\nBelow reviewer comments have been copy-pasted in black font, and our\nresponses are indicated in blue font. Line numbers in our responses refer to the\nrevised manuscript (except where indicated).\nReviewer #1 (Remarks to the Author):\nIt is nice that the authors have addressed most of my previous comments. I have some\nfollow-on suggestions, which are not substantial, mainly help the authors to correct\nsome minor errors and improve the presentation of the manuscript, so I’d classified my\nremaining comments as “minor revisions”. I hope my suggestions (sometimes very\ndetailed and specific) could help the authors to improve the paper and make it more\npowerful and impactful.\nWe thank Reviewer #1 for this favourable assessment of our paper, and for their\ninsightful and detailed comments in this round of review. We believe these comments\nhave significantly helped us to improve the paper. Thank you.\nReviewer comments:\n1) Line 87-138: this section can be shortened. Many of the texts within line 124-135 can\nbe moved into the Methods section when justifying your data/method.\nResponse: Agreed - thanks for the suggestions. We have now rephrased the text on\nincreasing ocean warming rate in this section (lines 109–157 of the revised manuscript),\nResponse to Reviewers\nand have moved the text on lines 125–131 into the Methods section and rephrased\nthese sentences as (lines 373–379):\n“Several factors can affect the estimate of ocean warming rate and should be\nfactored in, particularly over the pre-Argo period (Fig. 1, Supplementary Fig. 1). Such\nfactors include uncertainties due to data scarcity and data quality prior to the late\n1960s10, mapping artifacts from horizontal and vertical interpolation of the sparse\nhistorical data63, the increased spatial sampling coverage over time (particularly in the\ntropical and Southern Hemisphere oceans), and the shift of the observational network\nfrom a ship-based system to the Argo network during 2001–200333.”,\nand on lines 530–533 of the revised manuscript:\n“For linear trends within a period less than 15 years, decadal variability can be part\nof the signal as well as long-term trends due to anthropogenic warming16. Another factor\nis start- and end-point sensitivity in estimating linear trends and capturing nonlinearity in\nrates of OHC change16. Overall, the monotonic increase in ocean heat uptake remains\nclear since the 1980s.”.\n2) Line 90: “since the late 1970s”, should be “since the late 1960s”\nResponse: Thanks for the suggestion. Here in the manuscript, we want to highlight the\nperiod when XBTs were widely used for observing the upper ocean. We have now\nchanged the text to “…, which have been widely used since the 1970s” because the\ndeployment of XBTs began in the mid to late 1960s, and XBTs thereafter became the\nlargest source of data for the upper ocean thermal record during the 1970s−1990s\n(Goni et al. 2019), but were only widely used from the 1970s onwards.\n3) Line 99: “reanalysis”, should change to “objective analysis”\nResponse: Done – thanks.\nResponse to Reviewers\n4) Line 105: “discrepancies” better for “differences”\nResponse: Agreed and changed. Thanks for the suggestion.\n5) Line 93-95: ref.5 did not do this, you should cite ref.6 here.\nResponse: Done – thanks.\n6) Line 107: here ref.8 does not provide any analysis for rate change, you should\nreplace ref.8 by ref.15 here, for the very first claim of warming acceleration.\nResponse: Agreed - done. Thanks\n7) Fig.1, a very important figure of this study, can be improved. 1), the terminology\n“ocean heat accumulation” should be changed or clearly explained (please check\nCheng et al. 2022 NREE paper for terminology related to ocean warming, this term of\n“ocean heat accumulation” is rarely used and not clear). 2) what are the blue rectangles?\nWhy the central points of the blue bars are located at the top of the rectangles? These are\nnot clear. Please make it a self- explanatory figure, which is very helpful for readers to\nquickly get your points. 3) in panel-a, the start and end points of the three periods are not\napparent (one needs to read the figure caption to know this), probably you can add\nlight (maybe different level of grey) background shadings to better illustrate the three\nperiods.\nResponse: Agreed - good suggestions, thanks. We have now changed the terminology\n“ocean heat accumulation” to “ocean heat uptake” in Fig. 1b and updated all related text\nResponse to Reviewers\nin the manuscript accordingly. The blue rectangle bars in Fig. 1b represent the global\nocean heat uptake for every 11-year period (e.g. 1990–2000, 2000–2010, and 2010–\n2020), which is calculated as the linear trend of annual mean OHC multiplied by the\nlength of the period; superimposed blue error bars in Fig. 1b indicate the uncertainty\nrange of global ocean heat uptake (blue rectangle bars) from various products. We have\nnow reworded text on this point in the caption of Fig. 1 to be clearer on this:\n“b, Blue rectangle bars indicate the ensemble averaged global ocean heat uptake\n(1022 ) for every 11-year period across the measurement era (Methods). Superimposed\nerror bars indicate the 2 ensemble standard deviation uncertainty range ( 2 ) of\nJ\nglobal ocean heat uptake across various datasets.”.\n± ± σ\nWe have also added delineating arrows into Fig. 1a to better indicate the start and end\nyears of the pre-Argo period, WOCE era, and Argo era, respectively, as per the\nreviewer suggestion.\n8) Fig. 1b, Cheng et al. 2022 J. Climate indicates that linear trend is only valid for ~15\nyears segments, thus please consider change your segments accordingly to make the\nplot more meaningful.\nResponse: Thanks for the suggestion. We agree that the OHC trends over each 11-\nyear period might include both decadal variability and the long-term change. Overall, the\nincrease in OHC remains clear across the measurement era, particularly since the\n1980s. We did not change the 11-year ocean heat uptake in Fig. 1b to over a 15-year\nsegment. Instead, we have added text on this point at lines 530–533:\n“For linear trends within a period less than 15 years, decadal variability can be part\nof the signal as well as long-term trends due to anthropogenic warming16. Another factor\nis start- and end-point sensitivity in estimating linear trends and capturing nonlinearity in\nrates of OHC change16. Overall, the monotonic increase in ocean heat uptake remains\nclear since the 1980s.”\nResponse to Reviewers\n9) I strongly object using EN4, as explicitly stated regarding the EN4 data paper Good\net al. 2013: “It is important to note that the analyses will relax to climatology in the\nabsence of any observations. Care must therefore be taken if using the analyses for\napplications such as identifying trends in temperature or salinity, because a trend may\nbe unrealistic if analyzing periods when there were no observations.”. It is not\nrecommended to use EN4 for trend analysis as the current study. Furthermore, if you\nlook at Supplementary Fig.1, the big shift around 2000 is not physical, and the changes\nbefore 2000 are too little to be true.\nResponse: Thanks for pointing this out. We agree that there are large inconsistencies\nwithin the four ensembles of the EN4 data over the 1990s and a large shift of global\nOHC around 2000 (Supplementary Fig.1), which gives us limited confidence in any\nanalysis from the EN4 data over the pre-Argo period. Therefore, we have now excluded\nthe EN4 data over the pre-Argo period 1955–2004 and have only included the EN4 data\nover the Argo period 2005–2020 in this study. We have now updated Fig. 1 and Fig. 5\n(Fig. 2 of the revised manuscript), and all related text in the manuscript. We have also\nadded a statement of this point in the Methods section (lines 398–401):\n“Due to the large inconsistency within the four EN4.2.2 ensembles over the 1990s,\nand due to a large shift of global OHC around 2000–2004 in the EN4.2.2 ensemble\nmean, we have only included the EN4 data over the Argo period 2005–2020 in this\nstudy.”.\n10) Fig. 4: can the red bars be more transparent so one can see the dark blue bars\nbehind?\nResponse: Agreed and changed. Thanks for the suggestion.\nResponse to Reviewers\n11) Line 172-183 and Fig.5, I actually don’t find these discussions needed in this paper.\nThe difference of the trends for the two periods are unclear, as you correctly stated in\nline 187-188 “Note that trends within this relatively short period can indicate both the\nforcing signal due to anthropogenic warming as well as signals superimposed due to\ninternal variability”, the “regional intensification” of trend for 2010-2020 compared with\n2000-2010 period is not very meaningful and out of the scope of this study (stick to the\nOHC changes in different water masses would be interesting enough, and stick to 2005-\n2020 would be sufficient).\nResponse: Thanks for pointing this out. We have now moved Fig. 5 and Fig. 8 to the\nsection “Accelerating global ocean warming” to Fig. 2 and Fig. 3 in the revised\nmanuscript, respectively. We have also moved the text on lines 172–183 to the section\n“Accelerating global ocean warming”, and reworded this text to better indicate the\nspatial variability of ocean warming in the past few decades (lines 134–150 of the\nrevised manuscript):\n“The increased ocean warming is non-uniformly distributed across ocean basins.\nOverall, in each ocean basin, an increase in OHC is observed (values indicated in Figs.\n2a, b), with stronger warming in the mid-latitude Atlantic Ocean and the Southern\nOcean compared with other basins6,11. Total warming in the Southern Ocean is\nestimated to account for ~31% of the global upper 2000-m OHC increase from 1980–\n2000 to 2000–2010 (Fig. 2a), and almost half of the global OHC increase from 2000–\n2010 to 2010–2020 (values indicated in parentheses of Fig. 2b). Hence the Southern\nOcean has seen the largest increase in heat storage over the past two decades, holding\nalmost the same excess anthropogenic heat as the Atlantic, Pacific, and Indian Oceans\nnorth of 30°S combined (Fig. 2d). The most striking warming in the Southern Ocean is\nconcentrated on the northern flank of the Antarctic Circumpolar Current, the location of\ndeep mixed layers and subduction hotspots for Subantarctic Mode Water and Antarctic\nIntermediate Water, as well as the location of subtropical mode waters formation further\nequatorward (Fig. 3). The well-ventilated regions near western boundary current\nextensions in the North Atlantic and North Pacific also reveal large warming over the\nResponse to Reviewers\npast two decades. These hotspots of ocean warming are likely linked to enhanced\nuptake, subduction, and lateral spreading of heat associated with mode and\nintermediate waters that warrants further investigation.”.\n12) Fig. 8, probably you can merge this figure with Fig.6, to make each plot more\ninformative.\nResponse: Thanks for the suggestion. We have now moved Fig. 8 to the section\n“Accelerating global ocean warming” as Fig. 3 in the revised manuscript, to better\naddress the spatial variability of ocean warming. In this case, we did not merge Fig. 8\nwith Fig. 6, in part due to the new figure location, and also to avoid overcrowding the\ndiagram.\n13) Fig. 8, are you able to better separate SP-ESTMW and SP-WSTMW, NP-CMW and\nNP-STMW in the plot? Where are the edges of these water masses?\nResponse: Thank you for the suggestion. We have now updated Fig. 8 (Fig. 3 of the\nrevised manuscript) by adding the thickness contours of all Subtropical Mode Waters,\nNP-CMW, NP-IW, NA-SPMW, SAMW, and AAIW, and into the plot. The thickness of\nNA-MMW is not included in Fig. 3 for simplicity and clarity of this plot. Superimposed\ndashed contours in the new Fig. 3 indicate the 150-m, 200-m, 250-m, 550-m, 650-m,\nand 1000-m thickness of water masses based on the density constraints in\nSupplementary Table 3, and are represented by different colours. Now the geographic\nlocations of all these water masses, including SP-ESTMW, SP-WSTMW, NP-CMW, and\nNP-STMW, can be better indicated within the plot. Text on this point has also been\nadded to the caption of Fig. 3:\n“Superimposed contours represent the 150- , 200- , and 250- thickness of all\nSTMWs and NP-CMW (blue), 550- thickness of NP-IW (red), 650- thickness of\nm m m\nm m\nResponse to Reviewers\nSAMW and NA-SPMW (green), and 1000- thickness of AAIW (yellow), to indicate\ntheir geographic locations. Note that the water-mass thickness is estimated based on\nm\nthe density constraints given in Supplementary Table 3, and for simplicity, only the core\nthickness of each water mass is shown. The geographic constraints for estimating OHC\nchange of water masses are referred to in Supplementary Table 3.”.\n14) Fig. 6, please consider moving panel-d into the SI, because it is not directly linked to\nthe main flow.\nResponse: Thanks for the suggestion. We considered moving this panel into the\nsupplementary material, but that file will be separated from the main published paper.\nWe feel that keeping it in the main manuscript would give readers a straightforward view\nof the effects of ocean thermal expansion and haline contraction in the ocean steric sea-\nlevel rise.\n15) Fig. 9 caption, you don’t need to repeat the full name of all water masses.\nResponse: Agreed and changed. Thanks for the suggestion. We have now only kept\nthe full name of water masses in the caption of Fig. 3 and only the acronym of water-\nmass names in the caption of Fig. 9.\n16) Line 243: not all water masses are warming and thickening, so please consider\nrevise the title.\nResponse: Thanks for the suggestion. We considered changing the title of this section,\nbut we feel that the original title summarizes the section in a nice succinct way as most\nof the selected subtropical mode waters and SAMW have been warming and thickening,\nwhile the warming of AAIW and NPIW has also been substantial across the Argo period.\nResponse to Reviewers\nWe therefore decided to keep the previous title for this section. Since we have reworded\nthe second paragraph of this section (lines 273-286) as per the comment from Reviewer\n#3, the revised text for this section now fits better to the previous title.\n17) Line 308: this line probably starts a new section of “Summary and Discussion”?\nResponse: Agreed and changed – thanks.\n18) Line 337-340: this has repeated what you have said before. This paragraph needs\nto be shortened.\nResponse: Thanks for raising this. We have now rephrased and shortened the\n“Summary and discussion” section, and combined the text on lines 337–340 with the\nfirst paragraph of this section (lines 319–336 of the revised manuscript) to:\n“The world’s oceans hold by far the largest excess anthropogenic heat in Earth’s\nclimate system1. Knowing exactly where and how much heat is stored in the ocean is\nthus fundamental to our ability to understand and predict future climate change and sea-\nlevel rise. In addition, this knowledge can help inform how to optimise observing\nnetworks for monitoring the state of the ocean and climate system. Our study has\ndocumented a strong acceleration in global ocean warming since the 1990s, amounting\nto >25% increase in OHC during 2010–2020 relative to 2000–2010, and nearly a\ntwofold increase during 2010–2020 relative to 1990–2000. Approximately 89% of the\nincrease in global ocean heat uptake, which includes warming and volume changes, is\nconfined to global mode and intermediate water density layers during the Argo era\n2005–2020. Mode and intermediate waters play a vital role in absorbing and\nredistributing tracers such as heat, carbon, and oxygen in the global ocean. The\nwarming of 13 of these regionally defined mode and intermediate waters at fixed\ngeographic locations and isopycnal depths (i.e. with volume effects removed) accounts\nfor ~48% of global ocean warming despite their volume occupying just 24% of the\nResponse to Reviewers\nocean. This includes Subtropical Mode Waters in both hemispheres, which have overall\nshown a striking increase in OHC dominated by increases in volume, and Subantarctic\nMode Water and Antarctic Intermediate Water in the Southern Ocean, which have\nwarmed substantially and are responsible for more than one-third of global ocean\nwarming across the 2005–2020 period.”.\n19) Line 347-351: this should be shown combined of the quantitively statements within\n308-318.\nResponse: Thanks for the suggestions. We have now combined the text on lines 347–\n351 with the previous paragraph and rephrased the text into a new paragraph (lines\n319–336 of the revised manuscript), which can be seen in our response to comment\n#18 above.\nWe have also rephrased the last paragraph of the “Summary and discussion” section to\n(lines 355–363 of the revised manuscript):\n“Our work reveals accelerated warming in the upper 2000 of the ocean over the\npast several decades and highlights increasing heat uptake by mode and intermediate\nm\nwaters, with these two water masses responsible for the majority of ocean warming over\nthe Argo era (2005-2020), despite a limited area of interaction with the atmosphere.\nExactly how this heat uptake plays out over the coming decades and beyond remains\nhighly uncertain. For example, climate change induced warming and freshening at the\nsurface is projected to stratify the upper ocean, which will reduce the overturning of\nthese water masses, in turn reducing their capacity to uptake heat. This would have\nprofound implications for the rate of future anthropogenic climate change.”.\n20) Line 351-353: in previous paragraph, you already stated something similar about\n“better measure of heat”. Please improve the organization of the final two paragraphs.\nResponse to Reviewers\nResponse: Thanks for pointing this out. We have now removed this sentence on lines\n351–353, and have combined this paragraph with the previous paragraph (lines 319–\n336 of the revised manuscript), which can be seen in our response to comment #18\nabove. We have also rephrased the last paragraph of the “Summary and discussion”\nsection (lines 355–363 of the revised manuscript), which can be seen in our response to\ncomment #19 above.\n21) Line 356-357: “a total of 8”: different members of EN4 should be regarded as 1\ndataset.\nResponse: Agreed and changed. Thanks for the suggestion.\n22) Line 363: “relatively minor”, please consider to change it to a more objective\npresentation, maybe “better consistency over XX than YY”.\nResponse: Thanks for pointing this out. We have now revised this statement to (lines\n371–373): “The estimate of global and basin-wide OHC trends among these selected\nobservational-based products were previously shown to exhibit good agreement in the\nupper 2000- ocean over the Argo era44 (Fig. 1).”.\nm\n23) Line 366-369: this choice means the global OHC estimate will be biased because of\nthe smaller area.\nResponse: Yes, agreed. Therefore we have added a statement on this point in lines\n385–387: “As such, the resultant global OHC trends from SIO RG-Argo in Fig. 1 will be\nunderestimated if these excluded polar ocean regions are also warming5.”.\nResponse to Reviewers\nIn Figs. 2–10 of the manuscript, where we show OHC changes within ocean basins and\nin mode and intermediate waters, we focus on the region between 65°S-65°N due to the\nlow sampling coverage at higher latitudes and for better consistency of the data’s spatial\ncoverage. Text on this point has now been added into the caption of Fig. 2, and can\nalso be seen in the caption of Figs. 4, 5 and Supplementary Table 1, and on lines 154–\n157 in the revised manuscript: “Thus, for the remainder of this study we focus on the\nArgo era 2005–2020, a period common to all data products analyzed. We also limit the\nwater-mass decomposition to be between 65°S and 65°N, because of sparse\nobservations in the polar regions even after the Argo array was established.”.\n24) Line 382-384: IAP salinity is “absolute salinity”, so you should not need to do\nconversion for this dataset.\nResponse: Thanks for pointing this out. We have now corrected this statement as\n“The Gibbs-Sea Water (GSW) Oceanographic Toolbox from the TEOS-10 software64,65\nis applied to convert the SIO RG-Argo data and EN4.2.2 data to Conservative\nTemperature66 and Absolute Salinity67, as well as the IAP in-situ temperature to\nConservative Temperature, ...”.\n25) Line 419-421: just to check: is this following TEOS-10 standard?\nResponse: Yes, all variables in Eq. (1) were calculated following TEOS-10 standards,\nwhere c is a constant of 3992 J kg−1 K−1, ρ is potential density in kg m−3, and θ is the\np\npotential temperature with respect to 0 °C and is measured in degrees Celsius.\nResponse to Reviewers\n26) Line 507-512: how did you consider the reduction in degree of freedom when\ncalculating least-square linear trend?\nResponse: Thanks for pointing this out. Actually, we separately calculate the linear\ntrend over each selected period (e.g., 1955–2004, 2005–2020, 2010–2020), with the\ndegree of freedom chosen as 1 for the least-squares regression. We have now removed\n“using a linear least-squares regression” to avoid any confusion and revised that\nsentence as (lines 523–527 of the revised manuscript):\n“Linear trends of OHC are calculated within each selected time period for spatially\nintegrated annual mean OHC time series, to indicate the increase or decrease in OHC\nacross each selected period. For the trend maps, the fields are first vertically integrated\nover depth or isopycnal layers, and then the linear trends over the selected period are\ncomputed for each grid point.”.\n27) Line 540-542, IAP data link is: http://www.ocean.iap.ac.cn/\nResponse: Thanks and changed.\n28) Table 1. Because you aim to show the “accelerating global ocean warming”,\nprobably you can re-arrange the order of rows, starting from longest-period trend and\nend up with most-recent-period trend, in this way, one can nicely observe the increase\nin trends from up to bottom.\nResponse: Yes. Done – good suggestion, thanks.\n29) Figure 10: Please consider to make it similar to the style of Fig. 9. At least, you can\nmake bar plots at top and then spatial maps below, so one can better compare Fig.9\nResponse to Reviewers\nand Fig.10. To further improve the illustration, Fig.9 and 10 can be combined into a\npowerful/informative plot (deserve to do so for this high-level journal).\nResponse: Agreed - thanks for the suggestions. We have now updated Fig. 9 and\nmade Figs. 9 and 10 in the same style, i.e., both with spatial maps at the top and bar\nplots at the bottom. We also carefully considered combining Fig. 9 with Fig. 10,\nhowever, it came out too crowded for the combined new figure. So we have decided to\nkeep them as two separate plots in the revised manuscript.\n30) Fig. 3, the x-axis is not consistent between the first four panels and the last two,\nplease make it consistent, so one can better compare upper four panels with the lowest\ntwo.\nResponse: Agreed and changed. Thanks for the suggestion.\n31) Fig. 4: Actually I think you can combine Fig.4 into Fig.3, by make Fig.4 small and\nput it to the right side the Fig.3f. Actually in Fig.4, doing the statistics for every 0.1 kg\nm^-3 is not necessary, you can do it for 0.5 (or 0.25 at least), to simplify the plot without\nlosing the key information you want to convey.\nResponse: Thanks for the suggestions. We looked at combining these two figures,\nhowever, it came out too crowded for the combined new figure. In this way, we have\ndecided to keep Figs. 3 and 4 as two separate plots (Figs. 5 and 6 in the revised\nmanuscript), and Fig. 4 (now Fig. 6) shown for every 0.1 kg m^-3 density bin for\nconsistency with Figs. 2 and 3 (now Figs. 4 and 5).\nResponse to Reviewers\nReviewer #2 (Remarks to the Author):\nGeneral Comments.\nAs I said in the first review “I have enjoyed reading this paper. It has a wide appeal to\noceanographers and to a wider communications audience. This appeal is because the\nmanuscript documents in a more thorough way than previously the changes in a key\nwater-mass type that is found widely across the global oceans. The global perspective\nof this paper means that it details a coherent set of changes in both northern and\nsouthern hemispheres. “ and that “My belief is that this is a very relevant and timely\npaper for nature communications.“\nThe reviewers have responded to my comments very thoroughly. They have adopted an\nensemble approach through out all of the figures, provided estimates of error bars in\nfigures and tables.\nFrom my perspective they have thoroughly addressed my original comments.\nWe thank Reviewer #2 for this favourable assessment of our paper, and for their\ninsightful comments in these two rounds of review. We believe these comments have\nsignificantly helped us to improve the paper. Thank you.\nThe only quibble I now have is the use of the word acceleration in heat content. The\nacceleration that is presented is often based on differing periods of time, and trends\ncalculated over shorter and long periods will always be different. Thus this assertion of\nan acceleration can really only be made when the trend periods are presented for the\nsame length of time. The text should be altered to be clearer on this point.\nResponse to Reviewers\nResponse: Agreed - thanks for pointing this out. To address this, when indicating the\nincrease in ocean warming over different time periods in Table 1 and over the pre-Argo\nand Argo periods in Fig. 1a, we have converted the units to J/year to take into account\nthe different lengths of the different periods. The quantification and description of the\nincreased ocean warming rate over the past three decades in Fig. 1b are based on the\nsame length of periods over 1990–2000, 2000–2010, and 2010–2020. Therefore, the\nterm “acceleration” is appropriate when describing the increase in ocean warming. We\nhave now adjusted the statements on lines 109–133 and rearranged Table 1 to better\naddress this point.\nWe further agree that the choice of different lengths of time periods and start- and end-\npoints can also affect the estimate of linear trends due to the nonlinearity in ocean\nwarming and internal climate variability. For linear trends within a shorter record,\ninterannual and decadal variability of the climate system can have a larger influence.\nWe have now added text on this point at lines 530–533:\n“For linear trends within a period less than 15 years, decadal variability can be part\nof the signal as well as long-term trends due to anthropogenic warming16. Another factor\nis start- and end-point sensitivity in estimating linear trends and capturing nonlinearity in\nrates of OHC change16. Overall, the monotonic increase in ocean heat uptake remains\nclear since the 1980s.”.\nResponse to Reviewers\nReviewer #3 (Remarks to the Author):\nThis manuscript analyzes ocean heat content using eight (slightly) different g time-\nseries of global ocean heat content and the mean of three (slightly) different gridded\nmaps for ocean temperature and salinity fields. It focuses on contributions of mode and\nintermediate water warming to ocean warming as a whole. The revision is improved\nfrom the original, but the one of main issues with the manuscript brought up in the\nprevious review remain largely unaddressed, and another is not adequately addressed\nin the revision.\nWe thank and appreciate Reviewer #3 for the extensive and detailed comments in these\ntwo rounds of review, which have helped us to significantly improve the paper.\n1. The potential impact of increases in ocean sampling on the acceleration of ocean\nwarming between the 1990s and 2010-2020 remains under-acknowledged. For\ninstance, the claim of doubling rate over those time periods is made in the abstract\nwithout any mention of possible mapping artifacts. After a full paragraph (L105-121) on\nthe acceleration, only one oblique sentence (l124-126) briefly mentions this issue.\nLinking to acceleration of sea level rise might provide a bit more context.\nResponse: Thanks for pointing this out. We agree that the increase in ocean sampling\nand the shift of the observing network from a ship-based system to the Argo network in\nthe early 2000s can potentially influence the estimate of ocean warming rate in the past\nthree decades. We have now included a statement on this point to better acknowledge\nthese possible uncertainties, and also added text on the point of acceleration in global\nmean sea-level rise at lines 123–133:\nResponse to Reviewers\n“Note that there has been both increased ocean sampling and a shift of the\nobservational network from a ship-based system to the Argo network since the initiation\nof the global Argo array (2001–2003)33. This may impact the estimated increase in\nglobal ocean warming over the past three decades (Fig. 1). However, the rate of global\nmean sea-level rise has also been increasing since 1993 based on an independent\nestimate from satellite altimeter34, providing confidence in our results given that sea\nsurface height increases have been dominated by temperature effects since altimeter\nmeasurements began. Significant ocean warming and accelerating OHC changes are\nalso consistent with the increase in net radiative energy absorbed by Earth detected in\nsatellite observations5,35, something that is likely to continue throughout the 21st\nCentury6,15 in the absence of substantial greenhouse gas emissions reductions.”,\nand lines 373–379 of the Methods section:\n“Several factors can affect the estimate of ocean warming rate and should be\nfactored in, particularly over the pre-Argo period (Fig. 1, Supplementary Fig. 1). Such\nfactors include uncertainties due to data scarcity and data quality prior to the late\n1960s10, mapping artifacts from horizontal and vertical interpolation of the sparse\nhistorical data63, the increased spatial sampling coverage over time (particularly in the\ntropical and Southern Hemisphere oceans), and the shift of the observational network\nfrom a ship-based system to the Argo network during 2001–200333.”.\n2. Second, the manuscript still mis-interprets the volume changes in isopycnal classes\nFor instance, the long and elaborate description of possible reasons for changes in\nL261-293 glosses over the simplest explanation, that warmer surface conditions\npropagate down, and density horizons shift for the various water masses. What was\nonce a subtropical mode water density class has shifted towards what is now a subpolar\nmode water density class and what was once a subpolar mode water density class\nshifts to an intermediate water density class. This interpretation is very clear in the\nfigures, and pretty clear in the previous literature on the subject. For instance, look as\nthe quadrapoles at the boundaries between those water masses in Figure 2f and you\ncan see the all the water masses shifting to lighter density classes. The same effect can\nResponse to Reviewers\nbe seen even more clearly in Figure 3f. Figures 9d and 9e also show these changes\nreally clearly, with the SMPW supposedly shifting southward and eastward in the Indian\nand Southwest Pacific Ocean, and intermediate water supposedly\ndisappearing in the in the Southeast Pacific Ocean. It seems much more likely (and supported by\nother analyses in the literature) that the density of intermediate water has lessened over time,\nas has the density of SPMW, and STMW, as the surface conditions have warmed (and\nin some cases freshened). This is reflected in the change of volumes in fixed density\nlayers, as the water masses themselves shift to lighter layers.\nResponse: Thanks for pointing this out - agreed. We have now rephrased the text on\nlines 261–278 (now lines 273–286 of the revised manuscript), to address this point:\n“Among the mode and intermediate waters defined in this study, we find that OHC\ntrends of mode waters are dominated by their volumetric changes. Instead warming\n(fixed volume) exhibits a larger signal in the heat content trend for intermediate waters,\nincluding both Antarctic and North Pacific Intermediate Water (Figs. 3, 6 and 10;\nSupplementary Table 4). This is likely because while the warming of water masses\npropagates downward, the density-defined boundary of water masses shifts both\nvertically and laterally as a response to ocean warming (Supplementary Fig. 2). In\nparticular, the subtropical oceans and the Southern Ocean reveal a significant\ndeepening and poleward shift of mode water isopycnals. This indicates that the\nSubtropical/Subantarctic Mode Water has been expanding and shifting towards\nlocations that once were classified as Subpolar Mode Water and/or intermediate water.\nNote that although the mode water has shifted, this does not necessarily indicate a net\nwarming of the water mass. Instead, the intermediate water layer could potentially be\nlightened and warmed due to mixing with the warmed and deepened mode waters over\nthe subtropical oceans and the Southern Ocean18 (Fig. 10). In reality, a combination of\nthese two effects has likely played out.”.\nResponse to Reviewers\nWe have also added a plot of linear trends in zonally averaged temperature and density\nover the Argo era as Supplementary Fig. 2, to better indicate the vertical and lateral shift\nof density classes due to warming over the subtropical and Southern Oceans.\n3. This a more minor point, but on L218-220 and elsewhere the manuscript seem to\ngrapple with how heat builds up in the mode and intermediate waters. These waters\nwarm mostly because surface waters warm, even in winter, and thus the mode and\nintermediate waters warm when wintertime mixing deepens the reach of that warmer\nlayer.\nResponse: Thanks for the suggestion and apologies this was not clearer in the revised\nmanuscript. Text on this point has now been added (lines 226–230) to the revised\nmanuscript:\n“Therefore, a component of this accelerated OHC increase must be due to\nincreased heat uptake by mode and intermediate waters via wintertime ventilation, in\ncombination with deepening of the warmed water by mixing. In this way warming from\nthe surface can propagate into the ocean interior and spread laterally and vertically via\ncirculation and mixing with surrounding water masses.”.\n______________________________________________________________________\nREVIEWER COMMENTS\nReviewer #1 (Remarks to the Author):\nThe paper is ready to publish, with only one very minor language edifing suggesfion: Line 101, please\nchange “insufficient” to “limited” or something similar but less subjecfive.\nReviewer #3 (Remarks to the Author):\nThis revised manuscript is improved over the previous version. However, there are sfill a few points\nthat probably should be addressed by minor revision prior to publicafion. They follow, indexed by\nLine number (L).\n1. L46-47. Is the phrase \"and an increase of warming rates at both intermediate depths of 700–2000\n47 m, and in the deeper ocean\" supported by the references that follow? Is there a primary\nreference for warming below 2000 m amongst them?\n2. L104-106. The gridded products all use largely the same database, with variafions in quality\ncontrol, instrument bias correcfions, and mapping techniques, correlafions length and fime scales,\netc. So how is the ensemble average \"more reliable\"? Perhaps the most useful quanfity to emerge is\nreally the standard deviafion for the ensemble.\n3. L126-130 and 207-214. What about the addifion of water from land (mostly ice melt)? Doesn't\nthat also accelerate over the fime period in quesfion? How much of the accelerafion in sea level rise\nis due to accelerafion in movement of water from land to the ocean?\n4. L279-282. These statements may be true in terms of stafic water mass definifions based on fixed\ndensity (or temperature, or salinity) ranges. However, when a mode water is defined as a low PV\nregion of the water column at a certain locafion (for instance just north of the SubAntarcfic Front for\nthe SubAntarcfic Mode Water), or an intermediate water as being centered around a feature (such as\na salinity minimum north of the Antarcfic Circumpolar Current for Antarcfic Intermediate Water),\nwhat we see is that the density, temperature, and perhaps salinity of these water masses have\nshifted. So phenomenologically, the conclusions here seem misleading, and may confuse some\nreaders. They need to be revised.\n5. L432-435. Subsurface temperature and salinity changes are not solely governed by venfilafion fime\nscales. Isopycnal heave, owing to changes communicated by planetary (e.g., Rossby and Kelvin)\nwaves (e.g., Masada et al., 2010, Science, doi: 10.1126/science.1188703) can also effect these\nchanges.\nRESPONSE TO REVIEWER COMMENTS:\nBelow reviewer comments have been copy-pasted in black font, and our\nresponses are indicated in blue font. Line numbers in our responses refer to\nthe revised manuscript (except where indicated).\nReviewer #1 (Remarks to the Author):\nThe paper is ready to publish, with only one very minor language editing suggestion:\nLine 101, please change “insufficient” to “limited” or something similar but less\nsubjective.\nResponse: We have now changed the word “insufficient” to “limited”, as per the\nreviewer's suggestion. Thank you. And thanks to Reviewer #1 for their detailed and\ninsightful comments across the previous two rounds of review, which have helped to\nimprove the paper.\nReviewer #3 (Remarks to the Author):\nThis revised manuscript is improved over the previous version. However, there are\nstill a few points that probably should be addressed by minor revision prior to\npublication. They follow, indexed by Line number (L).\nResponse: We thank Reviewer #3 for their detailed and insightful comments in these\nthree rounds of review. We believe these comments have significantly helped us to\nimprove the paper.\n1. L46-47. Is the phrase \"and an increase of warming rates at both intermediate\ndepths of 700–2000 m, and in the deeper ocean\" supported by the references that\nfollow? Is there a primary reference for warming below 2000 m amongst them?\nResponse: Thanks for pointing this out. Yes, references 5, 7, and 9–11 all support\nthe statement of warming at intermediate depths of the ocean, and references 5, 7,\nand 10 further support the statement of warming in the deep ocean below 2000 m.\nWe have now separately cited these supporting references for each statement, and\nhave relocated the previous reference 43 (Bagnell and DeVries 2021) to this\nsentence, now cited as new reference 11, at lines 46–47 in the manuscript text:\n“..., and an increase of warming rates at both intermediate depths of 700–2000\nm5,7,9–12, and in the deeper ocean below 2000 m5,7,10,11.”.\n2. L104-106. The gridded products all use largely the same database, with variations\nin quality control, instrument bias corrections, and mapping techniques, correlations\nlength and time scales, etc. So how is the ensemble average \"more reliable\"?\nPerhaps the most useful quantity to emerge is really the standard deviation for the\nensemble.\nResponse: Agreed - thanks. We have now changed the word “reliable” to “robust” to\nmake this statement less subjective. A statement regarding the standard deviation of\nthe ensemble follows in the subsequent sentence (lines 106–108):\n“The uncertainty of this estimate is taken to be ±2 standard deviations (±2σ), to\nindicate the degree of agreement across the products.”.\n3. L126-130 and 207-214. What about the addition of water from land (mostly ice\nmelt)? Doesn't that also accelerate over the time period in question? How much of\nthe acceleration in sea level rise is due to acceleration in movement of water from\nland to the ocean?\nResponse: Thanks for pointing this out. The ice melt from glaciers, the Greenland\nand Antarctic ice sheets, and any land water storage changes has been estimated to\naccount for 19.4%, 15.2%, 8.6%, and 10.9% of global sea-level rise over the period\n1993–2018, respectively, while the remaining 45.9% is due to thermal expansion\neffects of the ocean (Table 9.5 of the latest assessment report of IPCC Working\nGroup I, ref. 1). Both the total sea-level rise and the thermosteric sea-level rise have\nbeen increasing across the altimeter measurement era 1993–2018 (compare the rate\nof sea-level rise over 1993–2018 and 2006–2018 in Table 9.5 of ref. 1, and see Fig.\n3 of Cazenave and Moreira (2022, ref. 42)). We have now revised the sentence on\nlines 126–130 as:\n“This may impact the estimated increase in global ocean warming over the past\nthree decades (Fig. 1). However, the rate of global mean sea-level rise has also been\nincreasing since 1993 based on an independent estimate from satellite altimeter\ndata1,35, providing confidence in our results given that half of the global sea surface\nheight increase is due to thermal expansion of the ocean since altimeter\nmeasurements began1.”,\nand have included clarifications of the ocean mass component of sea-level rise on\nlines 207–209:\n“An independent measure of OHC changes can be found in altimeter\nmeasurements of sea-level if salinity adjustments remain small and the ocean mass\ncomponent of sea-level change is known and accounted for.”,\nand on lines 215–218:\n“The mass component of sea-level change, due to land ice melt (including glaciers,\nand the Greenland and Antarctic ice sheets) and land water storage change, has\nalso significantly increased across the altimeter era, but display a more uniform\npattern of sea-level rise compared to the thermal expansion component41,42.”.\n4. L279-282. These statements may be true in terms of static water mass definitions\nbased on fixed density (or temperature, or salinity) ranges. However, when a mode\nwater is defined as a low PV region of the water column at a certain location (for\ninstance just north of the SubAntarctic Front for the SubAntarctic Mode Water), or an\nintermediate water as being centered around a feature (such as a salinity minimum\nnorth of the Antarctic Circumpolar Current for Antarctic Intermediate Water), what we\nsee is that the density, temperature, and perhaps salinity of these water masses have\nshifted. So phenomenologically, the conclusions here seem misleading, and may\nconfuse some readers. They need to be revised.\nResponse: Thanks for the suggestions.To better clarify what we can see from either\nthe density-based definition or the stratification- or salinity-based definition, we have\nnow included a statement of this point on lines 291–296:\n“In addition to density-based constraints, alternative criteria like stratification or\nsalinity within specific locations are often used to define mode and intermediate\nwaters26. These criteria also allow examination of changes in the density and\ntemperature of water masses. Within this framework, the poleward and downward\nmigration of isopycnals over the subtropical oceans and the Southern Ocean can\nserve as an indication of the warming and lightening of mode and intermediate\nwaters (Supplementary Fig. 2).”.\n5. L432-435. Subsurface temperature and salinity changes are not solely governed\nby ventilation time scales. Isopycnal heave, owing to changes communicated by\nplanetary (e.g., Rossby and Kelvin) waves (e.g., Masada et al., 2010, Science, doi:\n10.1126/science.1188703) can also effect these changes.\nResponse: Thanks for pointing this out. We have now included a citation to Li et al.\n(2023), Gunn et al. (2023), Purkey and Johnson (2010), and Masuda et al. (2010),\nand have added a clarification of this point in the manuscript text, lines 442–446:\n“Waters deeper than 2000 m are generally not ventilated rapidly enough to have\nexperienced the last 50 years of global warming, but note that both the abyssal\noverturning slowdown72,73 and oceanic internal waves74 can result in bottom ocean\nwarming over these time-scales, as seen in observations75.”."}]

Decision letter