First Assessment of the Earth Heat Inventory Within CMIP5 Historical Simulations

Historical Simulations Francisco José Cuesta-Valero1,2, Almudena García-García1,3, Hugo Beltrami1, and Joel Finnis4 1Climate & Atmospheric Sciences Institute, St. Francis Xavier University, Antigonish, NS, Canada. 2Environmental Sciences Program, Memorial University of Newfoundland, St. John’s, NL, Canada. 3Department of Remote Sensing, Helmholtz Centre for Environmental Research-UFZ, Leipzig, Germany. 4Department of Geography, Memorial University of Newfoundland, St. John’s, NL, Canada. Correspondence: Francisco José Cuesta-Valero (fcuestav@stfx.ca)

analyses included in von Schuckmann et al. (2020). Changes in sea ice volume in Church et al. (2011) are obtained from Levitus et al. (2005), and from the Pan-Arctic Ice Ocean Modeling and Assimilation System (PIOMAS, Zhang and Rothrock, 2003;Schweiger et al., 2019) in the case of von Schuckmann et al. (2020). All changes in ice mass are multiplied by the latent heat of fusion in order to obtain the corresponding estimate of cryosphere heat content.
Global averages of Ocean Heat Content (OHC), the heat content within the continental subsurface (ground heat content, 95 GHC), Atmosphere Heat Content (AHC) and heat uptake by ice masses (cryosphere heat content, CHC) were derived from the CMIP5 Historical experiments. The OHC values were estimated using the formulation for potential enthalpy described in McDougall (2003) and Griffies (2004) from simulated seawater potential temperature and salinity profiles ( Table 1 contains the list of variables employed for estimating each term of the EHI). Once the potential enthalpy has been determined, estimates of seawater density  and pressure profiles (Smith et al., 2010) allowed simulated heat content in the 100 ocean to be calculated as: where Q Ocean is the ocean heat per surface unit (in J m −2 ), S is salinity (in psu), θ is potential temperature (in • C), p is pressure (in bar), and z i , ρ i , H • i and ∆z i are depth (in m), density (in kg m −3 , potential enthalpy (in J kg −1 ) and thickness (in m) of the i-th ocean layer, respectively. This approach is based on the availability of both temperature and salinity profiles in 105 CMIP5 simulations, which allows to integrate changes in water density. Estimates of OHC from observational methods only consider temperature profiles, as salinity profiles are not routinely measured at the global scale. However, CMIP5 simulations yield similar changes in OHC from both methods ( Figure S1). Thus, we use the method described by Equation 1 to estimate OHC from simulations, since this approach includes simulated salinity profiles in the analysis, maximizing the information considered to estimate heat content. 110 The GHC series were estimated as in Cuesta-Valero et al. (2016) for all terrestrial grid cells. Subsurface thermal properties were computed taking into account spatial variations in soil composition (% of sand, clay and bedrock) and simulated subsurface water and ice amounts (Van Wijk et al., 1963;Oleson et al., 2010). The subsurface temperature profile was then integrated 115 where Q Ground is the subsurface heat storage per surface unit (in J m −2 ), and ρC i , T i and ∆z i are the volumetric heat capacity (in J m −3 K −1 ), the temperature (in K) and the thickness (in m) of the i-th soil layer, respectively. All CMIP5 GCMs present outputs for subsurface temperature, but not all models provide outputs for subsurface water and ice content in the same format (Table 1), hampering the estimate of thermal properties (ρC) in Equation 2. Indeed, two thirds of the GCMs provide the joint content of water and ice for each soil layer (mrlsl variable in CMIP5 notation), while the remaining third provide the total 120 water and ice content in the entire soil column (mrso variable). As in Cuesta-Valero et al. (2016), we considered water to be frozen in layers with temperatures below 0 • C and liquid water otherwise for models providing the mrlsl variable. For models providing the mrso variable, we distributed the water and ice content among the soil layers proportionally with layer thickness, considering ice in soil layers with temperature below 0 • C and liquid water otherwise.
The AHC series from CMIP5 simulations were estimated using the theoretical foundations of Trenberth (1997) and Previdi 125 et al. (2015). The simulated air temperature profile was integrated for all atmospheric grid cells together with estimates of wind kinetic energy, latent heat of vaporization and surface geopotential, which was determined as in Dutton (2002). Vertical atmospheric profiles were integrated in pressure coordinates as: where Q Atmosphere is atmospheric heat per surface unit (in J m −2 ), g is apparent acceleration due to gravity (in m s −2 ), p s is 130 surface pressure (in Pa), c p = 1000 J kg −1 K −1 is the specific heat of air at constant pressure, L = 2260 J kg −1 is the latent heat of vaporization, Φ s is the surface geopotential estimated from orography (in m 2 s −2 ), and T i , k i , q i and ∆p i are the air temperature (in K), specific kinetic energy (in J kg), specific humidity (in kg kg −1 ) and thickness (in Pa) of the i-th atmospheric layer, respectively.
For estimating the CHC series, the simulated cryosphere was divided into three terms: sea ice, subsurface ice and glaciers.

135
Variations in the mass of simulated sea ice and subsurface ice were multiplied by the latent heat of fusion (L f = 3.34 × 10 5 J kg −1 , Rhein et al., 2013) to obtain the heat absorbed in the melting process. The same method was applied to the change in snow mass in grid cells containing land ice within each CMIP5 GCM (glaciers or ice sheets, sftgif variable in the CMIP5 archive). Although this is not a satisfactory approach given the differences between snow and land ice, it is the only available approximation since CMIP5 GCMs do not typically represent terrestrial ice masses (Flato et al., 2013). Therefore, 140 the cryosphere heat content was estimated as where Q Cryosphere is absorbed heat per surface unit (in J m −2 ), ρ = 920 kg m −3 is ice density (Rhein et al., 2013), ∆ω is the change in subsurface ice mass per surface unit (in kg m −2 ), ∆p is the change in the proportion of sea ice at each ocean grid cell, ∆z is the change in thickness of sea ice at each ocean grid cell (in m), and ∆Ω is the change in snow amount at each cell 145 containing land ice (in kg m −2 ). It is important to note that nine of the CMIP5 GCMs did not provide outputs for the subsurface ice amount (mrfso variable) and that three of the models did not provide outputs for snow amount (snw variable, see  (Rhein et al., 155 2013;Palmer and McNeall, 2014;Trenberth et al., 2014;von Schuckmann et al., 2016). That is, if a model does not produce artificial sources or leakages of energy or mass (i.e., if the model conserves the total heat content in the system), the change in N and in EHC should be almost identical (Hobbs et al., 2016). Nevertheless, CMIP5 GCM simulations are prone to drift, particularly the ocean component due to incomplete model spin-up procedures (Sen Gupta et al., 2013;Séférian et al., 2016).
For this reason, potential drifts in estimates of heat content and the components of the radiative budget at the top of the atmo-160 sphere were removed by subtracting the linear trend of the corresponding preindustrial control simulation from the Historical simulations, which should correct artificial drifts in the simulated heat content within each climate subsystem (Hobbs et al., 2016). N estimates from the CESM1-CAM5 GCM constitute a particular case, since an unrealistic trend remained in the Historical experiment in comparison with other CMIP5 GCMs after removing the drift using data from the corresponding control simulation ( Figure S2). The rest of variables from this GCM were dedrifted using the trend estimated from the preindustrial 165 control simulation as in the other CMIP5 simulations, but the drift in the outgoing shortwave radiation and the outgoing longwave radiation at the top of the atmosphere could not be removed. Therefore, we used the trend estimated from the first five decades of the Historical simulation (1861( -1911 to remove the drift in N estimates, achieving a better comparison with the other CMIP5 GCMs ( Figure S2).
As a complement to the estimates of the EHI detailed above, we also estimated the partition of the simulated total heat content 170 among the ocean, the continental subsurface, the atmosphere and the cryosphere. A linear regression analysis was performed between the evolution of the simulated heat storage within each climate subsystem and the estimates of total heat content in the entire climate system to determine the partition of heat within the four climate subsystems (Figure 1). The slope of the linear fit was assumed to represent the simulated proportion of heat in the corresponding subsystem, thus providing estimates of OHC/N and OHC/EHC for the simulated proportion of heat in the ocean, GHC/N and GHC/EHC for the simulated proportion of heat 175 in the continental subsurface, AHC/N and AHC/EHC for the simulated proportion of heat in the atmosphere, and CHC/N and CHC/EHC for the simulated proportion of heat absorbed by the cryosphere.

Earth Heat Inventory
The CMIP5 ensemble mean overestimates the observed ocean heat content for the period 1972-2005 CE and underestimates 180 the observations for the continental subsurface and the cryosphere ( Figure 2). Additionally, the multimodel mean yields higher total heat in the climate system than observations, as expected due to the high OHC values reached by these simulations (Figure 2a). Indeed, the CMIP5 multimodel mean yields an OHC increase of 247 ± 172 ZJ (mean ± two standard deviations,    Table 2). However, the spread in the CMIP5 results is still large, with the difference between the highest and the lowest estimates of heat storage due to sea ice melting being more than double the value of the ensemble mean (5 ZJ).
Heat uptake by subsurface ice is the second contributor to the cryosphere heat content in all models after sea ice melting. between the multimodel mean and observations, the inter-model spread is large, with the difference between the maximum and minimum AHC from CMIP5 models reaching 5 ZJ, more than double the value of the observational estimate.

Heat Partition Within Climate Subsystems
The simulated heat storage within each climate subsystem has been assessed in the previous section, displaying a large intermodel spread among CMIP5 GCMs. This wide range of results hampers the assessment of the simulated Earth heat inventory, 240 particularly the evaluation of the represented ocean heat content and total heat in the climate system. Nevertheless, models may be distributing the total heat content among the four climate subsystems similarly. This section evaluates the partition of heat among climate subsystems within each CMIP5 GCM, testing whether models simulating higher values of N and EHC distribute this energy in the same proportion among climate subsystems as models simulating lower values of total heat content.
The simulated heat partitions by the thirty CMIP5 GCMs achieve a lower inter-model spread in comparison with the simu-  Table 2). Nevertheless, the ensemble mean presents a partition of heat in each climate subsystem similar to the results for the EHI. That is, the simulated proportion of energy in the ocean is larger than observations, the proportion of heat in the continental subsurface and in the cryosphere is lower than observations, and the proportion of heat in the atmosphere is in agreement with observations. Additionally, results vary depending on the metric used to characterize total heat content in the system, particularly for the ocean.  Table   2). The spread of OHC/EHC estimates is small, with values ranging from 91 ± 2 % (MIROC5) to 100 ± 1 % (MRI-CGCM3).
Nevertheless, the simulated proportion of heat in the ocean presents different results for some models when considering the 255 integration of the radiative imbalance at the top of the atmosphere as metric for total heat in the climate system (OHC/N, black dots in Figure 4a). The model spread is much larger for OHC/N estimates than for OHC/EHC estimates, ranging from EHC values displayed in Figure 2a. That is, some CMIP5 models yield excessively different values of N and EHC, suggesting the presence of non-conservation terms in the simulated energy budget (see Section 3.1 and Discussion section). Six models 260 obtain OHC/N estimates above 100%, which indicates that the simulated N in those models is much lower than EHC estimates (the BCC-CSM1.1-M, CANESM2, CMCC-CMS, MIROC-ESM, NOR-ESM-M and NOR-ESM-ME models in Figure 4a).
The opposite behavior occurs in other five models that simulate OHC/N values below 80% (the CESM1-CAM5, CMCC-CM, GFDL-CM3, IPSL-CM5A-LR, and IPSL-CM5B-LR models in Figure 4a), which is probably a excessively small proportion of heat stored in the ocean in comparison with observations (Hansen et al., 2011;Palmer et al., 2011;Rhein et al., 2013;Palmer 265 and McNeall, 2014;Trenberth et al., 2014;Hobbs et al., 2016;von Schuckmann et al., 2016von Schuckmann et al., , 2020. Estimates of the proportion of heat in the ground from CMIP5 GCMs show smaller differences between GHC/N and GHC/EHC than the retrieved proportion of heat in the ocean (Figure 4b). Both GHC/N and GHC/EHC estimates have a multimodel mean and 95% confidence interval of 2 ± 3 %, which is in agreement with estimates derived from Church et al.  Table 2).

Discussion
The thirty CMIP5 GCMs analyzed here simulate markedly different total heat contents within the climate system, independently of the analyzed metric (N, EHC and OHC values in Figure 2a), which may be caused by the different response from each model to the common Historical forcing. That is, different models simulate distinct responses to the common external forcing, as seen in the broad range of simulated equilibrium climate sensitivities in the literature (e.g. Knutti et al., 2017).

295
Indeed, Forster et al. (2013) assessed the response to the common forcing of a large ensemble of CMIP5 GCMs in terms of climate sensitivity, feedbacks and adjusted radiative forcing, showing that these models yielded a broad range of responses.
To test the potential relationship between total heat storage and model response, we performed a covariance analysis between some of the metrics used by Forster et al. (2013) to characterize the response of CMIP5 models and the estimated Earth heat inventory here ( Figure 6). The eighteen CMIP5 models in common with those analyzed in Forster et al. (2013) do not show 300 covariance between the heat storage within the different climate subsystems and equilibrium climate sensitivity nor with the transient climate response. However, the adjusted forcing during the last part of the Historical experiment (2001( -2005 presents significant correlation coefficients with N, EHC and OHC (red triangles in Figure 6). This is a reasonable result, as different adjusted forcings result from a spread of radiative imbalances at the top of the atmosphere and climate sensitivities, from which different N values arise -and therefore distinct heat storage within the ocean (Palmer and McNeall, 2014). The 305 relationship between adjusted forcing and heat storage, nevertheless, should be considered just as a potential line of research, since the estimates of radiative forcing from transient climate simulations depend on the method employed in the analysis (Forster et al., 2016), meaning that further work is needed to evaluate the robustness of this relationship.
The simulated proportion of heat in the ocean for some models shows markedly different results depending on the used metric for total heat content in the climate system ( Figure 4a). The different heat partition is caused by the discrepancies 310 between estimates of N and EHC within each GCM simulation (Figure 2a), which are probably related to non-conservation terms in the simulated energy budget by each GCM as discussed in Hobbs et al. (2016). That is, small numerical inconsistencies, insufficient spin up time, or the amount of water leaving the LSM component at the bottom of the soil column, among others, may prevent the conservation of energy in GCM simulations (Sen Gupta et al., 2013;Hobbs et al., 2016;Séférian et al., 2016;Trenberth et al., 2016). We applied a simple drift-removal technique to each variable considered in this study in order to 315 minimize the effect of possible non-conservation terms in our results (see Section 2). This method has shown good results in previous analyses including several CMIP5 experiments, although no perfect solution is available yet (Hobbs et al., 2016).
The low ground heat content achieved by the shallow LSM components (Figure 2b) alters the distribution of heat within models, mainly causing a higher proportion of heat stored in the ocean if considering EHC as metric for total heat content.
This can be seen in a covariance analysis between OHC/EHC estimates and the depth of the LSM components in the CMIP5  . Surprisingly, the simulated CHC indicates significant covariance with the depth of the employed LSM component (red diamond in Figure 6), although this should be the result of the different subsurface volume within CMIP5 models. That is, deeper models tend to simulate more subsurface ice and GHC than shallower models, and therefore more heat can be used to thaw the larger mass of subsurface ice. This result suggests another limit to the representation of the EHI within GCM simulations, as the lack of a sufficient continental subsurface volume alters the simulated heat uptake by the subsurface ice 330 masses. Nevertheless, further work is required to clarify this point.
The simulated cryosphere heat content and heat proportion are in better agreement with observations when ignoring the heat absorbed by terrestrial ice masses from the assessment, that is, considering only sea ice as cryosphere component (see results labelled as "only sea ice" in Table 2). The same can be said about the simulated proportion of heat in the ocean, which shows a reduction of 2% in the difference with observations if considering only sea ice as cryosphere (Table 2). This is caused by the 335 lack of a representation of land ice in CMIP5 simulations, as only the simulated heat uptake by sea ice can be directly compared with observations, and the method used here to approximate the melting of land ice in the models is not accurate enough. Our approach considers snow changes in grid cells indicated as land ice by the models, but results show that this method markedly underestimates heat uptake in comparison with observations ( Figure 3). Furthermore, the observed proportion of heat in the ocean yields different results if considering the whole cryosphere for estimating EHC or if considering only the change in 340 sea ice volume (Table 2). Therefore, heat uptake by terrestrial ice sheets and glaciers is important to improve the simulated EHI and the partition of heat within the four climate subsystems. CMIP5 GCMs currently include modules representing ice sheets, but such model components were not activated for generating the CMIP5 simulations analyzed here, probably due to atmosphere-ocean-ice-sheet simulations (Nowicki et al., 2016). Although these experiments are focused on understanding the contribution of ice sheets to sea-level rise, these simulations could be also useful to test if including land ice masses enhances the representation of the Earth heat inventory within GCMs, particularly the coupled experiments.

350
The ensemble of CMIP5 GCMs analyzed here overestimates the amount of heat stored in the ocean and underestimates the heat uptake by the cryosphere and the continental subsurface, while representing changes in atmosphere heat storage similar to observations. Models present a large inter-model spread of ocean heat content and total heat content in the system, probably related to the wide range of simulated responses to external forcing in these GCMs. The lack of an adequate representation of terrestrial ice masses and continental subsurface volume within CMIP5 models limits the amount of heat allocated within the 355 cryosphere and the continental subsurface. The issue of heat conservation within complex numerical simulations also affects the represented Earth heat inventory in the CMIP5 ensemble. Nevertheless, there is good agreement between simulated and observed atmosphere heat storage and heat uptake by changes in sea ice volume.
There are two main issues hindering the assessment of the EHI in CMIP5 models in comparison with observations, the non-conservation of energy in models and the markedly different amounts of simulated total heat content in the Earth system. Radiative Forcing (ERF) is similar in CMIP5 and CMIP6 (Smith et al., 2020), and the simulated radiative imbalance seems to present smaller inter-model variability in the CMIP6 ensemble than in the CMIP5 ensemble (Wild, 2020). Therefore, a future assessment of the simulated EHI within the CMIP6 ensemble is required to determine the performance of the new generation of models. Regarding the non-conservation of energy within the models, Irving et al. (2020)  heat inventory, as well as the associated phenomena relevant to society such as sea level rise or permafrost evolution. These issues will probably be present within the CMIP6 simulations, together with non-conservation of energy and drifts, but the comparison with observational references may help to mitigate this limitations in future generations of GCMs. For example, an extended sampling of the deepest part of the ocean will improve the observational estimate of OHC, and will provide a reference to evaluate deep heat uptake in GCMs, probably reducing the drift in these models (Irving et al., 2020;von Schuckmann 380 et al., 2020). Local and regional measurements of the state of glaciers and ice sheets may help to parameterize the evolution of ice masses in individual grid cells. That is, a simplified parameterization of land ice masses based on tiling (Essery et al., 2003;Best et al., 2004) could be implemented. This strategy has been successful in representing vegetation functional types at sub-grid scales (Melton and Arora, 2014), and it has been proposed to improve the representation of permafrost in land surface model components (Beer, 2016). Additionally, expanding the global network of subsurface temperature profiles will Author contributions. FJCV designed the study, analyzed the CMIP5 simulations, and produced all results and figures. All authors contributed to the interpretation and discussion of results. FJCV wrote the manuscript with continuous feedback from all authors.
Computational Excellence Network (ACENET-Compute Canada). We thank two anonymous reviewers for their constructive feedback. We are grateful to Fiammetta Straneo and Susheel Adusumilli for their help to obtain and interpret the sea ice volume data. This work was