First comprehensive assessment of industrial-era land heat uptake from multiple sources

The anthropogenically intensified greenhouse effect has caused a radiative imbalance at the top of the atmosphere during the industrial period. This, in turn, has led to an energy surplus in various components of the Earth system, with the ocean storing the largest part. The land contribution ranks second with the latest observational estimates based on borehole temperature profiles, which quantify the terrestrial energy surplus to be 6 % in the last 5 decades, whereas studies based on state-of-the-art climate models scale it down to 2 %. This underestimation stems from land surface models (LSMs) having a subsurface that is too shallow, which severely constrains the land heat uptake simulated by Earth system models (ESMs). A forced simulation of the last 2000 years with the Max Planck Institute ESM (MPI-ESM) using a deep LSM captures 4 times more heat than the standard shallow MPI-ESM simulations in the historical period, well above the estimates provided by other ESMs. However, deepening the LSM does not remarkably affect the simulated surface temperature. It is shown that the heat stored during the historical period by an ESM using a deep LSM component can be accurately estimated by considering the surface temperatures simulated by the ESM using a shallow LSM and propagating them with a standalone forward model. This result is used to derive estimates of land heat uptake using all available observational datasets, reanalysis products, and state-of-the-art ESM experiments. This approach yields values of 10.5–16.0 ZJ for 1971–2018, which are 12 %–42 % smaller than the latest borehole-based estimates (18.2 ZJ).

1 Introduction

The climate system has been experiencing a net heat gain in all its components (von Schuckmann et al., 2020) during the industrial period due to a radiative imbalance at the top of the atmosphere (Wild, 2020; Forster et al., 2021). The radiative imbalance is a consequence of the anthropogenically intensified greenhouse effect (Gulev et al., 2021). Accurately estimating the energy surplus partitioning among the different components helps in understanding their interactions and climate change influences on near-surface climate (Forster et al., 2021).

As a result of the energy imbalance at the surface, the land stored around 6 % (4.3 %–6.6 %, 21±2 ZJ) of the terrestrial energy surplus in the last 5 decades, as derived from observational estimates (von Schuckmann et al., 2020; Cuesta-Valero et al., 2021a, 2023), being the second-largest contributor after the ocean (ca. 90 %, 324±8 ZJ; Levitus et al., 2012; Abraham et al., 2013; von Schuckmann et al., 2020) to the total Earth system energy gain. Regardless of whether the amount associated with land heat uptake storage is comparatively small, climate models should realistically represent the partitioning of energy. State-of-the-art climate models quantify the land contribution to be roughly 2 % (Cuesta-Valero et al., 2021a), underestimating observational results. This underestimation of land heat uptake in climate models can alter the surface energy balance (Mottaghy and Rath, 2006; García-García et al., 2023), with impacts on soil hydrology (Krakauer et al., 2013), particularly in permafrost regions (Andresen et al., 2020). The amount of energy absorbed by the ground is also relevant for other soil processes, such as respiration or productivity (Pries et al., 2017), and for its biogeochemical activity (Soong et al., 2021).

Observational land heat uptake estimates have been derived from borehole temperature profile (BTP) log collections (Mareschal and Beltrami, 1992; Huang et al., 2000; Beltrami et al., 2002; Cuesta-Valero et al., 2021b, 2023). BTPs record transient temperature perturbations relative to the steady geothermal gradient. Assuming that the subsurface acts as a half-infinite fully conductive medium (Carslaw and Jaeger, 1959), BTPs provide information about the land heat uptake up to their logging date (Cuesta-Valero et al., 2021b). BTP logs are scarce and show uneven distributions of sampling dates and locations, as well as variable vertical resolution and depths. These factors may potentially have some impact on resulting estimates of land heat uptake (Beltrami et al., 2015; Melo-Aguilar et al., 2018, 2020). Nevertheless, estimates derived from BTPs are the only observational source of information to quantify the land contribution to the terrestrial energy budget so far (Forster et al., 2021; von Schuckmann et al., 2023).

State-of-the-art Earth system models (ESMs; Eyring et al., 2016) have been incapable of properly reproducing observational land heat uptake estimates. Their limited volume to store energy stems from imposing a bottom boundary condition placement (Smerdon and Stieglitz, 2006; Stevens et al., 2007) in their land surface models (LSMs), from ca. 3 to 10 m in the majority of Coupled Model Intercomparison Project Phase 5 (CMIP5; Taylor et al., 2012) and Phase 6 (CMIP6; Eyring et al., 2016) models. These depths are insufficient to account for the heat penetrated into the subsurface in the industrial period, rendering them inadequate for estimating long-term heat storage (Alexeev et al., 2007; MacDougall et al., 2008; González-Rouco et al., 2009; Cuesta-Valero et al., 2016; Steinert et al., 2021a). Only a small number of CMIP6 ESMs impose a deeper LSM depth (ca. 40 m), most of them using the Community Land Model (CLM; Lawrence et al., 2019) as the LSM (see Table 2). Furthermore, longer-term simulations spanning multi-centennial or millennial timescales, like the ones developed within the Paleoclimate Modelling Intercomparison Project (PMIP; Kageyama et al., 2018), should include deeper LSM depths to accommodate longer-lasting and deeper subsurface energy transfer (Smerdon and Stieglitz, 2006; Alexeev et al., 2007; Steinert et al., 2021a).

Some simulation-based efforts have been devoted to deriving unbiased estimates of land heat uptake using deepened LSMs. Standalone one-dimensional half-infinite space heat conduction forward models (González-Rouco et al., 2006, 2009) were used to yield artificial global mean surface temperature histories and heat uptake estimates, forced by either observational-based surface air temperature (SAT) (Huang, 2006) or ESM-simulated ground surface temperature (GST) data (Stevens et al., 2007; MacDougall et al., 2008). More recently, modifications of the LSM depth (Hermoso de Mendoza et al., 2020; González-Rouco et al., 2021) and soil hydrology (Steinert et al., 2021b) have been introduced into LSMs in standalone mode, and their impact on heat uptake and subsurface temperature variability has been extensively discussed. However, the impact of these LSM changes has not yet been explored within a fully coupled ESM.

This work analyzes for the first time a long-term ESM simulation including a deep LSM. A version of the MPI-ESM (Mauritsen et al., 2019), modified to include a deep zero-flux bottom boundary condition, is used to produce a 2000-year-long simulation including a PMIP4 representation of natural and external forcings (Jungclaus et al., 2017). This simulation is compared with another experiment that was made with the same setup, albeit with the standard shallow LSM version and different initial conditions. The land energy storage is computed in the industrial period for the deep and shallow simulations considering the progressive subsurface warming down to their zero-flux bottom boundary condition. In addition, land energy estimates were derived using an offline forward model driven by the GSTs of the shallow and deep simulations as boundary conditions. This constitutes a first assessment of the impacts of deepening the LSM on land energy storage, and it is extended to all available instrumental and modeling datasets. CMIP6 simulations and hybrid model–observational products (i.e., reanalysis) allow us, for the first time, to have a variety of observational and model-based land heat uptake estimates that enable an exhaustive comparison with BTP-based estimates.

Slivinski et al. (2019)Laloyaux et al. (2018)Poli et al. (2016)Muñoz Sabater et al. (2021)Kalnay et al. (1996)Ebita et al. (2011)Tardif et al. (2019)Rohde and Hausfather (2020)Osborn et al. (2021)Lenssen et al. (2019)Zhang et al. (2019)Willmott and Matsuura (2018)

Table 1 Characteristics of the observational-type global surface temperature datasets used in this work and their trends for various periods. The first seven rows account for reanalysis data, while the last five are for instrumental data. GST data were used, except in the cases marked with (¹), for which surface air temperature (SAT) data were used. The soil temperature at the first layer was taken as the GST in every case. Glacier areas (Antarctica and Greenland) were excluded from the calculations. Temperature trends in the three different periods assessed in Fig. 3 (1950–2000, 1960–2020, and 1971–2018) for every dataset and the means for reanalyses and observational subsets are also provided. (²) indicates cases in which the interval of trend evaluation is shorter due to limited data availability.

¹ SAT used instead of GST. ² Temperature trend value is given for a shorter period due to limited data availability.

Bi et al. (2020)Ziehn et al. (2020)Wu et al. (2019)Chen et al. (2019)Swart et al. (2019)Zhang et al. (2020)Danabasoglu et al. (2020)Gettelman et al. (2019)Cherchi et al. (2019)Lovato et al. (2022)Voldoire et al. (2019)Voldoire et al. (2019)Séférian et al. (2019)Golaz et al. (2019)Golaz et al. (2019)Döscher et al. (2022)Döscher et al. (2022)He et al. (2020)Li et al. (2020)Held et al. (2019)Dunne et al. (2020)Kelley et al. (2020)Kelley et al. (2020)Rind et al. (2020)Williams et al. (2018)Williams et al. (2018)Boucher et al. (2020)Lee et al. (2019)Tatebe et al. (2019)Hajima et al. (2020)Mauritsen et al. (2019)González-Rouco et al. (2021)Yukimoto et al. (2019)Seland et al. (2020)Seland et al. (2020)Wang et al. (2021)Sellar et al. (2019)

Table 2 CMIP6 experiments used in this work (Col. 1) and their corresponding horizontal resolution (Col. 2) and LSM component (Col. 3). The number of layers (Col. 4) in the LSM and the LSM depth (Col. 5), as well as temperature trends in the three different periods assessed in Fig. 3 (1950–2000, 1960–2020, and 1971–2018; Cols. 6–8) for every CMIP6 model and the multi-model mean, are also provided. Information on the deep LSM configuration of the MPI-ESM1.2-LR is also included and shown in bold.

1.1 MPI-ESM simulations

The first part of this work focuses on two simulations performed with version 1.2 of the Max Planck Institute for Meteorology ESM (MPI-ESM; Mauritsen et al., 2019), spanning the time interval 0–1850 CE and extended over the historical period and the SSP585 scenario. The interval encompassing 0–2100 CE is noted hereafter as P2k+. One of the experiments includes a standard shallow LSM while the other one uses a deep LSM (P2k+s and P2k+d, respectively). The shallow experiment used hereafter is the MPI-ESM contribution to PMIP4 (van Dijk et al., 2022). MPI-ESM comprises the ocean model MPIOM1.6 and the atmosphere model ECHAM6.3 (Stevens et al., 2013). The latter is directly coupled to the LSM JSBACH3.2 (JSBACH; Reick et al., 2021) through the surface exchange of mass, momentum, and heat. The low-resolution configuration (MPI-ESM1.2-LR) has an atmospheric horizontal resolution truncated to T63, corresponding approximately to a 200 km grid cell size. This resolution is shared by the LSM. The CMIP6 version of JSBACH uses a five-soil-layer discretization, with a zero-flux bottom boundary condition imposed at 10 m (P2k+s), which is too shallow to properly account for the conductive propagation of decadal and centennial temperature perturbations (Steinert et al., 2021a). The deeper version of JSBACH (González-Rouco et al., 2021; Steinert et al., 2021b) has a 12-soil-layer scheme down to 1417 m and is used here in the P2k+d experiment. The P2k+ simulations allow us to evaluate the effect of having a deep LSM on surface temperatures at long timescales. Also, having a long millennial P2k+d experiment provides realistic, i.e., consistent with external forcing, preindustrial initial conditions for subsurface temperatures down to the bottom boundary condition in 1850.

In addition to the P2k+ experiments, an ensemble of 30 historical and SSP585 simulations (hereafter denoted as 30ENS) run with the CMIP6 version of the MPI-ESM, MPI-ESM1.2-LR, are considered. The experiments use the same model version as P2k+s and serve the purpose of tackling the uncertainties in historical warming due to internal variability.

1.2 Reanalysis, observational, and CMIP6 model products

Global mean temperatures at the ground surface over land coming from seven global reanalyses and five gridded observational databases have also been used (Table 1) to derive land heat uptake estimates. In all cases, glacier areas were excluded. GST data were used from six of the sources (see Table 1), whilst SAT was taken as a surrogate of the GST evolution when GST was not available (six data sources; see Table 1). This substitution leans on the assumption that SAT and GST are strongly coupled at long-term scales (Melo-Aguilar et al., 2018). All sources cover the whole of the 20th century, except for ERA5-Land, NCEP1, and JRA55. Furthermore, all of them except LMRv2.1 cover the first decade of the 21st century. These datasets will allow us to obtain land heat uptake estimates complementary to BTP-based values.

Additionally, GST and subsurface temperature simulations of the historical period and SSP585 scenario from the available CMIP6 ESMs (Eyring et al., 2016) were used (see Table 2). This added up to 36 simulations stemming from 22 different ESMs and some model variants, 38 if accounting for P2k+ simulations since 1850. Both for reanalyses and CMIP6 ESMs, the temperature at the first soil layer was used as the GST. Only the first ensemble member was selected for every specific ESM. In the case of the CMIP6 historical and the SSP585 simulation ensemble with the MPI-ESM1.2-LR (30ENS), only the first member was considered for the analyses hereafter and is referred to as MPI-ESM1.2-LR. Some features of the LSM of every CMIP6 ESM, such as the number of soil layers, LSM depth, and horizontal resolution, can be found in Table 2. The ensemble of CMIP6 simulations complements land heat uptake estimates derived from the reanalysis and observational data with ESM experiments that considered changes in external forcing since 1850 CE, allowing us to tackle the system energy input regardless of the different realization of the internal variability. It is noteworthy that most CMIP6 ESMs include LSM depths lower than 14 m, which is not sufficient to properly represent the propagation of the annual temperature cycle into the soil column (Smerdon and Stieglitz, 2006). Only 14 out of 38 members impose an LSM depth below 40 m, still shallower than the 200 m LSM depth recommended to correctly simulate the propagation of the centennial climate change signal (González-Rouco et al., 2021; Steinert et al., 2021a). Thus, land heat uptake estimates for CMIP6 in this work will be derived from the integration of the LSM subsurface temperatures from every simulation, and corrected estimates will also be provided by using a standalone forward modeling approach (see Sect. 3).

2 Methods

2.1 Subsurface temperature generation

A half-infinite one-dimensional heat conduction forward model (Mareschal and Beltrami, 1992; Cuesta-Valero et al., 2022, 2023) is used to propagate yearly surface temperature data from the sources in Tables 1 and 2 and the P2k+ simulations (both d and s) and to produce forward temperature anomaly profiles (FTP hereafter):

$$\begin{array}{}\text{(1)}& \mathrm{FTP}(t_N,z)=\sum _{i=\mathrm{1}}^{N}T\left(t_i\right)\left[\mathrm{erfc}\left({\displaystyle \frac{z}{\mathrm{2}\sqrt{\mathit{\kappa }{t}_i}}}\right)-\mathrm{erfc}\left({\displaystyle \frac{z}{\mathrm{2}\sqrt{\mathit{\kappa }{t}_{i-\mathrm{1}}}}}\right)\right],\end{array}$$

where T(t_i) is the annual mean surface temperature anomaly at time step t_i, κ is thermal diffusivity, and z is depth. To compute FTPs at a certain time step, surface temperature data are trimmed to the interval spanning from the initial year (t₁; e.g., 1850, 1950, 1960, and 1971) to t_i. Temperature anomalies are calculated by computing annual anomalies with respect to the temperature value of the first year of the trimmed interval, t₁. This allows us to depart from equilibrium initial conditions. These FTPs are equivalent to the BTPs used to yield observational land heat uptake estimates, which are derived by subtracting the geothermal gradient (mean state) from the absolute log temperature values (Mareschal and Beltrami, 1992; Huang et al., 2000). An ensemble of three different FTPs per time step was obtained using different thermal diffusivity values resulting from a Monte Carlo random sampling method of local properties (Davison and Hinkley, 1997), which accounts for the uncertainty in temporal and spatial soil heterogeneity (Cuesta-Valero et al., 2022). The Monte Carlo method consisted of averaging 100 different global maps of thermal diffusivity with values randomly sampled from a realistic range of 0.5 to $\mathrm{1.5}\times {\mathrm{10}}^{-\mathrm{6}}$ m² s⁻¹. The sample size was the original number of grid points over land for every data source in Tables 1 and 2. Then, percentiles 2.5, 50, and 97.5 of the 100 values were taken as the three thermal diffusivity values for computing the global FTPs. This is the approach for every data source except for the P2k+ simulations, for which $\mathrm{0.79}\times {\mathrm{10}}^{-\mathrm{6}}$ m² s⁻¹ is taken as the thermal diffusivity to produce the FTP, the actual value used by JSBACH to resolve the temperature vertical scheme (Reick et al., 2021). All the FTPs were resolved down to a common depth of 1417 m so that results could be comparable to the estimates derived from the P2k+d simulation, limited by a zero-flux bottom boundary condition at 1417 m. This comparison is possible since this LSM depth ensures the thermal decoupling between the ground surface and the bottom layer, which is equivalent to the subsurface acting as a half-infinite conductive medium (González-Rouco et al., 2021; Steinert et al., 2021a). As the uncertainties associated with local thermal properties and subsurface depth are tackled when calculating the FTPs, global differences between mean FTPs and land heat uptake estimates from different sources cannot be due to these factors.

In addition, subsurface temperature data coming from both P2k+ and CMIP6 simulations are used. Temperatures at different subsurface levels for every time step are aggregated to yield subsurface temperature profiles. Subsurface temperature anomaly profiles (STPs hereafter) are subsequently computed by subtracting the initial temperature profile from this T(t_i,z), as follows:

$$\begin{array}{}\text{(2)}& \mathrm{STP}(t_i,z)=T(t_i,z)-T(t_{\mathrm{1}},z),\end{array}$$

where T(t_i,z) is the subsurface temperature at depth z and time step t_i (i.e., the subsurface temperature profile at time step t_i) and T(t₁,z) is the initial subsurface temperature profile.

Figure 1 Temperature evolution in MPI-ESM. Past 2k, historical, and SSP585 (P2k+) global means for (a) SAT, (b) GST, and (c) ST5 anomalies over land areas (excluding Greenland and Antarctica) in 0–2000 CE relative to 1850–1900 for a simulation run with the MPI-ESM1.2-LR with a shallow (subscript s, red lines) and a deep (subscript d, blue lines) version of its LSM, JSBACH. Annual time series are shown in light colors, while 31-year moving averages are portrayed in dark colors. (d) 31-year moving averages in panels a to c (see legend for colors) are represented here together with the ensemble of 30 MPI-ESM1.2-LR CMIP6 standard shallow historical and SSP585 simulations (30ENS, gray lines) in 1850–2100 CE for comparison. (e) SAT (black), GST (purple), and ST5 (pink) differences between P2k+s and P2k+d, filtered with a 31-year moving average. For panels a, b, c, d, and e, solid vertical lines indicate the 6 strongest volcanic events of 0–1850 CE, while dashed vertical lines indicate the 19 weakest events. (f, g, h) Range of temperature variability in annual (A) and 31 year moving average (MA31) time series in the preindustrial period (0–1850 CE) for SAT (f), GST (g), and ST5 (h). Whiskers represent percentiles 10 and 90 and boxes 25 and 75, respectively. (i) Temperature anomaly at 2100 with respect to 0–1850 CE (K) in P2k+ (see legend for colors). In this case, whiskers show the range of variability (percentiles 10 and 90) of the total warming (K) in 30ENS.

2.2 Land heat uptake estimation

The annual time evolution of the land heat uptake, Q_s(t), can be estimated by integrating the FTP anomalies along the whole soil column at every yearly time step as follows:

$$\begin{array}{}\text{(3)}& Q_{\mathrm{s}}\left(t\right)=A\sum _{j=\mathrm{1}}^n\mathrm{\Delta }z\left(j\right)C_{\mathrm{v}}\left(j\right){\displaystyle \frac{\mathrm{FTP}(t,j)+\mathrm{FTP}(t,j+\mathrm{1})}{\mathrm{2}}},\end{array}$$

where Δz(j), C_v(j), and FTP(t,j) are the thickness, volumetric heat capacity, and global mean FTP of layer j, respectively, and A is the Earth's land surface area excluding glaciers (1.34×10¹⁴ m²; Cuesta-Valero et al., 2023). FTPs were calculated using an evenly spaced layering of 1 m thickness. Land heat uptake from both P2k+ and CMIP6 (Table 2) was also calculated from the simulated model subsurface temperatures down to the zero-flux bottom boundary condition, i.e., STPs, by making use of Eq. (2), albeit with unevenly spaced temperature layers in these cases. Soil layering thickness for both shallow and deep configurations of JSBACH is given by González-Rouco et al. (2021), whilst it can be found in the references in Table 2 for CMIP6 models. Analogous to what is done to compute thermal diffusivity (see Sect. 3.1), three values of volumetric heat capacity were derived using a Monte Carlo approach (Davison and Hinkley, 1997; Cuesta-Valero et al., 2022) for every source in Tables 1 and 2 using randomly sampled values between 2.5 and 3.5×10⁶ Jm⁻³ K⁻¹. The combination of three values for volumetric heat capacity and three for thermal diffusivity yields nine land heat uptake estimates, which are represented here in terms of mean and variability (2 standard deviations).

Figure 2 Land heat uptake in the MPI-ESM. (a) FTPs (continuous lines) and MPI-ESM STPs (symbols; see legend) for P2k+d (blue) and P2k+s (red) and for the 30-member ensemble mean of MPI-ESM historical and SSP585 simulations (30ENS, black) in the years 2014 and 2100. The y axis is logarithmic. (b) Land heat uptake derived from the yearly time-step-wise vertical integration of FTPs and STPs for P2k+d (blue), P2k+s (red), and 30ENS (black). The time x axis is unevenly spaced to enhance land heat gain since 1850. For the 30ENS, the confidence interval is also shown in gray (p<0.05).

3 Results

3.1 Influence of a realistically deep LSM

The comparison of the MPI-ESM P2k+s and P2k+d simulated SAT, GST, and soil temperature at layer 5 (ST5) allows us to assess the thermal impact of increasing the LSM depth (Fig. 1). ST5 is selected because it is the deepest level for P2k+s, reaching 9.83 m, and it can be compared with the corresponding layer in P2k+d. Figure 1a to d shows the time evolution of SAT, GST, and ST5 for the two simulations, both in annual averages and with 31-year running means to emphasize the low-frequency response. The response to the common external forcing during the last 2 millennia (Jungclaus et al., 2017) generates similar interannual (e.g., volcanic events), multidecadal, and centennial variability in P2k+s and P2k+d for SAT and GST (Fig. 1a, b, and d), indicating no systematic changes in preindustrial mean state (Fig. 1f) and variability (Fig. 1g and h) from using the deep or shallow versions of the model. The resulting SAT and GST warming of about 5.5 °C at the end of the 21st century agrees well with the range of warming in the ensemble of MPI-ESM runs (Fig. 1e) and with a climate sensitivity of 3.2 °C (Mauritsen and Roeckner, 2020). Note that for each simulation SAT and GST virtually overlap, showing long-term SAT–GST coupling (Melo-Aguilar et al., 2018). However, the impacts of deepening the LSM are noticeable in ST5, consistent with González-Rouco et al. (2021), with P2k+d showing lower high-frequency volcanic cooling (e.g., 1257, 1457, and 1815; Fig. 1c), preindustrial variability (Fig. 1f, g, and h), and warming (0.8 °C by 2100; Fig. 1e and i) than P2k+s. The reduced warming is significant if compared to the range of 21st-century warming values produced by the 30ENS of MPI-ESM simulations.

Figure 3 Decadal multi-source surface temperature and heat uptake variability. (a) Global land mean temperature anomaly in 1850–2020 with respect to 1950 for reanalyses, observational, and a 38-member CMIP6 ensemble. The SAT is plotted both for observational and CMIP6 sources and LMRv2.1 (tags in bold), whereas the GST is given for the remaining reanalyses and CMIP6 simulations (see Tables 1, 2). GST and SAT spread ranges for the CMIP6 ensemble are also portrayed (p<0.05). (b) Decadal global land heat uptake in the industrial period (1850–2020) from the different data sources in (a). Red and blue tones in the heatmap represent decadal energy gains and losses, respectively. The horizontal black line separates CMIP6 models (38 members), including MPI-ESM P2k+ simulations, from reanalysis and observational data sources (12 members). CMIP6 sources are sorted downwards alphabetically, while the rest are ordered by length of decreasing periods. The last row represents the multi-source decadal mean.

Therefore, increasing the depth of the LSM produces significant changes in the subsurface temperatures, yet it does not change surface temperatures. The GSTs of the P2k+s and P2k+d can be considered equal in terms of their range of variability and long-term trends, indicating that the surface temperatures simulated by an LSM are not affected by its depth. This result offers the potential to use GSTs from different sources, e.g., reanalysis or CMIP6 simulations, as a boundary condition for forcing an infinite half-space standalone forward model and deriving realistic land energy storage estimates. This is first verified in Fig. 2 using the P2k+ runs. Figure 2a shows STPs from both P2k+s and P2k+d, the ensemble mean STP from 30ENS, and their respective FTPs by the end of the historical period (2014) and SSP585 (2100). Both STPs and FTPs can capture the different degrees of warming between 2014 and 2100 and agree to portray their penetration down to a depth of about 150 m (Fig. 2a). However, there is a remarkable difference when comparing the FTPs and STPs for P2k+s, P2k+d, and the 30ENS. Whilst the P2k+d FTP is capable of reproducing its analogous STP, FTPs and STPs for P2k+s diverge in depth. This is because the forward model is run with no zero-flux bottom boundary condition imposed, which can be considered analogous to the deep configuration of JSBACH (1417 m). Consequently, FTP warm anomalies are smaller in depth than the ones corresponding to the STP, which complies with a zero-flux bottom boundary condition at 10 m. Imposing a shallow boundary condition generates a warm bias with depth since the ground heat flux is halted from penetrating downwards due to the existence of a zero-flux condition (Smerdon and Stieglitz, 2006; Alexeev et al., 2007; Steinert et al., 2021a). Thus, FTPs derived from P2k+s and 30ENS reproduce P2k+d STP and FTP, realistically representing the half-infinite space thermal conduction that governs heat propagation in observational BTPs (Mareschal and Beltrami, 1992; Cuesta-Valero et al., 2021c).

Figure 2b shows the land heat uptake since 0 CE resulting from the stepwise integration of FTPs and STPs for the P2k+ simulations and FTPs for 30ENS. Both STP- and FTP-derived land heat uptake estimates for P2k+d and FTP-derived estimates for P2k+s depict a preindustrial land heat loss, which becomes more acute in the Little Ice Age (from 1400 to 1850; Miller et al., 2012), and a noticeable heat uptake since 1900. The STP-based heat uptake evolution for P2k+s misses the simulation of heat loss and largely underestimates the magnitude of the subsequent heat uptake shown by P2k+d, yielding a heat uptake of around 18 ZJ by the end of the 21st century, which is 4 times smaller than the values derived from deep FTPs and the P2k+d STP (about 75 ZJ). The agreement of P2k+s and P2k+d on the land heat uptake when the forward modeling approach is used evidences the viability of using surface temperatures as boundary conditions to derive FTPs and to estimate heat uptake when subsurface thermal schemes are too shallow.

3.2 Multi-source land heat uptake estimates

The previous analysis opens the possibility of estimating land heat uptake from a variety of sources. Global mean land surface temperature anomalies with respect to 1950 for the seven reanalyses, five observational databases, and the CMIP6 model ensemble are shown in Fig. 3a. By averaging 38 CMIP6 members, the effect of internal variability is minimized, thus emphasizing the response to forced variability. It can be seen that both SAT and GST CMIP6 multi-model means (thick lines) are fully coupled up to 2000, depicting sudden cooling events in the 1880 and 1990 decades, associated with the volcanic activity of Krakatoa or later El Chichón and Pinatubo (Self and Rampino, 1981; Mccormick et al., 1995). They also capture the cooling during the 1950s and 1960s and the subsequent warming. However, the warming trend is slightly higher for SAT than for GST, which produces a small temperature offset of about 0.1 K by 2020. This offset is negligible compared to the range in temperature variability of the 38 members (SAT, shaded gray; GST, scratched gold), so SAT and GST can be considered as coupled; therefore, SAT can be used as alternative boundary information to feed the forward model when GST data are not available. Regarding the observational and reanalysis datasets, surface temperatures lie within the range of CMIP6 multi-model variability. Nevertheless, observational data render higher temperature anomaly values than all the reanalyses, showing higher temperature trends in 1950–2000 and 1960–2020 (Table 1). Moreover, NCEP1 and JRA55 depict the coldest GST evolutions of all observational-based sources, even lying beneath the multi-model spread of CMIP6 models in some instances in the mid-70s.

Figure 4 Land heat uptake estimates from different sources and periods. Land heat uptake during 1950–2000 (a), 1960–2020 (b), and 1971–2018 (c) derived from FTPs (hollow points) of reanalysis and observational databases (Table 1), CMIP6 FTPs (right-hand hollow points) and direct integration of CMIP6 STPs (left-hand solid points; Table 2), and previous estimates (orange crosses): B02 (Beltrami, 2002), B02b (Beltrami et al., 2002), H06 (Huang, 2006), vS20 (von Schuckmann et al., 2020), and CV23 (Cuesta-Valero et al., 2023). The vS20 and CV23 estimates for 1950–2000 refer in fact to 1960–2000 and are marked with a star. The LSM depth for every CMIP6 model used in this work is given by gray shading (see legend above), while the multi-model mean for both the STP and FTP-based values is given by gold crosses. MPI-ESM P2k+d estimates are plotted in dark blue, and MPI-ESM P2k+s estimates are plotted in red.

The decadal evolution of the different surface temperature datasets impacts their land heat uptake, as shown in Fig. 3b, where decadal heat gains and losses are assembled for the time span of each of the 50 datasets collected in this work. All the sources agree on depicting a progressive land warming, only interrupted by transient coolings in the 1880s, 1900s, and 1960s, which coincide with the volcanic eruptions of Krakatoa (Self and Rampino, 1981), Péleé (Gueugneau et al., 2020), and Agung (Hansen et al., 1978). Moreover, heat gain is enhanced in 1980–2020, rendering a land uptake of around 14 ZJ. Even though land heat uptake intensifies in these decades, its increase is not linear, with 1990–2000 and 2010–2020 showing a deceleration in the heat uptake trend (0.4 and 0.8 ZJ, respectively) with respect to the 1970–2020 values (increasing land heat uptake of 1.2 ZJ per decade on average). The uncertainty associated with variability in soil properties is negligible and accounts for less than 0.1 ZJ decadally. Differences between sources in the intensity or the time occurrence of decadal heat gains or losses can result from differences in the representation of internal variability, as suggested by the higher consistency among estimates from observations compared to those from CMIP6 models. Products showing the greatest decadal temperature trends in Tables 1 and 2, such as ACCESS-ESM1-5 or GISTEMPv4, generate the largest storage of energy, with 16.4 and 17.2 ZJ during 1980–2020, respectively.

The estimates of land heat uptake derived from reanalyses, observational data, and CMIP6 models (Fig. 3) are compared to values previously reported in the literature derived from logged BTPs (Beltrami et al., 2002; Beltrami, 2002; von Schuckmann et al., 2020; Cuesta-Valero et al., 2023) and CRUTEM2 SAT (Huang, 2006) for the periods 1950–2000, 1960–2020, and 1971–2018 (Fig. 4). The selection of these time intervals enables the comparison with previous works. In the case of CMIP6 models, the heat estimates obtained through direct integration of both the STPs and FTPs are depicted. The MPI-ESM P2k+ estimates (Fig. 2) are also included. All the estimates from this work and the previous literature show an increasing trend in land heat uptake towards the present. The 1950–2000 interval (Fig. 4a) shows less warming in the reanalyses than in the observational derived estimates, which is particularly noteworthy for LMRv2.1 and JRA55. These two reanalyses show no heat uptake in 1950–2000, which stems from a long-lasting negative temperature anomaly in 1950–1980 that is not shown by the rest of the sources (Fig. 3a). Furthermore, CMIP6 depth-corrected estimates using FTPs are in range with the reanalysis and slightly below the observational estimates. Hence, all estimates based on FTPs agree overall with the previous estimates, with von Schuckmann et al. (2020) giving the highest land heat uptake values across all datasets. Figure 4b and c portray similar results for the periods 1960–2020 and 1971–2018, respectively. In both cases, the estimates derived from the different FTPs are consistent, with observational estimates at the upper end and some reanalyses showing the smallest land heat uptake values. For instance, land heat uptake estimates for the period 1971–2018 (Fig. 4c) range from 10.5 to 16.0 ZJ, yielding an ensemble mean of 13.9 ZJ for the CMIP6 FTPs, 12.2 ZJ for the reanalyses, and 15.0 ZJ for the observational-driven estimates. These differences are due to the smaller temperature trends shown by reanalyses and CMIP6 models when compared to observational datasets from the 1970s to the present, as shown by Fig. 3a (see Tables 1, 2 for the exact trend values).

Furthermore, Fig. 4 shows the dependency of land heat uptake on LSM depth for CMIP6 STP-derived estimates (note gray shading on the left side), evidencing a significant increase when the forward model is applied to correct estimates stemming from shallow ESMs (from 6.9 to 12.6 ZJ in 1960–2020, Fig. 4b, and from 7.0 to 13.9 ZJ in 1971–2018, Fig. 4c). Since 1950–2000 warming is very small for most of the CMIP6 models and very sensitive to the internal variability in each one (see Table 2), correcting the estimates in this case does not make a remarkable difference (Fig. 4a). Even after correcting the LSM depth bias on land heat uptake, the range in CMIP6 multi-model variability is still very large, with a standard deviation of 5.4 ZJ in 1960–2020 and 5.1 ZJ in 1971–2018. This large variability is explained by some outlier estimates (e.g., CAMS-CSM1-0 and GISS-E2-2-G) which yield very small or even negative land heat uptakes, linked to the small climate sensitivity of these models in the periods of study, as shown in Table 2. By contrast, the CMIP6 simulations showing the greatest temperature trends also show the largest corrected FTP values, some of them yielding greater estimates than the state-of-the-art estimates from von Schuckmann et al. (2020) and Cuesta-Valero et al. (2023). In all periods shown in Fig. 4, the multi-model CMIP6 mean is consistent with estimates yielded by reanalyses and observational datasets. Despite the consistency among the estimates obtained from the different datasets in this work, all of them yield results in the range of past literature values for the 1950–2000 period (Beltrami et al., 2002; Beltrami, 2002; Huang, 2006) but are lower than those in recent literature for 1960–2020 and 1971–2018, which quantified land heat uptake to be between 18.6 and 25 ZJ in 1971–2018 (von Schuckmann et al., 2020; Forster et al., 2021) and subsequently narrowed and lowered this range to be 18.0–18.4 ZJ (Cuesta-Valero et al., 2023).

4 Conclusions and discussion

We compare land surface and soil temperatures, as well as land heat uptake, from two simulations of the period 0–2100 CE made with two versions of the MPI-ESM: one including a standard shallow and the other one including a more realistic deep LSM. The deep LSM included in MPI-ESM generates about 4 times more heat uptake than the shallow version by the end of the historical period. This underestimation by the shallow LSM version only grows up to the end of the 21st century, reaching 18 ZJ for shallow and 75 ZJ for deep LSM versions. Also, while deepening the LSM produces relatively colder near-surface soil temperatures after 1850 CE, GST low-frequency variability and trends remain unchanged. This finding suggests the potential utility of employing GSTs to drive a forward model and of obtaining deep FTPs and subsequently unbiased estimates of land heat uptake. The concurrence of estimates from both STPs of P2k+d and FTPs of P2k+s and P2k+d supports the application of the forward modeling method to correct insufficient land heat uptake values stemming from shallow LSMs. The results portray a land heat gain in the industrial period which has intensified in the last decades of the 20th century and the first 2 decades of the 21st century. The FTP strategy is also applied to estimate land heat uptake from both observation-based products and CMIP6 simulations. Results are compared with previous estimates based on observational BTPs, thus providing additional information to assess consistency between model and observational-based estimates. Land heat uptake in 1971–2018 is between 10 and 16 ZJ for all sources, being around 14 ZJ for CMIP6 models. This figure doubles the 7 ZJ that CMIP6 models capture when their heat uptake capability is constrained by imposing a shallow zero-flux bottom boundary condition, which shows the relevance of the forward model in correcting model-based estimates. Reanalysis-based estimates are slightly smaller than those obtained from applying the forward model driven by observational gridded datasets or corrected CMIP6 simulations, consistent with the generally smaller reanalyses warming since 1950 relative to other datasets. SAT observation-based estimates render the highest values, consistent with SAT trends being slightly larger than those of GST. SAT-related values can thus be regarded as an upper-limit land heat uptake estimate.

Nevertheless, all observational-based and corrected CMIP6 estimates underestimate the most recent values of land heat uptake based on BTPs. This discrepancy may stem from different sources of uncertainty. One possibility is that since borehole-derived surface temperature and ground heat flux stories have a low time resolution (e.g., 25-year time step in Cuesta-Valero et al., 2021b), they may not consider the decadal cold anomalies occurring in 1960–1980 (Fig. 3). This discrepancy may play a major role since decadal surface temperature differences might alter inter-decadal land heat imbalance due to subsurface thermal inertia (Beltrami et al., 2015). Moreover, BTP logs are scarce in the last decades and non-existent after the 2000s, so borehole-based estimates for the last 2 decades stemming from von Schuckmann et al. (2020) and Cuesta-Valero et al. (2023) for 1960–2020 and 1971–2018 are yielded by extrapolating the ground heat flux temperature reconstructed in the last decades of the 20th century to the 21st century (2000–2018 and 2000–2020, respectively). This could be overlooking the global warming hiatus (1998–2012; Karl et al., 2015) effect on a relative slowdown in land heat uptake. Furthermore, the different methodological approaches used to yield observational and simulation-based (STP and FTP) heat uptake estimates might play a role. While STP- and FTP-based heat uptake values in a certain time period are rendered by integrating the differences in the temperature profiles between the initial and the final time steps of that period, previous observational land heat uptake values are yielded by estimating ground heat flux stories based on single-time BTPs at different locations to then reconstruct the global heat uptake evolution (Beltrami et al., 2002; Cuesta-Valero et al., 2023). This entails BTP-based estimates not representing net land heat variations in a certain time span but rather representing a continuous variation in land heat uptake by aggregating logs taken at different dates and locations. This subsequently implies that this methodology preserves land energy changes due to past climate signals (Beltrami et al., 2017), which may contribute to comparatively greater heat uptake estimates in recent decades. With these differences reconciled, this work contributes to bringing observational and model-based land heat uptake estimates to the same frame of comparison, by removing the bias associated with imposing LSM depths that are too shallow. The comparison between CMIP6-based and observational data-driven land heat uptake estimates can also serve as an additional metric to assess the realism of ESMs in reproducing current global warming.