2.1 Terminology for spatial scales

Throughout this paper we used a consistent terminology to refer to different characteristic length scales of landscape features and processes. This terminology is summarized in Table 1.

2.2 Model description

2.2.1 Multi-scale tiling to represent spatial heterogeneity

We used the concept of laterally coupled tiles (Langer et al., 2016; Aas et al., 2019; Nitzbon et al., 2019) to represent subgrid-scale spatial heterogeneities of permafrost terrain. In general, the tiling concept involves the partitioning of real-world landscapes into a certain number of characteristic units, which are associated with the major surface (or subsurface) heterogeneities found in the landscape. Each of these units is then represented by a single tile in a permafrost model, and multiple tiles can interact through lateral exchange processes. The tiling approach thus allows for simulating subgrid-scale heterogeneities and lateral fluxes in macro-scale models like LSMs without discretizing extensive landscape domains on a high-resolution mesh and hence keeps computational costs at a reasonable level.

For the present study, we applied the tiling concept at multiple scales to represent common spatial heterogeneities of ice-rich continuous permafrost lowlands (see Fig. 1). At the micro-scale, these landscapes are typically characterized by ice-wedge polygons which give rise to a regular patterning of the landscape (see Appendix A for an example site from northeastern Siberia). Here, we adopted the three-tile approach by Nitzbon et al. (2019), which partitions polygonal tundra lowlands into polygon centres, polygon rims, and inter-polygonal troughs. At the meso-scale, tundra lowlands are characterized by gently sloped terrain and often feature abundant thaw lakes that formed in the past due to thermokarst activity (see Appendix A). Here, we used three meso-scale tiles to represent a low-gradient slope which is efficiently drained at its lowest elevation (“drainage point”).

Figure 1Schematic illustration of ice-rich permafrost lowlands featuring spatial heterogeneity on different scales. At the meso-scale, the terrain is gently sloped and features larger landforms like thaw lakes. At the micro-scale, ice-wedge polygons entail a periodic patterning of the terrain which can have a low-centred (LCPs) or high-centred (HCPs) topography, depending on the grade of degradation. In this study, a tile-based modelling approach is pursued to represent these heterogeneities and investigate their effect on projections of landscape evolution and permafrost degradation.

Download

Figure 2 provides an overview of the four model setups which were investigated in this study. These setups differ with respect to the number of micro-scale tiles (N^μ) and meso-scale tiles (N^m), the product of which amounts to the total number of tiles ($N=N^{\mathit{\mu }}\cdot N^{\text{m}}$).

Figure 2Overview of the different tile-based model setups used to represent heterogeneity at micro- and meso-scales. An external reservoir (blue triangle) at a fixed elevation (e_res) prescribes the hydrological boundary conditions. (a) The single-tile setup used only one tile (homogeneous, H) and does not reflect any subgrid-scale heterogeneity. (b) The polygon setup reflects micro-scale heterogeneity of topography and ground ice distribution associated with ice-wedge polygons via three tiles: a polygon centre (C), polygon rim (R), and inter-polygonal trough (T). The parameters e_R and e_T indicate the (initial) elevation of rims and troughs relative to the centre. The micro-scale tiles interact through lateral fluxes of heat, water, snow, and sediment. (c) The low-gradient slope reflects meso-scale heterogeneity of an ice-rich lowland via three tiles (H^o,m,i), of which the outer one (H^o) is connected to a draining reservoir. The intermediate (H^m) and inner tiles (Hⁱ) are located at distances D^m along a low-gradient slope S^m. The meso-scale tiles interact through lateral surface and subsurface water fluxes. (d) The low-gradient polygon slope incorporates heterogeneities at both micro- and meso-scales via a total of nine tiles.

Download

• a.

  Single-tile. This is the simplest case (N^μ=1, N^m=1; Fig. 2a), which reflects homogeneous surface and subsurface conditions across all spatial scales via only one tile (H). Water can drain laterally into an external reservoir at a fixed elevation (e_res). This setup emulates the one-dimensional representation of permafrost in LSMs and corresponds to the setup used for the simulations conducted by Westermann et al. (2016).

• b.

  Polygon. This setting (N^μ=3, N^m=1; Fig. 2b) reflects the micro-scale heterogeneity associated with ice-wedge polygons via three tiles: polygon centres (C), polygon rims (R), and inter-polygonal troughs (T). The heterogeneity of the surface topography is expressed in different initial elevations of the soil surface of the tiles. The subsurface stratigraphies of the tiles differ with respect to the depth (d_x) and amount (θ_x) of excess ice, reflecting the subsurface ice-wedge network which is linked to the polygonal pattern at the surface. Lateral fluxes of heat, water, snow, and sediment are enabled among the tiles, and the trough tile is connected to an external water reservoir. This model setup has previously been used by Nitzbon et al. (2019, 2020a).

• c.

  Low-gradient slope. This setup (N^μ=1, N^m=3; Fig. 2c) reflects a meso-scale gradient of slope S^m, for which the outer (i.e. downstream) tile (H^o) is well-drained into an external reservoir. The intermediate (H^m) and inner (Hⁱ) tiles represent the landscape upstream of the outer tile at constant distances (D^m). The setting assumes a translational symmetry perpendicular to the direction of the slope; i.e. each tile represents the same areal proportion at the meso-scale. Lateral surface and subsurface water fluxes are enabled among the meso-scale tiles, while lateral heat fluxes as treated by Langer et al. (2016) were not considered at this scale.

• d.

  Low-gradient polygon slope. This setup (N^μ=3, N^m=3; Fig. 2d) reflects a meso-scale slope featuring ice-wedge polygons at the micro-scale. Each meso-scale tile includes three micro-scale tiles (C, R, and T), corresponding to the polygon setup described above. However, only the trough tile of the outer polygon (T^o) is connected to an external reservoir (well-drained), while the intermediate and inner trough tiles are hydrologically connected along the meso-scale slope.

In the initial configurations, the landscapes were assumed to be undegraded (i.e. the ice-wedge polygons were low-centred, LCPs), and no thaw lakes (water bodies, WBs) were present along the meso-scale slope. However, the model allowed thermokarst features like high-centred polygons (HCPs) and thaw lakes to develop dynamically as a consequence of excess ice melt at different scales.

2.3 Model settings and simulations

2.3.2 Soil stratigraphies and ground ice distribution

The subsurface composition of all tiles was represented via a generic soil stratigraphy (Table 2). The stratigraphy was based on previous studies that applied CryoGrid 3 to the same study area (Nitzbon et al., 2019, 2020a). It consists of two highly porous layers of 0.1 m thickness each which reflect the surface vegetation and an organic-rich peat layer. Below those, a mineral layer with silty texture follows. An excess ice layer of variable total ice content θ_i extends from a variable depth d_x down to a depth of 10.0 m. Between the variable excess ice layer and the mineral layer, we assumed an ice-rich intermediate layer of 0.2 m thickness and a total ice content of 0.65. Below the variable excess ice layer follow an ice-poor layer and bedrock down to the end of the model domain. The default ice content assumed for the excess ice layer of homogeneous tiles (no micro-scale heterogeneity) was ${\mathit{\theta }}_{\text{i}}^{\text{H}}=\mathrm{0.75}$, and the default depth of this layer was d_x=0.9 m. To reflect the heterogeneous excess ice distribution associated with ice-wedge polygons, different ice contents were assumed for the excess ice layers of polygon centres (${\mathit{\theta }}_{\text{i}}^{\text{C}}=\mathrm{0.65}$), rims (${\mathit{\theta }}_{\text{i}}^{\text{R}}=\mathrm{0.75}$), and troughs (${\mathit{\theta }}_{\text{i}}^{\text{T}}=\mathrm{0.95}$). However, the area-weighted mean ice content of ice-wedge polygon terrain is identical to the default value assumed for homogeneous tiles (Table 3).

Table 2Generic soil stratigraphy used to represent the subsurface composition of all tiles. An ice-rich layer of variable ice content (θ_i) is located at variable depth (d_x) from the surface. Note that the effective excess ice content (θ_x) is linked to the total ice content (θ_i) and the natural porosity (ϕ_nat) via the relation given in Eq. (2). The soil texture is used to parametrize the freezing-characteristic curve of the respective layer (see Westermann et al., 2016, for details).

Download Print Version | Download XLSX

Table 3Overview of the model parameters for different representations of the micro-scale. Note that on average the polygon tiles (C, R, and T) feature the same excess ice content and depth as the homogeneous tile (H). Setting the area of the homogeneous tile to A^μ=1 m² is an arbitrary choice, as it does not affect the magnitude of the lateral fluxes due to the assumption of translational symmetry.

Download Print Version | Download XLSX

2.3.5 Simulations

Table 4 provides an overview of all simulations which were conducted for the four different model setups introduced in Sect. 2.2.1. For both the single-tile and the polygon setups we conducted four model runs in which we varied the elevation of the external water reservoir (e_res) in order to reflect a broad range of hydrological conditions. For the low-gradient slope and the low-gradient polygon slope setups we conducted two model runs in which we varied the gradient of the slope (S^m) in order to reflect essentially flat (S^m=0.001) as well as gently sloped terrain (S^m=0.01). In all slope simulations, one tile (H^o and T^o, respectively) was connected to an external reservoir (see Fig. 2c and d; $e_{\text{res}}=-\mathrm{10.0}$ m).

Table 4Overview of the parameter variations in the simulations which were conducted for the four different model setups (see Fig. 2). Note that the last value for e_res in the single-tile and polygon setups – which reflects poorly drained conditions – was chosen to be 0.1 m less than the mean initial elevation of the respective setting (see Table 3).

Download Print Version | Download XLSX

The subsurface temperatures of each model run were initialized in October 1999 with a typical temperature profile for that time of year which was based on long-term borehole measurements from the study area (Boike et al., 2019). Using multi-year spin-up periods did not result in significant changes to the near-surface processes investigated in this study so that biases related to the initial subsurface temperature profile can be excluded. The analysed simulation period was the twenty-first century from January 2000 until December 2099.

3.1 Single-tile

Figure 3 shows the landscape configuration for selected years throughout the twenty-first century under RCP8.5 for the single-tile setup (see Fig. 2a) under four different hydrological boundary conditions. Until the middle of the simulation period (2050) no excess ice melt and associated ground subsidence occur, irrespective of the hydrological conditions (Fig. 3a–c, all columns; Fig. 8a). However, the maximum thaw depths are steadily increasing during that period (see Fig. 7a), with negligible differences between the four runs under different hydrological conditions. Between 2050 and 2075 excess ice melt occurs in all single-tile simulations, leading to ground subsidence of about 0.2 m in the two well-drained cases (Fig. 3d, first and second column), about 0.1 m in the intermediate case (Fig. 3d, third column) and about 0.5 m in the poorly drained case, where excess ice melt leads to the formation of a shallow surface water body (Fig. 3d, fourth column). By 2100 total ground subsidence amounts to about 0.8 m in the well-drained simulations (Fig. 4e, first and second column), and about 1.0 m in the simulation with $e_{\text{res}}=-\mathrm{0.5}$ m, where it is accompanied by the formation of a shallow surface water body (Fig. 4e, third column). In the poorly drained simulation ($e_{\text{res}}=-\mathrm{0.1}$ m), excess ice melt proceeds fastest, causing the surface water body to deepen and to reach a depth of about 2 m by 2100; a talik of about 1.5 m thickness has formed underneath the water body by the end of the simulation period (Fig. 3e, fifth column).

Figure 3Landscape configurations in selected years (a–e) simulated with the single-tile setup (see Fig. 2a) under four different hydrological boundary conditions (left to right). e_res denotes the elevation of the external water reservoir so that lower values correspond to better drainage.

Download

Figure 4Landscape configurations in selected years (a–e) simulated with the polygon setup (see Fig. 2b) under four different hydrological boundary conditions (left to right). e_res denotes the elevation of the external water reservoir so that lower values correspond to better drainage. The width of the three tiles corresponds to their areal fractions. Brighter colours reflect higher excess ice contents of the permafrost (see Table 3). Labels indicate the state of the micro-topography according to Eqs. (D1) to (D5)

Download

Overall, the simulation results indicate that permafrost degradation is strongest as soon as a limitation of water drainage results in the formation of a surface water body. The presence of surface water changes the energy transfer at the surface in different ways. First, it reduces the surface albedo (from 0.20 for barren ground to 0.07 for water), resulting in a higher portion of incoming shortwave radiation. Second, water bodies have a high heat capacity which slows down their freeze-back compared to soil. Third, the thawed saturated deposits beneath the surface water body have a higher thermal conductivity compared to unsaturated deposits, which allows heat to be transported more efficiently from the surface into deeper soil layers.

Figure 5Landscape configurations in selected years (a–e) simulated with the low-gradient slope setup (see Fig. 2c) for two different slope gradients (S^m). The outer tile is connected to an external reservoir ($e_{\text{res}}=-\mathrm{10.0}$ m) which allows for efficient drainage, while the intermediate and inner tiles can drain only to adjacent tiles.

Download

Figure 6Landscape configurations in selected years (a–e) simulated with the low-gradient polygon slope setup (see Fig. 2d) for two different slope gradients (S^m). The trough tile of the outer polygon is connected to an external reservoir ($e_{\text{res}}=-\mathrm{10.0}$ m) which allows for efficient drainage. The trough tiles of adjacent polygons are connected to each other and constitute the drainage channels for the entire model domain. Brighter colours reflect higher excess ice contents of the permafrost (see Table 3). Labels indicate the state of the micro-topography according to Eqs. (D1) to (D5).

Download

Figure 7Time series of the maximum thawed-ground thickness throughout the simulation period for all model setups (a–d) and parameter settings. e_res denotes the elevation of the external water reservoir so that lower values correspond to better drainage. S^m denotes the gradient of the meso-scale slope. To derive the maximum thawed-ground thickness, we first took the maximum annual thaw depth of each tile and then averaged these, weighted according to the areal fractions of the different tiles. We then took the 11-year running mean to obtain a smoothed time series. Dashed grey lines indicate the selected years for which the landscape configuration is explicitly shown in Figs. 3 to 6.

Download

Figure 8Time series of the accumulated ground subsidence throughout the simulation period for all model setups (a–d) and parameter settings. e_res denotes the elevation of the external water reservoir so that lower values correspond to better drainage. S^m denotes the gradient of the meso-scale slope. To obtain the accumulated ground subsidence, we first took the difference between the soil surface elevation in each year and the soil surface elevation at the start of the simulation period. We then averaged these differences, weighted according to the areal fractions of the different tiles. Dashed grey lines indicate the selected years for which the landscape configuration is explicitly shown in Figs. 3 to 6.

Download

Interestingly, during the initial phase of excess ice melt which occurs between 2050 and 2075, our simulations suggest a non-monotonous dependence of permafrost degradation on the drainage conditions. This is indicated by the fact that the lowest degradation both in terms of maximum thaw depth and accumulated subsidence is simulated for the intermediate case with $e_{\text{res}}=-\mathrm{0.5}$ m (see Figs. 7a and 8a). This can likely be attributed to contrasting effects of the hydrological regime in the active layer on thaw depths. When the near-surface ground is unsaturated (as in the simulations with $e_{\text{res}}=-\mathrm{1.0}$ m and $e_{\text{res}}=-\mathrm{10.0}$ m), the highly porous organic-rich surface layers (see Table 2) have an insulating effect on the ground below due their low thermal conductivity. At the same time, less heat is required to melt the ice contained in the mineral soil layers whose ice content corresponds to the field capacity than would be needed if their pore space was saturated with ice. In the intermediate case with $e_{\text{res}}=-\mathrm{0.5}$ m, the combination of dry, insulating near-surface layers and ice-saturated mineral layers beneath leads to the lowest thaw depths and hence the slowest initial permafrost degradation. However, as soon as a surface water body forms in that simulation (between 2075 and 2100), the positive feedback on thaw described above takes over, resulting in stronger degradation by 2100 compared to the well-drained settings ($e_{\text{res}}=-\mathrm{1.0}$ m and $e_{\text{res}}=-\mathrm{10.0}$ m), for which no surface water body forms during the simulation period.

In summary, the single-tile simulations illustrate the non-trivial relation between the ground hydrological regime and permafrost thaw and demonstrate the positive feedback associated with surface water body formation resulting from excess ice melt.

3.2 Polygon

Figure 4 illustrates the landscape evolution throughout the twenty-first century under RCP8.5 for the simulations with the polygon setup, i.e. with a representation of micro-scale heterogeneity typical for ice-wedge polygonal tundra. In addition, it shows the geomorphological state of the polygon micro-topography, i.e. whether it is low-centred (LCP), intermediate-centred (ICP), high-centred with inundated troughs (HCPi), high-centred with drained troughs (HCPd), or covered by a surface water body (WB), according to Eqs. (D1) to (D5) in Appendix D. Initially, all simulations feature undegraded ice-wedge polygons with a low-centred micro-topography (Fig. 4a).

Between 2000 and 2025 a shallow surface water body of about 0.5 m depth forms in the poorly drained simulation (e_res=0.0 m) as a result of excess ice melt in the centre, rim, and trough tiles (Fig. 4b, fourth column). The bottom of the water body has a high-centred topography. The landscape configuration does not change much until 2050 (Fig. 4c, fourth column), and the maximum thaw depths and accumulated ground subsidence increase more slowly than in the initial 2 to 3 decades (Figs. 7b and 8b, purple lines). We explain this interim stabilization with the subaqueous transport of sediment from the centres to the rims and from the rims to the troughs, where the additional sediment has an insulating effect on the ice wedge underneath. After 2050 excess ice melt proceeds faster, causing the surface water body to deepen, reaching a depth between about 1 m (centre tile) and more than 2 m (trough tile) by 2075 (Fig. 4 d, fourth column). We explain the acceleration of permafrost degradation at the beginning of the second half of the simulation period (Figs. 7b and 8b, purple lines) by a combination of additional warming from the meteorological forcing and positive feedbacks due to the surface water body (as explained for the single-tile simulation in Sect. 3.1). Until 2100 excess ice melt proceeds further, resulting in a water body depth of 2–4 m and the formation of an extended talik underneath (Fig. 4e, fifth column). By the end of the simulation period, the lake is not entirely bottom-freezing in winter. This constitutes another positive feedback on the degradation rate, since the heat from the ground can be transported much more efficiently through bottom-freezing water bodies than through those which do not bottom-freeze, due to the higher thermal conductivity of ice compared to that of unfrozen water.

The remaining three simulations show a similar landscape evolution during the first half of the simulation period (2000–2050), with excess ice melt occurring only in the trough tiles, resulting in a transition from LCP to ICP micro-topography (Fig. 4a–c, first to third column). In contrast to the single-tile setup, the polygon simulations show a monotonous relation between permafrost degradation and drainage conditions during this period. Higher elevations of the external water reservoir (i.e. poorer drainage) result in faster thaw-depth increase (Fig. 7b) and earlier onset of ground subsidence (Fig. 8b). Between 2050 and 2075, substantial excess ice melt occurs in the trough and rim tiles, involving a transition from the ICP to a pronounced HCP micro-topography (Fig. 4d, first to third column). In addition to the positive feedback through surface water formation, the degradation rate in the polygon simulations is accelerated by a positive feedback through lateral snow redistribution. Snow is deposited preferentially in micro-topographic depressions such as initially subsided troughs where it improves the insulation and traps heat in the subsurface during winter, enabling faster and deeper thaw in the subsequent summer. Note that this positive feedback is not represented in the single-tile simulations.

By 2075 the influence of the drainage conditions on the surface water coverage is most pronounced. While a water body extending over all tiles has formed in the poorly drained simulation (e_res=0.0 m), the polygon centres are (still) elevated above the water table in the intermediate case with $e_{\text{res}}=-\mathrm{0.5}$ m. For $e_{\text{res}}=-\mathrm{1.0}$ m, surface water is found only in the trough tile, while all surface water is drained in the simulation with $e_{\text{res}}=-\mathrm{10.0}$ m. Irrespective of the hydrological boundary conditions, permafrost degradation continues in the final decades until 2100. While the high-centred polygon in the well-drained simulation ($e_{\text{res}}=-\mathrm{10.0}$ m) remains free of surface water (Fig. 4e, first column), surface inundation increases in the runs with intermediate hydrological conditions (Fig. 4e, second and third column), resulting in the formation of a water body of 1–3 m depth underlain by a talik in the simulation with $e_{\text{res}}=-\mathrm{0.5}$ m. However, the water body is (still) bottom-freezing in 2100, in contrast to the water body with floating ice in the simulation with e_res=0.0 m.

Overall, the polygon simulations reveal a marked dependence of the pathways of landscape evolution and associated permafrost degradation on the hydrological conditions, consistent with previous results under recent climatic conditions (Nitzbon et al., 2019).

3.3 Low-gradient slope

While the single-tile and polygon simulations illustrated general feedbacks due to excess ice melt and micro-scale lateral fluxes, they cannot reflect heterogeneous landscape dynamics at the meso-scale. The simulations with the low-gradient slope setup illustrate how meso-scale lateral water fluxes in gently sloped terrain give rise to diverse pathways of landscape evolution and how these depend on the slope gradient (S^m).

Figure 5 shows the landscape evolution throughout the twenty-first century for the simulations with the low-gradient slope setup, i.e. with a meso-scale representation according to a slope but without micro-scale heterogeneities. Irrespective of the slope gradient, the simulated evolution of the outer tiles is very similar to that of the well-drained single-tile simulations ($e_{\text{res}}=-\mathrm{10.0}$ m and $e_{\text{res}}=-\mathrm{1.0}$ m) throughout the entire simulation period (compare Fig. 3, first and second columns, with Fig. 5, first and fourth column). The outer tiles are stable through the first half of the simulation period. During the second half, excess ice melt sets in, and the ground subsides at an increasing rate, reaching an accumulated subsidence of about 0.8 m by 2100. The similarity to the well-drained single-tile simulations can be explained by the fact that the outer tile is very efficiently drained ($e_{\text{res}}=-\mathrm{10.0}$ m) such that the lateral water input from the intermediate tile is directly routed further into the external reservoir. Hence, the “upstream” influence on the outer tile becomes negligible.

For both slope gradients, the evolution of the intermediate tile is similar to that of the inner tile. However, the evolution of the two tiles is different among the two simulations for different slope gradients. For the lower slope gradient (S^m=0.001), which reflects an essentially flat landscape, melting of excess ice and associated ground subsidence occur during the first half of the simulation period, and a shallow layer of surface water (about 0.2 m) forms in these tiles (Fig. 5a–c, second and third column). By 2050, the soil surface in the intermediate and inner tiles has subsided below the surface elevation of the outer tile such that the initial slope is dissipated. After 2050 excess ice melt proceeds faster, and the surface water body in the intermediate and inner tiles reaches a depth of almost 1 m by 2075 (Fig. 5 d, second and third column). The permafrost table in these tiles is lowered by about 2 m relative to its initial position. Between 2075 and 2100 a talik forms beneath the water body in the two inland tiles (Fig. 5 e, second and third column), which by 2100 reaches a depth of about 2.0 m. The ground subsidence in the outer tile during the second half of the simulation enables surface and subsurface water transport from the intermediate tile into the external reservoir and thus leads to a lowering of the water level of the water body in the intermediate and inner tiles (Fig. 5d–e). By 2100 the water body depth is about 1.2–1.5 m, which is significantly lower than the water body of about 2.0 m depth which forms in the poorly drained single-tile simulations (Fig. 3 e, fourth column). This difference in water body depths is despite the fact that excess ice melt and surface water formation set on several decades earlier in the low-gradient slope simulation than in the poorly drained single-tile simulation. These different pathways of water body and talik formation illustrate that there is no linear relationship between the water body depth and talik extent, since lateral interactions at the meso-scale can give rise to different transient dynamics.

For the higher slope gradient (S^m=0.01), which reflects moderately sloped terrain, initial excess ice melt occurs in the intermediate and inner tiles between 2000 and 2025, and by 2050 these have subsided by about 0.2 m (Fig. 5, fifth and sixth column), similar to the simulation with S^m=0.001. By 2075, ground subsidence has reached about 0.5 m in the intermediate and inner tiles. The overall shape of the slope does not change much because excess ice melt occurs also in the outer tile during that time. By 2100, the intermediate and inner tile have subsided by about 1.2 m, and the outer tile has subsided by about 0.8 m, resulting in a slightly concave slope. The presence of a thin layer of surface water in the intermediate and inner tiles in the years 2025, 2075, and 2100 indicates a wetter hydrological regime compared to the outer tile. However, as the moderate slope is sustained throughout the simulation, it facilitates drainage of surface water and precludes the formation of a surface water body as is the case in the simulation with the flatter topography (S^m=0.001). We suspect, however, that if the simulations would be prolonged beyond 2100, the excess ice melt in the intermediate and inner tiles would likely continue at a faster rate than in the outer tile and that this would result in the reversal of the initial slope and promote the formation of an inland surface water body.

While the overall dynamics of the three-tile slope simulations could also be reflected using a two-tile setup, we would like to point to the small but significant differences between the evolution of the intermediate tile and the inner tile. For both slope gradients, the intermediate tile shows a slightly stronger degradation rate than the inner tile. This can be explained due to the additional water input from the inner tile which sustains saturated conditions in the near-surface soil layers during the thawing season. The inner tile in turn lacks this lateral water input and is thus more likely to develop dry, thermally insulating conditions in the near-surface soil.

Overall, the slope gradient (S^m) has a similar influence on permafrost degradation in the low-gradient slope setup, as the elevation of the external reservoir (e_res) has in the single-tile simulations (see Figs. 7a and c and 8a and c). This can be attributed to the direct influence of the slope gradient on drainage efficiency. Due to the positive feedbacks related to the formation of a surface water body, the overall permafrost degradation in the setting with S^m=0.001 is almost twice as much as in the setting with S^m=0.01 for which no water body forms. We note that the maximum thawed ground and the accumulated ground subsidence by 2100 in both low-gradient slope simulations are within those simulated for the two extreme settings in the single-tile simulations.

3.4 Low-gradient polygon slope

Finally, we consider the most complex setup, which reflects gently sloped polygonal tundra via nine tiles, incorporating heterogeneities on both the micro- and the meso-scale. Figure 6 shows different stages of the simulated landscape evolution in the low-gradient polygon slope setup throughout the twenty-first century for two different slope gradients S^m.

Like the outer tiles in the low-gradient slope setup, the outer polygons show an evolution which is very similar to that of the well-drained (single-) polygon simulation with $e_{\text{res}}=-\mathrm{10.0}$ m, irrespective of the slope gradient (compare Fig. 6, first and fourth column, with Fig. 4, first column): until 2050 initial excess ice melt occurs in the trough tiles, entailing a transition from LCP to ICP topography; during the second half of the simulation period, permafrost degradation proceeds at an increasing rate and brings about a transition to a pronounced HCP topography with the troughs having subsided by about 2 m. Again, the similarity between the outer-polygon evolution and the well-drained single-polygon evolution can be explained by the efficient drainage of the outer-polygon troughs which leads to a negligible upstream influence.

For the lower slope gradient (S^m=0.001), the intermediate and inner polygons show an evolution which is similar to each other but different from each of the single-polygon simulations (compare Fig. 4 with Fig. 6, second and third column). Between 2000 and 2025 excess ice melt in the rim and trough tiles leads to ponding of surface water in the intermediate and inner polygons (Fig. 6b, second and third column) with the soil surface of the polygon centre being close to or elevated above the water table. Besides a lowering of the water table, the configuration of these polygons does not change much until 2050 (Fig. 6c, second and third column) when both have developed an HCPi topography, i.e. inundated rims and troughs which subsided below the level of the centres. Between 2050 and 2075, excess ice melt continues, leading to a pronounced high-centred topography (Fig. 6d, second and third column). Meanwhile, the surface water disappears from the rim and trough tiles, as a consequence of marked excess ice melt in the outer-polygon rim and troughs tiles during that period (Fig. 6d, first column). The soil surface elevation of the trough tiles increases from the outer towards the inner polygon, allowing for efficient drainage of the entire landscape. Between 2075 and 2100 all three polygons show a similar evolution, resulting in pronounced high-centred topographies with troughs about 2 m deep and rims that subsided by about 1 m in total. Note that between 2075 and 2100 the troughs of the intermediate and inner polygons subsided more than that of the outer polygon, resulting in an decreased drainage efficiency. Hence, surface water starts to pond again in the intermediate and inner troughs (Fig. 6e, second and third column). Overall, the low-gradient polygon slope simulation with S^m=0.001 shows a transient phase in which the formation of a water body is initiated, but its growth is impeded due to the formation of efficient drainage pathways through the subsiding troughs in the outer (“downstream”) polygon.

The simulated landscape evolution for the moderate slope gradient (S^m=0.01) differs from that for the lower slope gradient during the first half of the simulation period. Until 2025, excess ice melt is restricted to the trough tiles, leading to a transition from LCP to ICP topography in all three polygons (Fig. 6b, fifth and sixth column). In the intermediate and inner polygons, permafrost degradation proceeds faster than in the outer polygon such that these develop a HCPd topography by 2050. The faster degradation can be explained by the overall wetter hydrological regime compared to the very efficiently drained outer polygon. However, the meso-scale slope still allows for sufficient drainage of surface water such that – in contrast to the simulation with S^m=0.001 – the ponding of surface water is mostly precluded (Fig. 6c, fifth and sixth column). During the second half of the simulation period, the evolution of the polygons is very similar for both slope gradients. Ice-wedge degradation proceeds at an increasing pace between 2050 and 2100, resulting in polygons with pronounced HCP topography. While the troughs of the outer polygons are efficiently drained (HCPd, Fig. 6e, first and fourth column), some surface water pools in the troughs of the intermediate and inner polygons (HCPi, Fig. 6 second, third, fifth, and sixth column).

The overall similarity among the two low-gradient polygon slope simulations for different slope gradients is also reflected in the time series of maximum thaw depths (Fig. 7d) and accumulated ground subsidence (Fig. 8d). During a transient phase spanning roughly from 2010 until 2070, the permafrost degradation in the simulation with S^m=0.001 is slightly higher than in the simulation with S^m=0.01. Before and after this transient period, the maximum thaw depths and the accumulated ground subsidence are almost identical for the two settings. This is a qualitative difference to the low-gradient slope simulations without micro-scale heterogeneity, where the simulation with S^m=0.001 resulted in almost twice the amount of permafrost degradation compared to the simulation with S^m=0.01 (Figs. 7c and d and 8c and d). The projected permafrost degradation in both low-gradient polygon slope simulations is confined within the range spanned by the projections for the polygon setup (see Figs. 7b and d and 8b and d). However, by 2100 the projected maximum thaw depths and the accumulated ground subsidence are very close to the single-polygon simulation under well-drained conditions ($e_{\text{res}}=-\mathrm{10.0}$ m).

While the low-gradient polygon slope setup provides the most complex representation of landscape heterogeneities – and hence potentially captures most feedback processes – the projected pathways of landscape evolution and permafrost degradation showed little sensitivity to the initial gradient of the meso-scale slope, and, overall, the projected landscape response was more gradual than in the other setups which feature abrupt initiation or acceleration of permafrost degradation (Figs. 7 and 8). We also note that the low-gradient polygon slope setup is the only one for which the formation of a talik was not projected for any of the tiles and for any of the parameter settings (Fig. 6a–e).

4.1 Simulating permafrost landscape degradation pathways under consideration of micro- and meso-scale heterogeneities

The multi-scale tiling approach introduced in this study (see Sect. 2.2.1) facilitates the simulation of degradation pathways of ice-rich permafrost landscapes in response to a warming climate, under consideration of feedbacks emerging from heterogeneities and lateral fluxes on subgrid-scales. The following qualitative assessment based on observed landscape changes across the Arctic demonstrates that the model is able to realistically represent the dominant pathways of landscape evolution and important feedback processes induced by permafrost thaw and excess ice melt. The most idealized cases considered were single-tile simulations without heterogeneities on the micro- or meso-scale, corresponding to previous simulations with models that incorporate a representation of excess ice melt and ground subsidence (Lee et al., 2014; Westermann et al., 2016). The simulations with this setup reflect the fundamental modes of landscape change due to excess ice melt (i.e. thermokarst) (Kokelj and Jorgenson, 2013). These changes are ranging from the gradual subsidence of a well-drained “upland” setting to the formation of thaw lakes in a water-logged “lowland” setting. The marked sensitivity of the simulations to the prescribed hydrological conditions underlines that landscape evolution and permafrost degradation crucially depend on whether water from melted excess ice drains or ponds at the surface. On the one hand, the simulations reveal a positive feedback on permafrost degradation through ponding of water from melted excess ice at the surface. The presence of surface water alters the energy exchange with the atmosphere and increases the heat input into the ground (via albedo and thermal properties), which in turn allows for additional excess ice melt. This positive feedback is in agreement with previous observations (Connon et al., 2018; O'Neill et al., 2020) and modelling results (Rowland et al., 2011; Langer et al., 2016; Westermann et al., 2016). On the other hand, simulations under well-drained conditions favoured stabilization due to the improved insulation by dry organic surface layers, which is also in agreement with observations (Göckede et al., 2017, 2019) and modelling results (Martin et al., 2019; Nitzbon et al., 2019).

In the polygon setup, micro-scale heterogeneities associated with ice-wedge polygons and representing lateral fluxes of heat, water, snow, and sediment were considered. The simulated landscape evolution pathways under a warming climate correspond to those simulated and discussed by Nitzbon et al. (2020a) (for the stratigraphy of “Holocene deposits”). They range from the rapid formation of a deep surface water body to the development of high-centred polygons with a pronounced relief. While observations of initial ice-wedge degradation and the development of high-centred polygons are reported across the Arctic (Farquharson et al., 2019; Liljedahl et al., 2016), the likelihood of thaw lake formation in ice-wedge terrain is debated and likely depends on the overall wedge-ice content (Kanevskiy et al., 2017). Despite the same ground ice distribution in terms of excess ice content (θ_x) and depth (d_x) on average (see Table 3), the simulated permafrost degradation in the polygon settings exceeds that in the respective single-tile simulations, suggesting that positive feedbacks (e.g. preferential accumulation of snow and water in troughs) exceed the influence of stabilizing feedbacks (e.g. slumping of sediment from rims into troughs). Hence, the simulations suggest that the presence of micro-scale heterogeneities can crucially affect the timing and the rate of permafrost degradation in ice-rich lowlands. Ultimately, both the single-tile and the polygon setups are limited by the prescription of static hydrological boundary conditions which do not allow transient changes between different pathways but rather prescribe the possible landscape evolution a priori.

In the low-gradient slope setup, the hydrological conditions in the inland part of the model domain can develop dynamically, and the water level in a given part of the landscape adapts to topographic changes in adjacent parts due to the representation of meso-scale lateral water fluxes. From an almost flat initial topography (S^m=0.001, Fig. 5 left), the formation of a thaw lake in the inland has been simulated, succeeded by the gradual drainage of the lake in response to ground subsidence in the outer part. This setup thus captures in a stylized way the competing mechanisms of thaw lake formation and growth on the one hand and (gradual) lake drainage on the other hand. Gradual or partial drainage of mature thaw lakes has been observed (Morgenstern et al., 2013) and simulated (Kessler et al., 2012), and we note that also the thermokarst lakes in the southeastern and southwestern part of Samoylov Island are indicative of gradual drainage (Fig. A1a). This gradual lake drainage constitutes a potential negative feedback on permafrost degradation, since shallower water bodies are more likely to bottom-freeze, allowing for more efficient cooling of the ground during wintertime (Boike et al., 2015; Langer et al., 2016; O'Neill et al., 2020). Previous modelling studies have also demonstrated that the stability and the thermal regime of permafrost in the vicinity of thaw lakes is affected by meso-scale lateral heat fluxes from taliks forming underneath the lakes (Rowland et al., 2011; Langer et al., 2016). These effects have not been considered in this study but constitute another feedback due to lateral fluxes at the meso-scale.

Finally, the most complex pathways of landscape evolution are revealed when both micro- and meso-scale heterogeneities are taken into account. The simulation of a low-gradient polygon slope reflects the transition from low-centred polygonal tundra with inundated centres and troughs towards efficiently drained terrain with a high-centred polygon relief. This transient landscape evolution, which resembles the schematic evolution of polygonal tundra depicted by Liljedahl et al. (2016) well, has not been simulated in a numerical model before, as it involves an interaction between micro-scale (i.e. ice-wedge degradation) and meso-scale (i.e. drainage along slope) processes. Due to the inclusion of a wide range of feedback processes on multiple scales, we consider the simulations conducted with the low-gradient polygon slope setup to most realistically mirror the real-world dynamics of ice-rich permafrost lowlands. Interestingly, the projected thaw-depth increase and ground subsidence showed little sensitivity to the meso-scale slope gradient (Figs. 7d and 8d), suggesting a higher robustness of these projections against parameter variations compared to the more simple setups.

Despite the capacity of reflecting a wide range of landscape evolution pathways, certain forms of landscape change that have been observed in ice-rich permafrost landscapes cannot be reflected using the presented model framework due to a lack of necessary process parametrizations. For instance, real-world thaw lakes are known to drain “catastrophically” instead of gradually upon the incision of a lake basin by a thermo-erosional gully. This can happen as a result of either gully growth or lake expansion (Jones et al., 2011; Kessler et al., 2012), both of which are a consequence of lateral erosion. This process is, however, not represented in our model at the meso-scale. Similarly, the dynamic development of retrogressive thaw slumps and their interactions with other landforms would require additional parametrizations of mass-wasting processes as well as slope- and aspect-dependent modifications of incoming radiation (Lewkowicz, 1986). Beyond this, it should be noted that further potential feedbacks mechanisms are not represented in our model, such as lateral heat and sediment transport at the meso-scale or ecological processes such as vegetation succession, all of which potentially alter the ground thermal and hydrological regimes (Shur and Jorgenson, 2007; Kokelj and Jorgenson, 2013; Kanevskiy et al., 2017) and affect permafrost thaw.

4.2 Implications for simulating high-latitude land surface processes using macro-scale models

The presence of permafrost in high-latitude landscapes crucially affects hydrological and biogeochemical processes and makes their realistic representation within macro-scale LSM frameworks challenging (Chadburn et al., 2017; Andresen et al., 2020; Burke et al., 2020). Beyond the requirement to accurately simulate freeze–thaw processes near the surface, the presence of excess ice in the subsurface is an additional source of complexity. Indeed, its melting involves rapid changes in the topography which affect hydrological, ecological, and biogeochemical processes, and these in turn feed back on the thermal dynamics of the subsurface. Here, our results support previous findings which have suggested that micro- (Abolt et al., 2018, 2020; Nitzbon et al., 2019) and meso-scale (Langer et al., 2016; Jafarov et al., 2018) heterogeneities and lateral fluxes exert important controls on the simulated thermal, hydrological, and biogeochemical state of the subsurface in ice-rich terrain. Beyond the importance of subgrid-scale heterogeneity, our simulations suggest that interactions across scales can introduce further complexities and feedbacks in the land surface evolution in response to climate warming, which are not represented in LSMs. For example, the rapid subsidence of inter-polygon troughs due to melting of ice wedges can entail landscape-wide regime shifts of the soil moisture state and runoff generation. Similarly, the setups which reflect micro- and/or meso-scale heterogeneities enable simulating pathways of landscape evolution in which some parts of the landscape are inundated (e.g. deepening troughs), while other parts are subject to drained conditions. The spatial variability of the ground thermal and hydrological regimes directly affects other ecosystem processes such as the decomposition of soil organic carbon (Lara et al., 2015). In fact, several studies have suggested that small water bodies could constitute “hotspots” for carbon decomposition and greenhouse gas emissions (Laurion et al., 2010; Abnizova et al., 2012; Langer et al., 2015; Elder et al., 2020), which would be missed in projections of macro-scale models used to constrain the global permafrost carbon–climate feedback (Schuur et al., 2015).

The modelling approach of multi-scale landscape tiling pursued in our study constitutes a promising way to represent subgrid-scale heterogeneities of permafrost landscapes without the need to increase the spatial resolution of a macro-scale model. Our modelling suggests that key feedback processes acting at and across the micro- and meso-scales of ice-rich terrain could be represented in a nine-tile setup, corresponding to an increase in computational demand by about an order of magnitude. Despite mathematical and technical challenges that would need to be addressed before multi-scale tiling schemes could be incorporated within established LSMs (Fisher and Koven, 2020), we consider the concept of multi-scale tiling to bear considerable potential for a computationally affordable representation of subgrid heterogeneity and lateral interactions in high-latitude ecosystems within LSM frameworks.

4.3 Towards simulating the cyclic evolution of ice wedges and thaw lakes in thermokarst terrain

Beyond the discussed implications for simulating permafrost degradation in response to the expected climate warming within the course of the twenty-first century, our modelling work reveals novel potentials for investigating the geomorphological evolution of thermokarst terrain under both past and future climatic conditions. It has been suggested that thermokarst terrain is evolving in a cyclic manner on centennial-to-millennial timescales. For example, cycles of degradation and stabilization have been suggested to describe the evolution of ice wedges (Jorgenson et al., 2006, 2015; Kanevskiy et al., 2017) on these timescales. Similarly, but at a larger spatial scale, the cyclic evolution of thaw lakes has been hypothesized (Billings and Peterson, 1980; see Grosse et al., 2013, for a review of the hypothesis). While we do not want to discuss here to which extent ice-rich permafrost landscapes have been evolving or are evolving according to these theories, we would like to mention that the model framework presented in this study is able to capture a wide range of the processes involved in these cycles. Figure 9 provides a schematic of the cyclic evolution of thermokarst terrain and highlights the capacities of our model to simulate different pathways of (cyclic) landscape evolution.

Figure 9Schematic depiction of pathways of ice-rich permafrost landscape evolution as simulated within the presented model framework. The inner cycle reflects the cyclic evolution of ice-wedge polygons as described by Kanevskiy et al. (2017). The outer cycle involving the thaw lake stage reflects the thaw lake cycle as hypothesized by Billings and Peterson (1980). Formation of excess ice is lacking in the model such that the full cycles cannot be simulated. The expressions in brackets correspond to the simulated landscape configurations: a low-centred polygon (LCP), intermediate-centred polygon (ICP), high-centred polygon with inundated troughs (HCPi), high-centred polygon with drained troughs (HCPd), and water body (WB). Figure adapted from Jorgenson et al. (2015) and Nitzbon et al. (2019).

The initial and advanced degradation of ice wedges and the associated transition from low-centred to high-centred polygons was captured by the model setup by Nitzbon et al. (2019). Nitzbon et al. (2020a) complemented this by the stabilization of ice wedges due to sediment accumulation in the troughs and also showed that it is possible to simulate ice-wedge collapse and thaw lake formation in response to climate-warming scenarios. The simulations with representation of landscape heterogeneity and lateral water fluxes at the meso-scale, which are presented in this study, allow for capturing further feedbacks. First, the formation of efficient drainage pathways through degrading ice-wedge networks, as depicted by Liljedahl et al. (2016), was qualitatively captured by the simulation using the low-gradient polygon slope setup (Figs. 2d and 6). While the simulated landscape drainage did not completely interrupt the ice-wedge degradation (as indicated in Fig. 9), it did lead to a significantly slower degradation compared to the simulations without meso-scale heterogeneity (Fig. 7b and d). Second, representing meso-scale lateral water fluxes (Fig. 2c and d) enables our model in principle to simulate the drainage of thaw lakes, which would result in the exposure of unfrozen ground (talik) to the atmosphere. While it is not possible to simulate the catastrophic drainage of lakes as discussed above, subsidence of the terrain surrounding the thaw lake can lead to its gradual drainage on longer timescales than considered here. It is further likely that using a model setup with a more complex meso-scale landscape topology could result in more rapid drainage of water bodies.

As indicated in Fig. 9, several processes are not represented in our model such that it cannot capture the full cycles of ice-rich permafrost evolution. In order to represent the initial stabilization of ice wedges more realistically, further ecological processes like vegetation succession and organic matter formation would need to be represented (Jorgenson et al., 2015; Kanevskiy et al., 2017). However, lateral sediment transport into the troughs captures this initial stabilization stage to a certain extent. The advanced stabilization of ice wedges involves the formation of ice-rich layers consisting of segregation ice above the massive ice wedge. The process of ice-segregation would require a more sophisticated subsurface hydrology scheme, representing the migration of liquid pore water towards the freezing front, as well as vertical displacement of sediment resulting from frost heave (Ballantyne, 2018). Finally – and probably most importantly – the formation of wedge ice is not yet represented in the CryoGrid 3 model, which is necessary to form the initial stage of undegraded ice wedges from various other evolutionary stages (Fig. 9). In contrast to the formation of segregation ice, wedge-ice formation involves both vertical and horizontal processes. It presupposes the formation of frost cracks due to mechanical stress, which in turn is controlled by climatic conditions at the ground surface and at the top of permafrost as well as mechanical properties of the soil (Lachenbruch, 1962; Matsuoka et al., 2018). Despite the lack of parametrizations for these processes, our numerical model developments provide important progress towards an improved understanding of the long-term geomorphological evolution of periglacial landscapes under both past and future climatic conditions.