Projecting circum-Arctic excess-ground-ice melt with a sub-grid representation in the Community Land Model

To address the long-standing underrepresentation of the influences of highly variable ground ice content on the trajectory of permafrost conditions simulated in Earth system models under a warming climate, we implement a sub-grid representation of excess ground ice within permafrost soils using the latest version of the Community Land Model (CLM5). Based on the original CLM5 tiling hierarchy, we duplicate the natural vegetated land unit by building extra tiles for up to three cryostratigraphies with different amounts of excess ice for each grid cell. For the same total amount of excess ice, introducing sub-grid variability in excess-ice contents leads to different excess-ice melting rates at the grid level. In addition, there are impacts on permafrost thermal properties and local hydrology with sub-grid representation. We evaluate this new development with single-point simulations at the Lena River delta, Siberia, where three sub-regions with distinctively different excess-ice conditions are observed. A triple-land-unit case accounting for this spatial variability conforms well to previous model studies for the Lena River delta and displays markedly different dynamics of future excess-ice thaw compared to a single-land-unit case initialized with average excess-ice contents. For global simulations, we prescribed a tiling scheme combined with our sub-grid representation to the global permafrost region using presently available circum-Arctic ground ice data. The sub-grid-scale excess ice produces significant melting of excess ice under a warming climate and enhances the representation of sub-grid variability of surface subsidence on a global scale. Our model development makes it possible to portray more details on the permafrost degradation trajectory depending on the sub-grid soil thermal regime and excess-ice melting, which also shows a strong indication that accounting for excess ice is a prerequisite of a reasonable projection of permafrost thaw. The modeled permafrost degradation with sub-grid excess ice follows the pathway that continuous permafrost transforms into discontinuous permafrost before it disappears, including surface subsidence and talik formation, which are highly permafrost-relevant landscape changes excluded from most land models. Our development of sub-grid representation of excess ice demonstrates a way forward to improve the realism of excess-ice melt in global land models, but further developments require substantially improved global observational datasets on both the horizontal and vertical distributions of excess ground ice.

1 Introduction

Permafrost soils are often characterized by different types of ground ice that can exceed the pore space (Brown et al., 1997; Zhang et al., 1999). The presence of such “excess” ground ice can alter the permafrost thermal regime and landscape structure. Widespread thawing of permafrost is expected in a warmer future climate and modeling studies suggest large-scale degradation of near-surface permafrost at the end of the 21st century (Lawrence et al., 2008, 2011). Melting of ground ice due to active layer thickening releases water in the form of surface runoff, subsurface flow, or both, causing surface subsidence and modifying the local hydrological cycle (West and Plug, 2008; Grosse et al., 2011; Kokelj et al., 2013; Westermann et al., 2016). In addition to containing ground ice, some permafrost soils store massive amounts of carbon, which could be released to the atmosphere in the form of greenhouse gases upon thawing (Walter et al., 2006; Zimov et al., 2006; Schuur et al., 2008), possibly making a positive feedback to amplify future climate change (Koven et al., 2011; Schaefer et al., 2014; Burke et al., 2013). The existence of excess ice and its distribution in permafrost can significantly affect the rate of permafrost thawing (Westermann et al., 2016; Nitzbon et al., 2020) and in turn the rate of soil carbon release (Hugelius et al., 2014; Schuur et al., 2015; Turetsky et al., 2019). Therefore, better projections of excess-ice melt are critical to improve our understanding of the impacts of permafrost thaw on corresponding climatic impacts.

Previous studies address excess-ice modeling on the local or regional scale, in which the small study area makes it possible for detailed configurations of the cryostratigraphy of permafrost and excess ice based on observations. Simulations for the Lena River delta have retrieved the permafrost thermal dynamics fairly close to the observations with excess ice incorporated in the modeling (Westermann et al., 2016). A two-tile approach allowing lateral heat exchange between two land elements demonstrated that maintaining thermokarst ponds requires the heat loss from water to the surrounding land (Langer et al., 2016). A similar tiling approach has been applied to project the landscape changes due to permafrost thaw for ice-wedge polygons and peat plateaus with different features of ice melting and surface subsidence (Aas et al., 2019; Nitzbon et al., 2019).

On the global scale, the land components of Earth system models (ESMs) have significant capabilities of representing key permafrost physics. In the Community Land Model (CLM), for example, the representation of permafrost-associated processes has been continuously improved. By including key thermal and hydrological processes of permafrost, the CLM version 4 (CLM4) has reasonably reproduced the global distribution of permafrost (Lawrence et al., 2008, 2012; Slater and Lawrence, 2013). Projections based on the CLM4 under its highest warming scenario (RCP8.5) have shown over 50 % degradation of near-surface permafrost by 2100 (Lawrence et al., 2012). Moreover, the recently released CLM5 has more advanced representations of many biogeophysical and biogeochemical processes (Lawrence et al., 2019). A refined soil profile and upgraded snow accumulation and densification scheme in the CLM5 could contribute to simulating more realistic permafrost thermal regimes, whereas upgrades on biogeochemistry improve simulations of soil carbon release in response to permafrost thaw. In addition, an excess-ice physics scheme has been implemented in CLM4.5 (CLM4.5_EXICE) by Lee et al. (2014), which allowed for the first-order simulation of surface subsidence globally by modeling excess-ice melt under a warming climate.

The homogeneous distribution of excess ice throughout the grid cell in CLM4.5_EXICE (Lee et al., 2014) could cause biases in thaw trajectories in the warming climate. In nature, excess ice forms in a highly localized manner due to a variety of accumulation processes. For instance, segregated ice formed during frost heave differs substantially in excess-ice morphology from ice wedges that are formed from repeated frost cracking and freezing of penetrating water. Field measurements illustrate that the depth distribution of ground ice can vary substantially on the order of 10–50 m horizontally and 0–10 m vertically (De Pascale et al., 2008; Fritz et al., 2011). The horizontal grid spacing of ESMs, on the other hand, usually ranges from 1 to 2^∘ (∼100–200 km horizontal scale), which makes it impossible to represent localized excess ice. The mismatch in spatial scale between the model and the real world raises concerns for the reliability of excess-ice modeling in ESMs. Aside from the homogenously initialized excess ice in the grid cell, CLM4.5_EXICE initializes excess ice in the same soil depths globally (below 1 m), regardless of the varying active layer thickness in circum-Arctic permafrost areas (Lee et al., 2014). Such deficiencies in excess-ice parameterization hamper global projections of permafrost thaw including excess ice with ESMs.

To narrow the gap between the high spatial variability of excess ice and the coarse grid spacing in the ESMs, we applied a sub-grid approach in representing excess ice in permafrost soils within the CLM5 to investigate how presence and melting of excess ice affect land surface physics under a warming climate. We conducted idealized single-point simulations to examine the robustness of model development, and we furthermore conducted global simulations using a first-order estimate for the spatial distribution of excess ice and associated cryostratigraphies. Due to the lack of information in global excess-ice conditions, it is not the aim of this study to accurately project excess-ice melt and surface subsidence in the 21st century, but rather to develop and present a functioning process within a land surface model that can eventually bring permafrost thaw modeling towards a higher degree of accuracy on a global scale. The CLM5 with sub-grid excess-ice representation developed through this study would be ready to serve as a proper simulation tool on further advancing global excess-ice modeling once new datasets become available.

2 Methodology

2.1 Sub-grid representation of excess ice in CLM5

The CLM5 model utilizes a three-level tiling hierarchy to represent sub-grid heterogeneity of landscapes, which are (from top to bottom) land units, columns, and patches (Lawrence et al., 2019). There is only one column (the natural soil column) that is under the natural vegetated land unit, which represents soil including permafrost. In this study, we modify the CLM5 tiling hierarchy by duplicating the natural vegetated land unit, making extra land units for prescribing up to three different excess-ice conditions in permafrost (Fig. 1). The original natural vegetated land unit is considered “natural vegetated with no excess ice” (hereafter no-ice land unit), while we denote the additional land units as “natural vegetated with low content of excess ice” (hereafter the low-ice land unit), “natural vegetated with medium content of excess ice” (hereafter the mid-ice land unit), and “natural vegetated with high content of excess ice” (hereafter the high-ice land unit). The sub-grid initial conditions of excess ice are imported as part of the surface data, which include the variables of volumetric excess-ice contents, depths of the top and bottom soil layers of added excess ice, and the area weights of the four land units.

Figure 1 Modification of the CLM5 tiling hierarchy on the land unit level containing four natural vegetated land units for different excess-ice conditions.

We adopted the excess-ice physics from CLM4.5_EXICE (Lee et al., 2014), including thermodynamic and hydrological processes. The added excess ice is evenly distributed within each soil layer. Whereas the original CLM5 model already represents the dynamics of pore ice, our representation of excess-ice physics only addresses the ground ice bodies that exceed soil pore space. The volumetric excess-ice content in this study is defined as the ratio of the volume of excess ice in a soil layer to the volume of the whole soil layer. For example, a 50 % volumetric content of excess ice means the excess-ice body occupies 50 % volume of a soil layer, while the rest of the soil (and pore ice) occupies the other 50 % volume of the soil layer. If not otherwise notified, the parameter of volumetric ice content in this paper refers only to that of excess-ice bodies. After adding excess ice, the soil layer thickness increases accordingly. Because ice density is considered constant, the increase in soil layer thickness is linearly proportional to the volumetric content of excess ice. For example, adding an excess-ice body with a 50 % volumetric excess-ice content doubles the soil layer thickness of the corresponding soil layer. The revised algorithm for thermal conductivity and heat capacity of soil involves the effects of added excess ice, while the revised phase change energy equation allows excess ice to melt. The meltwater adds to soil liquid water in the same soil layer, and it can move to the above layer if the original layer is saturated. Such numerical implementation replicates how the melted excess ice eventually converts to runoff and discharges from the soil in the case of well-drained conditions. As excess ice melts, soil layer thickness decreases, which corresponds to surface subsidence due to excess-ice melt. In our model parameterization, excess ice only melts and does not re-form since the applied excess-ice physics does not account for the different ice formation processes.

Aside from sub-grid tiles for excess ice, we acknowledge that the version upgrade from CLM4.5 to CLM5 as the base model modifies the results of excess-ice melt compared to the results from Lee et al. (2014). By default, CLM5 represents soil with a 25-layer profile, for which the top 20 hydrologically active layers cover 8.5 m of soil. There are an additional 10 soil layers and it is 4.7 m deeper compared to the default hydrologically active soil layer profile in CLM4.5, not to mention the substantially more complex biogeophysical processes (Lawrence et al., 2019). Therefore, we developed the sub-grid representation of excess ice within the framework of the latest version of CLM. The duplicated land units prolong computation time by roughly 10 % compared to the original CLM5. We are, therefore, confident that our model development is highly efficient in addressing the sub-grid excess ice and subsequent permafrost thaw.

We examined the sensitivity of sub-grid excess-ice initialization by conducting idealized experiments (see the Supplement). Overall, for the same amount of excess ice in one grid cell located in the same depth, a higher volumetric excess-ice content along with a smaller area weight result in a later start of excess-ice melt and a lower melting rate. The different melting features from different sub-grid distributions of excess ice then lead to different hydrological impacts on the permafrost soil. The results of the idealized experiments show the necessity of introducing sub-grid configuration of excess ice to the CLM that has a typical horizontal grid spacing of 1–2^∘. More details are available in the Supplement.

2.2 Single-point simulations for the Lena River delta, Siberia

We conduct single-point simulations for the Lena River delta and compare the CLM5 model results to reference simulations with the CryoGrid3 model for the same location (Westermann et al., 2016). Abundant background information is available on the soil and ground ice dynamics from both observation and modeling, making the Lena River delta a suitable location to further evaluate our model development. The Lena River delta can be broadly categorized into three different geomorphological units that have distinctively different subsurface cryostratigraphies of excess ice (Schneider et al., 2009; Ulrich et al., 2009). In the eastern and central part of the river delta, ground ice has been accumulated in the comparatively warm Holocene climate. The subsurface sediments (hereafter denoted as “Holocene ground ice terrain”) are generally supersaturated with wedge ice that can extend up to 9 m underground with the volumetric contents of total ground ice (pore ice and excess ice) ranging from 60 % to 80 % (Schwamborn et al., 2002; Langer et al., 2013). On the other hand, higher excess-ice contents are found in Pleistocene sediments in the Lena River delta (hereafter the “Yedoma ice complex”), which are characterized by Yedoma-type ground ice (Schirrmeister et al., 2013), which can reach depths of up to 20–25 m and volumetric contents of total ground ice as high as 90 % (Schwamborn et al., 2002; Schirrmeister et al., 2003, 2011). Finally, the northwestern part of the delta features sandy sediments and is characterized by low excess-ice contents (hereafter denoted the “no-excess-ice terrain”; Rachold and Grigoriev, 1999; Schwamborn et al., 2002).

We determine the area weights of excess-ice land units in one single point based on the spatial pattern of three subregions (Fedorova et al., 2015). The cryostratigraphy and the volumetric contents of excess ice strictly follow those in Westermann et al. (2016). Noting that the excess-ice initialization scenario in Westermann et al. (2016) does not necessarily represent the realistic excess-ice condition for the Lena River delta, the purpose of applying the same excess-ice cryostratigraphy as in Westermman et al. (2016) is to evaluate our model development by addressing intercomparisons between model results. Meanwhile, we did not customize soil properties for different land units as in Westermann et al. (2016), as our model development does not support varying soil properties for different sub-grid land units. We also directly apply the snow accumulation physics in the CLM rather than customizing the snow density. By default, the current model does not form thermokarst lakes as the meltwater from excess-ice melt becomes surface runoff and is removed from the grid cell. To apply the sub-grid representation, we initialize the case with three land units (the triple-land-unit case) that respectively represent the three terraces in the Lena River delta. We also initialize an “average-ice single-land-unit” case without the sub-grid representation of excess ice. The excess-ice amount for each soil layer in the average-ice single-land-unit case is initially the same as that in the triple-land-unit case. The volumetric content of excess ice is determined by spatially averaging those for three excess-ice land units in the triple-land-unit case. Detailed information on the applied excess-ice conditions for both cases is listed in Table 1.

Table 1 The excess-ice initialization scenario in each of the three terraces (land units) for the Lena River delta, as well as that for the single-land-unit excess-ice initialization case.

n/a: not applicable.

We employed the single-point forcing data from Westermann et al. (2016) for the Lena River delta from 1901 to 2100, which is based on the CRU-NCEP (http://dods.extra.cea.fr/data/p529viov/cruncep/, last access: 3 March 2020) dataset for the historical period (1901–2005) and the CCSM4 model output under the RCP4.5 scenario for the projected period (2006–2100), but downscaled with in situ observations. We run 100-year spin-up simulations in order to stabilize the permafrost thermal regime after adding excess ice. Spin-up simulations are produced by running the model with cycled 1901–1920 climatological data. The purpose of spin-up simulations is to stabilize ground temperatures and volumes of excess-ice bodies. The 100-year length for spin-up is sufficient, as the model is run in satellite phenology (SP) mode that does not involve slowly evolving biogeochemical processes such as soil carbon accumulation. Moreover, we address idealized single-point simulations for additional permafrost locations with both continental and maritime climate that showcase the difference to Lee et al. (2014), the results of which are included in the Supplement.

2.3 Global simulations of excess-ice melt

The information available for the spatial distribution of excess ice and associated cryostratigraphies on the global scale is generally not as detailed as in the Lena River delta due to the lack of observations. For our global simulations we employ the widely used “Circum-Arctic Map of Permafrost and Ground-Ice Conditions” (hereafter the CAPS data; Brown et al., 1997) as a data source, while we translate the ground ice condition in the CAPS data to different excess-ice stratigraphies as model input data. The CAPS permafrost map categorizes the global permafrost area into classes coded by three factors: (i) permafrost extent (c: continuous; d: discontinuous; s: sporadic; i: isolated), (ii) visible ground ice content (h: high; m: medium; l: low), and (iii) terrain and overburden (f: lowlands, highlands, and intra- and intermontane depressions characterized by thick overburden cover; r: mountains, highland ridges, and plateaus characterized by thin overburden cover and exposed bedrock), resulting in more than 20 different varieties in permafrost characteristics (Fig. 2). For the simulations, we only use the CAPS distinction between the three classes: high, medium, and low ice contents. We qualitatively categorize excess-ice types with typical cryostratigraphies for which observations are available, recognizing that this is a crude first guess of the global distribution of ground ice which needs to be improved in future studies.

Figure 2 Spatial distribution of excess ground ice in the Northern Hemisphere modified from Brown et al. (1997). Compared to the original data, permafrost extents and ground ice contents are converted to definite numbers (percentages) for model computation.

The high-ice CAPS classes (e.g., chf, chr, and dhf) in central and eastern Siberia, as well as in Alaska, partly coincide with Yedoma regions (Kanevskiy et al., 2011; Grosse et al., 2013). The cryostratigraphy of the high-ice land unit is therefore broadly oriented at the excess-ice contents and distribution in intact Yedoma, which is characterized by massive ice wedges leading to typical average volumetric content of total ground ice in the range from 60 % to 90 % (Schwamborn et al., 2002; Kanevskiy et al., 2011). We therefore set the volumetric content of excess ice in the high-ice land unit to 70 %, and we put excess ice in all the soil layers between 0.2 m below the active layer and the bottom of the hydrologically active soil layer (8.5 m). The onset depth of the excess ice just below the active layer is based on the assumption of active ice aggradation which occurs at or below the permafrost table, e.g., the formation of wedge or segregation ice. Initializing high excess-ice content throughout the whole soil layer imitates the cryostratigraphy of Yedoma-type ice, while roughly 65 % of the high-ice land unit is located outside of the observed Yedoma regions (Schuur et al., 2015). The effects, limitations, and potential improvements of this initialization scenario will be mentioned in the discussion section. For the low-ice land unit, we assume both a significantly lower volumetric excess-ice content and a smaller vertical extent of the excess-ice body. The volumetric excess-ice content is set to 25 %, and we add excess ice at soil layers within 0.2 to 1.2 m below the active layer, which in particular represents sediments with segregated ice (e.g., Cable et al., 2018) but also accounts for a wide range of different excess-ice conditions found throughout the permafrost domain. For the mid-ice land unit, we set the volumetric excess-ice content to 45 % and put excess ice within 0.2 to 2.2 m below the active layer, making the volumetric excess-ice content and its vertical extent in between those for the low- and high-ice land units. The cryostratigraphies determine that excess-ice melt in the low-ice land unit can result in a maximum of 0.36 m of surface subsidence, while excess-ice melt in the medium-ice land unit can result in a maximum of 1.78 m of surface subsidence. For the high-ice land unit, the surface subsidence can be more than 10 m if all excess ice melts, which is expected to vary in space because of the different active layer thickness. For all three land units, the active layer thickness is determined by the soil temperature profile by the end of the spin-up in a no-ice case, which is the simulation by the original CLM5 model without excess ice incorporated. Non-permafrost regions in the CAPS data are assigned the no-ice land unit for 100 % of their area. We emphasize that the prescribed cryostratigraphies are a first-order approximation that can by no means represent the wide variety of true ground ice conditions found in the permafrost domain. Nevertheless, this makes it possible to gauge the effect of excess-ice melt on future projections of the permafrost thermal regime, when compared to “traditional” reference simulations without excess ice.

We design a tiling scheme prescribing the assignment of land units for each CAPS class based on previous observations and empirical estimates (Table 2). All CAPS classes in this study are categorized into three levels of volumetric ice content (5 %, 15 %, and 25 %) that are converted from the ranges (<10 %, 10 %–20 %, and >20 %) in the original CAPS data. The goal of our tiling scheme is to determine a combination of area weights of three excess-ice land units for each CAPS class, making the spatially averaged volumetric content of excess ice the same as that for the CAPS class. We assume that all CAPS classes have the same area fraction (20 %) of the low-ice land unit, and the CAPS classes with a higher ice content are due to the existence of the land units with a higher content of excess ice. We make this assumption based on previous studies that the segregated ice is widely distributed in permafrost. Observational studies have found segregated ice bodies in various continuous permafrost regions across the circum-Arctic including west central Alaska (Kanevskiy et al., 2014); Nunavik, Canada (Calmels and Allard, 2008); and Svalbard (Cable et al., 2018). In discontinuous-permafrost regions, segregated ice bodies also commonly exist underneath palsas and lithasas, including Fennoscandia (Seppälä, 2011); Altai and Sayan, Russia (Iwanhana et al., 2012); the Himalayas (Wünnemann et al., 2008); and Mongolia (Sharkhuu, 1999). The volumetric content of visible segregated ice bodies mentioned above ranges widely from 10 % to 50 % (Gilbert et al., 2016).

Table 2 The tiling scheme prescribing area weights of land units for each CAPS class. The detailed CAPS classes are shown in Fig. 2.

Note: for each class, the first letter is for the permafrost extent, the second for the excess-ice content, and the third for the terrain and overburden, following Brown et al. (1997).

Given the tiling scheme prescribed above, all CAPS classes are assigned a 20 % area of low-ice land unit. Correspondingly, the CAPS classes with 15 % volumetric ice content are assigned another 14 % area weight for mid-ice land unit on top of the CAPS classes with 5 % volumetric ice content, while the CAPS classes with 25 % volumetric ice are assigned another 22 % area for high-ice land unit on top of the CAPS classes with 15 % volumetric ice content. The classes of “chf” and “chr” are the exceptions as their corresponding regions are typically with the landscape of Yedoma or ice-wedge polygonal tundra or both (Kanevskiy et al., 2011; Grosse et al., 2013). We therefore assign only the low- and high-ice land units for these two CAPS classes. Summing up the land unit fractions for all the CAPS grid cells within each CLM grid cell obtains the area weights on the grid level that is stored in the surface data file. Figure 3 shows a schematic plot for the initialization scenario and the area covered by different excess-ice land units as the result of sub-grid excess-ice initialization in the global simulation case. Note that excess ice for some regions (e.g., southern Norway and the Alps) can completely melt out during the spin-up period since the CLM initial condition prescribes an overly warm (non-permafrost) soil temperature for these regions.

Figure 3 Schematic representation of the sub-grid excess-ice initialization scenario and maps showing the area weight (%) occupied by different excess-ice land units, i.e., the initial condition of excess ice in the global simulation.

In this study, we define the grid cells or land units with permafrost as the ones having at least one hydrologically active soil layer that has been frozen in the last consecutive 24 months. In this case, we define fully degraded permafrost when all land units in one grid cell have an active layer thickness of more than 6.5 m, recognizing that in reality permafrost at many localities may continue to exist to greater depths. We also prepare a “grid-average ice case” by applying the same total amount of excess ice as in the sub-grid ice case in each soil layer, but using only one land unit instead of three that account for the sub-grid variability of excess ice. The volumetric content of excess ice in the single land unit is calculated as the spatial average of those in the three land units in the triple-land-unit case. This grid-average ice case provides a reference to evaluate the effects of the sub-grid excess-ice representation on the global scale. Finally, we simulate a reference case without excess ice, denoted the “no-ice case” in the following. Details on the three cases for the global simulations are listed in Table 3. All global cases are forced by the third version of Global Soil Wetness Project forcing data (GSWP3; Kim et al., 2012), running in the satellite phenology (SP) mode. The International Land Atmosphere Model Benchmarking (ILAMB; Collier et al., 2018) project has indicated the superior performance of GSWP3 data forcing the CLM5 in the SP-only mode (https://webext.cgd.ucar.edu/I20TR/_build_090817_CLM50SPONLY_CRUNCEP_GSWP3_WFDEI/index.html, last access: 3 March 2020). We conducted a 100-year spin-up using the 1901–1920 climatology before conducting historical period simulations covering 1901–2005. The anomaly forcing under the RCP8.5 scenario on top of the 1982–2005 climatology forces simulations in the projected period.

Table 3 List of simulations conducted for this study.

3 Results

3.1 Excess-ice melt simulations for Lena River delta cryostratigraphies

By the end of the spin-up in the triple-land-unit case, the active layer thickness is 0.85, 0.55, and 0.45 m for the ice-poor terrain, the Holocene ice-wedge terrain, and the Yedoma ice complex, respectively. On the other hand, the active layer thickness for the average-ice single-land-unit case is 0.85 m, which is the same as in the no-excess-ice terrain in the triple-land-unit case. For the average-ice single-land-unit case, a small amount of excess ice (24 kg m⁻²) melts during the spin-up period, resulting in 2.6 cm surface subsidence throughout the grid.

For the Yedoma ice complex, there is very little excess-ice melt in the 1950s, and it stabilizes afterwards until the late 2000s when substantial ice melt and surface subsidence starts to occur. For the Holocene ground ice terrain, there is no excess-ice melt before the late 2010s. By the year 2100, the Yedoma ice complex has exhibited nearly 4 m of surface subsidence, while the Holocene ground ice terrain has about 0.6 m of surface subsidence (Fig. 4). For the average-ice single-land-unit case, the noticeable excess-ice melt and surface subsidence start in the late 2010s, which creates about 0.5 m of surface subsidence by 2100. The magnitude of surface subsidence in the average-ice single-land-unit case is lower than both the Holocene ground ice terrain and the Yedoma ice complex in the triple-land-unit case.

Figure 4 Annual freeze–thaw state for the three terraces for the triple-land-unit case, as well as for the average-ice single-land-unit case.

On the grid scale, the total excess-ice melt is higher in the average-ice single-land-unit case than in the triple-land-unit case (Fig. 5). By the year 2100, the average-ice single-land-unit case has about 30 kg m⁻² more excess-ice melt than the triple-land-unit case. The difference in excess ice on the grid level results from the different volumetric content of excess ice caused by the spatial averaging. In this way, the sub-grid representation of excess ice can potentially also provide more detailed and realistic representation of model variables on the grid level. This is particularly important for the CLM5, which serves as the land component in Earth system models, which requires the coupling between interacting components on the grid level.

Figure 5 Grid-mean excess-ice melt since 1900 for the single-point cases over the Lena River delta with and without the sub-grid excess-ice initialization.

Compared to Westermann et al. (2016), the CLM5 with sub-grid excess ice simulates slightly less (∼20 % less) surface subsidence by 2100 for both the central delta and ice complex. We consider this a good agreement as we do not expect a closer fit of the model results due to substantial differences in the model physics (for example, the Cryogrid3 simulations in Westermann et al., 2106, lack a representation of the subsurface water cycle). What these two studies have in common is the earlier start of excess-ice melt and more surface subsidence in the ice complex than in the central delta. The CLM5 with sub-grid excess ice also exhibits the varying active layer thickness with different excess-ice conditions as Cryogrid3 does. These results suggest that the new model development enables small-scale variability in excess-ice melt and subsequent impacts in agreement with previously published modeling efforts.

3.2 Global projection of permafrost thaw and excess-ice melt

Single-point simulations have shown that the varying excess-ice cryostratigraphies for different land units result in sub-grid variabilities of excess-ice melt and surface subsidence under the warming climate. The same features remain in the sub-grid ice case within the global simulations that excess ice in the low-ice land unit can completely melt out throughout the circum-Arctic permafrost region by the end of the 21st century (Fig. 6). The modeled magnitude of surface subsidence is similar to the ∼10 cm surface subsidence observed in Utqiaġvik and West Dock in the early 21st century (Shiklomanov et al., 2013; Streleskiy et al., 2017). The magnitude of surface subsidence is also comparable to the surface subsidence rate of 1–4 cm per decade on average over the North Slope of Alaska observed by satellite measurements since the 1990s (Liu et al., 2010). In comparison, the absence of surface subsidence for Arctic Alaska modeled by Lee et al. (2014) is due to an overly deep (1 m deep) excess-ice initialization depth. By the year 2100, most ice in the medium-ice land unit melts away in the sub-Arctic region, while there is less ice melt in the colder regions such as the North Slope of Alaska and central Siberia. The high-ice land unit has the greatest surface subsidence among the three because of its high excess-ice content, leading to 2–5 m of surface subsidence by the year 2100.

Figure 6 Maps showing sub-grid surface subsidence (m) in 2000, 2050, and 2100 in the low-, middle-, and high-excess-ice land units in the sub-grid ice case.

The existence of excess ice modulates the thermal regime of permafrost soil and is a major control on permafrost degradation trajectories in a warming climate. Permafrost with excess ice consistently exhibits delayed permafrost degradation compared to the no-ice case (Fig. 7). For the no-ice case modeled by the original CLM5, more than half of the permafrost area undergoes degradation by the end of the 21st century. By 2100, the only areas where permafrost remains are the North Slope of Alaska, northern Canada, and the majority of the land area in northern Siberia. The areas with remaining permafrost in the year 2100 under the RCP8.5 scenarios are substantially larger compared to the CLM4 simulations, in which nearly all permafrost in Eurasia becomes degraded (Lawrence et al., 2012). For the grid-average ice case, the presence of excess ice stabilizes the permafrost thermal regime and thus sustains a larger permafrost area on a global scale in the simulation. For example, permafrost areas in some subarctic regions in eastern and western Siberia, as well as part of the Arctic coastal regions in Yukon Territory, Canada, remain in the grid-average ice case by 2100. Compared to the grid-average ice case, even more permafrost areas are sustained in the sub-grid ice case, most of which are located in southern Siberia. In the subarctic regions in Alaska and northwest Canada as well as part of the central Siberia, permafrost degradation is delayed from the 2040s in the grid ice case to the 2080s in the sub-grid ice case. We emphasize that permafrost is only sustained according to the accepted temperature-based definition (ground material at temperature below zero for 2 consecutive years), but excess ice continuously melts in this process, which energetically is a different mode of permafrost degradation, similar to a negative mass balance of glaciers and ice sheets.

Figure 7 Maps showing the year of completed permafrost degradation (upper set of three maps), as well as the differences between cases (lower set of two maps). The purple color indicates the existence of permafrost in these grid cells by 2100. The difference in years is provided only for grid cell with completed permafrost degradation before 2100.

In the sub-grid ice case, the land units with high-excess-ice contents lead to more grid cells for which permafrost conditions remain in the year 2100 compared to the grid-average ice case. On the other hand, permafrost with excess ice only covers a fraction of a grid cell. Among the permafrost degradation trajectories in the three global simulation cases (Fig. 8), the sub-grid ice case can provide a more detailed picture on the timing of permafrost degradation. Grid cells become “partially degraded permafrost” if land units with excess ice still contain permafrost, which phenomenologically is a more realistic representation that also makes it possible to represent the permafrost distribution in the discontinuous and sporadic permafrost zones. On the other hand, only “fully degraded permafrost” and “remaining permafrost” can be distinguished for the no-ice and grid-average ice cases. Under the warming climate in the 21st century, the existence of excess ice, especially the high content of excess ice, has a stabilizing effect on soil temperature that delays the disappearance of permafrost on the sub-grid level. Therefore, by the year 2100, there are regions with partially degraded permafrost in between intact and degraded permafrost (Fig. 8). For example, in western Siberia, the Pacific coastal area of eastern Siberia, northwestern Canada, and along the Brooks Range in Alaska, taliks form for land units with low excess-ice contents, which leads to partially degraded permafrost regions. Therefore, permafrost degradation exhibits a gradual transition from continuous to discontinuous permafrost and to non-permafrost regions. Some of these regions also encounter substantial surface subsidence in the high-ice land unit (>5 m) (Fig. 6).

Figure 8 Maps of different stages of permafrost degradation diagnosed from the model output by the year 2100.

We further compare the total permafrost area (defined as land units with active layer thickness <6.5 m) in the three cases throughout time. The differences in permafrost area increase from the grid-average ice case and sub-grid ice case to the no-ice case at a rate of 1000 km² yr⁻¹ until 2050 (Fig. 9). After 2050, the area difference of permafrost in the grid-average ice case and no-ice cases rapidly increases, which reaches nearly 1 million km² by 2100. In the sub-grid ice case, the rate of increase remains relatively unchanged after 2050, resulting in an about 0.2 million km² larger permafrost area than that in the no-ice case.

Figure 9 Difference in modeled permafrost area versus time between the sub-grid ice case and no-ice case, as well as between the grid-average ice case and no-ice case.

4 Discussion

The aim of the sub-grid excess-ice representation in the CLM5 is to facilitate long-term global projection of excess-ice melt and surface subsidence in the permafrost regions. Results from our idealized sensitivity experiments (see the Supplement) imply that overly low volumetric content of excess ice, such as the grid-average ice case in this study and that in Lee et al. (2014), produces an overly early start of excess-ice melt and an overly high melting rate. This is because a higher content of excess ice covering a smaller area takes longer to absorb enough heat from the atmosphere to satisfy the latent heat of fusion requirements and start melting. Consequently, a good model performance relies not only on the updated sub-grid representation of excess ice in the global land model, but also on retrieving accurate initial conditions of excess ice. However, the corresponding observational data for both background excess-ice conditions and model evaluation are sparse, especially considering that drastic excess-ice melt as modeled until 2100 is only observed in a few locations today (e.g., Günther et al., 2015). In the following, we discuss the challenges and limitations of the sub-grid excess-ice framework and how this sub-grid representation can potentially help the development of other CLM components. Both single-point and global test simulations in this study have shown that excess ice melt under a warming climate is sensitive to its initialization depth. The active-layer-dependent excess-ice initialization in this study in the global simulation (sub-grid excess-ice case) yields excess-ice melt and surface subsidence rates in the early 2000s that are comparable to observations. The lower depths of the assumed excess-ice body control the termination of excess-ice melt which at the same time determines the onset of talik formation in many permafrost areas. Due to the scarcity of observational data, it is unclear to what extent the cryostratigraphies assumed in our tiling scheme can reproduce the true vertical extent of excess-ice bodies at least in a statistical sense. Even so, we manage to make the prescribed excess-ice condition as close to the previous results as possible. Firstly, our tiling scheme on the large scale strictly follows the CAPS data (Brown et al., 1997) in terms of the volumetric excess-ice content. Furthermore, statistics by Zhang et al. (2000) suggest the ranges of the vertical extent of ice-rich permafrost of 0–2 and 2–4 m respectively for the CAPS classes with low (5 %) and medium (15 %) ice content. Comparatively, the vertical extents of permafrost with excess ice prescribed by our tiling scheme are respectively 1.36 and 3.78 m for the same CAPS classes, both of which lie within the ranges in Zhang et al. (2000). The vertical extent of ice-rich permafrost for the high-ice land unit is much higher than that (4–6 m) in Zhang et al. (2000), but the unmelted part of the ice bodies does not strongly affect the overall rate of excess-ice melt, although the remaining ice can slightly change soil temperature and moisture of the surrounding permafrost. We therefore imply that our high-ice land unit initialization would not induce a strong bias in excess-ice melt projection in the 21st century.

Due to the lack of excess-ice datasets and observational evidence, our projections of excess-ice melt and surface subsidence likely have biases that arise from the need to make empirical estimates and simplifications for the excess-ice initialization scenarios in the global simulation cases. For example, as the CAPS data are mostly based on visible ice bodies (i.e., not pore ice) (Heginbottom et al., 1995), we used the reported volumetric ground ice content in the CAPS data to approximate the volumetric content of excess ice during model initialization. Further, the determination of volumetric contents of excess ice for three land units also results from sparse observations and empirical estimates. The prescribed excess-ice cryostratigraphies ignore ice morphology and the variation in volumetric content of excess ice with soil depth, regarding excess ice as homogeneous within each assigned sub-grid ice content type (low, mid, or high) (Fig. 3, upper panel). For the high-ice land unit, we simplify the cryostratigraphy initialization to Yedoma-type ice, which prescribes overly thick excess-ice bodies out of the Yedoma regions (Schurr et al., 2015). A deficiency in the current version of source code prevents us from initializing non-Yedoma wedged ice for the high-ice land unit where it occurs outside of the Yedoma region. Future versions of our model development will have more freedom in the stratigraphic configuration of excess ice, which will make it possible to prescribe different cryostratigraphies of the same land unit (e.g., the high-ice land unit) for different locations. Because of the above shortcomings in the excess-ice initialization, we do not expect the modeled excess-ice melt in this study to be an adequate representation of reality. Direct ingestion of new or improved observational datasets of excess-ice contents and cyostratigraphies would likely yield more accurate results; however, a spatially distributed global dataset with quantitative information on excess-ice stratigraphies does not exist at present. We emphasize that for a better projection of excess-ice melt, more observational data of excess-ice distribution and surface subsidence are required to further evaluate and validate the new model implementation of excess ice. On the regional scale, Jorgenson et al. (2008) presented a permafrost map of total ground ice volume for the uppermost 5 m of permafrost based on both observations and estimates for Alaska. In addition, O'Neill et al. (2019) compiled permafrost maps for northern Canada by paleographic modeling, mapping the abundances of three types of excess ice. Further improvements of model results depend on additional observationally constrained datasets of excess-ice conditions on the global scale.

The area weights of the excess-ice land units (Table 2) in the global simulation are obtained from the higher-resolution CAPS points located within a CLM grid cell. However, complex landscape development, such as thermokarst ponds, requires knowledge of the meter-scale distribution, for example the extent and geometry of individual ice wedges (Langer et al., 2016; Nitzbon et al., 2019), which cannot be represented with the still-coarse-scale excess-ice classes from the CAPS map. One possible solution to represent this could be to include another layer of sub-grid tiles below the CLM land unit level, where the individual tiles can interact laterally. This would allow for the representation of small-scale permafrost features within a large-scale land unit with a given excess-ice content. An example of how this could work is given by Aas et al. (2019), who simulated both polygonal tundra and peat plateaus with a two-tile interactive setup. This is also similar to the recent representation of hillslope hydrology by Swenson et al. (2019), where sub-grid tiles (on the column level in CLM) were used to represent different elements in a representative hillslope. In the future development of CLM, this could be part of a more generic tiling system where lateral heat and mass fluxes could be switched on and off to represent a wide range of land surface processes that are currently ignored or parameterized in land surface models (LSMs). Fisher and Koven (2020) have discussed the challenges and opportunities in such an adaptive and generic tiling system. We would also advocate for enhancing current tiling schemes in such a direction, which could substantially improve the realism in the representation of permafrost landscapes in LSMs. However, the success of such a tiling approach will rely heavily on the availability of adequate observational data, further highlighting the need for observational efforts and close collaboration between field scientists and modelers.

The more detailed simulation of permafrost degradation trajectory with a sub-grid representation of excess ice also builds more potential for better modeling the permafrost–carbon feedback with biogeochemistry activated (CLM5BGC). Excess ice stabilizes the permafrost thermal regime, therefore altering the rate of carbon release from the permafrost (Schuur et al., 2008). Improved projections of permafrost warming could also enhance modeling of vegetation type changes (e.g., shrub expansion) that determine the nitrogen uptake to the atmosphere (Loranty and Goetz, 2012). On the other hand, the possibility of simulating surface subsidence and excess-ice meltwater formation also opens the possibility of a more accurate representation of wetland formation. The increase in the area of wetland and soil moisture has an impact on the balance of CH₄ and CO₂ release from the permafrost as more organic matter could decompose in an anaerobic pathway (Lawrence et al., 2015; Treat et al., 2015). Compared to the parameterized inundated area simulation in the CLM5 (Ekici et al., 2019), a process-based wetland physics scheme together with the sub-grid representation of excess ice in this study would substantially contribute to the biogeochemical modeling over the circum-Arctic area.

5 Conclusion

This study develops a sub-grid representation of excess ice in CLM5 and examines the impacts of the existence and melting of excess ice on the sub-grid scale in a warming climate. Extra land units duplicated from the natural vegetated land unit in the CLM sub-grid hierarchy make it possible to prescribe up to three different excess-ice conditions in each grid cell with permafrost.

A test over the Lena River delta showcases that the sub-grid representation of excess ice can retrieve the sub-grid variability of annual thaw–freeze state and the excess-ice melt and surface subsidence through time. On the other hand, initializing excess ice homogeneously throughout the grid cell produces a smaller stabilization effect of excess ice to the permafrost thermal regime and the local surface subsidence under a warming climate. With a tiling scheme ingesting a global dataset of excess-ice condition into the CLM surface data, our model development shows the capability of portraying more details on simulating permafrost degradation trajectories. As excess ice thermally stabilizes the permafrost on the sub-grid scale, permafrost degrades with a trajectory from continuous permafrost to discontinuous permafrost and finally to a permafrost-free area. The modeled global pattern of permafrost therefore exhibits regions of discontinuous permafrost as the transition zone between the continuous permafrost and degraded permafrost.

This study, for the first time, used an ESM to project excess-ice melt and surface subsidence and permafrost degradation with sub-grid variability. The approach of duplicating tiles at the land unit level instead of the column level allows more freedom for further developments in this direction. Furthermore, the new CLM tiling hierarchy has much more potential than representing more accurate excess-ice physics as examined in this study. The accuracies of predicted excess-ice melt and surface subsidence trends are limited at present by the available global-scale dataset and studies on excess ground ice conditions; thus further advancement of excess-ice modeling will rely on new or improved observational studies or datasets of the excess ground ice conditions at the global scale. The model development in our study, therefore, lays the foundation for further advances focusing on excess-ice modeling and other processes in the CLM framework that could benefit from an improved sub-grid representation.