# Robust acceleration of Earth system heating observed over the past six decades

Audrey Minière, Karina von Schuckmann, Jean-Baptiste Sallée, Linus Vogt

<figure><img src="image_1.png" /></figure>

## To cite this version:

Audrey Minière, Karina von Schuckmann, Jean-Baptiste Sallée, Linus Vogt. Robust acceleration of Earth system heating observed over the past six decades. Scientific Reports, 2023, 13 (1), pp.22975. 10.1038/s41598-023-49353-1 . hal-04454855

HAL Id: hal-04454855

https://hal.science/hal-04454855v1

Submitted on 15 Feb 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

Check for updates

OPEN

# Robust acceleration of Earth system heating observed over the past six decades

Audrey Minière ID<sup>1,2✉</sup>, Karina von Schuckmann ID<sup>2</sup>, Jean-Baptiste Sallée ID<sup>3</sup> & Linus Vogt ID<sup>3</sup>

Global heating of the Earth system is unequivocal. However, detecting an acceleration of Earth heating has remained elusive to date, despite suggestive evidence of a potential increase in heating rates. In this study, we demonstrate that since 1960, the warming of the world ocean has accelerated at a relatively consistent pace of 0.15 ± 0.05 (W/m²)/decade, while the land, cryosphere, and atmosphere have exhibited an accelerated pace of 0.013 ± 0.003 (W/m²)/decade. This has led to a substantial increase in ocean warming, with a magnitude of 0.91 ± 0.80 W/m² between the decades 1960–1970 and 2010–2020, which overlies substantial decadal-scale variability in ocean warming of up to 0.6 W/m². Our findings withstand a wide range of sensitivity analyses and are consistent across different observation-based datasets. The long-term acceleration of Earth warming aligns qualitatively with the rise in CO₂ concentrations and the decline in aerosol concentration during the same period, but further investigations are necessary to properly attribute these changes.

In the past 150 years, Earth's climate has been warming at a rate that is unprecedented in at least the last 2000 years<fnref n="1" />. This human-caused warming has caused widespread adverse impacts and related losses and damages to nature and people, which will continue in the future as global climate continues to warm². Detecting changes in the rate of warming is crucial for informed decision-making in international climate negotiations, with the aim of limiting global warming to specific levels. However, it remains a significant challenge to detect such changes due to the substantial internal variability of the climate system on a decadal scale (e.g., ref. ³). In this paper, we address this challenge by examining the global heat accumulation rate across the entire climate system, including the ocean, atmosphere, cryosphere, and land. By focusing on this integrated view, rather than solely relying on changes in global mean surface temperature, we can mitigate the impact of variability and gain a more comprehensive understanding⁴,⁵.

Global heat accumulation in the climate system, resulting from the current positive Earth's Energy Imbalance (EEI) at the top of the atmosphere, is primarily dominated by changes in Global Ocean Heat Content (GOHC)⁴. GOHC changes account for approximately 90% of the total heat increase in the past fifty years, while land heating, ice melting, and atmospheric warming contribute around 5%, 3%, and 1% respectively⁶⁻⁸. Several studies have indicated an increase in the global heat accumulation rate in recent decades, with values rising from 0.50 [0.32 to 0.69] W/m² during the period 1971–2006 to 0.79 [0.52 to 1.06] W/m² (90% confidence interval) for the period 2006–2018 (ref. ⁴,⁶⁻²² and Fig. 1). Some studies have even suggested a potential doubling of EEI in the last decade compared to the previous one⁶,¹⁷.

Despite this suggestive body of evidence, no study has conducted an analysis of heat accumulation acceleration since 1960 to date. While the results presented in Fig. 1 provide insights, they represent accumulation rates computed over varying time spans, with higher rates calculated over decadal periods and lower rates calculated over multi-decadal periods. This variation in time spans makes it challenging to make definitive and quantified statements about acceleration. Additionally, the use of diverse datasets and methodologies can significantly impact the calculated accumulation rates. The only notable climate variable where acceleration has previously been detected in past decades is Global Mean Sea Level (GMSL)²³⁻³¹. This GMSL acceleration has been attributed to factors such as increasing GOHC leading to thermal expansion of seawater, declining land water storage, or increasing land ice melt²⁵⁻²⁷,³¹,³².

In this paper, we present the first observation-based quantification of the acceleration of Earth's system heating. Our study adopts a systematic approach, incorporating multiple datasets and employing various methods. We estimate the rate of change and acceleration of Earth's heat content using a collection of GOHC time series

<footnote n="1">¹Université Toulouse III - Paul Sabatier, Toulouse, France. ²Mercator Ocean International, Toulouse, France. ³LOCEAN-IPSL, Laboratoire d'Océanographie et du Climat: Expérimentation et Approches Numériques, CNRS/IRD/MNHN, Sorbonne Université, Paris, France. ✉email: aminiere@mercator-ocean.fr</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Assessment of observation-based Earth's energy imbalance (EEI) absolute values as available in the literature, considering various approaches and time periods (refer to Table S1 for references). The black circles highlight the EEI values derived from an international assessment conducted within the Global Climate Observing System⁷ (GCOS) framework. Notably, these values show a significant agreement with the EEI values estimated in the latest IPCC report⁴, which are represented by the purple circles.</figcaption></figure>

derived from in-situ temperature data spanning from 1960 to 2020. Additionally, we utilise ocean reanalyses data from 2005 to 2020 and satellite altimetry and gravimetry data from 2002 to 2020. To complement our analysis, we also incorporate non-ocean heat content time series covering the period from 1960 to 2020, and compare our findings to observations of the net radiative flux at the top of the atmosphere (TOA) spanning from 2001 to 2020.

## Heat content rate of change

Several research groups have developed four-dimensional global ocean temperature datasets, which enable the estimation of GOHC. In this study, we use an ensemble of ten freely accessible products (refer to Table S2 for the exhaustive list and associated references) and take a systematic approach to assess their consistency and discrepancies (see “Methods” section). We compare this ensemble mean to three other estimates of GOHC based on the ensemble mean of three ocean reanalyses, one satellite-derived estimate²⁰, and a composite ensemble⁷ of sixteen products developed within the framework of the Global Climate Observing System (referred to as GCOS heat content). It should be noted that not all products cover the same time period (see Fig. S1 and Table S2), and we have aimed to maximise the number of products used for each discussed time period throughout this paper.

None of these GOHC estimates can be considered flawless. The process of producing these estimates is accompanied by significant challenges stemming from observational gaps, historical changes in observational coverage, and potential sensor errors³³⁻³⁵. Consequently, each research group must make important assumptions regarding data quality control, data correction, and strategies for filling spatio-temporal gaps. These assumptions collectively contribute to the uncertainty associated with the reconstruction of the GOHC³⁶⁻³⁸. Unfortunately, producers do not always provide the GOHC uncertainty associated with their methodological choices (see Table S2), commonly referred to as internal GOHC uncertainty³⁹. Alternatively, one can compute an a posteriori estimate of GOHC uncertainty by determining the structural uncertainty, which is informed by the ensemble spread of a set of products³⁹⁻⁴². In this study, we aim to investigate how these different estimates of uncertainty, as well as our chosen statistical methodology for inferring time-derivatives, can impact the computation of GOHC rates of change.

The IAP product⁴³ (see Table S2) is one of the few that provides an estimation of its internal uncertainty. We took this opportunity to compare the internal uncertainty estimate of the IAP product with the structural uncertainty of the GCOS heat content⁶,⁷. Additionally, we calculated our own structural uncertainty based on two times the standard deviation of our set of products (see Methods). Both structural uncertainties encompass the IAP internal uncertainty in the GOHC anomaly time series (Fig. 2a). When propagated to determine the uncertainty of the GOHC rate of change, the structural uncertainty also provides the largest uncertainty estimate
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 2.** Global ocean heat content (GOHC) rate of change and its corresponding uncertainty estimates. (a) Comparison between structural and internal GOHC uncertainties. The GOHC from the IAP product⁴³ is represented by the plain blue line, surrounded by its associated internal uncertainty, depicted by blue shading. The dashed blue lines indicate the structural uncertainties derived from this study, while the red dashed lines represent the structural uncertainties from the GCOS product⁷. (b) Sensitivity test on GOHC rate of change uncertainties calculations. The bars display two times the WLS regression standard errors using the internal IAP uncertainty (shown in blue) and the structural GOHC uncertainty from this study (indicated by dashed blue), as well as the GCOS product (depicted by dashed red). (c) Sensitivity test on GOHC trends computation for different decades. Four methods are evaluated on the GOHC time series for the IAP product: WLS (medium blue), OLS (purple), LOWESS (red), QUADRATIC (yellow). For each method, the trend uncertainties accounting for sample autocorrelation are also shown (dashed lines). Further information on these methods can be found in the Methods section. All uncertainties are shown at the 95% confidence level.</figcaption></figure>

(Fig. 2b). Consequently, for the remainder of this study, we will employ the GOHC structural uncertainty in our calculations of GOHC rates of change. We note that our ensemble spread is only indicative of the true structural uncertainty, and might only provide an underestimate of the true uncertainty because we use an ensemble of opportunity. Furthermore, we tested four different methods for computing GOHC rates of change (see “Meth-ods” section). Although the choice of method has minimal impact on the computed rate of change itself, it does influence the associated uncertainty (Fig. 2c). Among the tested methods, the Weighted Least Squares regression (WLS) suggests the largest uncertainty. As a precautionary measure, we have selected to utilise this methodology for the remainder of this paper.

Regardless of the methodological choice (Fig. 2), the specific product used, and the time period considered within the past sixty years, all GOHC time series utilised in this study (Fig. 3), along with their associated rates of change (Fig. 4 and Fig. S2), exhibit consistency within their respective uncertainty ranges. This remarkable consistency instils high confidence in the finding that the global ocean has experienced a warming rate of
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 3.** Time evolution of global ocean heat content (GOHC) anomalies. This study's in-situ estimates of GOHC are represented by the blue dashed curve for individual products and the bold blue curve for the ensemble mean. These estimates are compared to the GOHC from GCOS⁷ (red curve), as well as reanalyses (pink curve) and satellite estimates²⁰ (MOHeaCAN; orange curve). The shadings on the graph indicate the structural GOHC uncertainties from this study (blue shading), GCOS (red shading), and reanalyses (pink shading). Additionally, the orange shading represents the internal uncertainty corresponding to the satellite-based GOHC estimate²⁰. All uncertainties are depicted at the 95% confidence level. The anomalies are presented relative to a baseline period of 2005–2020 (refer to the Methods section for detailed information on the GOHC processing). Refer to Table S2 for product references and additional details.</figcaption></figure>

0.66 ± 0.23 W/m², 0.69 ± 0.27 W/m², or 0.70 ± 0.24 W/m² during the period of 2006–2020, as indicated by the ensemble constructed in this study, the GCOS ensemble, or the ensemble of ocean reanalysis, respectively (Fig. 4) (all uncertainties mentioned here and hereafter are given at the 95% confidence level, and take into account serial autocorrelation and model-generated internal climate variability uncertainty, see Methods, Eq. 12). Only the indirect satellite-derived GOHC estimate suggests a slightly higher rate of change during the period 2006–2020, reaching a value of 0.87 ± 0.35 W/m² (Fig. 4). Nonetheless, all of these warming rates for the period 2006–2020 are greater than rates computed over longer time periods, particularly surpassing the rates for the period of 1993–2020, which stand at 0.61 ± 0.22 W/m² (or 0.61 ± 0.23 W/m² as estimated by GCOS), and significantly exceeding the rates for the period of 1971–2020, which amount to 0.46 ± 0.23 W/m² (or 0.48 ± 0.25 W/m² as estimated by GCOS).

Although the rates of heat content change for the land, atmosphere, and cryosphere are an order of magnitude smaller than the rates for the ocean, all components exhibit higher rates when focusing on more recent decades (Fig. 4). Importantly, this increased warming rate in shorter and more recent time periods is observed to occur at a comparable pace in terms of the percentage of increase across the different Earth system components. Compared to the period of 1971–2020, the in-situ GOHC rate was higher by 34 ± 18% (or 28 ± 24% based on the GCOS estimate) in the period of 1993–2020 and by 42 ± 16% (or 44 ± 34% based on the GCOS estimate) in the period of 2006–2020. The rate of change in heat content for the land component increased by a similar proportion, with a
<figure><img src="image_10.png" /><figcaption>**Figure 4.** Multi-decadal heating rates for the different components of the Earth System. The heating rates are provided for three main periods of the observing system: historical (1971–2020), satellite (1993–2020) and Argo era (2006–2020). This study's in-situ estimates of ocean warming rates (blue dots) are compared to GCOS ocean heating rates⁷ (red dots). The non-ocean heating rates (green dots) are computed from GCOS heat content time series⁷, and are the sum of heating rates for atmosphere (purple), land (yellow) and cryosphere (gray). The heating rates and their associated uncertainties are computed using WLS regression and are relative to the Earth's surface at the top-of-atmosphere (as described in the “Methods”). The uncertainties are displayed at the 95% confidence level and take into account serial autocorrelation (see “Methods”) (solid lines). For the ocean heating rates, the uncertainties taking into account model-generated internal variability are also displayed (dashed lines).</figcaption></figure>

percentage increase of 38 ± 78% and 51 ± 142% for the respective periods. The percentage increase in the rate of change for the cryosphere and atmosphere is also comparable, albeit slightly larger, at 60 ± 22% and 117 ± 49% for the atmosphere, and 56 ± 96% and 64 ± 184% for the cryosphere (see also ref. ⁷).

The increase in the rate of GOHC as we focus on more recent periods aligns with the wide range of estimates from individual published studies and international literature assessments (Fig. 1). However, the precise time evolution of this increase and the impact of comparing periods of different lengths, potentially affected by different processes, remain less clear. To address this, we calculate the heat content rates over consistent 10-year periods using a moving window spanning from 1960 to 2020 (Fig. 5). Despite a large error range in the early years of the time-series and significant decadal variability of up to 0.6 W/m², we observe a clear and steady low-frequency increase in the decadal GOHC rate from 1960 to 2020. The decadal GOHC rate has been consistently rising since the 1960s, with an increase of +0.91 ± 0.80 W/m² between the first decade (1960–1970) and the last decade (2010–2020) (average across this study and GCOS estimates, Fig. 5).

## Heat content acceleration

Over the past 20 years, in addition to products based on ocean in-situ measurements, we have gained access to other sources of evidence, such as ocean reanalysis and satellite data, which provide estimates of ocean heat content and energy flux at the top of the atmosphere⁴,⁴⁵ (Fig. 5). The increase in the Earth's decadal heat content rate estimated from these mostly independent sources of observations is notably consistent and exhibits a clear upward trend (Figs. 4, 5). This leads us to the question: can we formally detect an acceleration of the total Earth's heat content using the existing global climate observing system? To address this question, we employ a WLS
<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 5.** 1960–2020 time evolution of decadal heating rates of the Earth. The in-situ estimates of ocean heating rates are represented for each individual product (blue dashed curve, refer to Table S2 for references) and for the ensemble mean (bold blue curve). These estimates are compared to GCOS⁷ (bold red curve), reanalyses (pink curve) and satellite²⁰ (MOHeaCAN; orange curve) ocean heating rates. Non-ocean heating rates (green curve) are computed from GCOS heat content time series⁷, and equal to the sum of atmosphere, land and cryosphere heating rates. The 10-year means of the top-of-atmosphere (TOA) net radiative flux (black curve) are anchored on the 2005–2020 Earth Heat Inventory (EHI) trend of 0.75 ± 0.24 W/m² (refer to the Methods section for detailed information on TOA net flux anchoring). Heating rates and associated uncertainties are computed from WLS regression and are relative to the Earth's surface at the top-of-atmosphere (as described in the “Methods” section). The uncertainties take into account autocorrelation and are shown at the 95% confidence level for all estimates (see “Methods”) except TOA net radiative flux, where uncertainties of ± 0.1 W/m² have been reported⁴¹. The black triangles indicate major volcanic eruptions that have occurred since 1960. In the context of the heating rates, a positive value indicates that the Earth system is experiencing warming, while a negative value is associated with a cooling.</figcaption></figure>

methodology to estimate acceleration, consistent with our computation of the rate of change (see “Methods”). We calculate acceleration using a quadratic fit to the annual GOHC estimates shown in Fig. 3. We apply this approach over five multi-decadal periods (1960–2020, 1970–2020, 1980–2020, 1990–2020 and 2002–2020), to the ensemble of GOHC constructed in this study and to the GCOS ensemble, resulting in two acceleration estimates and associated uncertainties for each time period. All of these estimates robustly indicate a significant acceleration of GOHC since 1960, with an average rate of 0.15 ± 0.05 (W/m²)/decade. This rate is obtained from the average of the two accelerations and associated uncertainty estimated from this study and the GCOS ensembles while accounting for autocorrelation and the model-generated internal climate variability (see “Methods”, Eq. 12). Importantly, this GOHC acceleration remains remarkably consistent when computed over different multi-decadal time periods (Fig. 6). While not entirely independent, the consistency of the estimated acceleration when using different time periods and observation-based ensembles enhances our confidence in the robustness
<figure><img src="image_13.png" /><figcaption>**Figure 6.** Earth system constant heating acceleration. The acceleration rates are estimated from a quadratic in time WLS regression of the annual GOHC time-series (see “Methods” section). The in-situ ocean estimates from this study (blue bars) are compared to GCOS ocean estimates⁷ (red bars), and satellite ocean estimates²⁰ (MOHeaCAN; orange bars). The non-ocean components (green bars) are estimated by summing the atmosphere, land and cryosphere GCOS⁷ heat content time series. The top-of-atmosphere (TOA) estimates of warming acceleration (bars with black contours) are computed from a linear OLS regression accounting for autocorrelation¹⁷,⁴⁷, over global TOA net radiative flux (white bar) and near-global TOA net radiative flux (which excludes latitudes higher than 60°, represented by the hatched white bar). Uncertainties for all estimates are shown at the 95% confidence level and take into account serial autocorrelation (solid lines) (see “Methods”). For the ocean heat content accelerations, the uncertainties taking into account model-generated internal variability are also displayed (dashed lines). A positive value indicates that the heat content is accelerating, while a negative value suggests deceleration.</figcaption></figure>

of this multi-decadal GOHC acceleration. Similarly, acceleration of non-ocean heat content is also significantly detected over these multi-decadal time periods at a rate of 0.013 ± 0.003 (W/m²)/decade (green bars in Fig. 6).

The shorter the time period, the more sensitive the quadratic fit is to noise. As a result, the uncertainties associated with the computed acceleration over the past two decades (2002–2020) are much larger compared to those computed over a longer timespan (Fig. 6). Interestingly, however, the computed GOHC accelerations over these two decades remains significant (i.e., the lower bounds of the confidence intervals are above zero), and is notably larger than when computed over a longer timespan. The 2002–2020 in-situ GOHC acceleration is estimated at 0.30 ± 0.28 (W/m²)/decade (same as above: average across both this study and GCOS estimates, and accounting for autocorrelation and the model-generated internal climate variability uncertainty based on Eq. 12 in “Methods” section), double the value compared to the 1960–2020 in-situ GOHC acceleration estimate. This substantial in-situ GOHC acceleration estimate in 2002–2020 is supported by two independent estimates: one based on satellite-derived GOHC (0.42 ± 0.36 (W/m²)/decade) and the other on satellite-based energy flux at TOA (0.42 ± 0.21 (W/m²)/decade). However, our results also indicate that the acceleration estimates over two decades, in contrast to estimates over longer periods, are sensitive to methodological choices. This sensitivity prevents us from drawing firm conclusions regarding increased rate of acceleration over the past two decades (Figs. S4b, c, S5b, c, S6).

Discussion

Several recent studies have examined changes in GOHC rates over time¹⁶,¹⁷,²⁰,⁴⁶⁻⁴⁸. However, these studies have focused only on the most recent two decades and have not considered changes over a longer period of more than half a century, as we have done in this paper. Using different observational products and methodologies, these studies have estimated an increase in the energy flux at TOA at a rate of 0.42 ± 0.23 (W/m²)/decade for the period 2000–2020 (ref. ⁴⁷), or a rate of 0.38 ± 0.24 (W/m²)/decade for the period 2001–2020 (ref. ⁴⁶), and an increase in GOHC rates of 0.43 ± 0.40 (W/m²)/decade for the period 2005–2019 (ref. ¹⁷). These previous estimates align well with our multi-product estimate of the acceleration of in-situ GOHC from 2002 to 2020, which is 0.30 ± 0.28 (W/
m²)/decade. However, it is worth noting that due to the relatively short time span, the acceleration estimate and its associated uncertainty are sensitive to methodological choices. In contrast, our methodology in this study allows us to present compelling evidence that the acceleration of heat content in the Earth system began in the 1960s. GOHC has been steadily accelerating at a rate of approximately 0.15 ± 0.05 (W/m²)/decade since then, while other components of the climate system have been accelerating at a rate of 0.013 ± 0.003 (W/m²)/decade.

At the multidecadal scale, our findings indicate that the acceleration of Earth's heat content since the 1960s is robust to methodological choices and has remained relatively stable over a span of forty years. This provides empirical evidence supporting the notion that the acceleration is a result of long-term changes in the climate system. The observed multidecadal acceleration in heat content accumulation in the Earth system is qualitatively consistent with the documented likely increase in the rate of total anthropogenic effective radiative forcing since the 1970s, as estimated in the most recent IPCC report⁴⁹. This increase is primarily attributed to the growing concentrations of CO₂ and the declining concentrations of aerosols⁴⁹.

The rate at which heat accumulates in the Earth system is influenced by three main factors: radiative forcing, physical climate feedback, and land or sea surface temperature, which can modulate the intensity of feedback processes. Radiative forcing encompasses both natural factors (such as solar radiation and volcanic activity) and human-induced factors (such as greenhouse gas and aerosols emissions). The internal variability of the climate system can also affect this heat budget by influencing land or sea surface temperatures⁴⁶,⁵⁰. The extent to which the observed increase in Earth's heat content rates over the past two decades can be attributed to internal variability or forced by human activities has been a subject of debate¹⁷,⁴⁶. The in-situ ocean observation products presented in this study show variability in heat accumulation at a decadal scale, reaching up to 0.6 W/m², but the causality behind these variations remains unclear. The role of internal variability⁵⁰,⁵¹, changes in anthropogenic forcing¹⁷,⁴⁶, and the presence of uncertainties or undetected biases in the observing system³⁶,⁵²,⁵³ in explaining these changes still require further investigation.

In the past two decades, there is some indication that the acceleration of GOHC has increased compared to the previous sixty years, although this finding is sensitive to methodological choices. Raghuraman et al. (2021)¹⁷,⁴⁶, using climate model experiments, suggested that it is highly unlikely for the observed 2001–2020 trend in TOA net radiative flux to be solely explained by internal variability. However, they did propose that internal variability could contribute up to ±0.19 (W/m²)/decade over a 20-year period. Therefore, the possible increased acceleration over the past twenty years may result from the combination of internal variability in the climate system superimposed on the lower-frequency acceleration induced by human activities since the 1960s. However, there are also alternative possible explanations that suggest that the acceleration of the past twenty years could be linked to factors such as a significant rise in radiative forcing due to decreased aerosol concentration⁴⁸,⁵⁴⁻⁵⁶, and changes in clouds and sea-ice leading to increased warming due to change in climate feedback¹⁷,⁵⁷. Combined, these effects could induce a recent increase in rate of acceleration. To better quantify, understand and attribute this potential recent increase in acceleration, further investigation and quantification are needed.

Although consistent within their uncertainty ranges, it is worth noting a noticeable difference of 0.12 (W/m²)/decade between the central estimate of acceleration derived from in-situ ocean observation products (0.30 ± 0.28 (W/m²)/decade) and those obtained from remote sensing, such as satellite-derived ocean heat content (0.42 ± 0.36 (W/m²)/decade) or satellite-based energy net flux at the top of the atmosphere (0.42 ± 0.21 (W/m²)/decade). We should interpret this difference cautiously, considering that it is smaller than the uncertainty associated with each individual estimate. Nonetheless, it does raise questions and emphasises the distinction among the various products used in this study. An important factor contributing to the difference between these products is their spatial coverage. While the satellite-based energy flux estimate encompasses the entire globe, satellite-based ocean heat content excludes latitudes greater than 66°, and ocean in-situ products exclude regions beyond latitudes greater than 60° and the ocean below 2000 m depth. Consequently, we need to consider the potential impact of these different spatial coverages on our results. When we remove high latitudes poleward of 60° from our estimate of energy flux at TOA, the acceleration estimate is reduced by approximately 15%, bringing it closer to the ocean in-situ estimate. However, it is worth noting that satellite-based ocean heat content, which also excludes polar regions, still produces an acceleration estimate consistent with the energy flux at TOA. Therefore, a more plausible factor contributing to the difference may lie in the contribution of the deep ocean below 2000 m. Due to the lack of in-situ ocean coverage below 2000 m, we are unable to quantify the acceleration in this part of the ocean globally. However, one study has reported acceleration of deep ocean warming below 2000 m in the South Pacific Ocean⁵⁸. Additionally, Bagnell and DeVries (2021)¹⁸ attempted to reconstruct global deep ocean temperature change over the past century and demonstrated a significant increase in the rate of deep ocean warming from the 1990s-2000s, following a cooling phase that may have delayed the acceleration of full-depth GOHC.

Our findings are based on a comprehensive set of products derived from complex datasets that have inherent limitations in their coverage of vast ocean areas. Dealing with errors and uncertainties presents a significant challenge, especially when detecting trends and acceleration³. In this study, we have addressed these challenges by testing our results using various approaches to represent and propagate uncertainties in trends and acceleration (see Supplementary Information). Furthermore, we have included a diverse range of products, each employing different methodologies to construct their datasets. Through extensive sensitivity analyses, our results have consistently shown robustness, thereby increasing confidence in their validity. However, it is important to acknowledge that uncertainties persist, particularly concerning the limitations of the observing system during the early years of the analysed period³⁶,⁵⁹.

We reveal a previously uncharted acceleration of the GOHC since the 1960s. These observations serve as crucial indicators of climate change and play a vital role in enhancing our understanding of the Earth's response to human activities. In addition to quantifying the acceleration, our findings highlight the consistent insights provided by the current global climate observing system into past changes in Earth's heat content. It is imperative to prioritise the maintenance and improvement of the global climate observing system to ensure its continued
effectiveness in monitoring climate change in the future⁷. Furthermore, expanding the coverage of the global climate observing system to currently undersampled ocean regions and addressing data gaps in non-ocean components⁶,⁷ would enable more refined analyses of acceleration in ocean warming and reduce uncertainties in detecting and attributing global climate change.

## Methods

### In-situ GOHC timeseries and structural uncertainty

The global ocean heat content (GOHC) constitutes the major pillar (~90%) of the Earth Heat Inventory⁶,⁷, and currently, its rates of change provide the most accurate estimate of the absolute value of the Earth's energy imbalance.

When the four-dimensional gridded temperature datasets were available (see Table S2), we computed the GOHC (in Joules) at each month *t* by integrating the temperature *T* between 0 and 2000 m over the global ocean surface, as follows:

$$GOHC(t) = \rho * C * \sum_{x} \sum_{y} \sum_{z} T(t, x, y, z) * h(x, y, z) * A(x, y), \quad (1)$$

where *h* is the layer thickness in metres, *ρ* = 1030 kg/m³, the reference water density, *C* = 3980 J/°C/kg, the heat capacity of water, and *A*, the grid cell area at longitude and latitude (*x*,*y*) in m². We used a common mask for all gridded products, i.e., the most restrictive mask between all products (here the IPRC product is the one which has the smaller ocean domain), after masking the polar regions (i.e., poleward of 60° latitude) and shallow ocean areas (i.e., where bathymetry is less than 300 m). To account for deep ocean contribution (i.e., ocean below 2000 m), we added a linear trend on GOHC time series of 0.97 ± 0.48 ZJ/year (0.06 ± 0.3 W/m²) from 1992 to 2020 (ref. ⁶,⁷). We then computed annual averages of GOHC and calculated the GOHC anomalies relative to the 2005–2020 mean. The structural uncertainty is then computed by taking twice the ensemble standard deviation of GOHC anomalies. Note that in cases where structural uncertainty displays long-term variance (e.g., surpassing a 15-year period), the act of baselining each product over a 16-year period (2005–2020) may artificially reduce the inter-product spread over the baseline period.

### TOA net flux anchoring

The net radiative flux at the top of the atmosphere provides one of the most accurate estimates of the time evolution of the Earth's Energy Imbalance to-date, which can be determined to within 0.3 (W/m²)/decade (ref. ⁶⁰). However, its absolute value is more uncertain. For example, uncertainty resulting from calibration alone is 2 W/m² (ref. ⁶¹). There are also other sources of uncertainties associated with radiance-to-flux conversion and time interpolation (~0.2 W/m² for each)⁶⁰⁻⁶², or in assuming a 20 km reference level (0.1 W/m²)⁶³. Currently, the net imbalance from the standard CERES data products is ~4.3 W/m² (ref. ⁴⁴) which is much larger than the expected EEI to be 0.5-1W/m² (ref. ⁵). Therefore, to overcome this issue in its absolute value, the TOA net flux is commonly adjusted to be consistent with an estimate from ocean in-situ temperature change¹²,¹⁷. We chose to offset the TOA net radiative flux time-series to match the Earth Heat Inventory rate of change over 2005–2020 estimated from the GCOS ensemble at 0.75 ± 0.24 W/m² (ref. ⁷), such that TOA net radiative flux mean value over the 2005–2020 period is equal with the GCOS trend value. Applying this offset allows us to plot the TOA net radiative flux time-series on the same axis as other estimates, and has no implication on the calculation of TOA net flux trends.

### Heat content trend evaluation

There are many ways of estimating trends in a time series in the field of climate research (see for example the ref. ⁶⁴). We focus here on the most classical and often used techniques for estimating trends in the field of GOHC research. In Fig. 2, we tested the sensitivity of GOHC rates of change and its uncertainties to four methods which can be grouped into two types of calculations. One is based on a delta approach³⁹, and another one is based on a linear least squares approach (e.g., ref. ³⁷,⁶⁵⁻⁶⁸). For the uncertainty on the trend analysis, we have to consider four different aspects, i.e., statistical uncertainty (computed from the variance of the fit residuals), data uncertainty (structural or internal), uncertainty for serial autocorrelation, and uncertainty from internal climate variability. Both delta and least square approaches are applied on annual heat content timeseries and are described below, together with relevant uncertainty considerations.

### Delta approach

For the delta approach, the change in heat content series *y*(*t*) over a specific period, Δ*y*, is calculated by subtracting the first value from the last value over a specific period. We then computed the linear trend, *y*<sub>l</sub>, over the same period by dividing the change (in Joules) by the length of the period (in seconds). This method is widely used in the literature for estimating GOHC linear trends (e.g., ref. ⁴,¹¹,⁶⁹). To reduce the effect of high-frequency variability, data noise or changes in the observing system, before computing the trend, we first smoothed the time series *y*(*t*), using a weighted scatterplot smoothing approach⁷,⁵⁹ (named LOWESS in Fig. 2), or a quadratic fit (named QUADRATIC in Fig. 2). For the LOWESS fit, we use a span width of 25 years following Cheng et al. (2022) (ref. ⁷,⁵⁹).

To obtain an uncertainty range on our estimate of the rate of change, and to take into account the sensitivity of the calculation to interannual variability, we implemented a Monte-Carlo bootstrap to generate 1000 surrogates' series, *y*<sub>random</sub>(*t*), under the assumption of a given mean (our fitted time series, *y*<sub>fit</sub>(*t*))<sup>7,59</sup>. Each surrogate *y*<sub>random</sub>(*t*) consists of the fitted time series *y*<sub>fit</sub>(*t*) plus a randomly generated residual which follows a normal (Gaussian) distribution of standard deviation equal to the uncertainty associated to the time series *y*(*t*). We have then
repeated this approach, now randomly generating surrogates, $y_{random\_autocorr}(t)$, following a Gaussian distribution with a standard deviation equal to the residual variance, and including an autocorrelated component structured according to an AR(1) process (see Fig. 2b of ref. ⁵⁹). For both cases, the surrogates are smoothed with a LOW-ESS or QUADRATIC fit, and the trends are estimated from it, noted respectively $y_{random, t}$ and $y_{random\_autocorr, t}$. To associate an uncertainty to the trend, $y_t$, we compute the standard deviation of all the 1,000 trends $y_{random, t}$ noted as $\sigma_{random}$, and the standard deviation of all the 1,000 trends $y_{random\_autocorr, t}$ noted as $\sigma_{random\_autocorr}$. The 95% confidence interval for the linear trend, $y_t$, is then calculated as follows:

$$CI = \pm 2 * \sqrt{\sigma_{random}^2 + \sigma_{random\_autocorr}^2} \qquad (2)$$

With this method, we consider the data uncertainty, and we account for serial autocorrelation.

## Least squares approach

The Ordinary Least Square (OLS) approach is a classical method for estimating trends in key climate variables such as global mean surface temperature (GMST), global mean sea level (GMSL) or GOHC (e.g., ref. ⁸,³⁷). The standard error of OLS regression can be adjusted to consider the serial autocorrelation that can be very strong in the time series of climate variables, such as GMST⁷⁰⁻⁷² or radiative fluxes at the top of the atmosphere¹⁷,⁴⁷. However, the OLS regression does not consider the uncertainty associated with the variable for which we aim to estimate the trend, which is why some studies use other methods such as the Weighted Least Square (WLS) regression (e.g., ref. ³⁷,⁶⁶⁻⁶⁸). In this study, we consider these types of regressions together, with the aim of choosing the most suitable method for our case study in terms of estimating uncertainties and trends.

We regressed the equation $y = \beta t + \varepsilon$, where $t$ is the observed quantity (here the GOHC series), $t$ is the time vector, $\beta$ are unknown regression coefficients and $\varepsilon$ are the associated errors which are assumed to be Gaussian with mean zero. The regression was performed using either an ordinary least squares (named OLS in Fig. 2), or a weighted least squares (named WLS in Fig. 2) regressions.

The OLS regression estimates $\hat{\beta}_{OLS}$ and their associated variances are given by the following equations (ref. ⁶⁵):

$$\hat{\beta}_{OLS} = (X^T X)^{-1} X^T Y, \qquad (3)$$

and:

$$var(\hat{\beta}_{OLS}) = \hat{\sigma}^2 (X^T X)^{-1}, \qquad (4)$$

where $X$ is the design matrix (with ones in the first column and time values in the second column), $Y = y^T$. In Eq. (4), $\hat{\sigma}$ is the standard error of the regression, computed as:

$$\hat{\sigma}^2 = \frac{1}{N-2} \times \sum_{i=1}^{N} e_i^2, \qquad (5)$$

where $N$ is the sample size, and $e$ are the residuals of the regression.

To take into account the data uncertainty in the regression, we use a WLS fit, which estimates $\hat{\beta}_{WLS}$ and their associated variances, as follows (ref. ⁶⁵,⁶⁷):

$$\hat{\beta}_{WLS} = (X'^T X')^{-1} X'^T Y', \qquad (6)$$

and:

$$Var(\hat{\beta}_{WLS}) = (X^T WX)^{-1}, \qquad (7)$$

with:

$$W = diag(\frac{1}{w_{ii}^2}), \qquad (8)$$

where $W$ is a weighting matrix in which $w_{ii}$ are chosen to be the uncertainties associated to $y$ (for example the structural uncertainty of in-situ GOHC time series), $Y' = W^{\frac{1}{2}} Y$ and $X' = W^{\frac{1}{2}} X$. The advantage of WLS over OLS lies in its capacity to account for data uncertainty in the regression parameters and their associated standard errors. Note that for GOHC trend and acceleration evaluation, we only take into account the structural/internal uncertainty in the weighted matrix, and we do not consider the uncertainty raising for internal climate variability (denoted by $\sigma_{ICV}$). Assuming that $\sigma_{ICV}$ is constant over time: if $\sigma_{ICV}$ significantly surpasses the structural/ internal uncertainty, the $w_{ii}$ values would remain constant, and trend (or acceleration) values would closely resemble those from an OLS regression. Our evaluation shows that both OLS and WLS methods yield similar results within the ranges of uncertainty (Fig. S4 and S5). We adopt a conservative approach by preferring WLS, which is associated with uncertainties larger than OLS.

Both OLS and WLS methods assume white noise for fit residuals, i.e. they suppose no serial autocorrelation. However, we know that annual ocean heat content series exhibit autocorrelation (ref. ⁵⁹). We show that the residual autocorrelation structure of the heat content series is well captured—or even overestimated—by an autoregressive AR(1) process (see Supplementary Information Fig. S7). We then follow the methodology of Santer et al. (2008) (ref. ⁷¹) to compute a correction factor in order to adjust the standard error of the OLS and WLS parameters (computed in Eqs. 4 and 7 respectively). This correction factor is obtained via:
$$v = \frac{1 + \rho}{1 - \rho}, \qquad (9)$$

with $\rho$ the lag-1 (for 1-year) temporal autocorrelation coefficient of the regression residuals. Additionally, we test the influence of significant volcanic eruptions on the autocorrelation pattern of residuals. We show that the lag-1 coefficient for both the GCOS dataset and the dataset in our study exhibits minimal changes when discarding years after significant volcanic eruptions (as depicted in Fig. S7 of the Supplementary Information). As a result, we decide to retain the entire GOHC time series without excluding years following major volcanic eruptions. The corrected standard error, $\sigma_c$, is then computed as:

$$\sigma_c = \sigma \sqrt{v}, \qquad (10)$$

To compute the TOA net radiative flux trends, we use an OLS fit and adjust the standard error based on Eq. (9), as in Loeb et al. (2021) (ref. 17).

The 95% confidence interval for the trend is calculated based on $\pm 2$ times the corrected standard error of the regression as following:

$$CI = \pm 2\sigma_c, \qquad (11)$$

## Heat content acceleration evaluation

To compute heat content accelerations shown in Fig. 6, we regress a quadratic fit from the yearly heat content time series using a second-order WLS regression (i.e. we add a third column including a quadratic term $t^2$ in the design matrix X of Eqs. 5 and 6), using the GOHC structural uncertainty as weighting matrix. As for the linear least squares approach, the 95% confidence interval for the acceleration is calculated based on $\pm 2$ times the corrected standard error of the regression (Eq. 10). We also discuss another approach to detect acceleration on Supplementary Information (see Discussion S3 and Figs. S4, S5, S6, S10).

## Sensitivity of GOHC trend and acceleration estimates to internal variability

To assess the contribution of internal variability to trend values (Figs. 4 and 5) and accelerations (Fig. 6), we use pre-industrial simulations of GOHC from an ensemble of climate models of the CMIP6 exercise (46 models, each with one realisation), following the methodology outlined by Raghuraman et al. (2021, 2023) (refs. 46,73). For each of the X-length periods of interest, we randomly selected 3,000 segments of GOHC of length X from the model ensemble. From these segments, we calculated trends and accelerations using linear and quadratic fits (see Figs. S8 and S9 in the Supplementary Information). Subsequently, we computed the standard deviation of the 3000 trend and acceleration values, denoted as $\sigma_{ICV}$ (see Figs. S8 and S9 in the Supplementary Information). We use $\sigma_{ICV}$ as a proxy to quantify the portion of uncertainty in trends and acceleration introduced by internal variability, which we incorporate into confidence intervals as follows:

$$CI_{ICV} = 2 * \sqrt{\sigma_c^2 + \sigma_{ICV}^2}, \qquad (12)$$

where $\sigma_c$ represents the adjusted standard error of the linear and quadratic fit as described in Eq. (11). This new confidence interval is depicted by the dashed error bars in Figs. 4, 5, and 6.

## Reference surface

To ensure consistency with TOA net radiative flux estimate, all the heat content trend and acceleration values are given relative to the Earth's surface at the top-of-atmosphere, $S_{TOA}$, computed as follows: $S_{TOA} = 4\pi(R_T + z_{TOA})^2$ with $R_T$ the Earth's radius equals 6371 km, and $z_{TOA}$ the altitude of the top-of-atmosphere equals 20 km$^{20,63}$.

## Data availability

All datasets used in this study are freely available and can be downloaded from websites listed in Table S2.

Received: 13 June 2023; Accepted: 7 December 2023

Published online: 27 December 2023

## References

1. Summary for policymakers. In *Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change* (eds Masson-Delmotte, V. et al.) 3–32 https://doi.org/10.1017/9781009157896.001 (Cambridge University Press, 2021).

2. Pörtner, H.-O., Roberts, D.C., Poloczanska, E.S., Mintenbeck, K., Tignor, M., Alegría, A., Craig, M., Langsdorf, S., Löschke, S., Möller, V. & Okem A. (eds) *IPCC AR6 Working Group 2: Summary for Policymakers*. https://www.ipcc.ch/report/ar6/wg2/chapt er/summary-for-policymakers/.

3. Richardson, M. T. Prospects for detecting accelerated global warming. *Geophys. Res. Lett.* **49**, e2021GL095782 (2022).

4. Forster, P. et al. The Earth's energy budget, climate feedbacks, and climate sensitivity. In *Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change* (eds. Masson-Delmotte, V. et al.). 923–1054 https://doi.org/10.1017/9781009157896.001 (Cambridge University Press, 2021).

5. von Schuckmann, K. et al. An imperative to monitor Earth's energy imbalance. *Nat. Clim. Change* **6**, 138–144 (2016).

6. von Schuckmann, K. et al. Heat stored in the Earth system: Where does the energy go?. *Earth Syst. Sci. Data* **12**, 2013–2041 (2020).

7. von Schuckmann, K. et al. Heat stored in the Earth system 1960–2020: Where does the energy go?. *Earth Syst. Sci. Data* **15**, 1675–1709 (2023).
8. IPCC. *Climate Change 2021: The Physical Science Basis*. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change (2021) (in press).

9. Hansen, J. et al. Earth's energy imbalance: Confirmation and implications. *Science* **308**, 1431–1435 (2005).

10. IPCC. *Climate Change 2013: The Physical Science Basis*. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change. Vol. 1535 https://doi.org/10.1017/CBO9781107415324 (2013).

11. Trenberth, K. E., Fasullo, J. T. & Balmaseda, M. A. Earth's energy imbalance. *J. Clim.* **27**, 3129–3144 (2014).

12. Johnson, G. C., Lyman, J. M. & Loeb, N. G. Improving estimates of Earth's energy imbalance. *Nat. Clim. Change* **6**, 639–640 (2016).

13. Cheng, L., Abraham, J., Hausfather, Z. & Trenberth, K. E. How fast are the oceans warming?. *Science* **363**, 128–129 (2019).

14. Meyssignac, B. et al. Measuring global ocean heat content to estimate the earth energy imbalance. *Front. Mar. Sci.* **6**, 133 (2019).

15. Johnson, G. C. & Lyman, J. M. Warming trends increasingly dominate global ocean. *Nat. Clim. Change* **10**, 757–761 (2020).

16. Hakuba, M. Z., Frederikse, T. & Landerer, F. W. Earth's energy imbalance from the ocean perspective (2005–2019). *Geophys. Res. Lett.* **48**, e2021GL093624 (2021).

17. Loeb, N. G. et al. Satellite and ocean data reveal marked increase in Earth's heating rate. *Geophys. Res. Lett.* **48**, e2021GL093047 (2021).

18. Bagnell, A. & DeVries, T. 20th century cooling of the deep ocean contributed to delayed acceleration of Earth's energy imbalance. *Nat. Commun.* **12**, 4604 (2021).

19. Trenberth, K. E. & Cheng, L. A perspective on climate change from Earth's energy imbalance. *Environ. Res. Clim. 1*, 013001 (2022).

20. Marti, F. et al. Monitoring the ocean heat content change and the Earth energy imbalance from space altimetry and space gravimetry. *Earth Syst. Sci. Data* **14**, 229–249 (2022).

21. Su, H., Wei, Y., Lu, W., Yan, X.-H. & Zhang, H. Unabated global ocean warming revealed by ocean heat content from remote sensing reconstruction. *Remote Sens.* **15**, 566 (2023).

22. Cheng, L. et al. Past and future ocean warming. *Nat. Rev. Earth Environ.* **3**, 776–794 (2022).

23. Church, J. A. & White, N. J. A 20th century acceleration in global sea-level rise. *Geophys. Res. Lett.* **33**, 73 (2006).

24. Jevrejeva, S., Moore, J. C., Grinsted, A. & Woodworth, P. L. Recent global sea level acceleration started over 200 years ago?. *Geophys. Res. Lett.* **35**, 65 (2008).

25. Yi, S., Heki, K. & Qian, A. Acceleration in the global mean sea level rise: 2005–2015. *Geophys. Res. Lett.* **44**, 11905–11913 (2017).

26. Nerem, R. S. et al. Climate-change-driven accelerated sea-level rise detected in the altimeter era. *Proc. Natl. Acad. Sci.* **115**, 2022–2025 (2018).

27. Dangendorf, S. et al. Persistent acceleration in global sea-level rise since the 1960s. *Nat. Clim. Change* **9**, 705–710 (2019).

28. Ablain, M. et al. Uncertainty in satellite estimates of global mean sea-level changes, trend and acceleration. *Earth Syst. Sci. Data* **11**, 1189–1202 (2019).

29. Veng, T. & Andersen, O. B. Consolidating sea level acceleration estimates from satellite altimetry. *Adv. Space Res.* **68**, 496–503 (2021).

30. Guérou, A. et al. Current observed global mean sea level rise and acceleration estimated from satellite altimetry and the associated measurement uncertainty. *Ocean Sci.* **19**, 431–451 (2023).

31. Frederikse, T. et al. The causes of sea-level rise since 1900. *Nature* **584**, 393–397 (2020).

32. Chen, X. et al. The increasing rate of global mean sea-level rise during 1993–2014. *Nat. Clim. Change* **7**, 492–495 (2017).

33. Abraham, J. P. et al. A review of global ocean temperature observations: Implications for ocean heat content estimates and climate change. *Rev. Geophys.* **51**, 450–483 (2013).

34. Tan, Z. et al. A new automatic quality control system for ocean profile observations and impact on ocean warming estimate. *Deep Sea Res. Part Oceanogr. Res. Pap.* **194**, 103961 (2023).

35. Gouretski, V., Cheng, L. & Boyer, T. On the consistency of the bottle and CTD profile data. *J. Atmos. Ocean. Technol.* **39**, 1869–1887 (2022).

36. Boyer, T. et al. Sensitivity of global upper-ocean heat content estimates to mapping methods, XBT bias corrections, and baseline climatologies. *J. Clim.* **29**, 4817–4842 (2016).

37. Savita, A. et al. Quantifying spread in spatiotemporal changes of upper-ocean heat content estimates: An internationally coordinated comparison. *J. Clim.* **35**, 851–875 (2022).

38. Allison, L. C. et al. Towards quantifying uncertainty in ocean heat content changes using synthetic profiles. *Environ. Res. Lett.* **14**, 084037 (2019).

39. Palmer, M. D., Domingues, C. M., Slangen, A. B. A. & Dias, F. B. An ensemble approach to quantify global mean sea-level rise over the 20th century from tide gauge reconstructions. *Environ. Res. Lett.* **16**, 044043 (2021).

40. Thorne, P. W., Parker, D. E., Christy, J. R. & Mears, C. A. Uncertainties in climate trends: Lessons from upper-air temperature records. *Bull. Am. Meteorol. Soc.* **86**, 1437–1442 (2005).

41. MacIntosh, C. R., Merchant, C. J. & von Schuckmann, K. Uncertainties in steric sea level change estimation during the satellite altimeter era: Concepts and practices. *Surv. Geophys.* **38**, 59–87 (2017).

42. Camargo, C. M. L., Riva, R. E. M., Hermans, T. H. J. & Slangen, A. B. A. Trends and uncertainties of cam-driven sea-level change in the satellite altimetry era. *Earth Syst. Dyn.* **13**, 1351–1375 (2022).

43. Cheng, L. et al. Improved estimates of ocean heat content from 1960 to 2015. *Sci. Adv.* **3**, e1601545 (2017).

44. Loeb, N. G. et al. Clouds and the Earth's radiant energy system (CERES) energy balanced and filled (EBAF) top-of-atmosphere (TOA) edition-4.0 data product. *J. Clim.* **31**, 895–918 (2018).

45. Loeb, N. G., Manalo-Smith, N., Su, W., Shankar, M. & Thomas, S. CERES top-of-atmosphere Earth radiation budget climate data record: Accounting for in-orbit changes in instrument calibration. *Remote Sens.* **8**, 182 (2016).

46. Raghuraman, S. P., Paynter, D. & Ramaswamy, V. Anthropogenic forcing and response yield observed positive trend in Earth's energy imbalance. *Nat. Commun.* **12**, 4577 (2021).

47. Loeb, N. G. et al. Evaluating twenty-year trends in Earth's energy flows from observations and reanalyses. *J. Geophys. Res. Atmos.* **127**, e2022JD06686 (2022).

48. Kramer, R. J. et al. Observational evidence of increasing global radiative forcing. *Geophys. Res. Lett.* **48**, e2020GL091585 (2021).

49. Arias, P. A. et al. Technical summary. In *Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change* (eds. Masson-Delmotte, V. et al.). 33–144 https://doi.org/10.1017/9781009157896.001 (Cambridge University Press, 2021).

50. Baggensitos, D. et al. Earth's radiative imbalance from the last glacial maximum to the present. *Proc. Natl. Acad. Sci.* **116**, 14881–14886 (2019).

51. Church, J. A., White, N. J. & Arblaster, J. M. Significant decadal-scale impact of volcanic eruptions on sea level and ocean heat content. *Nature* **438**, 74–77 (2005).

52. Smith, D. M. et al. Earth's energy imbalance since 1960 in observations and CMIP5 models. *Geophys. Res. Lett.* **42**, 1205–1213 (2015).

53. Cheng, L. & Zhu, J. Artifacts in variations of ocean heat content induced by the observation system changes. *Geophys. Res. Lett.* **41**, 7276–7283 (2014).

54. Jenkins, S. et al. Is Anthropogenic global warming accelerating?. *J. Clim.* **35**, 4273–4290 (2022).

55. Quaas, J. et al. Robust evidence for reversal of the trend in aerosol effective climate forcing. *Atmos. Chem. Phys.* **22**, 12221–12239 (2022).
56. Stephens, G. L. et al. The changing nature of Earth's reflected sunlight. *Proc. R. Soc. Math. Phys. Eng. Sci.* **478**, 20220053 (2022).

57. Meyssignac, B., Chenal, J., Loeb, N., Guillaume-Castel, R. & Ribes, A. Time-variations of the climate feedback parameter λ are associated with the Pacific Decadal Oscillation. *Commun. Earth Environ.* **4**, 1–10 (2023).

58. Purkey, S. G. et al. Unabated bottom water warming and freshening in the South Pacific Ocean. *J. Geophys. Res. Oceans* **124**, 1778–1794 (2019).

59. Cheng, L., Foster, G., Hausfather, Z., Trenberth, K. E. & Abraham, J. Improved quantification of the rate of ocean warming. *J. Clim.* **35**, 4827–4840 (2022).

60. Loeb, N. G. et al. Toward optimal closure of the Earth's top-of-atmosphere radiation budget. *J. Clim.* **22**, 748–766 (2009).

61. Loeb, N. G., Kato, S., Loukachine, K. & Manalo-Smith, N. Angular distribution models for top-of-atmosphere radiative flux estimation from the clouds and the Earth's radiant energy system instrument on the Terra satellite. Part I: Methodology. *J. Atmos. Ocean. Technol.* **22**, 338–351 (2005).

62. Doelling, D. R. et al. Geostationary enhanced temporal interpolation for CERES flux products. *J. Atmos. Ocean. Technol.* **30**, 1072–1090 (2013).

63. Loeb, N. G., Kato, S. & Wielicki, B. A. Defining top-of-the-atmosphere flux reference level for Earth radiation budget studies. *J. Clim.* **15**, 3301–3309 (2002).

64. Visser, H., Dangendorf, S. & Petersen, A. C. A review of trend models applied to sea level data with reference to the "acceleration-deceleration debate". *J. Geophys. Res. Oceans* **120**, 3873–3895 (2015).

65. Wunsch, C. *The Ocean Circulation Inverse Problem*. https://doi.org/10.1017/CBO9780511629570 (Cambridge University Press, 1996).

66. Lyman, J. M. et al. Robust warming of the global upper ocean. *Nature* **465**, 334–337 (2010).

67. von Schuckmann, K. & Le Traon, P.-Y. How well can we derive global ocean indicators from Argo data?. *Ocean Sci.* **7**, 783–791 (2011).

68. Storto, A., Cheng, L. & Yang, C. Revisiting the 2003–18 deep ocean warming through multiplatform analysis of the global energy budget. *J. Clim.* **35**, 4701–4717 (2022).

69. Trenberth, K. E., Fasullo, J. T., von Schuckmann, K. & Cheng, L. Insights into Earth's energy imbalance from multiple sources. *J. Clim.* **29**, 7495–7505 (2016).

70. Santer, B. D. et al. Statistical significance of trends and trend differences in layer-average atmospheric temperature time series. *J. Geophys. Res. Atmos.* **105**, 7337–7356 (2000).

71. Santer, B. D. et al. Consistency of modelled and observed temperature trends in the tropical troposphere. *Int. J. Climatol.* **28**, 1703–1722 (2008).

72. Foster, G. & Rahmstorf, S. Global temperature evolution 1979–2010. *Environ. Res. Lett.* **6**, 044022 (2011).

73. Raghuraman, S. P., Paynter, D., Menzel, R. & Ramaswamy, V. Forcing, cloud feedbacks, cloud masking, and internal variability in the cloud radiative effect satellite record. *J. Clim.* **36**, 4151–4167 (2023).

## Acknowledgements

We thank F. Gues, M. Hakuba, A. Barnoud, M. Ablain, L. Cheng, F. Landerer, and T. Frederikse for helpful discussions. We also express our gratitude to the three anonymous reviewers for their valuable comments, which have significantly contributed to the improvement of the manuscript. JBS and LV have received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement no. 821001.

## Author contributions

All subsequent analysis for this paper was performed by A.M. and supervised by K.v.S and J.-B.S. L.V. helped estimate internal climate variability from CMIP6 pre-industrial simulations. All authors contributed to interpreting the results and writing the manuscript.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1038/s41598-023-49353-1.

**Correspondence** and requests for materials should be addressed to A.M.

**Reprints and permissions information** is available at www.nature.com/reprints.

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_14.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2023