High potential for loss of permafrost landforms in a changing climate

The presence of ground ice in Arctic soils exerts a major effect on permafrost hydrology and ecology, and factors prominently into geomorphic landform development. As most ground ice has accumulated in near-surface permafrost, it is sensitive to variations in atmospheric conditions. Typical and regionally widespread permafrost landforms such as pingos, ice-wedge polygons, and rock glaciers are closely tied to ground ice. However, under ongoing climate change, suitable environmental spaces for preserving landforms associated with ice-rich permafrost may be rapidly disappearing. We deploy a statistical ensemble approach to model, for the first time, the current and potential future environmental conditions of three typical permafrost landforms, pingos, ice-wedge polygons and rock glaciers across the Northern Hemisphere. We show that by midcentury, the landforms are projected to lose more than one-fifth of their suitable environments under a moderate climate scenario (RCP4.5) and on average around one-third under a very high baseline emission scenario (RCP8.5), even when projected new suitable areas for occurrence are considered. By 2061–2080, on average more than 50% of the recent suitable conditions can be lost (RCP8.5). In the case of pingos and ice-wedge polygons, geographical changes are mainly attributed to alterations in thawing-season precipitation and air temperatures. Rock glaciers show air temperature-induced regional changes in suitable conditions strongly constrained by topography and soil properties. The predicted losses could have important implications for Arctic hydrology, geo- and biodiversity, and to the global climate system through changes in biogeochemical cycles governed by the geomorphology of permafrost landscapes. Moreover, our projections provide insights into the circumpolar distribution of various ground ice types and help inventory permafrost landforms in unmapped regions.


Introduction
The presence of ground ice in the surface layers of permafrost is central to the geomorphic and hydroecological functioning of high-latitude and altitude permafrost landscapes (Rowland et al 2010, Kokelj et al 2014, Liljedahl et al 2016. Involved geomorphic processes give rise to permafrost landforms, such as ice-cored pingo mounds (Mackay 1988), ice wedges that form distinctive polygonal patterned ground (Washburn 1980), and slowly creeping rock glaciers in mountains (Barsch 1996). Ice-rich permafrost landforms are responsive to climatic perturbations as indicated by the intensified permafrost degradation during the past decades (Liljedahl et al 2016, Farquharson et al 2019. Relict permafrost landforms found in temperate climates (e.g. pingo remnants and ice-wedge casts in western and central Europe (Vandenberghe and Pissart 1993) imply that distributions of permafrost landforms have changed before and are likely to change in future climates (Washburn 1980, Grosse and Jones 2011, Liljedahl et al 2016. Projected increases in air temperatures and precipitation with Arctic amplification (Hoegh-Guldberg et al 2018) are anticipated to intensify permafrost warming (Biskaborn et al 2019) and associated thermokarst processes, i.e. thawing and settling of ice-rich permafrost (Jorgenson et al 2006, Grosse et al 2016, Douglas et al 2020. Landscapes with current thermokarst features or susceptibility to future thermokarst are estimated to cover up to 20% of circumpolar permafrost regions (Olefeldt et al 2016). Rising air and ground temperatures pose a potential constraint for pingo persistence across northern Siberia (Grosse and Jones 2011), Alaska (Jones et al 2012), and Canada (Mackay 1988), for example. Degradation of ice-wedge polygons has been observed across the Arctic, including continuous permafrost regions (Jorgenson et al 2006, 2015, Grosse et al 2016, Liljedahl et al 2016, Steedman et al 2017, Fraser et al 2018, Farquharson et al 2019. Moreover, downslope creeping of rock glaciers due to ground ice deformation has shown increased surface velocities related to increased air temperatures in some regions (Kääb et al 2007. The degradation of ground ice in permafrost landforms can induce thermokarst or mass-movement hazards (Kääb et al 2007) and thereby threaten infrastructure development (Rowland et (Lara et al 2018) significance, but also affect global systems through alterations in atmospheric (Schuur et al 2015) and aquatic (Serikova et al 2018) greenhouse gas fluxes from permafrost regions. These landforms also contribute to local biodiversity (Jorgenson et al 2006) by serving as highly specialized ecosystems for flora (Jorgenson et al 2015) and fauna (Fountain et al 2012). In some regions they also provide a freshwater supply for human use (Jones et al 2018).
At present, the information on current and potential future distributions of permafrost landforms at the circumpolar scale is sparse. This is despite an increasing amount of circumpolar to global permafrost extent simulations conducted both with processbased (e.g. Guo and Wang 2016) and empirical/statistical approaches (Chadburn et al 2017, Aalto et al 2018, Obu et al 2019. The lack of spatially detailed knowledge on geomorphic features poses a limitation to the assessments of climate change impacts on permafrost environments. We aim to reduce these knowledge gaps by applying statistical distribution modelling. Statistical modelling methods have been used in geomorphological contexts (Luoto et al 2010, Aalto et al 2014, Rudy et al 2017, Fewster et al 2020. In our previous work (Aalto et al 2018), we have demonstrated that a statistical approach produces similar permafrost extents with process-based models (e.g. Guo and Wang 2016) and other empirical models (Chadburn et al 2017, Obu et al 2019, and have relatively small prediction errors. The advantage of the approach over physically based models which require extensive parameterization of physical soil properties, for example, is that statistical modelling allows for predicting landform occurrence at high spatial resolution across the circumpolar region using readily available geospatial data on climate, soil and topography. Thus, our approach does not consider various processes affecting heat fluxes or ground ice accumulation in permafrost but is suitable for examining how different environmental factors affect the distribution of permafrost landforms across the circumpolar region. Our objective is to model, for the first time, the potential environmental spaces for the occurrence of pingos, ice-wedge polygons and rock glaciers across the Northern Hemisphere. These spaces encompass conditions where climate, topography and soil properties are suitable for landform presence. In the modelling we assumed the landforms to be in equilibrium with climate conditions for the considered period, that is, we did not make assumptions concerning whether the landforms are degrading or still growing. Based on the statistical modelling results, we present the circumpolar coverages of potential environmental spaces for the studied landforms during a recent period  and predict their regional changes for two future periods (2041-2060 and 2061-2080) using representative concentration pathway RCP4.5 (stabilization scenario) and RCP8.5 (rising radiative forcing pathway) climatechange scenarios. Moreover, the results offer insights into the fine-scale distribution of various ground ice types and help inventory permafrost landforms in unmapped regions.

Permafrost landform observations
Observations of pingos, ice-wedge polygons and rock glaciers across the Northern Hemisphere (i.e. presence data, figure 1(a)) made since late 20th century were compiled from available inventories, published studies and geomorphological maps (see section S1 (available online at stacks.iop.org/ERL/15/104065/mmedia) in supplementary material). Ice-cored pingo mounds (figure 1(b)) form in permafrost environments by injection and freezing of pressurized water in the soil, and results in massive-ice formation and upheaving of overlying surface material (Mackay 1988). Two primary pingo types are often distinguished: hydrostatic, predominantly occurring in low-lying continuous permafrost regions, and hydraulic, which are more typical to discontinuous permafrost and depend on water moving under a hydraulic gradient (Mackay 1988, Grosse and Jones 2011, Jones et al 2012. Distinct patterned ground features, ice-wedge polygons (figure 1(c)), occur nearly ubiquitously in unconsolidated deposits across the Northern Hemisphere permafrost domain and wedge ice is the most common form of ground ice in continuous permafrost (Bernard-Grand'Maison and Pollard 2018). The polygon shape results from the cyclical growth of ice wedges in thermal contraction cracks formed in susceptible soils during cold spells in wintertime (Washburn 1980, Kokelj et al 2014. Rock glaciers are bodies of poorly sorted debris and ice that move due to deformation of internal ice (Barsch 1996, Berthling 2011, figure 1(d)) in permafrost environments. They occupy most mountain regions in both continental and maritime climates including high-Arctic conditions (Evans 1993) and high elevations in warmer climates in, e.g. the Middle East (Gorbunov 2013). We included both active (flowing or creeping downward) and inactive (stagnant) rock glaciers that contain ice (Barsch 1996) regardless of whether they were primarily influenced by talus slope dynamics or by glacier dynamics (moraine-derived rock glaciers) (see section S1, supplementary material) but excluded relict landforms. In addition to presence data, the applied modelling required information on environments where the landforms are not present. Thus, we compiled datasets with 1349 absence grid cells for pingos, 494 for ice-wedge polygons and 1549 for rock glaciers (section S1, supplementary material).

Environmental data
In general, the formation of ground ice and associated landforms require permafrost conditions with efficient wintertime freezing and reasonably cool summers (Washburn 1980, Mackay 1988. Moreover, rainfall, topography, and properties of snow, vegetation and soil layers affect ground ice accumulation and preservation potential in the landforms (Barsch 1996, Boeckli et al 2012, Bernard-Grand'Maison and Pollard 2018, Knight et al 2019, Douglas et al 2020. We used geospatial data on physically relevant factors (climate, soil and topography predictors) at 30 arcsecond (~1 km) spatial resolution to examine their effects on landform occurrence.
To account for the assumed seasonal effects of climate we computed freezing and thawing degreedays (FDD and TDD, • C-days) and amounts of solid (Snowfall, mm) and liquid (Rainfall, mm) precipitation. Snowfall was defined as the sum of precipitation for months with average air temperature below 0 • C, and rainfall for those above 0 • C. The variables were computed separately for each considered time period and climate scenario using interpolated climate surfaces from the WorldClim database (Hijmans et al 2005). The database includes monthly air temperature averages and precipitation sums for the baseline period 1950-2000 (note that the majority of the permafrost landform observations were from this period) and for two future periods (2041-2060 and 2061-2080), for which we used RCP4.5 and RCP8. The GMTED2010 digital elevation model (Danielson and Gesch 2011) at 30 arc-second spatial resolution (~1 km) was employed to compute the topography predictors, which were used to represent the amount of incoming solar energy, and the transport of ground material and moisture across the landscape. Potential incident solar radiation (Mccune and Keon 2002) (PISR, MJ cm −2 a −1 ) was computed using slope, aspect and latitude, of which slope angle (Slope, • ) was used also as an independent predictor in the modelling. PISR was omitted from pingo and ice-wedge polygon modelling, because in flat areas it is argued to reflect the latitudinal gradient rather than local topography-mediated insolation conditions (Karjalainen et  Averaged proportions of coarse sediments (Coarse, grain size > 2 mm) and fine sediments (Fine, [sum of clay and silt] < 50 µm) as well as the content of soil organic carbon (SOC, g kg −1 ) within 2 m from the surface were derived from the Soil-Grids database (Hengl et al 2017). These factors were assumed to characterize the varying hydro-thermal conditions of permafrost soils. Finally, we determined water coverage (%) in the 30 arc-second grid cells using a 150 m resolution dataset on water bodies (Defourny 2016) to account for the proportion of thawed ground (lakes and taliks) in each grid cell. There also exist site-specific factors that we could not directly consider with geospatial data at a circumpolar extent. Pingo and ice-wedge polygon formation depends on groundwater flow, and the thickness and water flow permeability soil deposits (Jones et al 2012), while debris availability, lithology and meltwater dynamics are central for rock glaciers (Lilleøren andEtzelmüller 2011, Knight et al 2019).

Statistical modelling
We modelled the suitable environmental spaces for landform occurrence and factors affecting their distribution by relating the compiled observations of landform presence and absence to environmental conditions using multivariate statistical methods. We used four methods and model-averaging techniques (Aalto et al 2017) to reduce the prediction uncertainty related to the choice of a single modelling method (Thuiller et al 2009). Generalized linear modelling (GLM), generalized additive modelling (GAM), generalized boosting modelling (GBM) and random forest (RF), were implemented within BIOMOD2 environment (Thuiller et al 2009) in R (version 3.5.2, see model calibration details in section 2, supplementary material). Model evaluation was performed using 100-fold cross validation that at each run used a random sample of 70% of observations in the datasets to calibrate a model and to validate it with the remaining 30%. Model performance was evaluated at each round with the adjusted coefficient of determination (R 2 ) and two prevalence-independent statistical measures of classification accuracy; area under the receiving operating characteristic curve (AUC) and true skill statistic (TSS, Allouche et al 2006). Ensemble predictions were then formed by using the committee-averaging algorithm (Thuiller et al 2009, section 2, supplementary material). We evaluated the spatial reliability of the predictions by mapping the modelling agreement between the four methods (Luoto et al 2010).
To examine the environmental controls and sensitivities of permafrost landforms we estimated a modelling technique-independent variable importance score for each predictor (Thuiller et al 2009), and created response curves with each of the four methods to examine the shapes of relationships between each landform and predictors (Elith et al 2005) using response.plot2 function in BIOMOD2 environment (Thuiller et al 2009).
Projections for 2041-2060 and 2061-2080 were performed by substituting the climate parameters in the calibrated models with respective climate change scenario-based predictors (Hijmans et al 2005). Other predictors were not modified (Aalto et al 2017(Aalto et al , 2018 but used as they were in the baseline modelling in order to constrain suitable topography, surficial deposits and drainage conditions for landform occurrence in changing climatic conditions. These factors were considered to remain relatively unchanged within the study period, while vegetation properties that affect permafrost-climate responses (Shur and Jorgenson 2007) were assumed to change unpredictably, and thus not accounted for. Finally, areal coverages of projected environmental spaces were clipped and computed inside a recent permafrost zonation (Obu et al 2019). This delineation considers isolated permafrost areas in higher detail than the International Permafrost Association map (Brown et al 2002).

Results
Based on the predicted environmental spaces for current and future time periods, we classify three pathways of suitability for landform occurrence; 1) persisting conditions with suitable topographic and soil conditions for landform occurrence and climate remaining favourable, 2) environmental space loss, in which climatic conditions become unsuitable, and 3) new environmental spaces with potential for occurrences outside current distribution (i.e. suitable topography, soil and newly formed climatic conditions). Finally, we examine the factors affecting landform occurrence and evaluate the reliability of the predictions.
The largest losses of environmental spaces for pingos and ice-wedge polygons are concentrated on northwestern Russia and the middle reaches of major Arctic rivers (figures 2(a) and (b)). In the zone of discontinuous permafrost, suitable conditions for rock glacier occurrence show a distinctive shift towards higher elevations. As a result, areal losses are projected in most areas of current occurrences (figure 2(c)). In continuous permafrost regions, the projected losses are generally small.
New environmental spaces are projected to emerge where air temperatures and precipitation attain suitable levels for occurrence and other factors are favourable. For pingos and ice-wedge polygons, new suitable areas (circumpolar totals of 0.67 and 0.33 × 10 6 km 2 ) predominantly situate across Canada's Arctic Archipelago and in the Taymyr Peninsula (figures 2(a) and (b)). Pingo environments show a poleward shift, while suitable conditions for ice-wedge polygons are predicted in diverse landscapes amidst the recent occurrences. However, suitability for recent pingos, also, is predicted in many locations on Canadas's Arctic islands even though not visible in the small-scale map ( figure 2(a)). The largest areal increases in potentially suitable conditions for rock glacier occurrence (0.15 × 10 6 km 2 ) appear in northern Alaska, high-Arctic Canada, eastern Greenland and Far-Eastern Russia ( figure 2(c)). The general trajectory of the changes is towards colder regions of continuous permafrost. Based on the climate-change scenarios, both air temperature and rainfall show distinct increases at grid cells with observed landform occurrences (figure 3). Snowfall is projected to decrease at rock glacier sites, in contrast to the increase in rainfall (figures 3(c) and (d)).

Climatic controls and sensitivity of permafrost landforms
According to the analyses of variable importance and response shapes, the optimal conditions for pingos and ice-wedge polygons at a circumpolar scale are similar (figures 4(a)-(c), supplementary figure 3) and characterized by flat topography with moderately high soil moisture, low rain-and snowfall (<~300 mm), and TDD ( • C-days) of less than~2000. For ice-wedge polygons, FDD also have a notable contribution indicating that their occurrence requires at least~2500 • C-days. Notably, figure 3 suggests that FDD will drop below 2000 at many pingo and icewedge polygon sites regardless of the climate-change scenario.
Rock glacier occurrence increased on coarsesediment slopes with low wetness index, TDD less than~1200 • C-days, and FDD between~1000 and 4500 • C-days ( figure 4(d)). Increased probability of occurrence towards minimum FDD values in the dataset suggests that rock glaciers exist across a wide range of cold climates although most landforms in the dataset have relatively high FDD values ( figure  3(b)). Each landform shows a clear decline in probability of occurrence when TDD exceeds~1000 • Cdays (figures 4(b)-(d)). The response of rock glaciers to increasing TDD is especially sharp. An even more abrupt threshold is found between rainfall and both pingos and ice-wedge polygons.

Statistical and spatial evaluations
To assess model fit and predictive performance we used a split-sample approach; the models were calibrated using 70% of the data and evaluated against the remaining 30% of observations. Based on 100-fold cross-validation, our approach predicted well the permafrost landform occurrences (figure 5). Ensemble method yielded the highest classification accuracy (AUC and TSS) in all cases and thus superior performance. Rock glacier models had the highest R 2 , AUC and TSS values, yet respective ensemble averages in pingo and ice-wedge polygon models also indicated excellent predictive performance. Moderate between-model variability was observed. GBM and RF best explained the variation in each landform occurrence having R 2 values at or above 0.90 (figure 5).
In key regions of landform distribution, the models largely agreed (figure 6). High inter-model variability was encountered in poorly sampled regions, such as the Tibetan Plateau where GLM and GAM predicted notably larger ice-wedge polygon distributions than GBM and RF, which both yielded more constrained environmental spaces (figure 6(d)). The lower evaluation statistics of GLM and GAM in pingo modelling ( figure 5(a)) may explain the model disagreement in many regions (figure 6(a)). In cases like these, the ensemble method facilitates controlling for prediction uncertainty by allowing for a conservative model agreement-based assessment. Despite superior model performance, ensemble models failed to identify a few known regions for observed occurrence. For example, only GBM and RF predicted pingo occurrence potential in Svalbard and eastern Greenland (figure 6(a)). Evaluations of spatial uncertainty such as presented here help identifying regions and environmental conditions where additional investigations should be made and thereby improve future model performance.

Discussion
The predicted losses of environmental spaces are largely associated with regions of projected nearsurface permafrost thaw in the near future due to warming climate (e.g. Slater and Lawrence 2013, Chadburn et al 2017), but not exclusively, as areas of continuous permafrost in, e.g. northern Siberia may also become unsuitable for pingos and ice-wedge polygons. Based on our analyses, the environmental suitability for these landforms will be foremostly constrained by increasing TDD and rainfall. These . First studies considering thermokarst-inducing processes in a numerical model indeed project substantial surface permafrost degradation and widespread landscape collapse even for the currently still cold North-East Siberian permafrost region under a strongly warming climate scenario (RCP8.5) (Nitzbon et al 2020), which is largely in line with our findings. Alongside increasing air temperatures, increased rainfall has been shown to stimulate permafrost thaw (Douglas et al 2020).
Importantly, predicted losses do not imply complete thaw of permafrost at depth but rather that the permafrost landforms can be lost due to conditions becoming unfavourable for their occurrence. However, ground ice in the landforms can be still protected from thaw by ecosystem properties, such as insulation by thick layers of peat (Shur and Jorgenson 2007) or cooling effect of cold air circulation in blocky rock glacier surfaces (Jones et al 2019).
Local hydro-topographical and ecological properties can also determine whether prolonged degradation leads to drainage or wetting of landscape, and thus the potential release of greenhouse gases (Jorgenson et al 2015, Liljedahl et al 2016, Nitzbon et al 2020. Rainfall was the most important single variable delineating suitable conditions for pingos and ice-wedge polygons at the circumpolar scale, while TWI and FDD foremostly affected rock glacier occurrence potential. Even though these were central factors, the results implied that local topography and soil properties hold notable roles in constraining finer scale environmental suitability for the landforms. Based on the findings, it is suggested that increases in the amount of rainfall to above~200 mm (figures 4(b) and (c)) can severely limit pingo and ice-wedge polygon occurrence potential. In the case of ice-wedge polygons, TDD had nearly as an important of a contribution suggesting that air temperature is the central factor for their occurrence. Pingos also showed a secondary relation to TDD yet rainfall was more dominant in their case.
The effect of rainfall is due arguably to the increased conduction of heat into the ground   introduction of advective heat (Kane et al 2001). Rainfall also has been shown to induce abrupt thaw processes in ice-rich permafrost terrain (Kokelj et al 2015). Owing to the nonlinear response shape, the found effect also applies the other way; suitable conditions can emerge if currently very dry areas begin to receive higher rainfall (>~100 mm, figures 4(b) and (c)) in the future, provided that air temperature and other conditions are suitable. Some of rainfall's contribution could be attributed to its moderately strong correlation with TDD in the modelling data (supplementary figure 4). Notwithstanding, the analysis of variable importance suggested that rainfall had a greater relative effect than TDD, especially for pingos ( figure 4(a)).
The responses of pingos and ice-wedge polygons with rainfall, and rock glaciers with air temperature, were strong and supported by each method. The projections of future climate, however, include uncertainty, especially in mountainous areas (Hijmans et al  2005, Hoegh-Guldberg et al 2018). The interpolated climate surfaces cannot fully account for orographic precipitation patterns or wind-induced snow transport, for example. Ice-wedge polygons' response to snowfall indicated a lowered occurrence probability with increasing snowfall (supplementary figure 3) in line with the observations by Kokelj et al (2014). However, climatic predictors averaged over a series of decades cannot represent future dynamics of snow depth-controlled frost cracking under wintertime cold spells or other extreme air temperature, precipitation or disturbance (i.e. fire) events, which play a role in the formation and degradation of especially ice-wedge polygons (Jorgenson et al 2006, Jones et al 2015, Liljedahl et al 2016, Kanevskiy et al 2017. Based on experimental evidence (Allard andKasper 1998, Christiansen 2005), thermal contraction cracking occurs during prolonged cooling periods with ground temperatures below −20 to −15 • C. In this study, the moderate importance of FDD and snowfall for ice-wedge polygons suggests that higher temporal resolution of climate predictors (monthly to sub-monthly) and a more explicit representation of snow cover variation are needed to estimate thermal cracking potential.
Recent studies show that ice-wedge polygons formed in centennial to millennial time scales can be lost by melting of their uppermost portion in a few decades (Jorgenson et al 2006, Liljedahl et al 2016, Fraser et al 2018, Farquharson et al 2019.
Large pingos with thick insulating overburdens are to a degree resilient to short-term changes in climate (Mackay 1998), but evidence of collapsed and paleo-pingos in concurrent non-permafrost environments confirms their climate sensitivity over longer time scales (Vandenberghe andPissart 1993, Jones et al 2012). Degradation of ice-wedge polygons at local scales demonstrates large spatial variation owing to complex interactions between, e.g. ground ice, wetness conditions and microtopography (Kanevskiy et al 2017, Steedman et al 2017, Fraser et al 2018, Nitzbon et al 2019, 2020. Notwithstanding, some have argued that future thermokarst rates at broader scales mainly depend on the magnitude of regional climate change (Olefeldt et al 2016). Degradation of any of the studied permafrost landforms can also occur due to their natural lifecycles without external forcing (Washburn 1980, Mackay 1988. Few quantitative or spatial estimates exist to compare with the predicted coverages. Mackay (1972) estimated that ice-wedge polygons covered 2.6 × 10 6 km 2 of the Northern Hemisphere, which falls relatively close to our ensemble prediction for 1950-2000 (3.1 × 10 6 km 2 ). On a regional level, our ice-wedge polygon predictions for the main occurrence areas are in line with a recent Canada-wide modelling of wedge-ice content (O'Neill et al 2019), yet regional discrepancies are visible (section S3, supplementary material).
The discrepancies between the models (figures 4 and 5) are argued to have stemmed from different initial model assumptions and parameterizations (Thuiller et al 2009). For example, only GBM and RF captured a decline in rock glacier occurrence probability when FDD approached zero, i.e. milder freezing-season conditions (figure 4(d)). GLM and GAM overlooked this threshold and extended to too warm and flat areas, such as periglacial environments more characteristically affected by extensive solifluction processes, e.g. in northern Scandinavia (Aalto et al 2017). The wide distribution may also have been partly attributed to the inclusion of moraine-derived rock glaciers (section 1, supplementary material) which are easily mixed with debris-covered glaciers, and thus may extend below lower permafrost boundary limits. Moreover, some of the compiled observations might have been relict rock glaciers, and thus possibly lacked permafrost, as suggested by relatively high TDD and FDD for many objects in the dataset (figures 3(a) and (b)). The discrimination of intact rock glaciers from these landforms is often challenging, especially from optical remote sensing data from which a large part of observations was mapped.
Considering the projected new regions with suitable environmental spaces for landforms, we argue that the dependency of landform growth on local geomorphic and hydrologic conditions and glacial history means that they are likely to occur in more constrained environments than estimated in this study.
Prediction of pingo formation potential, for example, is difficult owing to its dependency on special groundwater pressure conditions (governed by underlying substrate properties) and dynamic processes such as lake drainage (Mackay 1998). Here, coarse sediments showed moderate importance and realistic response shapes with the landforms. Due to relatively few observations of soil properties within the permafrost domain available to produce the soil data layers (Hengl et al 2017), the data used here likely have a limited ability in depicting landform-scale variation in grain size proportions. Moreover, although we argue that the~1 km resolution is sufficient to delineate topo-climatic conditions for most rock glaciers, there remains a need to account for site-specific debris availability, lithology and melt-water dynamics, or groundwater conditions and sedimentary history in the cases of pingos and ice-wedge polygons. Not considering all constraining factors for landform development may cause modelling to overestimate the extent of new suitable conditions, and subsequently the presented proportions of lost environmental spaces are arguably conservative. Ultimately, we model suitable conditions for landform occurrences in recent and future climates, but cannot presently address the physical processes, such as ground-ice degradation or aggradation and associated heat fluxes, behind landform development using our approach. The vulnerability of permafrost to atmospheric warming is highly complex and depends on multiple dynamic processes acting in the surface water, snow, vegetation and soil layers and related positive and negative feedbacks Jorgenson 2007, Nitzbon et al 2019). Moreover, the thermal inertia of ice-rich ground encumbers permafrost warming (Grosse et al 2011). Ground ice and associated phase changes of water (i.e. the latent heat effect), as well as the content of unfrozen water in frozen soils (Hu et al 2020), are key constituents of the near-surface permafrost response to climate forcing that may not be uniform in different initial ground temperature conditions (Riseborough 1990). Not considering such transient effects lowers our ability to assess landform occurrence potential in the future, given that the used approach assumes equilibrium state between landforms and climate. Processbased models could be appropriate for assessing the physical responses and related feedbacks of climate change because of their explicit parameterization of heat fluxes and transient effects in the surface layers (e.g. Harris et al 2009). However, they are often conducted at coarse spatial resolutions, which limits their applicability in assessing permafrost landforms at a relevant scale in a hemispheric extent.
Statistical modelling allows for producing reasonably accurate information on the suitable environments for permafrost landforms at a high resolution. This information can be used to link locally resolved detailed land surface processes to globally relevant developments due to climate change. Predictions can also facilitate initial landform mapping in understudied regions although their accuracy at pixel resolution (~1 km 2 ) is to be taken cautiously due to the documented prediction errors (figure 5). In addition to likely profound geomorphic changes, ground ice degradation in the areas of projected lost environmental conditions may greatly alter hydrological and biogeochemical processes. In addition, the eventual disappearance of permafrost landforms represents a loss in local geodiversity, and thereby affects biodiversity (Antonelli et al 2018). From a human activity point of view, ground ice degradation and loss of permafrost landforms can affect freshwater availability (Jones et al 2018 and infrastructure integrity (Raynolds et al 2014, Kanevskiy et al 2017. Using permafrost landforms as proxies for massive ice, our modelling yielded novel information on fine-scale spatial variation of high ground ice content across permafrost landscapes. It should be noted that pingos occupy a very small proportion of landscapes relative to ice-wedge polygonal terrain or even large rock glaciers, and thus the potential amount of ground ice expected in a grid cell is different for each landform.

Conclusions
From this study the following conclusions could be drawn: • The suitable environmental conditions for permafrost landforms at the circumpolar scale are closely related to climate and prone to shift along projected regional changes in air temperature and precipitation. • On average, the landforms are projected to lose more than one-fifth of their suitable environments under a moderate climate scenario (RCP4.5) and around one-third under a very high baseline emission scenario (RCP8.5) by midcentury. • The sensitivity of pingos and ice-wedge polygons to the amount of rainfall was unambiguous, albeit TDD was almost as important for ice-wedge polygons. • Alongside air temperatures, precipitation regime is suggested to be a key constraining factor of geomorphic development outside mountainous environments. • For rock glaciers, air temperature and topography primarily controlled the occurrence potential. However, the approach had a limited ability to separate suitable environments for rock glaciers from those for other glacial or periglacial slope processes.
The predicted changes in environmental spaces suitable for the studied permafrost landforms will affect Arctic landscapes and their geomorphology across large areas. The potential implications should be considered when assessing impacts of climate change on natural and human systems.

Acknowledgments
This study was funded by the Academy of Finland (grants 315519, 307761). BMJ was supported by awards from the US National Science Foundation (OPP-1806213 and ICER-1927872) and GG was supported by ERC PETA-CARB (#338335), the HGF Impulse and Networking Fund (ERC-0013), and NASA NNX08AJ37G. The contribution by BE and KSL was funded by the University of Oslo, Norway. We wish to thank two anonymous reviewers for their highly valuable comments.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.

Author contributions
OK, ML and JH developed the original idea. OK led the compilation of observational data with contributions from BE, GG, BMJ and KSL, and geospatial data processing with JH and JA. OK performed the statistical analyses with contributions from JA, JH and ML. OK wrote the manuscript with contributions from all the authors.