Interactive comment on “ Determining the terrain characteristics related to the surface expression of subsurface water pressurization in permafrost landscapes using susceptibility modelling

Abstract. Warming of the Arctic in recent years has led to changes in the active layer and uppermost permafrost. In particular, thick active layer formation results in more frequent thaw of the ice-rich transient layer. This addition of moisture, as well as infiltration from late season precipitation, results in high pore-water pressures (PWPs) at the base of the active layer and can potentially result in landscape degradation. To predict areas that have the potential for subsurface pressurization, we use susceptibility maps generated using a generalized additive model (GAM). As model response variables, we used active layer detachments (ALDs) and mud ejections (MEs), both formed by high PWP conditions at the Cape Bounty Arctic Watershed Observatory, Melville Island, Canada. As explanatory variables, we used the terrain characteristics elevation, slope, distance to water, topographic position index (TPI), potential incoming solar radiation (PISR), distance to water, normalized difference vegetation index (NDVI; ME model only), geology, and topographic wetness index (TWI). ALDs and MEs were accurately modelled in terms of susceptibility to disturbance across the study area. The susceptibility models demonstrate that ALDs are most probable on hill slopes with gradual to steep slopes and relatively low PISR, whereas MEs are associated with higher elevation areas, lower slope angles, and areas relatively far from water. Based on these results, this method identifies areas that may be sensitive to high PWPs and helps improve our understanding of geomorphic sensitivity to permafrost degradation.


Introduction
Unusually warm conditions during recent years in the Arctic have led to changes in the thermal, hydrological, and geotechnical properties of the seasonal active layer and the uppermost permafrost (Kokelj et al., 2002;ACIA, 2005;IPCC, 2013).Water and ice enrichment at the base of the active layer has been well documented in the literature (Kokelj and Burn, 2005;Shur et al., 2005;Tarnocai, 2009;French and Shur, 2010).As precipitation or meltwater from ground ice percolates downwards through the active layer during the summer months, it is trapped above the impermeable permafrost table.During the fall freeze-back period this water undergoes refreezing, consequently developing an icerich transient layer at the base of the active layer (Hinkel et al., 2001;Kokelj and Burn, 2003;Shur et al., 2005).This transient layer undergoes episodic thaw during exceptionally warm years with thick active layer formation (Shur et al., 2005).This addition of moisture, as well as infiltration from late season precipitation, results in high pore-water pressures (PWP) at the base of the active layer (Zoltai and Woo, 1978;French, 2007;Yamamoto, 2014).PWP is the pressure that the water in the voids of saturated soil is under, and influences the shear strength of the soil (Mitchell, 1960;Morgenstern and Nixon, 1971;McRoberts, 1978;Harris, 1981).These pressures can lead to potentially hazardous forms of permafrost degradation and disturbance, and therefore it is important to understand how pore-water pressurization occurs across the landscape, particularly in relation to variable terrain characteristics.To predict areas at the landscape scale that have the potential for subsurface pressurization, we use susceptibility maps generated using predictive modelling approaches (Rudy et al., 2016a).Susceptibility mapping is based on the assertion that conditions which led to geomorphologic features in the past will also result in that same feature in the future (or present) (Varnes et al., 1984).Thus, areas identified as susceptible will have terrain characteristics similar to those in areas where this feature has already occurred.Recent landslide susceptibility studies in non-permafrost settings have begun to use nonlinear generalized additive models (GAMs; Hastie and Tibshirani, 1990;Brenning, 2008;Park and Chi, 2008;Goetz et al., 2011;Niu et al., 2014;Petschko et al., 2014).This research builds upon recent permafrost disturbance susceptibility modelling that used GAMs (Rudy et al., 2016a, b).For our susceptibility modelling we have monitored two landscape features: slope failures referred to as active layer detachments (ALDs) and expulsions of pressurized slurries referred to here as mud ejections (ME; Holloway et al., 2016) (Fig. 1).While these features are morphologically different, the processes causing their formation are similar, i.e., high PWPs caused by rainfall events and rapid thawing of ice lenses at the base of deep active layers (Washburn, 1956;Shilts, 1978;Zoltai and Woo, 1978;French, 2007;Lewkowicz, 2007).These latter features represent the surface manifestation of ephemeral high PWPs (Holloway et al., 2016) and are particularly important because MEs may act as an indicator of potentially hazardous pressures that may lead to slope failure.Hence, the objectives of this research are to (a) independently identify the modelled susceptibility regions for both ALD and ME features and (b) compare the landscape position of these features to predict areas susceptible to high PWP conditions and permafrost degradation landforms.

Study site
This study was undertaken at the Cape Bounty Arctic Watershed Observatory (CBAWO), on Melville Island, Nunavut, Canada (74 • 43 N, 109 • 35 W; Fig. 2).CBAWO is a multidisciplinary research station where monitoring of terrestrial and aquatic ecosystems has been ongoing since 2003.The landscape consists of rolling hills and broad valleys with relief generally < 100 m.The study site is underlain by thick continuous permafrost with a seasonally thawed active layer that can reach 0.7-0.9m by late July (Lewis et al., 2012;Rudy et al., 2013).CBAWO is climatically a polar semi-desert, with a mean annual air temperature (based on the nearest long-term weather station at Mould Bay, NWT, 300 km west) of −17.5 • C (1971 to 2000)  (Environment Canada, 2014).Precipitation primarily occurs as snow, which is extensively redistributed by wind and preferentially deposited on leeward slopes and in low-lying areas.Rainfall is infrequent but high-magnitude events do occur (Favaro and Lamoureux, 2014).Vegetation is composed of graminoid, prostrate dwarf shrub, and forb tundra.Vegetation cover is  (Rudy et al., 2013) and MEs (Holloway et al., 2016) at CBAWO.The contour interval is 10 m and was derived from a 1 m digital elevation model (DEM) obtained from stereo GEOEYE imagery acquired in (Collingwood, 2014).The inset map indicates the site location in the Canadian High Arctic.
www.the-cryosphere.net/11/1403/2017/The Cryosphere, 11, 1403-1415, heterogeneous and varies across the landscape, reflecting soil moisture and drainage conditions across a mesotopographic gradient (Atkinson and Treitz, 2012).The site is underlain by Devonian sandstone and siltstone bedrock comprising the Weatherall, Hecla Bay, Beverley Inlet, and Parry Island (Burnett Point Member) formations (Harrison, 1995), but outcrops are uncommon.The dominant surficial materials are late Quaternary glacial and marine sediments of unknown thickness and felsenmeer (Hodgson et al., 1984).
At CBAWO, ALDs and MEs occurred extensively in late July 2007 when exceptionally high summer temperatures and heavy rainfall resulted in deep active layer thaw (∼ 1 m) (Lamoureux and Lafrenière, 2009).Active MEs were also observed in 2011 and 2012 and corresponded to years with the highest mean July temperatures since 2003 at CBAWO (Holloway et al., 2016) and since 1948 when measurements began at Mould Bay, NWT (Environment Canada, 2014).ALDs and MEs occur when there is rapid thaw at the base of the active layer resulting in high PWPs and occur in similar ice-rich soil materials (Harris and Lewkowicz, 1993;Leibman, 1995;Lewkowicz and Harris, 2005;Lamoureux and Lafrenière, 2009).In the case of ALDs, these pressures lead to shear failure and downslope sliding of the active layer over the failure surface (Harris and Lewkowicz, 1993).ALDs are generally shallow, with a steep headwall or scarp and an un-vegetated slump scar often being ∼ 1 m deep (Fig. 1) (French, 2007).MEs form when high PWPs eject sediment slurries upward through pre-existing cracks or soil structures.They can occur as active (presently ejecting sediment) or inactive (dry and dormant) stratiform mounds on level terrain, and they naturally elongate downslope when occurring on slopes (Holloway et al., 2016).Previous research completed at CBAWO has found that these features appear to occur in distinct landscape settings; i.e., ALDs are commonly found on vegetated slopes, whereas MEs occur on high-elevation, flat, less-vegetated terrain (Holloway et al., 2016).

Data sources and processing
High-resolution (0.5 m) stereo panchromatic WorldView-2 data were collected on 15 July 2012.Using these data, a highresolution (1 m vertical resolution) digital elevation model (DEM) was derived using PCI Geomatica 10.3.2 (Collingwood, 2014).

Active layer detachments
An inventory of ALDs produced by Rudy et al. (2016b) was used in this study.The inventory of 131 ALDs was created by field mapping and through visual inspection of the WorldView-2 imagery (Fig. 2).To evaluate stable landscapes, an equal number (i.e., 131) of undisturbed points were randomly selected in ArcGIS with constraints defined for distance from a water source (> 10 m) and distance to an ALD (> 20 m).The location of mapped ALDs is available through the Polar Data Catalogue (Lamoureux, 2013).

Mud ejections
The dataset used in this study for ME locations was produced by Holloway et al. (2016).Locations were determined by field mapping in June-July 2012 and July 2013 and include a total of 228 MEs (Fig. 2).The data for the location of mapped MEs are available through the Polar Data Catalogue (Holloway, 2017).To remove large spatial clusters of MEs, a 10 m buffer zone was created around each mapped feature in ArcGIS 10.1 and areas where buffer zones intersected were treated as one large polygon with a single point representing a ME (location randomly determined) for every 10 clustered points (i.e., a cluster of 25 MEs would result in 3 points).A total of 6 clusters were removed, leaving a total of 78 MEs.To evaluate control areas where there were no MEs present, 78 sites were randomly generated in ArcGIS 10.1 to correspond to the number of MEs used in the model.

Model explanatory variables
Variables tested in the models were elevation, slope, distance to water, topographic position index (TPI), potential incoming solar radiation (PISR), distance to water, normalized difference vegetation index (NDVI; used only for the ME model), geology, and topographic wetness index (TWI).These variables were chosen as they all have the potential to contribute to areas having high PWPs.Elevation (m), slope angle ( • ), distance to water (m), TPI, and PISR were derived in ArcGIS 10.1 using the DEM.Elevation is used as a proxy for marine limit in the area, with more frostsusceptible soils and ground ice content being below marine limit (∼ 60 m a.s.l) (Barnett et al., 1977).Slope angle is considered an important factor in drainage and in gravitational movements like ALDs, which can occur on lowgradient slopes (Niu et al., 2005;van Westen et al., 2008).Similarly, distance to water is an indication of drainage and wetness of the landscape.TPI compares the elevation of each cell in a DEM to the mean elevation of a specified neighbourhood (50 m radius for this study) around that cell, which was used to evaluate drainage conditions for a location (Jenness, 2006;Guisan et al., 1999).PISR represents differences in intensity of solar radiation, which can control local temperature and evaporation and therefore soil moisture (van Westen et al., 2008).
An NDVI was used as a proxy for vegetation cover, derived from the multispectral WorldView-2 image acquired on 15 July 2012 (see Tucker, 1979).NDVI is a dimensionless radiometric measure that ranges from −1 (non-vegetated surfaces) to +1 (healthy, productive vegetation).NDVI was used only for the ME model, as ME location at the site was shown to be linked to vegetation cover, and vegetation cover determines patterns in soil moisture and ground ice conditions (Holloway et al., 2016).It was not used in the ALD model because original NDVI conditions are changed by ALDs as vegetation is removed (Rudy et al., 2013).
A geological map of the study area was digitized and georeferenced to the 15 July 2012 WorldView-2 image in Ar-cGIS 10.1 (Harrison, 1995).Additionally, a TWI (Beven and Kirkby, 1979) was calculated using Whitebox Geospatial Analysis Tools (Lindsay, 2014).TWI, a proxy for soil moisture (Beven and Kirkby, 1979), is used to quantify the factors controlling hydrological processes for a given area using elevation, slope, and the upstream area contributing to any given cell.While ground ice content is linked to high PWPs, it is not used as an input variable as ground ice maps were unavailable and impractical to attain.
To test for multicollinearity amongst the variables we used variance inflation factors (VIFs) and the Spearman's rank correlation coefficient (ρSp).VIFs estimate how much the variance of the regression coefficient is "inflated" because of linear dependence between the variables.A VIF can be calculated for each variable by deriving a linear regression of that variable on all the other variables.ρSp measures the statistical dependence between two ranked variables.

Generalized additive model
Modelling was performed as a case-control study with points for either ALDs or MEs as cases and randomly selected undisturbed points as controls.An equal number of undisturbed samples were collected to match the disturbed samples (ALDs or MEs) and resulted in an ALD dataset of 262 samples and a ME dataset of 456 samples (Table 1).This dataset was then separated into calibration and validation subsets, with 70 % disturbed and undisturbed samples for calibration and 30 % disturbed and undisturbed samples for validation.This resulted in 184 points for calibration of ALDs (92 each disturbed and undisturbed), 78 for validation (39 each disturbed and undisturbed) of ALDs, 320 (160 each disturbed and undisturbed) for calibration of MEs, and 136 (68 each disturbed and undisturbed) for validation of MEs.
To model the relationship between the response variable (either ALD or ME) and the terrain variables we used a GAM.GAMs are semi-parametric extensions of generalized linear models (GLMs) and provide the flexibility to represent the response's dependence on the predictor variable as either linear or nonlinear (Hastie and Tibshirani, 1990).This type of model is advantageous as nonlinear effects are known to exist in many geomorphologic studies (Goetz et al., 2011;Rudy et al., 2016b).To account for nonlinear variables, GAMs transform nonlinear predictor variables with a smoothing function.The GAMs in this study were fitted using a spline smoother with 4 degrees of freedom, allowing for the de-tection of complex nonlinear responses.The use of GAMs in susceptibility modelling has shown strong predictive performance and has been used in susceptibility modelling with positive results (Brenning, 2008;Jia et al., 2008;Goetz et al., 2011;Niu et al., 2014;Rudy et al., 2016a).The models were developed using R (version 2.15.3,R Development Core Team, 2013; see Rudy et al., 2016a, for details).Model selection was performed using the "dredge" function of the R package "MuMIn" in which the GAM was fitted through iterative evaluations of modified combinations of the terms in the global model (Barton, 2011).For both ALDs and MEs, model selection was based on two parameters, Akaike information criterion (AIC) and explained deviance.AIC measures the quality of each model in a set based on goodness of fit while penalizing for model complexity (Akaike, 1974).Explained deviance was calculated following Eq.( 1): where "null deviance" is the deviance of the model with only the intercept and "residual deviance" is the deviance that remains unexplained after all variables have been included.Output models were ranked by AIC and all models with AIC ≤ 10 were examined (Burnham and Anderson, 2004).This approach allows us to evaluate a wide range of possible models to ensure that each variable is informative and results in a model with the greatest explained deviance and lowest AIC (Hand, 1997;Hosmer and Lemeshow, 2000;Petschko et al., 2014).

Performance assessment
Two methods were used to assess model performance -the area under the receiver operating characteristic (AUROC) curve and a confusion matrix.The receiver operating characteristic (ROC) curve plots all possible combinations of sensitivities (i.e., percentage of correctly classified disturbance points) against the corresponding specificities (i.e., percentage of correctly classified undisturbed points) that can be achieved with a given classifier and is independent of the spatial density of disturbance (Goetz et al., 2011).Overall model performance is then determined by calculating the AUROC curve where the curve ranges from 0 to 1.A model that has an AUROC of 0.5 or less does not predict the occurrence of disturbance any better than chance, whereas a model with an AUROC of 1 represents a model with perfect prediction of the two classes.The quantitative-qualitative relationship between AUROC and prediction accuracy can be classified as follows: 0.9-1 is excellent; 0.8-0.9 is very good; 0.7-0.8 is good; 0.6-0.7 is average; and 0.5-0.6 is poor (Yesilnacar, 2005).To assess the performance of a presence-absence classifier the agreement between predictions and actual observations can be examined using a confusion matrix.For a disturbance to be predicted as present or absent, the predicted probability will be higher or lower than a pre-assigned prob- ability threshold.For this study, a threshold of 0.50 was selected to maximize the sensitivity-to-specificity ratio.

Permafrost disturbance susceptibility maps
For each GAM model, the results were interpolated (and extrapolated) with R packages "raster" and "rgdal" to produce a probability map that was classified into a permafrost disturbance susceptibility map.For this study, we classified our map into susceptibility zones, using the 50th, 75th, 90th, and 95th percentiles, representing very low (< 50), low (50-75), moderate (75-90), high (90-95), and very high (> 95) susceptibility to future disturbance.

Model fit and predictive power
ALDs and MEs were accurately modelled in terms of susceptibility to future disturbance across the study area (Table 1).Terrain variables included in the final ALD model were slope, elevation, PISR, TPI, distance to water, and TWI, resulting in an explained deviance of 45 %.In the ME model distance to water, NDVI, elevation, PISR, TPI, and TWI were included in the highest performing model with an explained deviance of 57 %.The ALD model calibrated with a sensitivity and specificity of 84 and 81 %, respectively, while the ME model had a sensitivity and specificity of 91 and 88 % (Table 1).These predictive metrics indicate that both GAM models are consistently identifying both disturbed ALDs and MEs and undisturbed points.AUROC values were consistently high for both: 0.91 for ALDs and 0.95 for MEs (Table 1).

Importance of predictor variables
AUROC was used to assess the discriminatory power of individual variables outside a statistical model using single variable models (Table 2).The strongest predictors for ALDs were slope and PISR (AUROC > 70 %), followed by elevation TPI and distance to water (AUROC > 60 %).Although TWI was weakly related to ALD occurrence (AU-ROC < 60 %), its interactions with the other model variables led to an increase in the performance of the full model.
For MEs, all variables with the exception of TWI had AU-ROCs > 70.Differences in values between undisturbed and disturbed ALD points were checked for significance using a Wilcoxon rank sum test (for continuous variables).In the ALD model, all continuous variables with the exception of TWI and distance to water are statistically significant at the 99.9 % level, whereas in the ME model all continuous variables with the exception of TWI and PISR are statistically significant at the 99.9 % level.
VIFs and ρSp were used to examine correlations between predictor variables in both models.For the ALD model, slope and PISR had the strongest correlation with ρSp equal to −0.52, while all other variables had weaker correlations with |ρSp| < 0.5.All variables in the ALD model had VIFs below two.For the ME model, elevation and distance to water had the strongest correlation with ρSp equal to 0.54, while all other variables had weaker correlations with |ρSp| < 0.5.All variables in the ME model had VIFs below two.
Bivariate plots were constructed to view the relationship between the probabilities of an ALD or a ME occurring and the terrain variables that had the largest influence on the models (Fig. 3).When the slope is low the probability of an ALD occurring is also low, but as the slope increases to ∼ 12 • the probability peaks and then starts to decline slightly as the slope increases.The opposite is the case for MEs, where the probability is highest at low slopes and decreases with steeper slopes.For PISR, the probability of ALDs is highest when PISR of 1100 MJ m −2 and decreases in a logistic pattern as PISR increases.For MEs, the probability is low when PISR ∼ 1150-1200 MJ m −2 and then peaks at 1250 MJ m −2 .For distance to water, the probability of ALDs is highest when the distance to water is small and decreases as the distance to water becomes greater (although slightly increasing again at 500 m).For MEs, the probability is lower with smaller distances to water and peaks at 600 m, after which it steadily decreases as the distances get higher.For TPI, the probability of ALDs is highest with a TPI of −2.5, while for MEs it is highest at a TPI of 1, indicating the ALDs occur in landscape concavities and MEs occurring in areas that are more convex to flat.For elevation, both ALDs and MEs have low probability at low elevations, albeit the prob-  (m) 308(192)/294(202) 287(239)/413(276) 310(258)/627(202) 339(277)/652(209) 346( 283)/650(279) 66/83 NDVI -/0.18(0.17)-/0.04(0.12)-/-0.02(0.15)-/−0.05(0.17ability of ALDs is higher than that of MEs.The probability of ALDs peaks at about 50 m, where it decreases as elevation increases.For MEs, the probability is highest at 80 m and decreases as the elevation increases.When considering TWI, the probability of ALDs is high when TWI is low and decreases as TWI increases.The probability of MEs is low when TWI is low, peaks between 5 and 10, and then declines again.For NDVI, the probability of MEs peaks at −0.1 and declines as NDVI increases.

Permafrost disturbance susceptibility maps
The spatial predictions of the GAM for ALDs and MEs (Fig. 4) indicate that ALDs have a high probability of occurring on moderate to steep slopes with relatively low PISR, whereas MEs are more likely to occur at high elevations with low slope angles and far from water sources (i.e., drier).For ALDs, the areas of high susceptibility on the landscape are found along river channels and on steeper upland slopes.
Areas of low susceptibility are found on plateaus, far from water sources with relatively higher PISR.Although highsusceptibility zones account for a small portion of the landscape, there is a greater density of disturbance at these loca- tions (Table 3).By contrast, the areas of high ME susceptibility are found on dry, barren uplands and plateaus.

Spatial extent of susceptibility zones
The spatial extent of moderate-to high-susceptibility zones was compared between the two models to identify key terrain attributes responsible for high PWPs (Fig. 5, Table 4).Very little overlap exists between the models (< 1 % of the study area, Table 3).Terrain characteristics are similar between the features in low susceptibility zones but differ in the modelled high-susceptibility zones (Tables 2, 4).ALDs are commonly found on sloped terrain while MEs tend to occur on plateaus.Slope is the main variable driving ALD initiation, while distance to water is the most important variable explaining ME formation.

Landscape distribution and terrain controls over features formed by pore-water pressurization
Landscapes composed of fine-grained surficial sediments are susceptible to a wide range of permafrost degradation processes, including the development of high PWP in the active layer (Harris and Lewkowicz, 1993;Lewkowicz and Harris, 2005).These PWPs are associated with ALDs (Lewkowicz and Harris, 2005) where soil liquid limits are exceeded.Similarly, MEs have been documented in settings in proximity to ALDs and the occurrence of MEs is associated with deep active layer thaw and potentially with summer rainfall (Edlund, 1989;Lewkowicz, 2007;Holloway et al., 2016).While soil PWP measurements are not available to confirm pressurization in these instances, the inferred mechanism is diapirization of sediment slurries from the base of the active layer caused by pore-water pressurization due to ice thaw (Lewkowicz, 2007).Similar climatic conditions and active layer processes are attributed to both ALD and ME formation www.the-cryosphere.net/11/1403/2017/The Cryosphere, 11, 1403-1415, 2017 (Holloway et al., 2016;Lamoureux and Lafrenière, 2009).However, when reaching some threshold of PWP, the landscape response varies depending on localized terrain characteristics, and high PWPs are expressed at the surface as ALDs or MEs.This analysis demonstrates that the distribution of ALDs and MEs that have developed since 2007 at CBAWO is largely distinct in terms of their spatial occurrence.The terrain associated with MEs include high-elevation sites, frequently on plateau or interfluve locations with low slope angles (Table 2).In contrast, ALDs are found at downslope landscape positions on low to moderate slope angles, often in areas of convergent slope drainage and topographic concavity, a pattern observed elsewhere (Lewkowicz and Harris, 2005;Rudy et al., 2016a).Hence, while surficial materials are broadly similar across CBAWO, the landscape zonation of these two features appears to follow a slope continuum.
The precise relationships between pressurization events in warm years and feature development remain speculative but provide explanation for the distinct distribution of ALDs and MEs.In 2007, the warmest year since regional records began in 1949, deep active layer development and late July rainfall triggered widespread ALD formation (Lamoureux and Lafrenière, 2009).At the time of ALD initiation, adjacent slopes were characterized by numerous MEs in soil that was minimally disturbed (Fig. 6), indicating both the presence of pressurized fluids in the disturbance and the formation of MEs in areas with limited or no slope failure.Similar conditions were observed with MEs associated with terminated active layer fractures in 2012, further suggesting the presence of fluid slurries in situations approaching those that generate ALDs (Figs. 1 and 6).These observations suggest that MEs, while clearly reflecting evidence for subsurface soil water pressurization (i.e., Lewkowicz, 2007), also likely play a stabilization role through pressure release to the surface (Urciuoli and Pirone, 2013).Holloway et al. (2016) noted that MEs were primarily located in areas where surface soil was relatively dry and inferred these soil conditions may contribute to soil structure conducive to slurry diapir formation and PWP release.
Collectively, these observations strongly suggest that MEs result from soil water depressurization where active layer failure and movement is not possible or necessary.By contrast, ALDs are associated with sufficient pressurization to induce slope fracturing and downslope movement.

Susceptibility modelling and landscape implications
Modelling results provide a means to evaluate spatial patterns of features formed by high PWPs across the entire landscape.The independently generated models show strong mutual exclusivity of the locations where ALD and ME susceptibility is high, representing ∼ 1 % of the study area (Figs. 4 and 5).
Although at the landscape scale this overlap appears min- imal, these areas are of particular interest since this is indicative of where we expect a transition in soil stability and hence, areas of key potential insight into fine-scale landscape controls over disturbance.Figure 6 shows a photo of one of these areas where both feature types are present on the westcentral part of the study area (Fig. 2).These locations will be monitored in the future to determine PWP compared to areas where models suggest low susceptibility to ALD and ME and thus lower PWP.Considerably more of the landscape is highly susceptible (in the high and very high zones) to disturbance in terms of ALDs rather than MEs (Table 3).Terrain does not vary as much between MEs and ALDs in the low-/moderatesusceptibility zones, whereas it varies substantially in the very high-susceptibility zone (Table 2).This landscape zonation of ME and ALD activity is consistent with field observations at CBAWO and constitutes the first recognition of a broader landscape pattern of soil pore-water pressurization landforms.While results indicate that many areas of the CBAWO landscape appear to be less susceptible to excess pressurization processes as expressed by MEs and ALDs, we can interpret these results to reflect controls over how pressurization affects landscapes and the locations and terrain types most susceptible.Further, the mutual exclusivity of these modelled high-susceptibility areas over space suggests differential responses to pressurization across different terrain factors.
The distinct zonation of ME in upland areas with low slopes indicates that soil water pressurization results in artesian fluid release at the surface in years exhibiting warm conditions and deep active layer thaw (Holloway et al., 2016).MEs also occur coincident with conditions that result in ALD formation (Fig. 6) and in cases where ALD formation was initiated but ultimately terminated.MEs appear to be caused by fluid release to the surface that reduced subsurface pressures, essentially acting as a fluid pressure release mechanism.While this process is observed in other geotechnical settings where drainage stabilization is used for landslide mitigation (Urciuoli and Pirone, 2013), it is not clear whether the presence of a ME will stabilize the soil surface and prevent soil fracturing or any slope movement.Given that ME zones occur where ALDs are not generally observed, we interpret the ME susceptibility zone to represent areas conducive to pore-water pressurization but where slope disturbance is inhibited, primarily by low slope angles.However, at this time it is unclear the extent to which MEs represent a soil stabilization mechanism.
The combined extent of modelled high-susceptibility areas provides a key set of terrain conditions that appear to be associated with soil water pressurization during warm summers at CBAWO.ALDs appear in many permafrost landscapes (French, 2007;Lewkowicz, 2007;Lewkowicz and Harris, 2005;McRoberts and Morgenstern, 1974;Rutter et al., 1973;Morgenstern and Nixon, 1971), while MEs have had less attention.The literature and our observations indicate that MEs occur on Ellesmere (Lewkowicz, 2007), Banks, and Melville islands in the Canadian High Arctic and further south in Nunavut and the Northwest Territories in very different vegetation and terrain conditions (Shilts, 1978).Despite these observations, the broader distribution of MEs remains poorly understood, and the terrain conditions that are associated with their occurrence at CBAWO may not apply to other settings.Hence, while MEs are indicative of subsurface fluid pressurization in permafrost settings, the function they play in dissipating potential slope instability requires further investigation, particularly with field observations to determine the localized processes of water accumulation within the subsurface.

Conclusion
Independent GAM models accurately captured the terrain controls that appear to affect the distribution of ALDs and MEs on the landscape.The susceptibility models demonstrate that ALDs are most probable on hillslopes with gradual to steep slopes and relatively low PISR, whereas MEs are associated with higher elevation areas, low slope angles, and areas relatively far from water (drier).Both features are thought to be formed when high soil PWPs are generated and dissipated, resulting in slope failure or liquefied sediment ejecting to the soil surface.This analysis reveals that these features are found in proximity, but in largely discrete areas at CBAWO.Hence, the joint zones of susceptibility appear to be sensitive to development of high PWPs and point to a possible larger-scale landscape pattern to active layer response to climate warming.While further research and field measurements is necessary to elucidate these patterns, these results provide incentive and potential to move towards a more systematic interpretation of the distribution of high soil PWP and related landforms to improve our understanding of geomorphic sensitivity to permafrost degradation and related geohazards in Arctic landscapes.

Figure 1 .
Figure 1.(a) Active layer detachment at Cape Bounty Arctic Watershed Observatory (CBAWO), Melville Island, NU, on 28 July 2007.Clay slurry is evident along scar track immediately post disturbance.(b) Active mud ejection occurring on a plateau on 13 July 2012.(c) Clay slurry pooling in a crack at the headwall of a recently initiated active layer detachment on 16 July 2012.(d) Field of inactive MEs on a plateau on 18 June 2012.

Figure 2 .
Figure2.Map of the study area showing locations of ALDs(Rudy et al., 2013) and MEs(Holloway et al., 2016) at CBAWO.The contour interval is 10 m and was derived from a 1 m digital elevation model (DEM) obtained from stereo GEOEYE imagery acquired in(Collingwood, 2014).The inset map indicates the site location in the Canadian High Arctic.

Figure 3 .
Figure 3. Bivariate plots of the probability of the presence of ALDs (black line) and MEs (grey line) in relations to terrain variables.The grey shading illustrates the 95 % confidence bands for each plot.

Figure 4 .
Figure 4. ALD (a) and ME (b) susceptibility maps for CBAWO showing areas of very low to very high susceptibility to disturbance.Marine limit is denoted by the 60 m contour line, and the random points for undisturbed locations are indicated.

Figure 5 .
Figure 5. Areas where the ALD and ME maps at CBAWO have overlapping susceptibilities (superimposed over the 1 m shaded relief DEM).The marine limit is denoted by the 60 m contour line.

Figure 6 .
Figure 6.MEs adjacent to an ALD at CBAWO on 28 July 2007.The photo was taken at the edge of the ALD looking out towards the MEs and the adjacent terrain.

Table 1 .
Model performance metrics for training (70 % samples) and testing (30 % samples).Sensitivity and specificity refer to the ability of the respective model to accurately identify the presence and absence of disturbance.AUROC is a measure of overall model fit without the addition of a user-specified threshold.

Table 2 .
Mean and standard deviation (in brackets)values for terrain variables in each susceptibility zone for ALDs and MEs and performance metrics for ALD and ME single variable models.

Table 3 .
Area and disturbance density of susceptibility zones.

Table 4 .
Terrain values for areas of overlap between ALD and ME models.See units for each terrain variable in Table2.