Heat stored in the Earth system 1960–2020: where does the energy go?

The Earth climate system is out of energy balance, and heat has accumulated continuously over the past decades, warming the ocean, the land, the cryosphere, and the atmosphere. According to the Sixth Assessment Report by Working Group I of the Intergovernmental Panel on Climate Change, this planetary warming over multiple decades is human-driven and results in unprecedented and committed changes to the Earth system, with adverse impacts for ecosystems and human systems. The Earth heat inventory provides a measure of the Earth energy imbalance (EEI) and allows for quantifying how much heat has accumulated in the Earth system, as well as where the heat is stored. Here we show that the Earth system has continued to accumulate heat, with 381±61 ZJ accumulated from 1971 to 2020. This is equivalent to a heating rate (i.e., the EEI) of 0.48±0.1 W m⁻². The majority, about 89 %, of this heat is stored in the ocean, followed by about 6 % on land, 1 % in the atmosphere, and about 4 % available for melting the cryosphere. Over the most recent period (2006–2020), the EEI amounts to 0.76±0.2 W m⁻². The Earth energy imbalance is the most fundamental global climate indicator that the scientific community and the public can use as the measure of how well the world is doing in the task of bringing anthropogenic climate change under control. Moreover, this indicator is highly complementary to other established ones like global mean surface temperature as it represents a robust measure of the rate of climate change and its future commitment. We call for an implementation of the Earth energy imbalance into the Paris Agreement's Global Stocktake based on best available science. The Earth heat inventory in this study, updated from von Schuckmann et al. (2020), is underpinned by worldwide multidisciplinary collaboration and demonstrates the critical importance of concerted international efforts for climate change monitoring and community-based recommendations and we also call for urgently needed actions for enabling continuity, archiving, rescuing, and calibrating efforts to assure improved and long-term monitoring capacity of the global climate observing system. The data for the Earth heat inventory are publicly available, and more details are provided in Table 4.

1 Introduction

The Earth energy imbalance (EEI) is the most fundamental indicator for climate change, as it tells us if, how much, how fast, and where the Earth's climate is warming, as well as how this warming evolves in the future (Hansen et al., 2011, 2005; von Schuckmann et al., 2016). The EEI is given by the difference between incoming solar radiation and outgoing radiation, which determines the net radiative flux at the top of the atmosphere (TOA). Today, the Earth climate system is out of energy balance; consequently, heat has accumulated continuously over the past decades, warming the ocean, the land, the cryosphere, and the atmosphere, determining the Earth heat inventory (Fig. 1, von Schuckmann et al., 2020). This planetary warming is human-driven and results in unprecedented and committed changes to the Earth system (Fig. 1) (IPCC, 2021), with adverse impacts for ecosystems and human systems (IPCC, 2022a). As long as this imbalance persists (or even increases) planet Earth will keep gaining energy, increasing planetary warming (Hansen et al., 2005, 2017). Today, the EEI can be best estimated from the quantification of the Earth heat inventory, complemented by direct measurements from space (von Schuckmann et al., 2016; Loeb et al., 2021). In addition, the Earth heat inventory as derived from multiple sources of measurements and models also allows researchers to unravel where the energy – mostly in the form of heat – is stored in the Earth system across all components (von Schuckmann et al., 2020). Results of the first internationally driven initiative on the Earth heat inventory (von Schuckmann et al., 2020) not only show how much and where heat has accumulated in the Earth system but also show for the first time that the Earth energy imbalance has increased over the recent decade. This increase is expected to have fundamental implications for the Earth's climate, and several potential drivers have been discussed recently (Loeb et al., 2021; Hakuba et al., 2021; Kramer et al., 2021).

The Earth system responds to an imposed radiative forcing through a number of feedback mechanisms, which operate on various timescales. Earth's radiative response is complex, comprising a variety of climate feedback mechanisms (e.g., water vapor feedback, cloud feedback, ice–albedo feedback) (Forster et al., 2021). Conceptually, the relationships between EEI, radiative forcing, and surface temperature change can be expressed as the following (Gregory and Andrews, 2016):

$$\begin{array}{}\text{(1)}& \mathrm{\Delta }{N}_{\mathrm{TOA}}=\mathrm{\Delta }{F}_{\mathrm{ERF}}-\left|{\mathit{\alpha }}_{\mathrm{FP}}\right|\mathrm{\Delta }{T}_{\mathrm{S}},\end{array}$$

where ΔN_TOA is the Earth's net energy imbalance at TOA (in Wm⁻²), ΔF_ERF is the effective radiative forcing (Wm⁻²), ΔT_S is the global surface temperature anomaly (K) relative to the equilibrium state, and α_FP is the net total feedback parameter (Wm⁻² K⁻¹), which represents the combined effect of the various climate feedback mechanisms. Essentially, α_FP in Eq. (1) can be viewed as a measure of how efficient the system is at restoring radiative equilibrium for a unit surface temperature rise. Thus, ΔN_TOA represents the difference between the applied radiative forcing and Earth's radiative response through climate feedback associated with surface temperature increase (e.g., Hansen et al., 2011). Observation-based estimates of ΔN_TOA are therefore crucial to our understanding of past climate change and for refining projections of future climate change (Gregory and Andrews, 2016; Kuhlbrodt and Gregory, 2012). The long atmospheric lifetime of carbon dioxide means that ΔN_TOA, ΔF_ERF, and ΔT_S will remain positive for centuries, even with substantial reductions in greenhouse gas emissions, and lead to substantial sea-level rise, ocean warming, and ice shelf loss (Cheng et al., 2019; Forster et al., 2021; Hansen et al., 2017; IPCC, 2021a; Nauels et al., 2017). In other words, warming will continue even if atmospheric greenhouse gas (GHG) amounts are stabilized at today's level, and the EEI defines additional global warming that will occur without further change in forcing (Hansen et al., 2017). The EEI is, in principle, less subject to decadal variations associated with internal climate variability than global surface temperature and therefore represents a robust measure of the rate of climate change and its future commitment (Cheng et al., 2017b; Forster et al., 2021; Loeb et al., 2018; Palmer and McNeall, 2014; von Schuckmann et al., 2016).

Figure 1 Schematic overview on the central role of the Earth heat inventory and its linkage to anthropogenic emissions, the Earth energy imbalance, change in the Earth system, and implications for ecosystems and human systems. The Earth heat inventory plays a central role for climate change monitoring as it provides information on the absolute value of the Earth energy imbalance, the total Earth system heat gain, and how much and where heat is stored in the different Earth system components. Examples of associated global-scale changes in the Earth system as assessed in Gulev et al. (2021) are drawn, together with major implications for the ecosystem and human systems (IPCC, 2022b). Upward arrows indicate increasing change, downward arrows indicate decreasing change, and turning arrows indicate change in both directions. The percentages for heat stored in the Earth system components are provided over the period 2006–2020 (see Sect. 6).

The heat gain in the Earth system from a positive EEI results in directly and indirectly triggered changes in the climate system, with a variety of implications for the environment and human systems (Fig. 1). One of the most direct implications from a positive EEI is the rise of global mean surface temperature. The accumulation and storage of surplus anthropogenic heat leads to ocean warming and thermal expansion of the water column, which together with terrestrial ice melt leads to sea-level rise (WCRP Global Sea Level Budget Group, 2018). Moreover, there are various facets of impacts from ocean warming such as on climate extremes, which are provided in more detail in a recent review (Cheng et al., 2022a). The heat accumulation in the Earth system also leads to warming of the atmosphere, particularly to a temperature increase in the troposphere, leading to water vapor increase and changes in atmospheric circulation (Gulev et al., 2021).

On land, the heat accumulation leads to an increase in ground heat storage, which in turn triggers an increase in ground surface temperatures that may increase soil respiration and may lead to a decrease in soil water, depending on the climatic and meteorological conditions and factors such as land cover and soil characteristics (Cuesta-Valero et al., 2022; Gulev et al., 2021). Moreover, inland water heat storage increases, leading to increases in lake water temperatures that may result in algal blooms and lake stratification, and typically leads to a decrease in lake ice cover. Heat gain in the Earth system also induces an increase in permafrost heat content, which in turn leads to disruptive changes in ground morphology, CH₄ and CO₂ emissions, and a decrease in permafrost extent and ground ice volume. More details are synthesized in Cuesta-Valero et al. (2023a). In the cryosphere, associated changes include a loss of glaciers, ice sheets, and Arctic sea ice (IPCC, 2019, 2021a). These human-induced changes have already impacted ecosystems and have adverse impacts on human systems (Fig. 1). Particularly, they have emerged for ecosystem structure and species ranges and phenology (timing of life cycles), and they include adverse impacts such as for water security and food production; health and wellbeing; and cities, settlements, and infrastructure (IPCC, 2022b; see their Fig. SPM.2).

Regularly assessing, quantifying, and evaluating the Earth heat inventory creates a unique opportunity to support the call to action and solution pathways as assessed during the sixth assessment cycle of the Intergovernmental Panel on Climate Change (IPCC). Moreover, the Earth heat inventory allows for a regular stocktaking of the implementation of the Paris Agreement¹ while monitoring progress towards achieving the purpose of the agreement and its long-term goals based on best available science. These assessment outcomes further emphasize the need to extend the Global Climate Observing System (GCOS) beyond the strict scientific observation of the climate state to also support policy and planning (GCOS, 2021). Science-driven studies driven by an Earth system view and backed by concerted multidisciplinary and international collaborations play a critical role to support these objectives (Crisp et al., 2022; Dorigo et al., 2021; von Schuckmann et al., 2020). With this second study, we aim to contribute to a more frequent and regular science-driven update of the state of the Earth heat inventory as an important indicator of climate change.

Based on the quantification of the Earth heat inventory published in 2020 (von Schuckmann et al., 2020), we present the updated results of the Earth heat inventory over the period 1960–2020, along with the long-term Earth system's heat gain over this period, and the partitions of where the heat goes for the ocean, atmosphere, land, and cryosphere. Section 2 provides the updates for ocean heat content, which is based on improved evaluations (e.g., trend evaluation method) and the addition of further international data products of subsurface temperature. Updated estimates and refinements for atmospheric heat content are discussed in Sect. 3. For the land component in Sect. 4, an improved uncertainty framework is proposed for the ground heat storage estimate, and new evaluations for inland freshwater heat storage and thawing of permafrost have been included (Cuesta-Valero et al., 2022). An update of the heat available to melt the cryosphere is described in Sect. 5 based on re-enforced international collaboration. In Sect. 6, the updated Earth heat inventory is established and discussed based on the results of Sects. 2–5. In the final section, challenges and recommendations for future improved estimates are discussed for each Earth system component, with associated recommendations for future evolution of the observing system.

2 Heat stored in the ocean

Global ocean heat content (OHC) can be estimated directly from subsurface temperature measurements, which is one of the variables of the in situ component of the Global Ocean Observing System (GOOS²) and which has continued to evolve during the past decades (Abraham et al., 2013; Gould et al., 2013; Moltmann et al., 2019). The evolution of the ocean observing system for subsurface temperature measurements is provided, for example, in Cheng et al. (2022a), leveraging the transition from historical measures to modern autonomous techniques, which achieved near-global coverage in the year 2006 (the so-called golden Argo era). Different research groups have developed gridded products of subsurface temperature fields and ocean heat content using different processing methodologies (Abraham et al., 2022; Boyer et al., 2016; Cheng et al., 2022b; Gulev et al., 2021; Li et al., 2022; Savita et al., 2022). Additionally, specific Argo-based products are listed on the Argo web page (http://www.argo.ucsd.edu/, last access: 12 July 2022). Near-global OHC can also be indirectly estimated from spatial geodetic measurements by combining sea surface height from altimetry and ocean mass from gravimetry to solve the sea-level budget equation (Meyssignac et al., 2019; Dieng et al., 2017; Llovel et al., 2014). Spatial geodetic OHC is available since 2002 and provides full-depth OHC variations (Marti et al., 2022; Hakuba et al., 2021). Ocean reanalysis systems have also been used to deliver estimates of near-global OHC (Trenberth et al., 2016; von Schuckmann et al., 2018), and their international assessments show increased agreement with increasing in situ data availability for the assimilation, particularly when the Argo project had achieved nearly global-scale data sampling (Fig. 2) (Palmer et al., 2017; Storto et al., 2018, 2019; Meyssignac et al., 2019).

This initiative relies on the availability of regular updates to data products, their temporal extensions, and direct interactions with the different research groups. A complete view of all subsurface ocean temperature products can only be achieved through a concerted international effort and over time, particularly accounting for the continued development of new or improved OHC products. In this study, we do not achieve a holistic view of all available products but present a starting point for future international regular assessments of global OHC. A first established international ensemble mean and standard deviation of near-global OHC up to 2018 was established in von Schuckmann et al. (2020), which has now been updated up to 2020 and further extended with the addition of five new products (Fig. 3). The ensemble spread gives an indication of the agreement among products and can be used as a proxy for uncertainty. Compared to the results in von Schuckmann et al. (2020), the spread has increased, which can be referred back to the additional use of data products, the impact of year-to-year variations, and the refined use of the ensemble spread approach (see below).

Although there has been a tremendous improvement in in situ subsurface temperature measurements over time, estimates of global OHC remain an area of active research to minimize the major effects from different data processing techniques of the irregular (in space and time) in situ database and associated sampling characteristics, followed by the choice of the climatology used in the mapping process and data bias corrections, which today induce discrepancies between the different estimates (Boyer et al., 2016; Cheng et al., 2014; Gouretski and Cheng, 2020; Cheng et al., 2018; Good, 2017; Savita et al., 2022; Allison et al., 2019). Concerns about common errors in the products remain. Accurate understanding of the uncertainties of the product is an essential element in their use. So far, a basic assumption is that the error distribution for the observations is Gaussian with a mean of zero, which has been approximated by an ensemble of various products. However, a more complete understanding of any apparent trends requires the determination of systematic errors (e.g., systematic calibration errors) (or the impacts of changing observation densities through a synthetic profile approach; Allison et al., 2019) and of instrument technologies (Wong et al., 2020). These elements can result in biases across the ensemble, or produce artificial changes in the energetics of the system (Wunsch, 2020). For example, Li et al. (2022) estimated that assuming linear vertical interpolation with sparse historical vertical profiles results is an underestimation of global ocean heat content (and ocean thermal expansion) trends since the 1950s on the order of 14 % compared with a more sophisticated vertical interpolation scheme (Barker and McDougall, 2020; Li et al., 2022), with the greatest systematic underestimates at latitudes 15–20^∘ N, and Li et al. (2022) also found that interannual differences between various eXpendable BathyThermograph (XBT) corrections were similar to the differences when only higher-quality hydrographic data were included, implying the need for improved time-dependent XBT corrections. The uncertainty can also be estimated in other ways, including some purely statistical methods (Levitus et al., 2012; MacIntosh et al., 2017; Cheng et al., 2019) or methods explicitly accounting for the error sources (Gaillard et al., 2016; Lyman and Johnson, 2014; von Schuckmann and Le Traon, 2011). Each method has its caveats; for example, the error covariances are mostly unknown and must be estimated a priori. For this study, adopting a straightforward method with a “data democracy” strategy (i.e., all OHC estimates have been given equal weight) has been chosen as a starting point, which is different from the ensemble approach adopted in the Sixth Assessment Report (AR6) (Forster et al., 2021).

Figure 2 Ensemble mean time series and ensemble standard deviation (95 %, shaded) of global ocean heat content (OHC) anomalies relative to the 2005–2020 climatology for the 0–300 m (gray), 0–700 m (blue), 0–2000 m (yellow), and 700–2000 m depth layer (green). The ensemble mean is an outcome of an international assessment initiative, and all products used are referenced in the legend of Fig. 3. The trends derived from the time series are given in Table 1. Note that values are given for the ocean surface area between 60^∘ S and 60^∘ N and are limited to the 300 m bathymetry of each product.

Figure 3 Trends of global ocean heat content (OHC) as derived from different products (colors) and using LOWESS (see text for more details). References are given in the figure legend, except for CMEMS (CORA, Copernicus Marine Ocean Monitoring Indicator, 2023), EN.4.2.2.c14 (Good et al., 2013) with Cheng et al. (2014), XBT (Gouretski and Cheng, 2020), and Mechanical Bathythermograph (MBT) bias corrections, as well as for the method by Palmer et al. (2007). CSIRO-GEOMAR-NOC (Argo) (Domingues et al., 2008; Wijffels et al., 2016; Roemmich et al., 2015), CSIRO-GEOMAR-NOC (hist) (Church et al., 2011; Domingues et al., 2008), NOC (National Oceanographic Centre) (Desbruyères et al., 2017), and the Argo dataset MOAA GPV (Hosoda et al., 2008) are also included. Results from the Copernicus Marine Service Global Reanalysis Ensemble Product have been added as well (Copernicus Marine Ocean Monitoring Indicator: global ocean heat content) for comparison but are not considered for the ensemble mean in Fig. 1. The ensemble mean and standard deviation (95 % confidence interval) are indicated in black. The shaded areas show trends from different depth layer integrations, i.e., 0–300 m (light turquoise), 0–700 m (light blue), 0–2000 m (purple), and 700–2000 m (light purple). For each integration depth layer, trends are evaluated over the three study periods, i.e., historical (1960–2020), altimeter era (1993–2020), and golden Argo era (2006–2020). See text for more details on the international assessment criteria. Note that values are given for the ocean surface area (see text for more details). References as indicated in the legend include the following: Cheng et al. (2017a); Gaillard et al. (2016); Good et al. (2013); Ishii et al. (2017); Kuusela and Giglio (2022); Levitus et al. (2012); Li et al. (2017, 2022); Lyman and Johnson (2014); Roemmich and Gilson (2009); and von Schuckmann and Le Traon (2011).

The continuity of this activity will help to further expand international collaboration and to unravel uncertainties due to the community's collective efforts on data quality as well as on detecting and reducing processing uncertainties. It also provides up-to-date scientific knowledge of ocean warming. Products used for this assessment are referenced in the caption of Fig. 3. Estimates of OHC have been provided by the different research groups under homogeneous criteria: all estimates use a coherent ocean volume limited by the 300 m isobath (700 m for Li et al., 2022) of each product and are limited to 60^∘ S–60^∘ N, since most observational products exclude high-latitude ocean areas because of the low observational coverage, and only annual averages have been used. The ocean areas within 60^∘ S–60^∘ N include 91 % of the global ocean surface area, and limiting to the 300 m isobath neglects the contributions from coastal and shallow waters, so the resultant OHC trends will be underestimated if these ocean regions are warming. For example, neglecting shallow waters is estimated to account for more than 10 % for 0–2000 m OHC trends (Savita et al., 2022; von Schuckmann et al., 2014) and about 4 % for the Arctic area (J. Mayer et al., 2021). The assessment is based on three distinct periods to account for the evolution of the observing system, i.e., 1960–2020 (i.e., “historical”), 1993–2020 (i.e., “altimeter era”), and 2006–2020 (i.e., “golden Argo era”). All time series go up to 2020 – which was one of the principal limitations for the inclusion of some products. Our final estimates of OHC for the 0–300, 0–700, 700–2000, and 0–2000 m depth layers are the ensemble average of all products, with the uncertainty range defined by the standard deviation (2σ, 95 % confidence interval) of the corresponding ensemble used (Fig. 2).

For the trend evaluation, we have followed the most recent study by Cheng et al. (2022b) and used a locally weighted scatterplot smoothing (LOWESS) approach to reduce the effect of high-frequency variability (e.g., year-to-year variability), data noise, or changes in the observing system as it relies on a weighted regression (Cleveland, 1979) within a prescribed span width of 25 years for the historical and altimeter era and 15 years for the recent period (2006–2020). The change in OHC(t) over a specific period, ΔOHC, is then calculated by subtracting the first value from the last value of the fitted time series, OHC_LOWESS(t), to obtain the trend while dividing by the considered period. To obtain an uncertainty range on the trend estimate and to take into account the sensitivity of the calculation to interannual variability, we implement a Monte Carlo simulation to generate 1000 surrogate series OHC_random(t), under the assumption of a given mean (our “true” time series OHC(t)) (Cheng et al., 2022b). Each surrogate OHC_random(t) consists of the fitted true time series OHC(t) plus a randomly generated residual which follows a normal (Gaussian) distribution and which is included in an envelope equal to 2 times the uncertainty associated with the time series. Then, a LOWESS fitted line is estimated for each of the 1000 surrogates. The 95 % confidence interval for the trend is then calculated based on ±2 times the standard deviation (±2σ) of all 1000 trends of the surrogates. However, the use of either trend estimates following a linear, LOWESS approach, or the approach discussed in Palmer et al. (2021) leads to consistent results within uncertainties (not shown).

In agreement with Cheng et al. (2019) and Gulev et al. (2021), our results confirm a continuous increase in ocean warming over the entire study period (Fig. 2). Moreover, rates of global ocean warming have increased over the three different study periods, i.e., historical up to the recent decadal change. The trend values are all given in Table 1. The major fraction of heat is stored in the upper ocean (0–300 and 0–700 m depth). However, heat storage at intermediate depth (700–2000 m) increases at a nearly comparable rate as reported for the 0–300 m depth layer (Table 1, Fig. 3). There is a general agreement among the 16 international OHC estimates (Fig. 3). However, for some periods and depth layers the standard deviation (95 % confidence level) reaches maxima at about 0.3 W m⁻². All products agree on the fact that global ocean warming rates have increased in the past decades and doubled since the beginning of the altimeter era (1993–2020 compared with 1960–2020) (Fig. 3). Moreover, there is a clear indication that heat sequestration took place in the 700–2000 m depth layer over the past six decades linked to an increase in OHC trends over time (Fig. 3). Ocean warming rates for the 0–2000 m depth layer reached record rates of 1.03 (0.62) ±0.2 W m⁻² over the period 2006–2020 for the ocean (global) area, consistent with what has been reported in Johnson et al. (2022).

Table 1 OHC trends using LOWESS (locally weighted scatterplot smoothing; see text for more details) as derived from the ensemble mean (Fig. 2) for different time intervals, as well as different integration depths. The regression was done for each time period (1960–2020, 1971–2020, 1993–2020, 2006–2020). A time window of 25 years was used for the periods that allowed it (1960–2020, 1971–2020, 1993–2020). For the period 2006–2020, a time window of 15 years was used. Note that values are given in watts per square meter (Wm⁻²) relative to the global surface (a factor of 0.61 for the ocean surface is used considering the area 60^∘ S–60^∘ N, >300 m bathymetry). See also text and Figs. 2–3 for more details. Additionally, values for satellite-derived estimates of OHC have been added for the most recent period, which are updated after Hakuba et al. (2021) and Marti et al. (2022).

For the deep OHC changes below 2000 m, we adapted an updated estimate from Purkey and Johnson (2010) (PG10 hereinafter) from 1992 to 2020, which is a constant linear trend estimate (0.97±0.48 ZJ yr⁻¹, 0.06±0.03 W m⁻²) derived from a global integration of OHC below 2000 m using basin-scale deep-ocean temperature trends from repeated hydrographic sections. Some recent studies strengthened the results in PG10 (Desbruyères et al., 2016; Zanna et al., 2019). Desbruyères et al. (2016) examined the decadal change of the deep and abyssal OHC trends below 2000 m in the 1990s and 2000s, suggesting that there has not been a significant change in the rate of decadal global deep/abyssal warming from the 1990s to the 2000s, and the overall deep-ocean warming rate is consistent with PG10. Using a Green's function method and Estimating the Circulation and Climate of the Ocean (ECCO) reanalysis data, Zanna et al. (2019) reported a deep-ocean warming rate of ∼0.06 W m⁻² during the 2000s, consistent with PG10 used in this study. Zanna et al. (2019) show a fairly weak global trend during the 1990s, which is different from observation-based estimates. This mismatch might come from how surface–deep connections are represented in ECCO reanalysis data and the use of time-mean Green's functions in Zanna et al. (2019), as well as from the sparse coverage of the observational network for relatively short time spans. Furthermore, combining hydrographic and deep-Argo floats, a recent study (Johnson et al., 2019) reported an accelerated warming in the South Pacific Ocean in recent years, but a global estimate of the OHC rate of change over time is not available yet, and the rates of warming may vary by ocean basin. Comparison of the results in Table 1 with OHC estimates derived from the space geodetic approach (Hakuba et al., 2021; Marti et al., 2022) shows overall agreement within uncertainties.

Before 1992, we assume zero OHC trend below 2000 m due to insufficient global observations below 2000 m, following the methodology in some studies (Cheng et al., 2017a, 2022b), IPCC AR5 (Rhein et al., 2013), and IPCC AR6 (Forster et al., 2021; Gulev et al., 2021). The deep warming is likely driven by decadal variability in deep water formation rates, which could have been in a non-steady-state mode prior to 1990, introducing additional uncertainty to the pre-1990 OHC estimates. Using surface temperature observations and assuming the heat is advected by mean circulation, Zanna et al. (2019) show a near-zero (small cooling trend) OHC trend below 2000 m from the 1960s to 1980s, suggesting the trend before 1992 might be small. The derived time series of PG10 after 1991 and zero-trend before 1992 is used for the Earth energy inventory in Sect. 5. A centralized (around the year 2006) uncertainty approach has been applied for the deep (>2000 m depth) OHC estimate, following the method by Cheng et al. (2017a), which allows us to extract an uncertainty range over the period 1993–2018 within the given (lower (0.96–0.48 ZJ  yr⁻¹), upper (0.96+0.48 ZJ yr⁻¹)) range of the deep OHC trend estimate. We then extend the obtained uncertainty estimate back from 1992 to 1960, with zero OHC anomaly.

3 Heat available to warm the atmosphere

The heat content of the atmosphere is small compared to those of the other Earth subsystems. Yet it is by no means negligible, since, in relative terms, the atmospheric heat gain is rapid over the recent decades and has a high impact on human life (Fig. 1) (IPCC, 2021a). Atmospheric observations show a warming of the troposphere and a cooling and contraction of the stratosphere since at least 1979 (Pisoft et al., 2021; Steiner et al., 2020a). In the tropics, the upper troposphere has warmed faster than the near-surface atmosphere since at least 2001, as seen with the new observation technique of GPS radio occultation (Gulev et al., 2021; Ladstädter et al., 2023; Steiner et al., 2020a, b), while observations based on microwave soundings have likely underestimated tropospheric temperature trends in the past (Santer et al., 2021; Zou et al., 2021).

Recently, a continuous rise in the height of the tropopause has been observed for 1980 to 2020 over the Northern Hemisphere (Meng et al., 2022). The increase is equally due to tropospheric warming and stratospheric cooling in the period 1980 to 2000, while the rise after 2000 resulted primarily from enhanced tropospheric heat gain. Moreover, indications exist of a widening of the tropical belt (Grise et al., 2019; Fu et al., 2019; Staten et al., 2020) as well as of changes in the seasonal cycle (Santer et al., 2022). However, changes in atmospheric circulation and related conditions for extreme weather are still subject to uncertainty (Cohen et al., 2020), while the occurrence of heat-related extreme weather events has clearly increased over recent decades (Cohen et al., 2020; IPCC, 2021b), with high risks for societies, economies, and the environment (Fischer et al., 2021).

A regular assessment of atmospheric heat content changes is hence critical for a complete overview of energy and mass exchanges with other climate components and for a complete energy budgeting of Earth's climate system.

3.1 Atmospheric heat content

In a globally averaged and vertically integrated sense, heat accumulation in the atmosphere arises from a small imbalance between net energy fluxes at the top of the atmosphere (TOA) and the surface (denoted s). The total energy budget of the vertically integrated and globally averaged atmosphere (indicated by the global averaging operator 〈.〉) reads as follows (Mayer et al., 2017):

$$\begin{array}{}\text{(2)}& \begin{array}{rl}〈{\displaystyle \frac{\partial E_{\mathrm{A}}}{\partial t}}〉& =\langle N_{\mathrm{TOA}}\rangle -\langle F_{\mathrm{s}}\rangle -\langle F_{\mathrm{snow}}\rangle \\ & -\langle F_{\mathrm{PE}}\rangle ,\end{array}\end{array}$$

where the vertically integrated atmospheric energy content E_A per unit surface area (J m⁻²) reads as

$$\begin{array}{}\text{(3)}& E_{\mathrm{A}}=\underset{z_{\mathrm{s}}}{\overset{z_{\mathrm{TOA}}}{\int }}\mathit{\rho }\left(c_{\mathrm{v}}T+g\left(z-z_{\mathrm{s}}\right)+L_{\mathrm{e}}q+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{V}^{\mathrm{2}}\right)\mathrm{d}z.\end{array}$$

In Eq. (2), N_TOA is the net radiation at top of the atmosphere, F_s is the net surface energy flux defined as the sum of net surface radiation and latent and sensible heat fluxes, F_snow denotes the latent heat flux associated with snowfall, and F_PE additionally accounts for sensible heat of precipitation. See Mayer et al. (2017) or von Schuckmann et al. (2020) for a discussion of the last two terms, which are small on a global scale and hence often neglected.

Equation (3), formulated in mean-sea-level altitude (z) coordinates and used here for integrating over observational data, provides a decomposition of E_A into sensible heat energy (sum of the first two terms, internal heat energy, and gravity potential energy), latent heat energy (third term), and kinetic energy (fourth term), where ρ is the air density, c_v the specific heat for moist air at constant volume, T the air temperature, g the acceleration of gravity, L_e the temperature-dependent effective latent heat of condensation L_v or sublimation L_s (the latter relevant below 0 ^∘C), q the specific humidity of moist air, and V the wind speed. We neglect atmospheric liquid water droplets and ice particles as separate species, as their amounts and especially their trends are small.

In computing E_A for the purpose of this update to the von Schuckmann et al. (2020) heat storage assessment, we continued to use the formulations described therein, including that we refer to (geographically aggregated) E_A as atmospheric heat content (AHC) in this context. This acknowledges the dominance of the heat-related terms in Eq. (3). Briefly, in deriving the AHC from observational datasets, we accounted for the intrinsic temperature dependence of the latent heat of water vapor in formulating L_e (for details, see Gorfer, 2022), while the reanalysis derivations approximated L_e by constant values of L_v, as this simplification is typically also made in the assimilating models (e.g., ECMWF-IFS, 2015). As another small difference, the observational estimations neglected the kinetic energy term in Eq. (3), while the reanalysis estimations accounted for it. However, the resulting differences in AHC anomalies from any of these differences are negligibly small, especially when considering trends over time.

3.2 Datasets and heat content estimation

Turning to the actual datasets used, the AHC and its changes and trends over time can be quantified using various data sources. Reassessing possible data sources, we extended the high-quality datasets that we used in the initial von Schuckmann et al. (2020) assessment. In particular, we updated the time period from 2018 to 2020 and improved the backward extension from 1980 to 1960. Specifically, the adopted datasets and the related AHC data record preparations can be summarized as follows.

Atmospheric reanalyses combine observational information from various sources (radiosondes, satellites, weather stations, etc.) and a dynamical model in a statistically optimal way. These data have reached a high level of maturity, thanks to continuous improvement work since the early 1990s (Hersbach et al., 2018). Especially reanalyzed thermodynamic state variables, like temperature and water vapor that are most relevant for AHC computation, are of high quality and suitable for climate studies, although temporal discontinuities introduced from changing observing systems continue to deserve due attention (Berrisford et al., 2011; Chiodo and Haimberger, 2010; Hersbach et al., 2020; M. Mayer et al., 2021).

We use the latest generation of reanalyses, including ECMWF's fifth-generation reanalysis ERA5 (Bell et al., 2021; Hersbach et al., 2020), Japan Meteorological Agency (JMA)'s reanalysis JRA55 (Kobayashi et al., 2015), and NASA's Modern-Era Retrospective analysis for Research and Applications version 2 (MERRA2) (Gelaro et al., 2017). ERA5 and JRA55 are both available over the full joint time frame of this heat storage assessment from 1960 to 2020, while MERRA2 complements these from 1980 to 2020. The additional JRA55C reanalysis variant of JRA55, included for initial intercomparison in von Schuckmann et al. (2020), is no longer used since it is available to 2012 only, and due to its similarity to JRA55 it is not adding appreciable complementary value.

In addition to these three reanalyses, the datasets from two climate-quality observation techniques are used for complementary observational AHC estimates. These include the Wegener Center (WEGC) multi-satellite radio occultation (RO) data record (WEGC OPSv5.6; Angerer et al., 2017; Steiner et al., 2020b), over 2002–2020, and a radiosonde (RS) data record derived from the high-quality Vaisala sondes RS80/RS92/VS41 (WEGC Vaisala; Ladstädter et al., 2015), covering 1996–2020. These RO and RS datasets provide atmospheric profiles of temperature, specific humidity, and density that are vertically completed by colocated ERA5 profiles in domains not fully covered by the data (e.g., in the lower troposphere for RO or at polar latitudes for RS). Similar to dropping the JRA55C reanalysis variant for no longer adding appreciable further value, the simplified AHC proxy data based on microwave sounding unit (MSU) observational data, intercompared in von Schuckmann et al. (2020), are no longer used.

From the observational data, the AHC is estimated by first evaluating Eq. (3) (using all terms for total and the third term only for latent AHC) at each available profile location and subsequently deriving it as volumetric heat content, for up to the global scale, from vertical integration, temporal averaging, and geographic aggregation according to the approach summarized in von Schuckmann et al. (2020) and described in detail by Gorfer (2022). For the reanalyses, the estimation is based on the full gridded fields. Applying the approach for cross-check with reanalysis profiles subsampled at observation locations only, confirms its validity as it accurately leads to the same AHC results as from the full gridded fields.

Overall, the ensemble spread of all the atmospheric datasets used is deemed a reasonable proxy for the uncertainty in the ensemble-mean annual AHC anomaly data, in particular since 1980 during the “satellite observation era” (e.g., Hersbach et al., 2020; Steiner et al., 2020a). However, the uncertainties of the trend estimates, i.e., of the AHC increase rates (“AHC gain”) obtained from linear fitting to the anomaly data over periods of interest (see Sect. 3.3), are weakly depending on these data uncertainties anyway, since the trend uncertainties are dominated by the interannual natural variability in the data, which is significantly larger than the data uncertainties expressed by the ensemble spread (see Fig. 4).

3.3 Atmospheric heat content change since 1960 and its amplification

Figure 4 shows the resulting global AHC change inventory over 1960 to 2020 (61-year record), in terms of total AHC anomalies for each data type (Fig. 4a) and for the ensemble mean with trends for selected periods and uncertainty estimates (Fig. 4c). The selected trend periods align with those for ocean data and with availability of atmospheric datasets (see Sect. 3.2 above) and represent a reference trend 1961–2000 plus recent trends of the last about 30, 20, and 15 years. Latent AHC anomalies, a key component of the AHC (Matthews et al., 2022), are also shown (Fig. 4b and d). Compared to von Schuckmann et al. (2020), the AHC data have the El Niño–Southern Oscillation (ENSO) signal removed (with ENSO regressed out via the Niño 3.4 Index; and cross-check with non-ENSO-corrected data is showing that trend differences are reasonably small). However, variability due to volcanic eruptions is still included and may somewhat influence the trends over 1993–2020, which start in the cold anomaly after the Pinatubo eruption (Santer et al., 2001).

The latent AHC (Fig. 4b and d), which accounts for about one-quarter of the total AHC, exhibits a qualitatively similar temporal evolution as total AHC, but with larger relative uncertainty compared to the total AHC. The RO and RS datasets in Fig. 4b show some differences, particularly the low latent AHC values in the 1990s and early 2000s from the RS WEGC Vaisala dataset likely stem from known dry biases of the RS80/RS90/RS92 humidity sensors (Wang et al., 2002; Verver et al., 2006; Vömel et al., 2007). Estimated trends based on these RS data are thus likely too high, although the overall increase in latent AHC is also substantial in the other datasets.

Figure 4 Annual-mean global AHC anomalies from 1960 to 2020 of total AHC (a, c) and latent-only AHC (b, d), respectively, of three different reanalyses and two different observational datasets shown together with their mean (a, b) and the mean AHC anomaly shown together with four representative AHC trends and ensemble spread measures of its underlying datasets (c, d). The in-panel legends identify the individual datasets (a, b) and the selected trend periods together with the associated trend values (plus 90 % confidence range) and ensemble spread measures (c, d), with the latter including the time-average standard deviation and minimum/maximum deviations of the individual datasets from the mean.

The results clearly show that the AHC trends have increased from the earlier decades, represented by the 1961–2000 trend of near 1.7 TW. We find the mean trend to be about 2.5 times higher over 1993–2020 (about 5.3 TW) and about 4 times higher in the most recent two decades (about 6–7 TW), a period that is already covered by the RO and RS records. Latent AHC trends in the most recent periods are about 3 times larger than the 1961–2000 reference period. Since 1971, the heat gain in the atmosphere amounts to 5±1 ZJ (see also Fig. 8).

The remarkable amplification of total AHC and latent AHC trends is highlighted in Fig. 5 and summarized in Table 2 for the representative recent periods vs. the 1961–2000 reference period. The 1961–2000 and 1993–2020 periods were covered by reanalysis only, while the WEGC Vaisala RS dataset additionally covers the 2001–2020 and 2006–2020 periods, and the RO dataset covers the most recent period (see dataset descriptions in Sect. 3.2). The larger diversity of recent datasets induces more spread; for example, the RS dataset shows an amplification factor of near 4.5 in the global total AHC gain for 2001–2020, while the amplification factors from the reanalyses range from 2.6 to 3.8. Amplifications are generally largest in the Southern Hemisphere extratropics, where the 1961–2000 reference gain is smallest, and weakest in the tropics. In the most recent period (2006–2020), the amplification factors are strongest, with the RS and RO datasets on the high end of the spread (near a factor of 5 in global total AHC) and somewhat smaller but still high from the reanalyses (around a factor of 4).

For the latent AHC amplification factors, we see moderate values in the 1993–2020 period in the global mean and tropics. In the tropics, the lower uncertainty bound for amplification is slightly below 1 during all three recent trend periods. The spread of the amplification factors increases for the most recent periods, which is, on the one hand, due to the shorter duration. The range increase is also related to the inclusion of the RS and RO datasets after 2000, which contribute the largest and smallest latent AHC gain amplification factors. For 2006–2020, the global mean amplification factor from RO is about 2, whereas from the RS dataset it is near 5. Regarding latitudinal bands, the amplification factors are again strongest in the extratropics, where the 1961–2000 reference gains are also smallest, exhibiting a large spread, especially in the southern extratropics. The relatively large amplification factors of the RS WEGC Vaisala dataset are likely exaggerated due to the well-documented dry bias of the early RS humidity sensors as noted above (Wang et al., 2002; Vömel et al., 2007; Verver et al., 2006).

Despite the uncertainties and spread described, the overall message from Fig. 5 and Table 2 is very clear and substantially reinforcing the evidence from the initial von Schuckmann et al. (2020) assessment: the trends in the AHC, including in its latent heat component, show that atmospheric heat gain has strongly increased over the recent decades.

Figure 5 Amplification of long-term trends in AHC anomalies (“AHC gain”) for total AHC (left) and latent-only AHC (right) in four geographic domains (global, Northern Hemisphere extratropics, tropics, Southern Hemisphere extratropics) for three recent time periods (legend upper left) expressed as a ratio of the trend of each period relative to the trend in the previous-century reference period (1961–2000) (noted below the reference line where the amplification factor equals 1). The amplification factor for each recent-trend case (for the four domains of both total and latent AHC) is depicted for the mean anomaly serving as best estimate (larger black circles) and the related recent trends in the individual-dataset anomalies (colored circles as per upper-right legend). The related 90 % uncertainty range (black error bar) is estimated from the spread (standard deviation) of the individual-dataset amplification factors. The trend in the mean anomaly over 1961–2000 is used as the reference AHC gain.

Table 2 Long-term trend values in mean AHC anomalies (AHC gains; in units of zettajoule per decade (ZJ per decade) and terrawatt (TW), the latter listed in parentheses) and amplification factors vs. the 1961–2000 reference gain (lines marked “Ref.”), for total AHC (left block) and latent-only AHC (right block) for the three recent time periods in four geographic domains as illustrated in Fig. 5. The AHC gain and amplification values are listed together with their 90 % confidence ranges.

4 Heat available to warm land

In previous studies, the land term of the Earth heat inventory was considered as the heat used to warm the continental subsurface (Hansen et al., 2011; Rhein et al., 2013; von Schuckmann et al., 2020). Temperature changes within the continental subsurface are typically retrieved by analyzing the global network of temperature–depth profiles, measured mostly in the Northern Hemisphere, southern Africa, and Australia. Each temperature profile records changes in subsurface temperatures caused by the heat propagated through the ground due to alterations in the surface energy balance (Cuesta-Valero et al., 2021b). Such perturbations in the subsurface temperature profiles can be analyzed to recover the changes in past surface conditions that generated the measured profile, allowing a reconstruction of the evolution of ground surface temperatures and ground heat fluxes at decadal to centennial timescales (Beltrami et al., 2002; Beltrami and Mareschal, 1992; Demezhko and Gornostaeva, 2015; Hartmann and Rath, 2005; Hopcroft et al., 2007; Jaume-Santero et al., 2016; Lane, 1923; Pickler et al., 2016; Shen et al., 1992). Although previous estimates only considered changes in ground temperatures for representing the heat storage by exposed land, ground heat storage has been found to be the second largest term of the Earth heat inventory, accounting for 4 % to 6 % of the total heat in the Earth system (von Schuckmann et al., 2020, Sect. 6).

The ground heat is, nevertheless, not the only energy component of the continental landmasses. Other processes with large thermodynamic coefficients, such as permafrost thawing and the warming of inland water bodies, occur across large areas, leading to the exchange of large amounts of heat with its surroundings over time. To account for those heat exchanges, a recent study (Cuesta-Valero et al., 2023a) has estimated the heat uptake by permafrost thawing and the warming of inland water bodies, as well as ground heat storage from subsurface temperature profiles, resulting in a comprehensive estimate of continental heat storage. Therefore, our estimate is different to terrestrial or land estimates, as we take into account the subsurface and water bodies of the continental landmasses and thus not the land surface. The authors used the same global network of subsurface temperature profiles as in von Schuckmann et al. (2020) to estimate ground heat storage but applied an improved inversion technique to analyze the profiles. This new technique is based on combining bootstrap sampling with a widely used singular value decomposition (SVD) algorithm (e.g., Beltrami and Mareschal, 1992) to retrieve past changes in surface temperatures and ground heat fluxes, which also resulted in smaller uncertainty estimates for global results (Cuesta-Valero et al., 2022). Heat uptake from permafrost thawing was estimated using a large ensemble of simulations performed with the CryoGridLite permafrost model (Nitzbon et al., 2022). Ground stratigraphies required for this purpose, including ground ice distributions, were generated using various global ground datasets. For soil properties, we used the datasets described in Masson et al. (2003) and Faroux et al. (2013); for soil organic carbon, we used the dataset described in Hugelius et al. (2013), and for excess ground ice content, we used Brown et al. (1997). Latent heat storage due to melting of ground ice is evaluated to a depth of 550 m over the Arctic region. Uncertainty ranges are evaluated using 100-parameter ensemble simulations with strongly varied soil properties and soil ice distributions. The climate forcing at the surface is based on a paleoclimate simulation performed by the Commonwealth Scientific and Industrial Research Organization (CSIRO), providing the initialization of the permafrost model, and data from the ERA-Interim reanalysis since 1979 onwards. Heat storage by inland water bodies was estimated by integrating water temperature anomalies in natural lakes and reservoirs from a set of Earth system model (ESM) simulations participating in the Inter-Sectoral Impact Model Intercomparison Project phase 2b (ISIMP2b) (Frieler et al., 2017; Grant et al., 2021; Golub et al., 2022). Heat storage is then computed using simulations with four global lake models following the methodology presented in Vanderkelen et al. (2020) but replacing the cylindrical lake assumption in that study for a more detailed lake morphometry, which leads to a more realistic representation of lake volume.

Figure 6 Continental heat storage from Beltrami et al. (2002) (black), von Schuckmann et al. (2020) (gray), and Cuesta-Valero et al. (2023a) (red). Gray and red shadows show the uncertainty range of the heat storage from von Schuckmann et al. (2020) and Cuesta-Valero et al. (2023a), respectively.

Figure 6 shows the three main estimates of heat gain by the continental landmasses since 1960. The first global estimate of continental heat storage was provided by Beltrami et al. (2002), consisting of changes in ground heat content for the period 1500–2000 as time steps of 50 years (black line in Fig. 6). These estimates were retrieved by inverting 616 subsurface temperature profiles constituting the global network of subsurface temperature profiles in 2002, yielding a heat gain of 9.1 ZJ during the second half of the 20th century. A comprehensive update was included in von Schuckmann et al. (2020) using the results of Cuesta-Valero et al. (2021b) (gray line in Fig. 6), with the main difference consisting of the use of a larger dataset with 1079 subsurface temperature profiles. Since many of these new profiles were measured at a later year than those in Beltrami et al. (2002), the inversions from this new dataset were able to include the recent warming of the continental subsurface, yielding higher ground heat content than those from Beltrami et al. (2002). Concretely, the estimates in von Schuckmann et al. (2020) showed a heat gain of 24±5 ZJ from 1960 to 2018.

Recently, a new estimate of continental heat gain including the heat used in permafrost thawing and in warming inland water bodies was presented in Cuesta-Valero et al. (2023a) (red line in Fig. 6), achieving heat gains of 24±2 ZJ since 1960 and 21±2 ZJ since 1971 (see also Fig. 8). Despite considering the heat stored in permafrost thawing, the warming of inland water bodies, and the warming of the ground, the retrieved continental heat storage is similar to the values from ground warming in von Schuckmann et al. (2020). There is a difference of ∼3 ZJ between the average ground heat storage in Cuesta-Valero et al. (2022) (21.6±0.2 ZJ) and in von Schuckmann et al. (2020) (24±5 ZJ), which is similar to the heat storage in inland water bodies and the heat storage due to permafrost thawing together (see below). That is, the decrease in ground heat storage in the new estimates is compensated by the heat storage in inland water bodies and permafrost degradation. Another important result is the narrower confidence interval in estimates from Cuesta-Valero et al. (2023a), which is directly related to the new bootstrap technique used to invert the subsurface temperature profiles (Cuesta-Valero et al., 2022). This new bootstrap technique offers a more adequate statistical framework than the technique used in von Schuckmann et al. (2020) as demonstrated in Cuesta-Valero et al. (2022); thus, we are confident in the robustness of the lower uncertainty estimate for ground heat storage presented here. Heat storage within inland water bodies has reached 0.2±0.4 ZJ since 1960, with permafrost thawing accounting for 2±2 ZJ. Therefore, ground heat storage is the main contributor to continental heat storage (90 %), with inland water bodies accounting for 0.7 % of the total heat and permafrost thawing accounting for 9 %. Despite the smaller proportion of heat stored in inland water bodies and permafrost thawing, several important processes affecting both society and ecosystems depend on the warming of lakes and reservoirs, as well as on the thawing of ground ice (Gädeke et al., 2021). Therefore, it is important to continue quantifying and monitoring the evolution of heat storage in all three components of the continental landmasses.

5 Heat utilized to melt ice

Changes in Earth's cryosphere affect almost all other elements of the environment including the global sea level, ocean currents, marine ecosystems, atmospheric circulation, weather patterns, freshwater resources, and the planetary albedo (Abram et al., 2019). The cryosphere includes frozen components of the Earth system that are at or below the land and ocean surface: snow, glaciers, ice sheets, ice shelves, icebergs, sea ice, inland water body ice (e.g., lake, river), permafrost, and seasonally frozen ground (IPCC, 2019). In this study, we estimate the heat uptake by the melting of ice sheets (including both floating and grounded ice), glaciers, and sea ice at the global scale (Fig. 7). Notwithstanding the important role that snow cover plays in the Earth's energy surface budget as a result of changes in the albedo (de Vrese et al., 2021; Qu and Hall, 2007; Weihs et al., 2021), its influence on the temperature of underlying permafrost (Jan and Painter, 2020; Park et al., 2015), or on sea ice in the Arctic (Perovich et al., 2017; Webster et al., 2021) and Antarctica (Eicken et al., 1995; Nicolaus et al., 2021; Shen et al., 2022), estimates of changes in global snow cover are still highly uncertain and not included in this inventory. However, they should be considered in future estimates. Similarly, changes in lake ice cover (Grant et al., 2021) are not taken into account here and warrant more attention in the future. Permafrost is accounted for in the land component (see Sect. 4).

Figure 7 Heat uptake (in ZJ) and mass loss (in trillions of tonnes) for the Antarctic Ice Sheet (grounded and floating ice, green), glaciers (orange), Arctic sea ice (purple), Greenland Ice Sheet (grounded and floating ice, red), and Antarctic sea ice (blue), together with the sum of the energy uptake within each one of its components (total, black). Uncertainties are 95 % confidence intervals provided as shaded areas. See the text for more details.

We equate the energy uptake by the cryosphere (glaciers, grounded and floating ice of the Antarctic and Greenland ice sheets, and sea ice) with the energy needed to drive the estimated mass loss. In doing so, we assume that the energy change associated with the temperature change of the remaining ice is negligible. As a result, the energy uptake by the cryosphere is directly proportional to the mass of melted ice:

$$\begin{array}{}\text{(4)}& E=\mathrm{\Delta }M\times (L+c\times \mathrm{\Delta }T),\end{array}$$

where, for any given component, ΔM is the mass of ice loss, L is the latent heat of fusion, c is the specific heat capacity of the ice, and ΔT is the rise in temperature needed to bring the ice to the melting point. For consistency with previous estimates (Ciais et al., 2014; Slater et al., 2021; von Schuckmann et al., 2020), we use a constant latent heat of fusion of 3.34×10⁵ J kg⁻¹, a specific heat capacity of 2.01×10³ J kg^−1 ∘C⁻¹, and a density of ice of 917 kg m⁻³. Estimating the energy used to warm the ice to its melting point requires knowledge of the mean ice temperature for each component. Here we assume a temperature of −15 ^∘C for floating ice in Greenland, −2 ^∘C for the floating ice in Antarctica, $-\mathrm{20}\pm \mathrm{10}$ ^∘C for grounded ice in Antarctica and Greenland, and 0 ^∘C for sea ice and glaciers. Although this assumption is poorly constrained, the energy required to melt ice is primarily associated with its phase transition, and the fractional energy required for warming is a small percentage (<1 % ^∘C⁻¹) of the total energy uptake (Slater et al., 2021). Nevertheless, we include an additional uncertainty of ±10 ^∘C on the assumed initial ice temperature within our estimate of the energy uptake. An overview of all datasets used and their availabilities are provided in Table 3 and are further described in the following.

Table 3 Overview on data used and their availability for the estimate of heat available to melt the cryosphere over the period 1979–2020. Backward extension to 1971 for the heat inventory is based on the assumption of negligible contribution. General specifications include constant values for latent heat of fusion of 3.34×10⁵ J kg⁻¹, specific heat capacity of 2.01×10³ J kg^−1 ∘C⁻¹, density of ice of 917 kg m⁻³ for first-year ice and 882 kg m⁻³ for multiyear ice; see also Ciais et al. (2014), Slater et al. (2021), and von Schuckmann et al. (2020). Other component specifications are provided in the table.

Grounded ice losses from the Greenland and Antarctic ice sheets from 1992 to 2020 are estimated from a combination of 50 satellite-based estimates of ice sheet mass balance produced from observations of changes in ice sheet volume, flow, and gravitational attraction, compiled by the Ice Sheet Mass Balance Inter-comparison Exercise (IMBIE³) (Shepherd et al., 2018, 2019). To extend those time series further back in time, we use ice sheet mass balance estimates produced using the input–output method, which combines estimates of solid ice discharge with surface mass balance estimates. Satellite estimates of ice velocity are available from the Landsat historical archive from 1972, allowing for the calculation of ice discharge before the 1990s, while surface mass balance is estimated from regional climate models. We extend the IMBIE mass balance time series backwards to 1979 for Greenland using Mouginot et al. (2019) and Mankoff et al. (2019) and for Antarctica from 1972 to 1991 using Rignot et al. (2019).

Changes in Antarctic floating ice shelves due to thinning between 1994 and 2017 are derived from satellite altimetry reconstructions (Adusumilli et al., 2020). There were no estimates of ice shelf thinning between 1979 and 1993; therefore, we assume zero mass loss from ice shelf thinning during that period. Changes in Antarctic ice shelves due to increased calving in the Antarctic Peninsula and the Amundsen Sea sector are derived from ERS-1 radar altimetry (Adusumilli et al., 2020) for 1994–2017. For the 1979–1994 period, we only have data for changes in the extent of the Antarctic Peninsula ice shelves from Cook and Vaughan (2010). These are converted to changes in mass using an ice shelf thickness of 140±110 m ice equivalent, which represents the range of ice thickness values for the portions of Antarctic Peninsula ice shelves that have collapsed since 1994 (Adusumilli et al., 2020). Once icebergs break off from large Antarctic floating ice shelves, the timescales of dissolution of the icebergs are largely unknown; therefore, we assumed a linear rate of energy uptake between 1979–2020. For icebergs, we use an initial temperature of −16 ^∘C, which was the mean ice temperature in the Ross Ice Shelf J-9 ice core (Clough and Hansen, 1979). There are no large-scale observations or manifestations of significant firn layer temperature change for the Antarctic Ice Shelf; for example, there is no significant trend in the observationally constrained model outputs of surface melt described in Smith et al. (2020). Therefore, the change in temperature of any ice that does not melt is assumed to be negligible.

Changes in the floating portions of the Greenland Ice Sheet include ice shelf collapse, ice shelf thinning, and tidewater glacier retreat. As in von Schuckmann et al. (2020), we assume no ice shelf mass loss before 1997 and estimate a loss of 13 Gt yr⁻¹ after 1997 based on studies of the Zacharie Isstrom, C. H. Ostenfeld, Petermann, Jakobshavn, 79N, and Ryder glaciers (Mouginot et al., 2015; Moon and Joughin, 2008; Münchow et al., 2014; Motyka et al., 2011; Wilson et al., 2017). We assign a generous uncertainty of 50 % to this value. For tidewater glacier retreat, we note a mean retreat rate of 37.6 m yr⁻¹ during 1992–2000 and 141.7 m yr⁻¹ during 2000–2010 (Carr et al., 2017). We assume that the former estimate is also valid for 1979–1991 and that the latter estimate is valid for 2011–2020. Assuming a mean glacier width of 4 km and thickness of 400 m, we estimate mass loss from glacier retreat to be 9.3 Gt yr⁻¹ during 1979–2000 and 35.1 Gt yr⁻¹ during 2000–2020. Based on firn modeling, we assessed that warming of Greenland's firn has not yet contributed significantly to its energy uptake (Ligtenberg et al., 2018).

The contributions from both the Antarctic and Greenland ice sheets to the EEI are obtained by summing the mass loss from the individual components (ice shelf mass, grounded ice mass, and ice shelf extent) for each ice sheet separately, and, given that the datasets used for each component are independent, the uncertainties were summed in quadrature. This is then converted to an energy uptake according to the equation above.

Glaciers are another part of the land-based ice, and we here include glaciers found in the periphery of Greenland and Antarctica (but distinct from the ice sheets) in our estimate. We build our estimate on the international efforts to compile and reconcile measurements of glacier mass balance, under the lead of the World Glacier Monitoring Service (WGMS⁴). Up to 2016, the results are based on Zemp et al. (2019), who combined geodetic mass balance observations from digital elevation model (DEM) differencing on long temporal and large spatial scales with in situ glaciological observations, which are spatially less representative but provide information of higher temporal resolution. Through this combination, they achieve coverage that is globally complete yet retains the interannual variability well. For 2017 to 2021, the numbers are based on the ad hoc method by Zemp et al. (2020), which corrects for the spatial bias of the limited number of recent in situ glaciological observations that are available with short delay (WGMS, 2021), to derive globally representative estimates. Error bars include uncertainties related to the in situ and spaceborne observations, extrapolation to unmeasured glaciers, density conversion, and to glacier area and its changes. For the conversion from mass loss to energy uptake, only the latent heat uptake is considered, which is based on the assumption of ice at the melting point, due to a lack of glacier temperature data at the global scale. Moreover, since the absolute mass change estimates are based on geodetic mass balances, mass loss of ice below floatation is neglected. While this is a reasonable approximation concerning the glacier contribution to sea-level rise, it implies a systematic underestimation of the glacier heat uptake. While to our knowledge there are no quantitative estimates available of glacier mass loss below sea level on the global scale, it is reasonable to assume that this effect is minor, based on the volume–altitude distribution of glacier mass (Farinotti et al., 2019; Millan et al., 2022). Further efforts are underway within the Glacier Mass Balance Intercomparison Exercise (GlaMBIE⁵), particularly to reconcile global glacier mass changes, including estimates from gravimetry and altimetry, and to further assess related sources of uncertainties (Zemp et al., 2019).

Sea ice, formed from freezing ocean water and further thickened by snow accumulation, is not only another important aspect of the albedo effect (Kashiwase et al., 2017; Zhang et al., 2019) and water formation processes (Moore et al., 2022) but also provides essential services for polar ecosystems and human systems in the Arctic (Abram et al., 2019). Observations of sea ice extent are available over the satellite era, i.e., since the 1970s, but ice thickness data – required to obtain changes in volume – have only recently become available through the launch of CryoSat-2 and ICESat-2. For the Arctic, we use a combination of sea ice thickness estimates from the Pan-Arctic Ice-Ocean Modeling and Assimilation System (PIOMAS) between 1979 and 2011 (Schweiger et al., 2019; Zhang and Rothrock, 2003) and CryoSat-2 satellite radar altimeter measurements between 2011 and 2020 when they are available (Tilling et al., 2018; Slater et al., 2021). PIOMAS assimilates ice concentration and sea surface temperature data and is validated with most available thickness data (from submarines, oceanographic moorings, and remote sensing) and against multidecadal records constructed from satellite (Labe et al., 2018; Laxon et al., 2013; Wang et al., 2016). We note that the PIOMAS domain does not extend sufficiently far south to include all regions covered by sea ice in winter (Perovich et al., 2017). Given that the entirety of the regions that are unaccounted for (e.g., the Sea of Okhotsk and the Gulf of St Lawrence) are only seasonally ice covered since the start of the record, this should not influence the results. We convert monthly estimates of sea ice volume from CryoSat-2 satellite altimetry to mass using densities of 882 and 916.7 kg m⁻³ in regions of multiyear and first-year ice, respectively (Tilling et al., 2018). During the summer months (May to September), the presence of melt ponds on Arctic sea ice makes it difficult to discriminate between radar returns from leads and sea ice floes, preventing the retrieval of summer sea ice thickness from radar altimetry (Tilling et al., 2018). As a result, we use the winter-mean (October to April) mass trend across the Arctic for both CryoSat-2 and PIOMAS estimates for consistency. According to PIOMAS, winter Arctic sea ice mass estimates are 19 Gt yr⁻¹ (6 %) smaller than the annual mass trend between 1979 and 2011 (−324 Gt yr⁻¹) and so are a conservative estimate of Arctic sea ice mass change (Slater et al., 2021). The uncertainty in monthly Arctic sea ice volume measurements from CryoSat-2 ranges from 14.5 % in October to 13 % in April (Tilling et al., 2018; Slater et al., 2021), and it is estimated as $\pm \mathrm{1.8}\times {\mathrm{10}}^{\mathrm{3}}$ km³ for PIOMAS (Schweiger et al., 2011).

Satellite radar altimeter retrievals of sea ice thickness in the Southern Ocean are complicated by the presence of thick snow layers with unknown radar backscatter properties on Antarctic sea ice floes. As a result, no remote sensing estimates are available for Antarctic sea ice, and we use sea ice volume anomalies from the Global Ice-Ocean Modeling and Assimilation System (GIOMAS; Zhang and Rothrock, 2003), which is the global equivalent to PIOMAS. GIOMAS output has been recently validated against in situ and satellite data by Liao et al. (2022). We compute Antarctic sea ice trends as annual averages between January and December. In the absence of a detailed characterization of uncertainties for these estimates, we use the uncertainty in GIOMAS sea ice thickness of 0.34 m (Liao et al., 2022) to estimate the uncertainty in GIOMAS sea ice volume to be $\pm \mathrm{4.0}\times {\mathrm{10}}^{\mathrm{3}}$ km³, using an annual mean sea ice extent of 11.9×10⁶ km² (Lavergne et al., 2019). One caveat to this is that the observational estimates have their own significant uncertainties (Kern et al., 2019; Liao et al., 2022). For future updates of the Earth heat inventory, we also aim to include observation-based (remote sensing) estimates in the Southern Ocean (Lavergne et al., 2019).

Our estimate of the total heat gain in the cryosphere amounts to 14±4 ZJ over the period 1971–2020 (see also Fig. 8 and Sect. 6) (assuming negligible contribution before 1979 according to the data availability limitation), which is consistent with the estimate obtained in von Schuckmann et al. (2020) within uncertainties. Approximately half of the cryosphere's energy uptake is associated with the melting of grounded ice, while the remaining half is associated with the melting of floating ice (ice shelves in Antarctica and Greenland, Arctic sea ice). Compared to earlier estimates and in particular the 8.83 ZJ estimate from Ciais et al. (2014), this larger estimate is a result of both the longer period of time considered and the improved estimates of ice loss across all components, especially the ice shelves in Antarctica. Contributions to the total cryosphere heat gain are dominated by the Antarctic Ice Sheet (including the floating and grounded ice, about 33 %) and Arctic sea ice (about 26 %), directly followed by the heat utilized to melt glaciers (about 25 %). The Greenland Ice Sheet amounts to about 17 %, whereas Antarctic sea ice is accounted for with a nonsignificant contribution of about 0.2 %.

6 The Earth heat inventory: where does the energy go?

Evaluations of the heat storage in the different Earth system components as performed in Sects. 2–5 now allow for the establishment of the Earth heat inventory. Estimates for all Earth system components cover a core period of 1971–2020, except for the cryosphere where negligible contribution is assumed before 1979. Our results reconfirm a continuous accumulation of heat in the Earth system since our estimate begins (Fig. 8). The total Earth system heat gain in this study amounts to 380±62 ZJ over the period 1971–2020. For comparison, IPCC AR6 obtained a total heat gain of 434.9 (324.5 to 545.5) ZJ for the period 1971–2018 and is hence consistent with our estimate within uncertainties (Forster et al., 2021). However, it is important to note that our estimate still excludes some aspects of Earth heat accumulation, such as, for example, the shallow areas of the ocean, which are challenging to be quantified with respect to gaps in the observing system.

Figure 8 Total Earth system heat gain in ZJ (1 ZJ =10²¹ J) relative to 1960 and from 1960 to 2020. The upper ocean (0–300 m, light blue line, and 0–700 m, light blue shading) accounts for the largest amount of heat gain, together with the intermediate ocean (700–2000 m, blue shading) and the deep ocean below 2000 m depth (dark blue shading). The second largest contributor is the storage of heat on land (orange shading), followed by the gain of heat to melt grounded and floating ice in the cryosphere (gray shading) and heating of the atmosphere (magenta shading). Uncertainty in the ocean estimate also dominates the total uncertainty (dot-dashed lines derived from the standard deviations (2σ) for the ocean, cryosphere, land, and atmosphere). See Sects. 2–5 for more details of the different estimates. The dataset for the Earth heat inventory is published at the German Climate Computation Centre (DKRZ; https://www.dkrz.de/, last access: 29 March 2023) (see Sect. 7). Consistent with von Schuckmann et al. (2020), we obtain a total heat gain of 381±61 ZJ over the period 1971–2020, which is equivalent to a heating rate (i.e., the EEI) of 0.48±0.1 W m⁻² applied continuously over the surface area of the Earth (5.10×10¹⁴ m²). The corresponding EEI over the period 2006–2020 amounts to 0.76±0.2 W m⁻². The LOWESS method and associated uncertainty evaluations have been used as described in Sect. 2.

The estimate of heat storage in all Earth system components not only allows for obtaining a measure of how much and where heat is available for inducing changes in the Earth system (Fig. 1) but also to improve the accuracy of the Earth system's total heat gain. In 1971–2020 and for the total heat gain, the ocean accounts for the largest contributor with an about 89 % fraction of the global inventory. The second largest component in the Earth heat inventory relies on heat stored in land with an about 6 % contribution. The cryosphere component accounts for about 4 % and the atmosphere for about 1 %. For the most recent era of best available GCOS data for the Earth heat inventory since the year 2006, the fractions amount to about 89 % for the ocean, about 5 % for land, about 4 % for the cryosphere, and about 2 % for the atmosphere.

The change of the Earth heat inventory over time allows for an estimate of the absolute value of the Earth energy imbalance. Our results of the total heat gain in the Earth system over the period 1971–2020 is equivalent to a heating rate of 0.48±0.1 W m⁻² and is applied continuously over the surface area of the Earth (5.10×10¹⁴ m²). For comparison, the heat gain obtained in IPCC AR5 amounts to 274±78 ZJ and 0.4 W m⁻² over the period 1971–2010 (Rhein et al., 2013). In IPCC AR6, the total heat rate has been assessed by 0.57 (0.43 to 0.72) W m⁻² for the period 1971–2018 and 0.79 (0.52 to 1.06) W m⁻² for the period 2006–2018 (Forster et al., 2021). Consistently, we further infer a total heating rate of 0.76±0.2 W m⁻² for the most recent era (2006–2020).

Thus, the rate of heat accumulation across the Earth system has increased during the most recent era as compared to the long-term estimate – an outcome which reconfirms the earlier finding in von Schuckmann et al. (2020) and which had then been concurrently and independently confirmed in Foster et al. (2021), Hakuba et al. (2021), Loeb et al. (2021), Liu et al. (2020), Raghuraman et al. (2021), and Kramer et al. (2021). The drivers of a larger EEI in the 2000s than in the long-term period since 1971 are still unclear, and several mechanisms are discussed in literature. For example, Loeb et al. (2021) argue for a decreased reflection of energy back into space by clouds (including aerosol cloud interactions) and sea ice and increases in well-mixed greenhouse gases (GHG) and water vapor to account for this increase in EEI. Kramer et al. (2021) refer to a combination of rising concentrations of well-mixed GHG and recent reductions in aerosol emissions to be accounting for the increase, and Liu et al. (2020) address changes in surface heat flux together with planetary heat redistribution and changes in ocean heat storage. Future studies are needed to further explain the drivers of this change, together with its implications for changes in the Earth system.

Besides heat, which is the focus of this study, Earth also stores energy chemically through photosynthesis in living and dead biomass with plant growth. Recent studies (Friedlingstein et al., 2022; Denning, 2022; Crisp et al., 2022) on the Global Carbon Budget and carbon cycle show that approximately 25 % of the added anthropogenic CO₂ is removed from the atmosphere by increased plant growth, which is a result of fertilization by rising atmospheric CO₂ and nitrogen inputs and of higher temperatures and longer growing seasons in northern temperate and boreal areas (Friedlingstein et al., 2022). This significant increase in carbon uptake by the biosphere indicates that more energy is stored inside biomass, together with the stored carbon. The quantification of the additional amount of energy stored inside the biosphere is outside the scope of this study.

7 Data availability

The time series of the Earth heat inventory are published at DKRZ (https://www.dkrz.de/, last access: 24 January 2023) under https://www.wdc-climate.de/ui/entry?acronym=GCOS_EHI_1960-2020. More details are given below.

• von Schuckmann et al. (2023): data for ocean heat content (Sect. 2) and the total heat inventory as presented in Sect. 6 are integrated.

• Kirchengast et al. (2022): data for the atmospheric heat content are distributed (Sect. 3).

• Cuesta-Valero et al. (2023b): data for the ground heat storage, together with the total continental heat gain, are provided (Sect. 4).

• Vanderkelen and Thiery (2022): data for inland freshwater heat storage are included (Sect. 4).

• Nitzbon et al. (2022b): data for permafrost are delivered (Sect. 4).

• Adusumilli et al. (2022): data for the cryosphere heat inventory are provided.

The Digital Object Identifiers (DOIs) for data access are provided in Table 4.

Table 4 Overview of Digital Object Identifiers (DOIs) for data access for the components of the Earth heat inventory and associated references. The results are presented in Fig. 8.

8 Conclusions

The Earth heat inventory is a global climate indicator integrating fundamental aspects of the Earth system under global warming. Particularly, the Earth heat inventory provides the best available current estimate of the absolute value of the Earth energy imbalance (Cheng et al., 2017a, 2019; Hakuba et al., 2021; Hansen et al., 2011; Loeb et al., 2012, 2022; Trenberth et al., 2016; von Schuckmann et al., 2020). Moreover, its evaluation enables an integrated view of the effective radiative climate forcing, Earth's surface temperature response, and the climate sensitivity (Forster et al., 2021; Hansen et al., 2011, 2005; Palmer and McNeall, 2014; Smith et al., 2015). Additionally, its quantification informs about the status of global warming in the Earth system as it integrates the heat “in the pipeline” that will ultimately warm the deep ocean and melt ice sheets in the long term (Hansen et al., 2011, 2005; IPCC, 2021b). The Earth heat inventory also reveals how much and where surplus anthropogenic heat is available for melting the cryosphere and warming the ocean, land, and atmosphere, which in turn allows for an evaluation of associated changes in the climate system, and it is essential to improve seasonal-to-decadal climate predictions and projections on century timescales to enable improved planning for and adaptation to climate change (Hansen et al., 2011; von Schuckmann et al., 2016, 2020). Regular international assessment on the Earth heat inventory enables concerted international and multidisciplinary collaboration and advancements in climate science, including contributing to the development of recommendations for the status and evolution of the global climate observing system (GCOS, 2021; von Schuckmann et al., 2020).

This study builds on the first internationally and multidisciplinary Earth heat inventory in 2020 (von Schuckmann et al., 2020) and provides an update on total Earth system heat accumulation, heat storage in all Earth system components (ocean, land, cryosphere, atmosphere), and the Earth energy imbalance up to the year 2020. Moreover, this study improved earlier estimates and further extended and fostered international collaboration, allowing researchers to move towards a more complete view of where and how much heat is stored in the Earth system through the addition of new estimates such as for permafrost thawing, inland freshwater (Sect. 4), and Antarctic sea ice (Sect. 5). Results obtained reveal a total Earth system heat gain of 381±61 ZJ over the period 1971–2020, with an associated total heating rate of 0.48±0.1 W m⁻². About 89 % of this heat is stored in the ocean, about 6 % on land, about 4 % in the cryosphere, and about 1 % in the atmosphere (Figs. 8, 9). The analysis additionally reconfirms an increased heating rate which amounts to 0.76±0.2 W m⁻² for the most recent era (2006–2020). The drivers for this change still need to be elucidated, and they most likely reflect the interplay between natural variability and anthropogenic change (Loeb et al., 2021; Kramer et al., 2021; Liu et al., 2020); their implications for changes in the Earth system are reflected in the many record levels of change in the 2000s reported elsewhere (e.g., Cheng et al., 2022b; Forster et al., 2021; Gulev et al., 2021).

The Paris Agreement builds upon the United Nations Framework Convention on Climate Change, and for the first time all nations agreed to undertake ambitious efforts to combat climate change, with the central aim to keep global temperature rise this century well below 2 ^∘C above pre industrial levels and to limit the temperature increase even further to 1.5 ^∘C. Article 14 of the Paris Agreement requires the Conference of the Parties serving as the meeting of the Parties to the Paris Agreement (CMA) to periodically take stock of the implementation of the Paris Agreement and to assess collective progress towards achieving the purpose of the agreement and its long-term goals through the so-called Global Stocktake of the Paris Agreement (GST)⁶ based on best available science. The Earth heat inventory provides information on how much and where heat is accumulated and stored in the Earth system. Moreover, it provides a measure of how much the Earth is out of energy balance, and when combined with directly measured net flux at the top of the atmosphere, it also enables us to understand the change of the EEI over time. This in turn allows for assessing the portion of the anthropogenic forcing that the Earth's climate system has not yet responded to (Hansen et al., 2005) and defines additional global warming that will occur without further change in human-induced forcing (Hansen et al., 2017). The Earth heat inventory is thus one of the key critical global climate change indicators defining the prospects for continued global warming and climate change (Hansen et al., 2011; von Schuckmann et al., 2016, 2020). Hence, we call for an implementation of the Earth heat inventory into the Global Stocktake.

Figure 9 Schematic presentation on the Earth heat inventory for the current anthropogenically driven positive Earth energy imbalance (EEI) at the top of the atmosphere (TOA). The relative partition (in %) of the Earth heat inventory presented in Fig. 8 for the different components is given for the ocean (upper: 0–700 m, intermediate: 700–2000 m, deep: >2000 m), land, cryosphere (grounded and floating ice), atmosphere, and EEI for the periods 2006–2020 and 1971–2020 (for the latter period, values are provided in parentheses). The total heat gain (in red) over the period 1971–2020 is obtained from the Earth heat inventory as presented in Fig. 8.

The quantifications presented in this study are the result of multidisciplinary global-scale collaboration and demonstrate the critical importance of concerted international efforts for climate change monitoring and community-based recommendations for the global climate observing system. For the ocean observing system, the core Argo sampling needs to be sustained – which includes the maintenance of shipboard collection of reference data for validation – and complemented by remote sensing data. Extensions such as into the deep-ocean layer need to be further fostered, and technical developments for the measurements under ice and in shallower areas need to be sustained and extended. Moreover, continued efforts are needed to further advance bias-correction methodologies, uncertainty evaluations, data recovery, and processing of the historical dataset. Spatial geodetic observations and the closure of the sea-level budget serve as a valuable constraint for the full-column OHC. Although the independent estimates agree within uncertainty, the geodetic approach suggest slightly larger OHC linear trends, especially since 2016. Though efforts are underway to investigate the emerging discrepancy (e.g., Barnoud et al., 2021), the causes are not yet fully understood and require further investigation.

For the ground heat storage, the estimate had been hampered by a lack of subsurface temperature profiles in the Southern Hemisphere, as well as by the fact that most of the profiles were measured before the 2000s. Subsurface temperature data are direct and independent (not proxy) measurements of temperature, yielding information on the temporal variation of the ground surface temperature and ground heat flux at the land surface. A larger spatial-scale dataset of the thermal state of the subsurface from the last millennium to the present will aid in the continued monitoring of continental heat storage, provide initial conditions for land surface model (LSM) components of Earth system models (ESMs) (Cuesta-Valero et al., 2019), and serve as a dataset for validation of climate models' simulations (Cuesta-Valero et al., 2021a, 2016). Progress in understanding climate variability through the last millennium must lean on additional data acquisition as the only way to reduce uncertainty in the paleoclimatic record and on changes to the current state of the continental energy reservoir. Remote sensing data are expected to be very valuable to retrieve recent, past, and future changes in ground heat flux at short timescales with near-global coverage. However, collecting subsurface temperature data is urgent as we must make a record of the present thermal state of the subsurface before the subsurface climate baseline is affected by the downward-propagating thermal signal from current climate heating. Furthermore, an international organization should take responsibility to gather and curate all measured subsurface temperature profiles currently available and those that will be measured in the future, as the current practices, in which individual researchers are responsible for measuring, storing, and distributing the data, have led to fragmented datasets, restrictions in the use of data, and loss of the original datasets. Support from GCOS for international data acquisition and curating efforts would be extremely important in this context.

For the permafrost estimates, the primary sources of uncertainty arise from lacking information about the amount and distribution of ground ice in permafrost regions, as well as measurements of liquid water content (Nitzbon et al., 2022a). Permafrost heat storage is defined as the required heat to change the mass of ground ice at a certain location; thus, monitoring changes in ground ice and water contents would be required to improve estimates of this component of the continental heat storage. Nevertheless, the current monitoring system for permafrost soils is focused on soil temperature, and the distribution of stations is still relatively scarce in comparison with the vast areas that need to be surveyed (Biskaborn et al., 2015). Due to the current limitations in the observational data, a permafrost model was used to estimate the heat uptake by thawing of ground ice. This approach retrieves latent heat fluxes in extensive areas and at depths relevant to analyze the long-term change in ground ice mass, but this is done at the cost of ignoring other relevant processes, such as ground subsidence, to balance model performance with computational resources. Including permafrost heat storage in the Tibetan Plateau is a priority for the next iteration of this work, as well as to explore new methods to evaluate model simulations using the available observations in permafrost areas.

For inland water heat storage, a better representation of lake and reservoir volume would be possible by better accounting for lake bathymetry using the GLOBathy (Khazaei et al., 2022) dataset and results from the upcoming Surface Water and Ocean Topography (SWOT) mission. These improvements in the representation of lake volume and an updated lake mask will be available in the upcoming ISIMIP3 simulation round, next to improved meteorological forcing data (Golub et al., 2022). In contrast to Vanderkelen et al. (2020), the heat storage in rivers is not included in this analysis due to the high uncertainties in simulated river water volume. To reduce the uncertainty in river heat storage, the estimation of river water storage should be improved, together with an explicit representation of water temperature in the global hydrological models (Wanders et al., 2019). These improvements will be incorporated into ISIMIP3 and will lead to better estimates of inland water heat storage, thus enhancing future estimates of continental heat storage. In the long run, these model-based estimates could be supplemented or replaced by observation-based estimates, which would however require a large, global-scale effort to monitor lake and river temperatures at high spatial resolution and over long time periods. Estimates for inland water heat storage and permafrost heat storage in this analysis depend heavily on model simulations, which is a particular challenge for analyzing and adding uncertainty ranges, as the sources of uncertainty in model simulations differ from those in observational records (Cuesta-Valero et al., 2023a). Future estimates should hence focus on a hybrid approach considering in situ measurements, reanalyses, remote sensing data, and model simulations, consistent with the methods employed for deriving cryosphere and atmosphere heat storage for the Earth heat inventory.

For the cryosphere, sustained remote sensing for all of the cryosphere components is critical for quantifying future changes over these vast and inaccessible regions; in situ observations are also needed for process understanding and in order to properly calibrate and validate them. For sea ice, observations of the albedo, area, and ice thickness are all essential – the continuation of satellite altimeter missions with high inclination, polar-focused orbits is critical for our ability to monitor sea ice thickness in particular. Observations of snow thickness with multifrequency altimeters and microwave radiometers are essential for further constraining sea ice thickness estimates. For ice sheets and glaciers, reliable gravimetric, geodetic, and ice velocity measurements; knowledge of ice thickness and extent; snow/firn thickness and density; and the continuation of the now three-decade long satellite altimeter records are essential for understanding changes in the mass balance of grounded and floating ice. The recent failure of Sentinel-1b, which in tandem with Sentinel-1a could be used to systematically measure ice speed changes every 6 d, means that images are now being acquired every 12 d and thus an earlier launch of Sentinel-1c should be encouraged to regain the ability to monitor ice speed changes over short timescales. The estimate of glacier heat uptake is particularly affected by lacking knowledge of ice melt below sea level and, to a lesser degree, lacking knowledge of firn and ice temperatures. This lack of observations is likely related to most studies on glaciers focusing on their contribution to sea-level rise or seasonal water availability, where melt below sea level and warming of ice do not matter much. However, it becomes obvious here that this gap introduces a systematic bias in the estimate of cryospheric energy uptake, which is presumably small compared to the other components but unconstrained. Although the Antarctic sea ice change and the warming of Greenland and Antarctic firn are poorly constrained or have not significantly contributed to this assessment, they may become increasingly important over the coming decades. Similarly, there exists the possibility for rapid change associated with positive ice dynamical feedback at the marine margins of the Antarctic Ice Sheet. Sustained monitoring of each of these components will, therefore, serve the dual purpose of furthering the understanding of the dynamics and quantifying the contribution to Earth's energy budget. In addition to data collection, open access to the data and data synthesis products, as well as coordinated international efforts, are key to the continued monitoring of the ice loss from the cryosphere and its related energy uptake.

For the atmosphere, there is a need to sustain and enhance a coherent operational long-term monitoring system for the provision of climate data records of essential climate variables. Observations from radiosonde stations within the GCOS Reference Upper-Air Network (GRUAN) and from satellite-based Global Navigation Satellite System (GNSS) radio occultation deliver thermodynamic profiling observations of benchmark quality and stability from surface to stratopause. For climate monitoring, it is of critical importance to ensure continuity of such observations with global coverage over all local times. This continuity of radio occultation observations in the future is not sufficiently guaranteed as we are facing an imminent observational gap in the middle to high latitudes for most local times (IROWG, 2021), which is a major concern. Thus, there is an urgent need for satellite missions in high-inclination orbits to provide full global and local-time coverage in order to ensure global climate monitoring. Operational radio occultation missions need to be maintained as support for a global climate observing system, and long-term availability and archiving of measurement data, metadata, and processing information need to be ensured.

In summary, we also call for urgently needed actions for enabling continuity, archiving, rescuing and calibrating efforts to assure improved and long-term monitoring capacity of the global climate observing system for the Earth heat inventory and to complement with measurements from space for assessing the changes in EEI (e.g., Loeb et al., 2021; von Schuckmann et al., 2016). Particularly, the summarized recommendations include the following:

• We need to sustain, reinforce, and even to establish data repositories for historical climate data (archiving).

• We need to reinforce efforts for recovery projects for historical data and associated metadata information (rescuing).

• We need to sustain and reinforce the global climate observing system for assuring the monitoring of the Earth heat inventory targets, such as for the polar, deep, and shallow ocean, as well as of top-of-the-atmosphere radiation fluxes (continuity).

• We need to foster calibration measurements (in situ) for assuring the quality and reliability of large-scale measurement techniques (e.g., remote sensing) and autonomous components (e.g., Argo) (calibrating).

A continuous effort to regularly update the Earth heat inventory is important as this global climate indicator crosses multidisciplinary boundaries and calls for the inclusion of new science knowledge from the different disciplines involved, including the evolution of climate observing systems and associated data products, uncertainty evaluations, and processing tools. The outcomes have further demonstrated how we are able to evolve our estimates for the Earth heat inventory while bringing together different expertise and major climate science advancements through a concerted international effort. All of these component estimates are at the leading edge of climate science. Their union has provided a new and unique insight into the inventory of heat in the Earth system, its evolution over time, and the absolute values. The data product of this effort is made available and thus can be used for climate model validation purposes. The results also demonstrate that further efforts are needed for uncertainty evaluations, such as, for example, the use of synthetic profile analyses. Indeed, improving the climate observing system will allow for reduced uncertainties for estimating the Earth heat inventory. However, further evaluations are needed to unravel uncertainties in the different components of the Earth heat inventory, which rely for example on nonhomogeneous data sampling and large data gaps, the use of different measurement types and statistical approaches, instrumental bias corrections, and their joint analysis of mode-based quantifications.

This study has demonstrated the unique value of such a concerted international effort, and we thus call for a regular evaluation of the Earth heat inventory. This updated attempt presented here has focused on the global area average only, and evolving into regional heat storage and redistribution, the inclusion of various timescales (e.g., seasonal, year to year), and other climate study tools (e.g., indirect methods, ocean reanalyses) would be an important asset of this much-needed regular international framework for the Earth heat inventory. This would also respond directly to the request of GCOS to establish the observational requirements needed to further monitor the Earth's cycles and the global energy budget (GCOS, 2021). The outcome of this study will therefore directly feed into GCOS assessments of the status of the global climate observing system, and the identified observation requirements will guide the development of the next generation of in situ and satellite global climate observations as specified by GCOS by all national meteorological services and space agencies and other oceanic and terrestrial networks.