First Quantification of the Permafrost Heat Sink in the Earth's Climate System

Due to an imbalance between incoming and outgoing radiation at the top of the atmosphere, excess heat has accumulated in Earth's climate system in recent decades, driving global warming and climatic changes. To date, it has not been quantified how much of this excess heat is used to melt ground ice in permafrost. Here, we diagnose changes in sensible and latent ground heat contents in the northern terrestrial permafrost region from ensemble‐simulations of a tailored land surface model. We find that between 1980 and 2018, about 3.9+1.4−1.6 $3.9\genfrac{}{}{0pt}{}{+1.4}{-1.6}$  ZJ of heat, of which 1.7+1.3−1.4 $1.7\genfrac{}{}{0pt}{}{+1.3}{-1.4}$  ZJ (44%) were used to melt ground ice, were absorbed by permafrost. Our estimate, which does not yet account for the potentially increased heat uptake due to thermokarst processes in ice‐rich terrain, suggests that permafrost is a persistent heat sink comparable in magnitude to other components of the cryosphere and must be explicitly considered when assessing Earth's energy imbalance.

we are not aware of comparable long-term records of changes in ground water/ice contents that would allow to estimate the overall latent heat uptake in a similar manner. Permafrost thaw and ground ice melt are being reported to occur across the Arctic (Farquharson et al., 2019;Jorgenson et al., 2006;Liljedahl et al., 2016), and projected to continue in the coming decades and centuries (Burke et al., 2020;Lawrence et al., 2012). Hence, a potentially substantial heat sink, namely the latent heat uptake through melting of ground ice in the permafrost region, has to date not been quantified and is therefore missing in global assessments of Earth's energy budget (Forster et al., 2021;von Schuckmann et al., 2020).
In absence of suitable long-term observations, numerical models provide an alternative possibility to constrain plausible ranges for heat contents in and fluxes between different compartments of the climate system Cuesta-Valero, García-García, Beltrami, and Finnis (2021). A land surface model (LSM) tailored to realistically assess permafrost heat uptake requires (a) representation of ground freezing, (b) a sufficiently deep subsurface domain (Hermoso de Mendoza et al., 2020;, and (c) a reasonable prescription of ground ice contents and excess ice.
Here, we use the permafrost model CryoGridLite (Langer, Nitzbon, Groenke, et al., 2022) to simulate changes in sensible and latent heat content in the northern terrestrial permafrost region, and provide a first estimate of the heat uptake through both warming and thawing of permafrost during the past almost four decades . We further investigate spatial patterns of heat uptake in the study region, and put our findings into context with previous assessments of Earth's ice and energy imbalances.

Sensible and Latent Ground Heat Contents
For a given location (x, y) and year (t), we consider the mean annual ground heat content h = h S + h L (J m −2 ) which is composed of a sensible part h S (x, y, t) linked to the subsurface temperature T(z) (K), and a latent part h L (x, y, t) linked to the subsurface volumetric liquid water content θ w (z) (m 3 m −3 ): where C(z) (J m −3 K −1 ) is the volumetric heat capacity, and L sl = 334 ⋅ 10 6 J m −3 the volumetric latent heat of fusion of water. In Equations 1 and 2 the integration is done over depth z down to a fix lower boundary z lb , and <.> t denotes the annual average of the integral. Warming of ground causes an increase in sensible heat content, while thawing of (perennially) frozen ground causes an increase in latent heat content.
Here, we used a numerical model which is run for the entire northern terrestrial permafrost region to simulate subsurface temperatures and liquid water contents. Similar to Cuesta-Valero et al. (2016) and Cuesta-Valero, García-García, Beltrami, and Finnis (2021) we applied Equations 1 and 2 to diagnose mean annual ground heat contents for each year (t) and grid cell (x, y) of the model domain (see Text S5 in Supporting Information S1). Subsequently, we consider temporal changes in the mean annual sensible and latent heat contents both at the level of the grid cells (h S/L (x, y)), as well as aggregated over the entire model domain (H S/L (t) = ∑ (x,y) h S/L (x, y, t) Ω(x, y) (J), with Ω (m 2 ) being the land area per grid cell): with t 1 − t 0 > 0 being the considered time period. Positive (negative) values of ΔH S/L indicate that the ground acts as a heat sink (source) to the embedding climate system.

Model Simulations
We conducted parameter-ensemble simulations with the transient one-dimensional permafrost model CryoGri-dLite (Langer, Nitzbon, Groenke, et al., 2022) which is a simplified fast version of CryoGrid3 (Westermann et al., 2016). CryoGridLite solves the one-dimensional heat equation under consideration of phase change of soil water. Unlike CryoGrid3, it does not solve the full surface energy balance but instead uses a temperature as the upper boundary condition similar to Jafarov et al. (2012) and Westermann et al. (2013). This and the use of an implicit backward integration scheme make CryoGridLite more than two orders of magnitude faster than CryoGrid3.
We closely followed Langer, Nitzbon, Groenke, et al. (2022) to conduct distributed model simulations for the northern terrestrial permafrost region on a 1°-latitude by 1°-longitude grid. For our model setup, we used a subsurface domain down to a depth of z lb = 550 m, allowing for a much more realistic representation of the longterm transient thermal regime and heat reservoirs than the shallow subsurface representation typical for LSMs (Hermoso de Mendoza et al., 2020;. At the bottom of the model domain, a geothermal heat flux (Q geo [W m −2 ]) according to Davies (2013) was prescribed as the lower boundary condition. At the surface, the model was driven by timeseries of daily mean near-surface air temperatures (T air ) and snowfall rates. To account for the modification of ground surface temperatures through snow, CryoGridLite comprises a scheme which simulates the dynamic build-up and melting of a seasonal snowpack. We used an integration timestep of 1 day for solving the heat equation, which was found to be sufficient for numerical stability and accuracy. Mathematical details are provided in Supporting Information S1 (Text S1); details on the physical processes and the numerical solution scheme in Langer, Nitzbon, Groenke, et al. (2022).

Ground Stratigraphies
The ground stratigraphies which determine the thermal properties of the subsurface were based on the Open-ECOCLIMAP global database of land surface parameters Masson et al. (2003) and Faroux et al. (2013), the soil organic carbon map by Hugelius et al. (2014), and the global map of soil thicknesses by Pelletier et al. (2016). CryoGridLite does not comprise a dynamic hydrology scheme such that the volumetric composition in terms of mineral, organic, and water/ice contents remains constant throughout the simulation. However, the uncertainty and variability with respect to the subsurface hydrological conditions and their effect on the heat transfer have been addressed by performing ensemble simulations with different water/ice contents and distributions (see Section 2.2.3).

Forcing Data and Model Spin-Up
The meteorological forcing data were based on the ERA-Interim reanalysis (Dee et al., 2011;Simmons & Poli, 2015) spanning the period from 1980 to 2018 which is also the main analysis period of our study. We performed a model spin-up to establish a transient equilibrium state of the ground thermal regime at the beginning of the analysis period. For this, we applied anomalies derived from paleo-climate simulations by the CSIRO Mk3L climate model (Bothe et al., 2013;PAGES 2k-PMIP3 group, 2015;Phipps et al., 2013) to the ERA-Interim data (see Text S2 in Supporting Information S1 for details). The spin-up consisted of two phases: First, a single simulation with default parameter values was run from year 500-1600 C.E. for each grid cell of the model domain. Second, an ensemble of 50 simulations with variations in the model parameters was run from 1,600 to 1,980 (see Text S3 in Supporting Information S1 for details). Our results are based on subsequent period from 1980 until 2018 for which the parameter-ensemble was continued using the original ERA-Interim forcing data.

Parameter Ensemble
The heat uptake of the ground and its partitioning in the permafrost region are affected by various environmental factors. The ground thermal regime and, thus, the sensible heat content are crucially affected by the seasonal snowpack (e.g., Jan & Painter, 2020; A. G. Slater et al., 2017;Zweigel et al., 2021). Potential changes in the latent heat content are primarily affected by the amount and the distribution of ground ice. To account for both the uncertainty and the spatial variability associated with these environmental factors, we conducted an ensemble of 50 independent simulations with a random variation of several parameters for each grid cell of the model domain. Following Langer, Nitzbon, Groenke, et al. (2022), we varied the model parameters prescribing the depths and pore water/ice contents of several subsurface layers, as well as the threshold for the maximum height of the snowpack.
In addition to the variation of pore ice/water, we also considered the presence of excess ice, based on the map by Brown et al. (1997). Specifically, we increased the water/ice content by a volumetric fraction (χ) within a domain from the ground surface (z = 0) down to a depth z χ and reduced the sediment fraction in this part of the stratigraphy by a factor 1 − χ. The depth (z χ ) and the amount (χ) of excess ice have been varied in the parameter ensemble using uniform prior probabilities based on the ground ice map (z χ ∈ [10, 20 m]; .3] for low/medium/high ground ice content). Note that we treated excess ice only statically such that it modifies the total ground ice content and the ground thermal properties, while dynamic stratigraphic changes due to excess ice formation (e.g., through ice segregation) or excess ice melt (e.g., through ground subsidence and surface water formation) were not considered.
With the ensemble simulations we captured a wide range of environmental conditions which determine ground heat storage and heat exchange with the atmosphere on a local scale. In addition, we used the ensemble simulations for model evaluation, and to estimate confidence intervals for our results.

Model Evaluation
We evaluated our simulation results against observations of ground temperatures and thaw depths (see Figures S1 and S2, and Text S4 in Supporting Information S1). Based on this evaluation, we consider our model setup as suited for providing a first-order estimate of permafrost heat uptake. For a more detailed discussion of the model evaluation we refer to Langer, Nitzbon, Groenke, et al. (2022).
We diagnosed that 2.2 +0.6 −0.6 ZJ (1980-2018; 56% of the total heat uptake) of sensible heat have accumulated in the permafrost region's landmass which is equivalent to an average sensible heat uptake rate of 0.10 +0.03 −0.03 W m 2 , and   Slater et al. (2021) (Sl21) amounts to about 11.3% of the global land heat uptake in that period estimated by Cuesta-Valero, García-García, Beltrami, González-Rouco, and García-Bustamante (2021)  A similar amount of heat taken up in the permafrost region is used for the melting of ground ice. We found that 1.7 +1.3 −1.4 ZJ (1980-2018; 44% of the total heat uptake) of latent heat were taken up by the subsurface of the model domain through melting of ground ice, equivalent to an average latent heat uptake rate of 0.08 +0.06 −0.07 W m 2 . This would correspond to about 10.9% of the heat accumulation in Earth's cryosphere, if the contribution of permafrost was added to the estimate by von Schuckmann et al. (2020) for the cryosphere without permafrost (Figure 1b). Converting the permafrost latent heat change into an average annual ground ice melt rate (using a specific latent heat of fusion of water of L sl = 330, 000 J kg −1 ) results in a value of m pi = 140.6 Gt yr −1 (1994-2017) which is about 9.0% of the combined melt rates of grounded ice (543 Gt yr −1 ) and floating ice (653 Gt yr −1 ) estimated for the same period (T. Slater et al., 2021) (Table 1, Sl21).
Our findings reveal a duality of the permafrost region: on the one hand, it is a component of Earth's continental landmass, which is warming at a similar rate as the global continental landmasses; on the other hand, it is a component of Earth's cryosphere, which constitutes a substantial heat sink when compared to different components of the remaining cryosphere.

Comparison to Other Components of the Cryosphere
While the latent heat uptake by permafrost has been steadily increasing from about 0.03 ZJ yr −1 in the 1980s to about 0.06 ZJ yr −1 in the 2010s (Figure 1c), the other cyrosphere components show an even more pronounced increase during this period, which is particularly strong for the ice sheets of Greenland and Antarctica from the 1990s to the 2000s. Consequently, the relative share of permafrost in the heat uptake by the entire cryosphere declined from about 16% during the 1980s (1980)(1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989) to about 11% during the 2010s (2010-2017), suggesting a decreasing importance of the permafrost heat sink relative to the other cryospheric heat sinks. We hypothesize that this is due to a slower response of ground ice (permafrost) to recent warming compared to that of floating ice (ice shelves, sea ice) and grounded ice (ice sheets, glaciers) (Figure 1d), which are directly exposed to the atmosphere. Permafrost in turn is primarily a subsurface phenomenon protected from the atmosphere through the active layer as well as isolating buffer layers such as vegetation (Shur & Jorgenson, 2007;Stuenzi et al., 2021), which can delay and dampen the response to warming at the surface. The considered period of about four decades is, however, too short to draw firm conclusions regarding the longer-term strength of the different cryospheric heat sinks.
All cryospheric components including permafrost ground ice have in common that their capacity as a heat sink is limited by the amount of ice stored in them. These sinks for latent heat are thus qualitatively different from the sinks for sensible heat provided by the ocean, the land, and the atmosphere. With the projected shrinking of Earth's grounded, floating and ground ice masses (Fox-Kemper et al., 2021), the cryospheric heat sinks will also eventually shrink, implying that a larger fraction of the EEI has to be stored in the oceans, land, and the atmosphere. Indeed, a decrease of the cryospheric heat storage and a concurrent increase in the shares of the ocean and atmosphere during the past decade has been reported (von Schuckmann et al., 2020(von Schuckmann et al., , 2023. This shift in the partitioning of the EEI happens along with a current trend toward larger absolute values of the EEI (Loeb et al., 2021).

Contrasting Patterns of Sensible and Latent Ground Heat Changes
Except for very few locations, the entire northern permafrost region has been simulated to constitute a net heat sink during the analysis period (1980-2018; Figure 2f), indicating that permafrost ground has been warming or thawing (or both) across the Arctic during the last decades. Our simulations further suggest contrasting spatial patterns of sensible versus latent heat uptake. The sensible heat change is most pronounced in areas featuring cold permafrost, that is, the northernmost latitudes as well as a vast portion of central eastern Siberia (Figure 2d), roughly coinciding with the zone of continuous permafrost. The latent heat change in turn is dominating in the southern part of the model domain, in a band from central and southern Alaska through Canada and the southern coast of the Hudson Bay to northern Quebec, as well as from northern Finnoscandia through western and southern Siberia, terminating at the Sea of Okhotsk, coinciding with the transitional region between the discontinuous and the continuous permafrost zones (Figures 2a and 2e). The overall tendency of higher sensible heat changes in more northerly regions is in line with ground warming trends derived from borehole measurements (Biskaborn et al., 2019;S. L. Smith et al., 2022).
To better understand the co-variation of changes in sensible and latent heat contents, Figure 3 shows the respective linear trends for all grid cells of the model domain. For the great majority (85.0%) of locations, both latent and sensible heat contents are showing positive trends during the analysis period. Sensible heat trends show a bimodal distribution with one mode being slightly positive, and the second mode being substantially positive between 0.15 and 0.20 W m −2 . Latent heat trends in turn follow a unimodal distribution which has its mode at slightly positive values, but with a long tail toward more positive values of up to about 0.6 W m −2 . In other words, while for the majority of locations the ground is warming but not thawing, there is also a substantial amount of locations, where permafrost is thawing without significant warming trends. Figures 3a and 3b further reveal that the locations with substantial warming are characterized by cold mean annual ground temperatures (MAGT < −5°C) and shallow maximum thaw depths (TD < 2 m). At the same time, those locations which show substantial thawing are characterized by MAGT close to 0°C, and a wide range of thaw depths mostly exceeding 2 m. Overall, considering changes and trends in sensible and latent ground heat contents complements the traditional approach of assessing the thermal state of permafrost in terms of ground temperatures and thaw depths (Biskaborn et al., 2019;S. L. Smith et al., 2022).

Limitations and Implications
An observation-based assessment of the heat uptake through terrestrial permafrost is not only limited by the sparsity of existing long-term monitoring sites across the permafrost region, but necessarily incomplete since changes in latent heat contents are not being monitored systematically, with few exceptions (e.g., Boike et al., 2018Boike et al., , 2019. This concerns (a) the melting of pore ice, which could be measured as changes in the liquid water content in the active and permafrost layers (Nicolsky & Romanovsky, 2018), and (b) the melting of excess ice, which is manifesting in ground subsidence, but technically challenging to be monitored at large scales (e.g., Liu & Larson, 2018).
Our model-based assessment compensates this lack of observations by providing consistent estimates of both sensible and latent heat changes across the northern terrestrial permafrost region. Nevertheless, it should be considered a first-order approximation, as it is subject to several sources of uncertainty. Besides the uncertainty of the climate forcing, particularly before the reanalysis period, there is only limited knowledge on the pan-Arctic distribution of ground ice at present, which constitutes the primary source of uncertainty for our estimate of latent heat changes. We addressed this factor of uncertainty by randomly varying the ground ice contents and distribution in our ensemble simulations. Consequently, the range between the 5th and 95th percentiles of estimated pan-Arctic latent heat uptake (2.7ZJ) is more than two-fold of the corresponding range for the sensible heat uptake (1.3ZJ) (Table 1, ERA-Interim; Figures 2g and 2h). Further shortcomings of our model are the lack of a dynamic hydrology scheme, no detailed treatment of vegetation effects, a static excess ice representation, and a simplified freezing characteristic. Similarly, detailed snow processes such as rain-on-snow events (Westermann Figure 2. (a-c) Permafrost zones (a; from dark to light brown: continuous, discontinuous, sporadic, isolated) and excess ground ice contents (b) according to Brown et al. (1997). Excess ice contents (χ, without massive ground ice) have been considered in the simulations, but total ice contents are usually dominated by pore ice in the soil's pore space (c). (d-i) For each 1°-by-1° grid cell of the model domain, the maps are showing the ensemble mean (d-f) and the difference between the 95th and 5th percentiles (g-i) of the simulated mean change rate in sensible (d, g), latent (e, h), and total (f, i) ground heat contents between 1980 and 2018. The heat changes during this period, integrated over the entire model domain are indicated in the lower left corner of panels (d-f), and the spread between the 5th and 95th percentile in panels (g-i). et al., 2011), depth hoar formation (Gouttevin et al., 2018), or wind compaction (Zweigel et al., 2021) were not considered explicitly, and assumed to be of secondary importance. Despite these limitations, our model captures the essential physical mechanisms to assess the exchange of heat between the atmosphere and the land surface, its transport and storage within the subsurface, and its partitioning during the phase change of ground water. Indeed, our model is capable of reproducing recent borehole temperatures well (see Text S4 in Supporting Information S1 and Langer, Nitzbon, Groenke, et al. (2022)).
Overall, we consider our model-based estimate of the permafrost heat sink to be valid but probably conservative for the following reasons: First, we did not take into account thaw processes acting in ice-rich permafrost terrain, so-called thermokarst. These processes are increasingly abundant in a warming Arctic and involve the thawing of permafrost containing excess ice as well as the melting of massive ground ice (Turetsky et al., 2019(Turetsky et al., , 2020. Thermokarst activity involves a considerable uptake of latent heat, and has been observed and projected to occur also in areas featuring cold permafrost (Farquharson et al., 2019;Nitzbon et al., 2020) which otherwise are primarily a sink for sensible heat. While there are recent advances in modeling thermokarst Lee et al., 2014;Westermann et al., 2016), global assessments are still limited by the lack of information on the distribution of ground ice. Second, our model domain only comprises the northern terrestrial permafrost region. We did not consider changes in latent heat in other areas featuring seasonally or perennially frozen ground such as the Tibetan plateau, Alpine regions, and Antarctica. The heat uptake through thawing of submarine permafrost (Overduin et al., 2019) was also not assessed in our analysis. While the net heat change in these areas is very likely outweighed by the northern terrestrial permafrost region, they probably also constitute net heat sinks in recent decades (S. L. Wilkenskjeld et al., 2022).
Our findings underline the importance of the northern permafrost region as a component of the terrestrial cryosphere that is relevant to the global climate beyond its dormant carbon pool. Key steps to an improved representation of permafrost biogeophyiscs in LSMs include, (a) increasing the depth of the subsurface , (b) representing excess ice melt (Lee et al., 2014) and formation (Aga et al., 2023), and (c) parameterizing subgrid-scale heterogeneity and lateral fluxes Nitzbon et al., 2019Nitzbon et al., , 2021N. D. Smith et al., 2022).

Conclusions
We have provided a first estimate of heat uptake by Arctic terrestrial permafrost in recent decades. We conclude that northern permafrost constitutes a substantial heat sink in the Earth's climate system and is special because 9 of 11 of its dual function as a part of the continental landmass and as a component of the cryosphere. Over the past four decades, the northern permafrost region has absorbed roughly equal amounts of sensible and latent heat, resulting in both warming and thawing of permafrost. While previous assessments of Earth's energy imbalance have considered only the warming of the continental landmass, our results suggests that the melt rate of ground ice is comparable in magnitude to that in other major components of the cryosphere, and thus must be explicitly considered in assessments of Earth's energy and ice imbalances.