Underestimated Land Heat Uptake Alters the Global Energy Distribution in CMIP6 Climate Models

Current global warming results in an uptake of heat by the Earth system, which is distributed among the different components of the climate system. However, current‐generation climate models deliver heat inventory and partitioning estimates of Earth system components that differ from recent observations. Here we investigate the global heat distribution under warming by using fully‐coupled CMIP6 Earth system model experiments, including a version of the MPI‐ESM with a deep land model component, accommodating the required space for more realistic terrestrial heat storage. The results show that sufficiently deep land models exert increased subsurface land heat uptake, leading to a heat uptake partitioning among the Earth system components that is closer to observational estimates. The results are relevant for the understanding of Earth's heat partitioning and highlight the importance of the land heat sink in the Earth heat inventory.


Introduction
"Where does the energy go?" is a long-standing question in climate research that addresses the amount of heat uptake by the different components of the Earth system as a response to a positive global energy imbalance (von Schuckmann et al., 2020(von Schuckmann et al., , 2023)).Currently, the Earth system is in a positive radiative imbalance due to anthropogenic emissions of greenhouse gases (Hansen et al., 2005(Hansen et al., , 2011)), which gradually increases its heat storage.For the period of 1971-2020, 381 ± 61 ZJ were observed to be stored in the Earth system, accounting for ocean (89%), land (6%) and atmosphere (1%) warming, and cryosphere melting (4%) (Cuesta-Valero et al., 2016, 2021, 2023;von Schuckmann et al., 2020von Schuckmann et al., , 2023)).
Determining the distribution of this heat among the climate system components is important because the heat storage of different Earth system components and the transfer between them largely dictates the state and variations of the global climate (Hansen et al., 2011).For example, heat uptake within these subsystems drives variations in thermodynamic and hydrological cycles (Rhein et al., 2013) and has an impact on the global carbon cycle, ecosystem functioning, agricultural production, environment, and society (Gornall et al., 2010).
Earth system models (ESMs) are employed to estimate the evolution of Earth's heat uptake and its distribution in the future.Cuesta-Valero et al. (2021) provide the first complete assessment of the heat distribution in ESMs participating in the Climate Model Intercomparison Project (CMIP) Phase 5.They find a fractional ocean heat uptake of 96% ± 4%, land warming of 2% ± 3%, and accumulation of heat in the atmosphere and cryosphere at 1% ± 1%, respectively, of the total system heat uptake.These results represent an overestimation of the amount of ocean heat storage and an underestimation of heat uptake by the cryosphere and land components in comparison to the observational estimates (Cuesta-Valero et al., 2021).Previous literature has suggested that land heat uptake under climate change conditions in ESMs is compromised by too shallow land surface model components (LSMs) because a subsurface bottom boundary condition too close to the surface biases the representation of subsurface heat propagation and heat distribution (Alexeev et al., 2007;Cuesta-Valero et al., 2016, 2021;González-Rouco et al., 2009, 2021;MacDougall et al., 2008;Smerdon & Stieglitz, 2006;Steinert et al., 2024;Steinert, González-Rouco, de Vrese, et al., 2021;Steinert, González-Rouco, Melo-Aguilar, et al., 2021;Stevens et al., 2007).Despite recent improvements in modeling land processes in ESMs (Fisher & Koven, 2020), only limited attention has been directed toward the effect of land model depth and its impact on the representation of terrestrial thermodynamics.Therefore, the bulk of the current generation of LSMs within ESMs still has depths that range between 3 and 14 m (Burke et al., 2020;Cuesta-Valero et al., 2016), allocating limited space for subsurface heat propagation and hence, land-climate feedbacks (González-Rouco et al., 2009).The exceptions are ESMs using the Community Land Model (Lawrence et al., 2019;Oleson et al., 2013) and the ORCHIDEE model (Cheruy et al., 2020), extending their bottom boundaries down to around 45 and 90 m, respectively, which leads to a higher and more realistic land heat storage (Cuesta-Valero et al., 2021).
While the impacts of realistically deep land components have been investigated with standalone LSMs and analytical models, they have not yet been considered in the frame of ESM coupled simulations (Phillips, 2020).Here, we hypothesize that ESMs with deeper land model components not only show larger land heat storage but are consequently adjusting the distribution of excess heat among the climate system components toward fractions closer to observational estimates.We therefore explore the influence of a deepened land model component considering an ensemble of existing CMIP6 climate model simulations and support our analysis by a currentgeneration ESM that we modified to encompasses both a deep and a shallow land model setup.We consider the land model depth of CMIP6 ESMs in the assessment of the Earth system heat inventory and investigate the effects of a more realistic land heat sink on the global heat distribution.

CMIP6 Ensemble
Simulations of 33 CMIP6 ESMs are considered for our analysis of the Earth energy inventory (Table S1 in Supporting Information S1).Their land model depth varies between 3 and 90 m, with 61% lying between 3 and 14 m.Some ESMs use an identical land model coupled with different ocean or atmosphere components.Only one ensemble member per model is considered.For our analysis, we use the historical (1850-2014) and the future (2015-2100) SSP5-85 forcing scenario periods (Gidden et al., 2019;Meinshausen et al., 2017).Additionally, the climate trends of all available years of the piControl simulations were used to detrend the transient CMIP6 climate projections.This is a necessary step to remove spurious drifts from incomplete spin-up processes in global climate models because ESMs require millennia under stable conditions to reach an equilibrium state, mostly due to the long timescales of the global ocean circulation (Gupta et al., 2013;Peacock & Maltrud, 2006;Séférian et al., 2016).Previous analyses have shown that this dedrifting technique mitigates the non-conservation terms in the simulated Earth heat inventory within transient climate simulations (Hobbs et al., 2016;Irving et al., 2020), which is crucial for our analysis.For the calculation of the Earth system heat uptake and distribution (see Section 2.4), we use a total of 20 output variables provided by the CMIP6 protocol that are summarized in Table S2 of Supporting Information S1.

MPI-ESM Ensemble
CMIP6 models provide outputs for the LSM component with a single depth configuration.The large inter-model variability in global heat uptake and non-conservation of energy within the models (Irving et al., 2020;Wilde & Mulholland, 2020) is thus challenging to attribute to differences in land model depth.Therefore, we support our multi-model ensemble with coupled simulations of the Max Planck Institute for Meteorology Earth System Model (MPI-ESM; Mauritsen et al., 2019) with two different LSM setups.The first setup equals the standard CMIP6 version of MPI-ESM, which uses the JSBACH LSM (Reick et al., 2021) with five layers (9.83 m; SHALLOW hereafter).For the second setup, we introduce modifications to the vertical structure of JSBACH, which uses an extended 18-layer discretization with a larger model depth (1,417 m; DEEP hereafter) that is well suited for the simulation of centennial to millennial timescales (González-Rouco et al., 2021;Steinert et al., 2024;Steinert, González-Rouco, de Vrese, et al., 2021;Steinert, González-Rouco, Melo-Aguilar, et al., 2021).The two simulations were prepared by restarting the ocean component MPIOM from an existing CMIP6 control simulation.Prior to the latter, a spin-up run was conducted over several centuries starting from the Polar Science Center ocean hydrographic climatology (Steele et al., 2001).Using an identical fully coupled model with two different depth configurations reduces the uncertainties associated with model differences in simulating the global energy budget.

Heat Conduction Forward Model
The dependence of the land heat uptake on the depth of the LSMs is evaluated using an offline one-dimensional numerical heat conduction forward model (HFMod; García-Pereira et al., 2024).This model numerically discretizes the vertical heat conduction equation (Carslaw & Jaeger, 1959): where T represents the temperature at a specific depth z and time step t, and κ denotes the thermal diffusivity.The numerical solution is attained using a Finite Time Centered Space finite difference scheme, consisting of an explicit Euler scheme for time resolution and a second-order central difference for the spatial dimension (Roache, 1998).The resulting numerical solution of the heat conduction equation is expressed as: with T n i representing the temperature value at soil level i during time step n.The model produces the subsurface vertical temperature structure at each time step, employing a thermal equilibrium boundary condition at the ground surface, that is, T[z = 0] = T 0 , and a zero-flux boundary condition at the bottom, that is, T[z = z max ] = 0. To ensure the stability of the numerical scheme, governed by the Courant-Friedrichs-Lewy condition Δt ≤ Δx 2 /2κ (Courant et al., 1928), a constant time step of 1 day and a soil layer thickness of 1 m was employed.The zero-flux bottom boundary condition was imposed at an infinite depth (HFMod infinite solution) and specific depths (HFMod limited solution).

Earth Energy Calculations
Earth system heat content estimates were obtained from CMIP6 outputs following the methods described in Cuesta-Valero et al. (2016, 2021).The data used for these calculations is summarized in Table S2 of Supporting Information S1.Ocean heat content (OHC) integrates the simulated temperature and salinity profiles using the algorithm for potential enthalpy (Griffies, 2004;McDougall, 2003).To convert the retrieved potential enthalpy into heat content, seawater density and pressure are estimated for each ocean layer as functions of depth, temperature, and salinity (McDougall et al., 2003;Smith et al., 2010).
Land heat content (LHC) is estimated by integrating the simulated subsurface temperature and water profiles for each model.Volumetric heat capacity is derived from sand, clay, water content, and ice content at each model layer (Oleson et al., 2010;Van Wijk et al., 1963), with sand and clay data retrieved from the ECOCLIMAP project (Champeaux et al., 2005).Simulated soil water and ice content at each layer are classified as water if the temperature of the layer is 0°C or above and ice otherwise.For simulations providing only column-integrated soil data, the water and ice contents are distributed proportionally within the soil layers above the bedrock and then classified as water or ice as in the layer-resolved cases (Cuesta-Valero et al., 2016, 2021).The thermal properties of bedrock are considered to start at 3.8 m and expand until the bottom layer of the model for all models to have a common depth reference, even for models that do not define bedrock (including the MPI-ESM simulations used herein).The volumetric heat capacity for bedrock is assumed to be ρ C = 3.0 • 10 6 Jm 3 K 1 (Cuesta-Valero et al., 2016, 2021;MacDougall et al., 2008MacDougall et al., , 2010)).
Several variables are used to derive atmospheric heat content (AHC).The simulated atmospheric profiles of temperature, wind, and specific humidity are integrated following the theory in Trenberth (1997) and the practical method described in Previdi et al. (2015), in combination with simulated values of surface pressure and geopotential (Dutton, 2002).The specific heat capacity and latent heat of vapourization are considered constant for the entire atmospheric profile, with values of c p = 1,000 Jkg 1 K 1 and L v = 2,260 Jkg 1 , respectively.
Regarding the cryosphere heat content (CHC), simulated changes in ice masses are considered, with model outputs allowing us to estimate ice melt in different components.All but three models provide enough variables to estimate changes in sea ice mass, with a varying number of models presenting enough variables to estimate changes in subsurface ice and snow mass (Table S2 in Supporting Information S1).None of the models explicitly simulate changes in glacier extent or the ice caps in Greenland and Antarctica, a known limitation of CMIP models (Cuesta-Valero et al., 2021;Flato et al., 2013).The heat content in the different simulated components of the cryosphere is estimated as the change of mass multiplied by the latent heat of fusion L f = 3.34 • 10 5 Jkg 1 (Rhein et al., 2013).

Land Heat Uptake Across CMIP6 Models
The land surface temperature (here defined as Ground Surface Temperature-the first model layer) in the CMIP6 ensemble shows a pronounced increase over the 20th and 21st centuries, following the imposed socio-economic pathway SSP5-85 (Figure 1a).The spread among the 33 participants ranges between about 4 and 10 K at the end  S1 in Supporting Information S1) for the period 1850-2100 relative to pre-industrial conditions.The blue and gray scales stand for model ECS in (a) and model depth in (b), respectively.(c) Global average cumulative land heat uptake (land heat uptake sensitivity; excluding glaciers) in relation to equilibrium climate sensitivity (bubble size; Table S1 in Supporting Information S1).Bubble colors in panels (c) and (d) refer to the respective model depths as indicated in the legend in panel (b).(d) Global average cumulative land (excluding glaciers) heat uptake per degree of global land surface temperature change (land uptake efficacy) with land model depth (Table S1 in Supporting Information S1) for the period of 1850-2100 relative to pre-industrial conditions.Shading indicate the range of results from the heat conduction forward model with the soil lines representing values of 1.0 • 10 6 m 2 s 1 for the thermal diffusivity, and 3.0 • 10 6 Jm 3 K 1 for the volumetric heat capacity (Figure S2 in Supporting Information S1).
of the 21st century.With the land heat uptake directly related to the increase in surface temperature, the former is expected to increase accordingly.However, the land energy uptake among the CMIP6 models show a surprisingly large spread of 6-143 ZJ by 2100 (Figure 1b).Higher land surface temperatures do not lead to a proportional increase in land heat uptake in many CMIP6 models (Figure 1c).This is caused by differences in land model depth (Figure 1b).Equilibrium climate sensitivity (ECS; Meehl et al., 2020;Schlund et al., 2020;Zelinka et al., 2020, Table S1 in Supporting Information S1) can be excluded to cause the range of land heat uptake for the different CMIP6 models because the spread persists when accounting for the ECS (Figure 1c; bubble size).The modeldepth dependency of the land heat uptake leads to a notable separation into two groups in the CMIP6 land heat content (Figures 1b and 1c).This separation motivates a split of the CMIP6 ensemble into shallow and deep model groups according to their respective land model depth (Table S1 in Supporting Information S1).Shallower models between 3 and 14 m depth (Figure 1c, orange) show no significant relationship between land heat uptake and surface warming (R 2 = 0.00).These models are limited to 31 ZJ (6-93 ZJ; indicating minimum and maximum values) heat uptake between 1850 and 2100, despite some of them simulating up to 9.2 K of global mean temperature increase.Contrarily, the group of deeper models between 42 and 90 m (Figure 1c, turquoise) correlate linearly with the amount of surface warming (R 2 = 0.86), with heat uptake of 101 ZJ .The two extra MPI-ESM setups fall well into these two depth categories (Figure 1c).As MPI-ESM has a small ECS (Meehl et al., 2020;Schlund et al., 2020;Zelinka et al., 2020), the amount of warming and respective energy uptake is on the lower end of the CMIP6 ensemble.
To remove the influence of differences in land surface warming across ESMs, we define land heat uptake efficacy (LHUe) to be the land heat uptake divided by the surface warming (Figure 1d).The LHUe shows large differences between the shallow and the deep CMIP6 model groups.Generally, deeper land models equal larger land energy uptake.However, even the deepest CMIP6 models are shallower than 170 m, which was identified as the recommended depth to fully resolve the propagation of the surface temperature signal in historical and climate change scenario experiments (Steinert, González-Rouco, Melo-Aguilar, et al., 2021).This discrepancy suggests that the deep CMIP6 models are not saturated in land heat uptake and would take up even more energy if they were deeper.We employ the HFMod heat conduction forward model (see Section 2.3) that allows us to define the relationship between LHUe and land model depth in limited and infinitely deep land model configurations (Figure 1d).Since there is no official information on soil parameters for all CMIP6 models, we run the HFMod with parameter values for volumetric heat capacity and thermal diffusivity (Figure S1 in Supporting Information S1) over a range of common values of 0.5-1.5 • 10 6 m 2 s 1 for the thermal diffusivity, and 2.5-3.5 • 10 6 Jm 3 K 1 for the volumetric heat capacity (e.g., Cuesta-Valero et al., 2023).The results (Figure 1d; shading) illustrate that deviations in the effectiveness of heat conduction across models can govern a range of LHUe that captures the values of the deep CMIP6 model group.With a larger LSM depth and thus increased land heat uptake capacity, MPI-ESM SHALLOW and DEEP simulations have the same LHUe as the deep CMIP6 model average.MPI-ESM DEEP has slightly lower LHUe than the HFMod solution, consistent with lower global average values of 0.79 • 10 6 m 2 s 1 for the thermal diffusivity and 2.0 • 10 6 Jm 3 K 1 for the volumetric heat capacity for MPI-ESM (Steinert, González-Rouco, Melo-Aguilar, et al., 2021).Any deviations of the CMIP6 model results from the solution of the heat conduction forward model might be additionally influenced by their consideration of soil moisture for the land heat uptake calculations, whereas the heat conduction forward model only considers thermal diffusion.The results confirm that LHUe saturates at 17.9 ZJ/K with a depth of about 200 m, where the diffusion of surface temperature anomalies with depth is attenuated completely.

Redistribution of Earth System Energy
Given the more realistic land heat uptake in simulations with a deeper land model, the question is how this affects the overall Earth energy budget and distribution among the climate system components.In the case where the energy is underestimated in the land system in the shallow LSM configurations, the assumption is that this energy is incorrectly distributed in other Earth system components, namely the ocean, atmosphere, and cryosphere.Across all climate system components, the spread of the subsystem cumulative heat uptake between 1850 and 2100 in the CMIP6 models is relatively large (Figure 2).As expected, the ocean heat uptake is largest with 2828 ZJ (1609-3621 ZJ).As stated above, the land heat storage shows the largest uncertainty as a consequence of the depth limitations of some models and stores 56 ZJ (6-143 ZJ) across all models.Various ESMs match historical ocean heat uptake well, while cross-model uncertainty remains large (Table S3 in Supporting Information S1).The differences between the shallow and deep model configurations are most prominent in the land heat uptake.MPI-ESM SHALLOW and DEEP align well with the CMIP6 averages, but the low ECS of MPI-ESM prevents more intense heat accumulation despite its LSM depth being sufficiently larger.Nonetheless, both the deep CMIP6 model group and MPI-ESM DEEP match observational estimates (added as relative changes to the CMIP6 ensemble mean in 1960) between 1960 and 2020 (Figure 2; von Schuckmann et al., 2023) with a difference of 6.6 and 2.9 ZJ at 2020, respectively, whereas the shallow models show considerably lower values with a difference of 17.4 and 16.9 ZJ at 2020, respectively.Here, the CMIP6 shallow model range of relative cumulative heat uptake by 2020 encompasses 0.2-20.4ZJ, while the CMIP6 deep models range 13.9-23.2ZJ-in comparison to the observed 23.9 ZJ.The atmosphere and cryosphere show the smallest heat uptake of 25 ZJ (15-37 ZJ) and 21 ZJ (9-28 ZJ), respectively, between 1850 and 2100.The spread of results for the cryosphere is partly caused by the lack of data availability (see Table S2 in Supporting Information S1), with some model output repositories missing variables to calculate the energy storage in some of the subcomponents of the cryosphere, likely causing the inconsistency with the observations.For example, the CMIP6 version of MPI-ESM does not include soil freezing and thus does not consider the contribution of energy from soil ice melt.Additionally, many models lack a dynamic representation of glaciers and ice caps, thus not including the heat exchange from their melting or growth.Although this can be considered a limitation of our analysis, the small contribution of the cryosphere heat uptake will likely not change the results considerably.
Solely considering absolute changes in the energy inventory can lead to inconsistencies in the effect of land model depth-related changes on the global heat distribution (Figures S3 and S4 in Supporting Information S1) due to model structural differences, internal variability, or model drift and the associated magnitude of energy conservation in many CMIP6 models (Irving et al., 2020).We therefore consider the heat distribution of climate subsystems as a fraction of the total Earth Heat Content (EHC) instead of the top of the atmosphere (TOA) energy imbalance (Cuesta-Valero et al., 2021).Due to the reported inconsistencies between the global energy budget and conservation (Irving et al., 2020;Wilde & Mulholland, 2020), we consider the EHC as the sum of heat contents in all climate system components.
For the shallow CMIP6 model group, our calculations find a fractional heat uptake distribution of 97.4% for the ocean, 1.1% for the land, and 0.8% and 0.6% for the atmosphere and the cryosphere, respectively (Figures 3a-3d).Note that only a subset of model experiments from the original CMIP6 ensemble is used here (17/20 shallow, and 9/13 deep models) based on data availability and excluding outliers that exceeded two standard deviations of heat  S1 in Supporting Information S1) and the two MPI-ESM setups DEEP and SHALLOW for (a) ocean, (b) land, (c) atmosphere, (d) cryosphere (sea ice, land ice, snow), (e) the Earth Heat Content (EHC), which is the sum of all components in (a-d), and (f) integrated top of the atmosphere (TOA) energy imbalance in simulations between 1850 and 2100 following SSP5-85 forcing.For panel (e), EHC is only shown for a subset of models that allowed the calculation of all heat uptake components.For comparison, observational estimates (purple) from von Schuckmann et al. (2023) were added between 1960 and 2020.RMSE values are provided between the model-derived and observational estimates.Note the differences in y-axis between the panels.uptake percentages from the mean.The heat uptake partitioning highlights the limitations of a compromised land heat sink deviating substantially from the expected 5%-6% (von Schuckmann et al., 2020(von Schuckmann et al., , 2023)).The comparison between the deep and shallow CMIP6 model groups shows a considerable increase in percentages of the land heat uptake from 1.1% to 3.6% for the period 1850-2100.While the energy budgets of the atmosphere and the cryosphere stay relatively the same, the ocean loses heat fraction in the deep CMIP6 model group by about 1.5%.Fractional ocean (land) heat uptake ranges 95.5%-98.5% (0.3%-2.7%) and 93.5%-98.9%(2.6%-5.4%)for the 5th-95th percentiles in the shallow and deep CMIP sub-ensemble, respectively.The mean differences between the sub-ensembles of shallow and deep ESMs are statistically significant for the land and ocean heat uptake (t-test with p-value <0.05).These results are consistent with the heat content changes obtained by comparison of the MPI-ESM SHALLOW and DEEP simulations.Land and ocean heat uptake changes show similar magnitude in the DEEP model by about 2.0% (changing from 1.0% to 3.0%) and 2.0% (changing from 98.0% to 96.0%), respectively.In contrast to the CMIP6 ensemble, the DEEP MPI-ESM shows no decrease of heat in the cryosphere, likely due to the limitation in the parameterization of cryospheric elements in MPI-ESM.Land heat uptake remains slightly lower than the observational estimates, likely due to the dominating ocean heat uptake in both the CMIP6 and MPI-ESM ensembles, but consistent with the uncertainty in the observational estimate (Cuesta-Valero et al., 2021;von Schuckmann et al., 2023).Compared to the EHC-based estimates, taking the TOA imbalance as a reference for the global heat uptake between the shallow and deep CMIP6 model groups provides very similar results for the direction and magnitude of change of the land, atmosphere and cryosphere heat uptake fractions, while the changes of the ocean heat uptake are smaller (Figure S2 in Supporting Information S1).For TOA-based estimates, only land heat uptake change is statistically significant.We note that these results are consistent with the analysis of CMIP5 models over the historical period (Cuesta-Valero et al., 2021), indicating that, with respect to the main driver of land heat uptake -the land model depth, there has been no apparent improvement between the last and the current iteration of CMIP models (Cuesta-Valero et al., 2016).
Considering the global heat distribution over time also allows the assessment of temporal changes of the relative Earth system heat distribution (Figure 3e).The results show heat content percentages that are relatively constant over time.Larger variations at the beginning of the period are expected as the absolute amount of heat exchange in the Earth system is relatively small and can generate larger fractional changes.As the absolute Earth system heat content increases in the 21st century under the SSP5-85 scenario, the relative heat uptake of the ocean and land stabilizes.Given this stabilization of the global heat fractions, it is expected that the present-day ratios of global  S1 in Supporting Information S1) for the period 1850-2100, and (e) the temporal evolution of the land and ocean heat uptake fractions [% of EHC] over the 21st century.

Geophysical Research Letters
10.1029/2023GL107613 heat uptake will not change substantially by the end of the 21st century.As most of the cryosphere is frozen water, its heat is expected to be transferred to the ocean and the atmosphere as part of a changing hydrological cycle under warming.

Summary and Conclusions
This study evaluates estimates of the land heat uptake from climate models with deep and shallow land model components under transient climate conditions and assesses their resulting influence on the global heat distribution.Since the average of the current generation of CMIP6 models does not realistically account for the land heat content in the Earth system (Cuesta-Valero et al., 2021;von Schuckmann et al., 2020), we separate them into two groups of shallow and deep models according to their land model depth and compare them to a set of simulations with MPI-ESM that uses two corresponding setups of a shallow and a deep land model component.
A clear influence of land model depth on the Earth heat inventory in Earth system models has been established.Deeper climate models generate Earth system heat distribution more consistent with observational and numerical estimates than the group of shallower models.The deeper group of CMIP6 models with depths between 42 and 90 m has a more realistic land heat uptake with improvements over the shallow model group.Furthermore, comparing shallow and deep CMIP6 land models with the standard shallow and a sufficiently deep MPI-ESM setups facilitates direct insight into the changes in the Earth system heat distribution.The results illustrate that the ocean can lose a fraction of the global heat budget to the land due to the increased land heat sink capacity.In both the CMIP6 models and MPI-ESM, the impacts of changes in the land heat uptake on the global heat inventory are relatively small because of the large natural variations of ocean heat content given its massive mass, thermal capacity and heat distribution capability through water advection (Levitus et al., 2000).The results for the ocean and land from climate models with deeper land model components are in closer agreement to the most recent observational estimates (von Schuckmann et al., 2023) than those with shallow components.Changes in the heat uptake fraction of the atmosphere and cryosphere are negligible and are partially sensitive to the representation of dynamic cryospheric elements in current-generation ESMs.Throughout the 21st century, the ratios of heat uptake between the climate subsystems stay relatively constant, even for a strong climate change scenario.Because energetic model biases in absolute terms tend to be larger between models than between intra-model climate subsystems (Cuesta-Valero et al., 2016, 2021;Irving et al., 2020;Wilde & Mulholland, 2020), individual impacts of a misrepresented energy distribution in the climate system may vary between models.Nevertheless, Earth energy distribution differences of a few percent of energy uptake may influence the terrestrial carbon cycle, hydro-thermodynamic coupling and biogeophysical properties of the land surface (Cuesta-Valero et al., 2023;Gulev et al., 2021).
Our results confirm that deeper models across CMIP6 are necessary to have a realistic representation of the land heat content and the corresponding subsurface thermodynamic state.We conclude that around two-thirds of CMIP6 models are compromised by an insufficiently deep land model depth.Further, we argue that even the deeper CMIP6 models are still not deep enough to capture the full potential for land heat storage since analytical solutions suggest deeper models of at least 170 m to be necessary (Steinert, González-Rouco, Melo-Aguilar, et al., 2021).Their land heat uptake efficacy comparable to the heat conduction model solutions likely stems from uncertainties in the consideration of soil thermal properties.Additional computational requirements from deeper LSMs may add up with the use of high-density layer discretization and increasing process fidelity, but generally remains small compared to the computational costs of other Earth system components (González-Rouco et al., 2021;Steinert, González-Rouco, Melo-Aguilar, et al., 2021).Overall, these results highlight the relevance of the land heat sink in simulating the global energy inventory under global warming.

Data Availability Statement
The calculation of the Earth heat inventories is based on CMIP6 data that is publicly available on servers provided by ESGF (https://esgf-node.llnl.gov/projects/cmip6/)and following the methods described in Cuesta-Valero et al. (2021).Observational estimates for the heat uptake of all Earth system component were taken from von Schuckmann et al. ( 2023) accessible under https://www.wdc-climate.de/ui/entry?acronym=GCOS_EHI_1960-2020.Data of the modified MPI-ESM experiments used herein are available from Zenodo (Steinert et al., 2024).

Figure 1 .
Figure 1.(a) Global average land surface temperature (excluding glaciers) and (b) land heat uptake in 33 ESMs participating in CMIP6 (TableS1in Supporting Information S1) for the period 1850-2100 relative to pre-industrial conditions.The blue and gray scales stand for model ECS in (a) and model depth in (b), respectively.(c) Global average cumulative land heat uptake (land heat uptake sensitivity; excluding glaciers) in relation to equilibrium climate sensitivity (bubble size; TableS1in Supporting Information S1).Bubble colors in panels (c) and (d) refer to the respective model depths as indicated in the legend in panel (b).(d) Global average cumulative land (excluding glaciers) heat uptake per degree of global land surface temperature change (land uptake efficacy) with land model depth (TableS1in Supporting Information S1) for the period of 1850-2100 relative to pre-industrial conditions.Shading indicate the range of results from the heat conduction forward model with the soil lines representing values of 1.0 • 10 6 m 2 s 1 for the thermal diffusivity, and 3.0 • 10 6 Jm 3 K 1 for the volumetric heat capacity (FigureS2in Supporting Information S1).

Figure 2 .
Figure 2. Cumulative heat uptake [ZJ] of 33 CMIP6 models (see TableS1in Supporting Information S1) and the two MPI-ESM setups DEEP and SHALLOW for (a) ocean, (b) land, (c) atmosphere, (d) cryosphere (sea ice, land ice, snow), (e) the Earth Heat Content (EHC), which is the sum of all components in (a-d), and (f) integrated top of the atmosphere (TOA) energy imbalance in simulations between 1850 and 2100 following SSP5-85 forcing.For panel (e), EHC is only shown for a subset of models that allowed the calculation of all heat uptake components.For comparison, observational estimates (purple) from vonSchuckmann et al. (2023) were added between 1960 and 2020.RMSE values are provided between the model-derived and observational estimates.Note the differences in y-axis between the panels.

Figure 3 .
Figure3.Cumulative heat uptake for (a-d) global heat uptake in reference to EHC [% of EHC] of land, ocean, atmosphere, cryosphere (sea ice, land ice, snow), respectively, separated into shallow and deep model groups according to their land model depth (TableS1in Supporting Information S1) for the period 1850-2100, and (e) the temporal evolution of the land and ocean heat uptake fractions [% of EHC] over the 21st century.