Topographic and Ground-Ice Controls on Shallow Landsliding in Thawing Arctic Permafrost

An increase in Arctic shallow landsliding is a potential consequence of climate warming. Warmer summer-air temperatures and larger rainfall events drive heat into the active layer, melting ice and decreasing soil shear stress. Topography has the potential to exacerbate landsliding by controlling the distribution of ground ice and the movement of water in the subsurface. We demonstrate that shallow Arctic landslides initiate


Introduction
Exceptional weather events are becoming more frequent in the Arctic causing permafrost to thaw rapidly (Dobricic et al., 2020;Lewkowicz & Way, 2019). Extreme summer air-temperatures and more frequent precipitation events drive thaw fronts deeper into ice-rich permafrost soils (Stocker et al., 2013), melting previously thermally stable ice. Melting ice in thawing permafrost increases pore-water pressures within the active layer that destabilize the overlying soil causing retrogressive thaw slumps (Balser et al., 2014;Lewkowicz & Way, 2019), solifluction (Harris et al., 2011), and active-layer detachment (ALD) failures (Lewkowicz, 1990). When these mass-wasting processes occur in proximity to populated areas they damage infrastructure (Hanna et al., 1998), cause economic loss (Allard et al., 2012), and threaten public safety (Larsson, 1982). By 2050, two-thirds of populated areas in permafrost regions will be affected by mass-wasting processes associated with permafrost degradation (Hjort et al., 2018). Long-term studies of permafrost degradation show that mass-wasting processes facilitate and enhance the mobilization of stored carbon stored in permafrost soils (Lamoureux & Lafreniere, 2009). Permafrost contains twice the amount of carbon than is currently in the atmosphere (Hugelius et al., 2014), making landsliding a potential contributor to the growing greenhouse gas emissions from degrading permafrost (Lamoureux & Lafrenière, 2018). A better understanding of the processes controlling hillslope instabilities will help quantify the contribution of mass-wasting processes to Arctic greenhouse gas emissions and evaluate their threat to infrastructure and public safety.
Abstract An increase in Arctic shallow landsliding is a potential consequence of climate warming.
Warmer summer-air temperatures and larger rainfall events drive heat into the active layer, melting ice and decreasing soil shear stress. Topography has the potential to exacerbate landsliding by controlling the distribution of ground ice and the movement of water in the subsurface. We demonstrate that shallow Arctic landslides initiate in zero-order drainage basins consistent with models of shallow landsliding in non-permafrost environments. However, the low average slopes and low concavity of Arctic hillslopes cannot create pore-water pressures high enough to generate landsliding. Instead, two-dimensional slope stability modeling suggests that the vertical distribution of ground-ice distributions controls landslide susceptibility. High ground-ice concentrations close to the potential failure plane act as a stronger control than high average ice volumes or rapid thawing. Our results demonstrate that landslide susceptibility is strongly affected by topographic controls on ground ice and hydrology.
Plain Language Summary With a rapidly warming climate, landslide activity is increasing in the Arctic. Warmer Arctic temperatures thaw the permanently frozen soil (permafrost) that underlies hillslopes. Thawing creates meltwater, decreasing the strength of the soil, and in some cases initiating landsliding. We mapped the location of shallow landslides in Alaska. These landslides demonstrate a clear spatial pattern, forming above the channel network in depressions called zero-order drainage basins. The similarity between the spatial patterns of permafrost landslides and rainfall-triggered landslides suggests that the topography, hydrology, and spatial distribution of ground ice act to enhance landslide potential in zero order drainage basins. Modeling of hydrology and ground ice effects on slope stability, suggests a weak hydrologic control and a strong ground ice control on landslide triggering. Here we argue that topography affects the distribution of ground ice, in particular its concentration at potential landslide failure planes. We conclude that the strong relationship between topography and landsliding in the Arctic could improve estimates of shallow landslide susceptibility.
Assessing the susceptibility of Arctic hillslopes to mass-wasting processes is particularly challenging due to the strong control of ground ice on pore-water pressures in the active layer (Harris et al., 2011). Of the mass-wasting processes, ALDs are particularly susceptible. ALDs are shallow translational landslides that fail along a shear-plane between the base of the active layer and the top of the permafrost (Lewkowicz & Harris, 2005). Prior to ALD failure, annual plug-like solifluction movements of active-layer soil (Harris & Lewkowicz, 2000;Harris et al., 2008;Lewkowicz & Harris, 2005) progressively reduces basal shear strengths to residual values Lewkowicz & Harris, 2005). When prolonged periods of high summer air-temperatures drive thaw fronts beyond the average active-layer depth and into the ice-rich layers above the permafrost melting ice lenses and collapsing pore spaces elevate excess pore-water pressures (Harris & Lewkowicz, 2000;Harris et al., 2011). The excess pore-water saturates the already weakened active-layer soil reducing shear-strengths further, creating a shearing plane and failure (Harris et al., 2001;Harris & Lewkowicz, 2000). The high excess pore-water pressures and residual shear strengths mean that ALDs can happen in low-relief topography, and on slopes as shallow as 3° (Rudy et al., 2016;Swanson, 2012). The importance of topography on ALD generation is typically discussed in the context of aspect controls on active layer thickness. North-facing slopes have a prevalence for higher rates of ALD failures due to thinner active layers that increase the likelihood of ground being ice closer to the surface, making them more susceptible to thaw (Rudy et al., 2016). Additionally, north-facing slopes have a greater amount of snow accumulation that remains on the slope later into the thawing season. This leads to a source of moisture for ground-ice aggradation during autumn freeze-back (Rudy et al., 2016;Young et al., 1997).
In non-permafrost settings, topography controls shallow landsliding through both the relationship between slope and shear stress and the hillslope hydrology, which controls pore-water pressure (Montgomery & Dietrich, 1994). Here, convergent topography drives throughflow that increases pore-water pressures and triggers shallow landslides because the infiltration of intense precipitation (Iverson, 1997). Given the complex interaction between ground-ice aggradation/degradation and topography, we still have a developing understanding of their relationship to landsliding in thawing permafrost (Harris & Lewkowicz, 1993;Mackay, 1981). In this study, we seek to understand the spatial distribution of ALD failures and whether this can provide insight into their triggering mechanisms. We perform a topographic landslide-analysis based on a mapped ALD inventory, comparing this to a model of slope stability under different ground-ice conditions to infer a mechanism for shallow landsliding on permafrost hillslopes.

Topographically Based Model
We developed an ALD inventory for a 100 km 2 region in the central Brooks Range (68.01°N, 155.85°W) (Figure 1a). It is a region of continuous cold-permafrost with low-relief topography that is overlain by till-derived colluvium and covered in Arctic tundra vegetation (Balser, 2015;Jorgenson, 2009;Trochim, Jorgenson et al., 2016). Permafrost extends to a depth of 295 m (GTN-P, 2021) and the active layer varies in thickness between 57 and 200 cm (Balser, 2015;CALM, 2021).Widespread ALD failures occurred in this region after a period of exceptionally high summer-air temperatures during the summer of 2004 (Swanson, 2010(Swanson, , 2014. We mapped 188 ALD failures ( Figure 1a) using Maxar Technologies 2008-derived imagery accessed via Google Earth with a spatial resolution of 50 cm. Each ALD has a minimum mapped unit of 25 m 2 . We defined an ALD as a soil-exposed runout zone with a compressional tongue of soil material at its foot ( Figure 1b).
We compared our mapped ALD polygons to topographic derivatives of a 5-m interferometric synthetic-aperture radar-derived digital elevation model (available at USGS, 2015). To establish the importance of convergent topography to ALD failures we created a series of stream networks with upstream drainage area thresholds ranging from 12,500 to 125,500 m 2 at intervals of 2,500 m 2 . Mapped ALD rasters were compared to these stream networks, which we buffered to account for the width of each ALD. Typically, stream networks with ranges of 25,000-50,000 m 2 mapped most closely to the streams visible on the imagery. Buffer widths varied from 10 to 100 m in increments of 10 m for each stream network and represent the hillslope distance above a given stream network. We assessed the spatial coincidence of ALDs and the stream network by comparing the area of ALD failures and non-failed terrain (Figure 1c).
We compared the distribution of ALD failures with a commonly used topographically derived landslide susceptibility model called SHALSTAB (Dietrich & Montgomery, 1998;Montgomery & Dietrich, 1994). SHALSTAB calculates landslide susceptibility using the infinite-slope assumption for soil-mantled landscapes (where a thin soil overlay bedrock). In the model, the degree of topographic convergence dictates the pore-water pressure generated by throughflow (Montgomery & Dietrich, 1994). The model was applied by assuming that the permafrost represented an impermeable surface (analogous to a bedrock surface in a non-permafrost landscape), allowing us to assess the importance of throughflow in the generation of ALDs. We implemented the model using peak and residual soil shear-strength parameters (Table A1 and Figure S1) derived from geotechnical observations of ALD failures of Harris and Lewkowicz (2000) and hydraulic conductivity from Anderson and Anderson (2010

Slope-Stability Model
We explored hypotheses for how topographic controls on ground ice might promote ALD failures by modeling slope stability. The hypotheses we tested were: (a) higher volumes of ice lenses along the permafrost-active layer boundary drive failure, (b) vertical distribution of ice-lens in the active-layer where ice lenses are concentrated at the base of the active layer are more likely to fail, and (c) excess pore-water dissipation rates caused by advancing thaw fronts strongly control failure potential. To design our hillslope and implement our slope-stability tests we used the geotechnical software GeoStudio (Krahn, 2004). For our stability analysis we applied the 2-D limit equilibrium-based Morgenstern-Price method of slices (Morgenstern & Price, 1965) which outputs a factor of safety (FoS) value to indicate the level of stability. The value derived from a FoS calculation is, in simple terms, the ratio between the maximum resisting forces and the destabilizing forces along a selected shear surface within a hillslope. The condition where the FoS = 1 means that the shearing resistance is fully mobilized and in equilibrium with the destabilizing forces, a FoS > 1 implies a stable slope, whilst a FoS < 1 implies failure. Where the FoS > 1 the slope is considered closer to failure the closer the value is to 1. As there is no field-based data to parameterize the model, we used the dimensions and soil geotechnical properties from well-constrained laboratory tests of ALD failures . By modeling a laboratory slope we can test our hypotheses regarding the role of ice lens distributions on ALD initiation while avoiding uncertainties in the material parameters and failure dimensions of a natural slope.
Our modeled hillslope was 15 × 7 m with a gradient of 24° and a maximum active-layer thickness of 1.5 m (normal to the surface). The active layer consisted of 18,060 cells with each cell being 5 × 2.5 cm. We simulated the passage of thawing from the surface to a depth of 1.5 m at increments normal to the hillslope surface. Between 0 and 1 m, thaw depths were considered at increments of 12.5 cm and below 1 m at increments of 5 cm (Figures 2a and S2). For each analysis, the thaw-front depth defined an interface between an upper layer of thawed soil and a lower layer of impermeable permafrost. We randomized, at any particular depth, the location of cells that represent ice lenses. When thawed, cells with no ice lenses were given a unit weight of 19.62 kN m −3 , cohesion of 3.5 kPa, an internal angle of friction of 31° and geostatic pore-water pressures . Melting ice lenses generate pore-water pressures 20% greater than the hydrostatic Kern-Luetschg & Harris, 2008), so as the thaw front passed each depth increment all of the cells identified as ice lenses were given elevated excess pore-water pressure values. Apart from our third test, excess pore-water pressures remained elevated throughout the test. In order to exclusively understand the role of ice lenses in ALD failures we omitted plug-like flow in our hillslope model by assigning peak values instead of residual values for our basal shear strengths.
To test our first hypothesis, we controlled the volume of ice lenses at a depth of 1.5 m in three experiments. First, we controlled the location of ice-lens cells, considering an even distribution (Figure 3a, blue line). Second, we randomly varied excess pore-water pressures by 0%-20% of the hydrostatic value in each ice lens cell (Figure 3a, orange line). Third, for each test we calculated the mean of five tests with a random distribution of ice-lens width (which is a function of the occurrence of collocated ice-lens cells) (Figure 3a, black line and gray range).
For the second hypothesis, we modeled the depth distribution of ice lenses based on three distributions of ground ice; uniform, linear, and power law (Figures 2b-2d). We varied the total volume of ice within the soil column, first considering scenarios with the same total volume of ice within the soil column (2,230 icelens cells) and second by varying the total volume of ice lenses to maintain a 60% ice volume at a depth of 1.5 m. Excess ice volumes greater than 60% of the soil mass are rare in the field (Allard et al., 1996;Kokelj & Burn, 2003;Morse et al., 2009), so we chose this as a maximum value for all tests. For our uniform distribution we randomly assigned 2,230 ice-lens cells evenly through the active layer to reach an average of 12% by volume (Figure 2b). For our linear distribution we allocated 2,230 ice-lens cells linearly with depth reaching 25% at 1.5 m (Figure 2c). Linear and uniform distribution of ice lenses represent end-member scenarios for ground ice distributions rather than any particular formation from two-sided freezing, these distributions help to put into context the results from our more representative power-law distribution. To account for concentrated ice lenses at the base of the active layer in permafrost regions with two-sided freezing, we fitted a power-law function of the form I vol = 11.841 z3.995 to field data from Kokelj and Burn (2003), (Figure 2d), where z is a negative depth below the ground surface (m) and I vol is the volume of ice relative to the soil mass (%).
For the third hypothesis we considered three tests designed to represent a zero, slow, and fast rate of excess pore-water dissipation (the rate of drainage of water from the thaw front) ( Figure S4). We represented the dissipation of elevated excess pore-water pressure by returning pore-water pressures to hydrostatic values once outside of the excess pore-water pressure layer. In test 1 we assumed that the dissipation of excess pore-water pressure was zero (rapid thaw advance). For test 2 excess pore-water pressures were limited to cells within 12.5 cm of the thaw front (a slower thaw advance). In test 3, the excess pore-water pressure layer thickness was decreased to 2.5 cm (the slowest thaw advance). To contextualize the excess pore-water pressure layer, we used data from Harris et al. (2008, Figure 19b). From this we calculated the average dissipation time of 43.97 h for excess pore-water pressures and combined this with the reported average thaw rate of 1.48 mm hr −1 . We calculate that their experiment had an excess pore-water pressure layer thickness of 6 cm.

Results and Discussion
ALD failures form within convergent topography in the uppermost parts of catchments (drainage areas of less than 75,000 m 2 ; Figures 1a, 1c, and S3). Figure 1c shows a correlation between the location of ALD failures and the channel network relative to hillslopes (represented by the different buffer widths). Here, 40% of the ALD failures in the landscape occur within 15 m of the channel network at drainage areas greater than 50,000 m 2 , this increases to a maximum of 76% at a drainage area of 12,500 m 2 , while only occupying 22% of the total land area (Figures 1c and S3). Increasing the buffer width also increases the area of representative ALD failures but requires a buffer that covers considerably greater land area. A 14 m buffer is representative of the average ALD failure widths (Figure 1c) in our study area (14.2 ± 3.8 m) and within the range of previously mapped ALD failures (10-30 m) in the Arctic (Swanson, 2014). Drainage densities in our study site are low, consistent with other permafrost catchments, where stream networks are ephemeral (Arp et al., 2012). While Arctic landscapes do not show the typical ridge and valley topography of non-permafrost landscapes (Grieve et al., 2018;Hack & Goodlett, 1960), the coincidence of landsliding and drainage area (Figure 1c) suggests that catchment topography influences landsliding.
In non-permafrost landscapes, the relationship between catchment topography and landsliding is usually associated with the hillslope hydrological controls on throughflow (Hack & Goodlett, 1960;Montgomery & Dietrich, 1994). We explored the potential for a similar hydrologic control using a landslide susceptibility map derived from SHALSTAB. The resulting map performed poorly because the slope angles upon which ALDs fail are much lower than shallow landslides in non-permafrost landscapes (Figures 1a and S1). SHALSTAB and similar models are driven by a steady state sub-surface hydrology, where throughflow saturates the soils and generates hydrostatic pore-water pressures. In our permafrost example, pore-water pressures created by the model are not sufficient to initiate failure at lower slope angles despite studies showing throughflow in this region in the thawed active layer (Evans et al., 2020;Kane et al., 2013;Voytek et al., 2016). We conclude, based on this experiment, that a simple hydrological mechanism cannot fully explain the propensity of ALDs to initiate in areas of convergent topography. With studies indicating that regions of convergent topography have higher ground ice content (Balser, 2015; we developed a series of tests to better understand how melting ground ice can trigger ALD failure. To do this we controlled the distribution and rate of ice lens melting in a 2-D slope-stability model.

MITHAN ET AL.
10.1029/2020GL092264 6 of 10 We tested the role of ground ice in generating excess pore-water pressures that can lead to lower FoS values by modifying the distribution and rate of ice-lens melt in a 2-D slope-stability model. Our experiment shows that the distribution of excess pore-water pressures and the size of ice lenses along a failure plane strongly influence the FoS (Figure 3a). We observe that increasing the proportion of ice lenses along the failure plane causes a linear reduction in the FoS, with the positioning of ice lenses on the plane having a small impact on the FoS (Figure 3a). For our test with a constant excess pore-water pressure of 17.6 kPa, the minimum FoS was 1.005, while variable excess pore-water pressures yielded a higher minimum FoS of 1.065. The random distribution of ice-lens width affects the FoS most when the proportions of ice lenses along the failure plane are between 20% and 80%. At low and high proportions, the variability in ice lens size is small, as all ice lenses are either very small (single pixels) or very large. For intermediate proportions of ice lenses, lower factors of safety are associated with more frequent wide ice lenses, suggesting that having wide lenses with larger gaps between them produces lower factors of safety than evenly spaced, narrow ice lenses. This implies that slow and prolonged freezing of this zone (which is liable to produce thicker, more continuous lenses) may precondition slopes to failure by increasing ice content and therefore raising pore-water pressure well above the hydrostatic. Of the three variables tested (width, ice lens size distribution, and magnitude of excess pore-water pressure), the modeled slope is most sensitive to the magnitude of excess pore-water pressure generated by each lens. For example, with an average excess pore-water pressure value of 16.1 kPa, our random scenario generates a minimum FoS of 1.063 (Figure 3a orange line), compared with a FoS of 1.005 for a consistent excess pore-water pressure of 17.6 kPa ( Figure 3a, blue line). ALD failures are therefore sensitive to the amount of ice along the failure plane and the excess pore-water pressures that it generates.
When we account for the depth distribution of ice lenses into our model, we found that the magnitude of excess pore-water pressures generated from melting ice lenses at depth controls the FoS. Linear and uniform ice lens distributions generate the most stable slopes with factors of safety 15%-20% higher than for a power distribution (Figure 3b). The difference in the FoS reflects the proportion of ice at the failure plane; with low ice-volume in the uniform (12%) and linear (25%) distributions relative to the power-law distribution (60%) (Figure 3b). This result is confirmed when we set the ice-lens concentration to 60% at the failure plane for every distribution scenario ( Figure S2). At these high ice-concentrations, the uniform and linear distributions produce lower factors of safety because they have higher ice-concentrations through the active layer. Therefore, it is not the total volume of ice lenses in the active layer that governs stability, instead it is the proportion of ice lenses found along the failure plane. This numerical experiment provides constraint on the potential ground ice conditions that would be necessary to trigger ALDs.
With most observed ALDs initiating shortly after extreme warming events (Lewkowicz & Harris, 2005), we examined if the rate at which a thaw front advances into the active layer will impact slope stability. Surprisingly, the rate made no measurable difference to the FoS in our model runs ( Figure S4), as the peak excess pore-water pressure is more important than the rate of pore-water pressure dissipation. Hence the number and size of ice lenses on that failure plane has a stronger effect on the failure potential. Our results demonstrate a strong spatial control on the initiation of ALDs, that could provide a framework for the development of landslide susceptibility maps in the Arctic. ALDs form within a narrow range of slopes and drainage areas within our study area. These areas tend to be convergent, and a few hundred meters upslope of the channel heads. Hence, it may be possible to model landslide susceptibility as a simple function of those parameters (as demonstrated in Figure 1). A greater understanding of the spatial and topographic controls on ground-ice formation is necessary to improve the physical basis for potential landslide susceptibility modeling. In particular, the variation of ground ice concentrations as a function of drainage area in a catchment may be a key variable in future estimates of landslide susceptibility.
While there is no explicit temporal component to our modeling, the observation that ALD failure depends on the development of large ice lenses at depth, has implications for the temporal distribution of landslides ( Figure 4). If we consider the mass of ice at a potential failure plane as the balance between ice generated at a particular depth by cryosuction and the total amount of thaw within the summer, then landslide susceptibility will vary, even within a steady state climate due to random weather events. We can explore the temporal evolution of landslide potential in a changing climate (Figure 4). At average thaw depths there is a low concentration of ice lenses because in any given year there is an equal likelihood of melting or expanding ice lenses. In a steady state climate, a warmer summer, where the active layer is deeper (possibly two or three standard deviations above the average) than the average, will increase the probability of a landsliding event. When average summer air-temperatures rise due to a warming climate, average thaw fronts deepen, meaning that a smaller deviation from the average is required to generate an ALD failure event. In the tenth year of our conceptual model (Figure 4), extended high summer air-temperatures push the maximum thaw depth beyond the base of the active layer and into a region of concentrated ground ice in the upper permafrost. Our conceptual model can help explain the observed concurrent nature of ALD failure events after periods of higher than average warm air temperatures observed in low-relief Arctic landscapes (Lamoureux & Lafrenière, 2018;Lewkowicz & Harris, 2005;Swanson, 2014). Our observations show that ALD failures have a propensity to initiate in zero-order drainage basins (convergent topography) and that this failure is likely related to a topographic control on ground-ice concentrations. This suggests that landslide susceptibility models for ALD failures can be developed based on simple topographic metrics (slope and drainage area), but that going beyond susceptibility to understand hazard requires an understanding of the temporal (and spatial) development of ground ice.

Conclusion
Our observations show that ALD failures have a propensity to initiate in zero-order drainage basins and that landslide triggering is likely related to a topographic control on ground-ice concentrations in the active layer at the permafrost table. We found that a non-linear distribution of ice lenses in the active layer produces a less stable slope at the failure plain. Conceptually, this has implications for the temporal evolution of shallow landslide potential. As average thaw fronts deepen, only a small deviation from the average thaw depth is required to push thaw fronts into zones of concentrated ice lenses in the upper permafrost. This suggests that landslide susceptibility models for ALD failures can be developed based on simple topographic metrics, but that going beyond susceptibility to understand hazard requires an understanding of the temporal and spatial development of ground ice in Arctic hillslopes.

Data Availability Statement
Datasets for this research are available in this in-text data citation references: Mithan et al. (2020), with a CC BY 4.0 license.