Massive Ice Control on Permafrost Coast Erosion and Sensitivity

High overall rates of permafrost cliff retreat, coupled with spatial variability, have been accompanied by increased uncertainty over future landscape dynamics. We map long‐term (>80 years) retreat of the shoreline and photogrammetrically analyze historic aerial imagery to quantify the processes at a permafrost coast site with massive ground ice. Retreat rates have been relatively constant, but topographic changes show that subsidence is a potentially critical but often ignored component of coastal sensitivity, exceeding landward recession by over three times during the last 24 years. We calibrate novel passive seismic surveys along clear and variable exposures of massive ground ice and then spatially map key subsurface layers. Combining decadal patterns of volumetric change with new ground ice variation maps enables past trends to be interpreted, future volumetric geomorphic behavior to be better constrained, and improves the assessment of permafrost coast sensitivity and the release of carbon‐bearing material.


Introduction
Arctic permafrost is found on 34% of Earth's coasts (Lantuit et al., 2012) and is estimated to release 14 teragrams of particulate organic carbon to the nearshore zone each year, equivalent to the annual Arctic riverine inputs or vertical net methane emissions from terrestrial permafrost . Pan-Arctic coasts on average recede by~0.5 m a −1 (Lantuit et al., 2012), but coasts consisting of ice-rich permafrost sediments have been shown to be retreating much faster in recent decades (Radosavljevic et al., 2016), exceeding 20 m a −1 in places (Solomon, 2005). Constraining erosion rates, responses, and trends is essential to mitigating the threats to Arctic communities and infrastructure (Radosavljevic et al., 2016), and establishing the resultant material fluxes from coastal erosion impacts on global carbon release and climate feedbacks (Couture et al., 2018).
Ice-rich permafrost terrain comprises a surficial "active layer," a shallow (generally <1 m) layer that thaws seasonally (Lyon et al., 2009); perennially frozen "permafrost" that can reach hundreds of meters thick (Lantuit et al., 2012) and areas of segregated ice known as "massive ice" (Mackay, 1972;Murton, 2009). Massive ice can have multiple forms and origins (Gilbert et al., 2016;Mackay, 1980) but can generally be defined as large bodies of ground ice where gravimetric moisture contents exceed 250% by weight (Pollard, 1990). This epigenetic volume expansion (formed within existing sediments as permafrost aggrades) typically results in uplift, altering local morphology and producing a landscape of domed mounds at a range of scales (Gilbert et al., 2016). Massive ice contains only trace concentrations of organic carbon but can account for over 90% of the eroding cliff face exposure where ice cored domes are exposed at the coast (Pollard, 1990).
Ground ice layers are sporadically and often unpredictably distributed (Mackay, 1972;Pollard, 1990), restricting the best estimates of volume and inferred carbon release to a standard factor to account for ice content in erosion rates (e.g., an average value of 19% used by Couture et al., 2018). The variation in massive ice presence and associated processes often remains restricted to visible exposures and boreholes (Dallimore et al., 1996), which are often limited in distribution and accessibility. Consequently, the role of massive ice in controlling erosion rates and failure processes is generally omitted from models projecting permafrost coast change (Smith et al., 2007), with notable exceptions (for example Dallimore et al., 1996). Here we analyze historic aerial imagery to establish past trends in a key massive ice type-site on the ice-rich permafrost coast of the Beaufort Sea, Arctic Canada. A sequence of volumetric changes shows that despite high rates of coastal erosion the area was dominated by the rapid collapse of a massive ice dome. Within this context we use a passive seismic survey approach to map massive ice layers for the first time, to better understand and constrain the future evolution of the site and to assess the implications for improving estimates of organic, potentially carbon-rich material release into the sea.

Regional Setting and Methods
The Tuktoyaktuk Coastlands, northwest Canada, are located within the southern Beaufort Sea area and are dominated by ice-rich permafrost landscapes (Rampton, 1988). This generally flat, deltaic landscape is punctuated by ice cored plateaus and domes across the tundra. Here, we focus on Peninsula Point (Figure 1a), 6 km southwest of Tuktoyaktuk, and within the Pingo Canadian Landmark (a national historic site managed by Parks Canada), it is a representative type-site for intrasedimental massive ice (Gilbert et al., 2016;Mackay & Dallimore, 1992;Murton, 2009). Compounding the dynamic erosion processes in the area (Murton, 2005), reduced sea-ice (Overeem et al., 2011) and frozen ground seasons (Laberge & Payette, 1995;Liljedahl et al., 2016), increasing storm intensity (Vermaire et al., 2013) and a relative 2.5 mm a −1 sea-level rise (Hill et al., 1993) have intensified the degradation of permafrost coasts in the region.

Change Detection: Shoreline Retreat and Photogrammetric Modeling
Aerial images from the Natural Resources Canada archive (https://www.nrcan.gc.ca/maps-tools-publications/satellite-imagery-air-photos/air-photos/national-air-photo-library/9265) have been georeferenced using common features to those identifiable within a 2018 orthoimage. Check points excluded from the georeferencing process have been used to produce the error bounds in the shoreline positions. The older positions have been georeferenced from the work of Mackay and Dallimore (1992), and as neither the raw images nor the calculated registration errors were available, the highest check point error (15 m in 1985, accounting for possible image distortion, plotting accuracy, and georeferencing errors) has been assumed for the shoreline positions in 1935;1971 and1950 (Figure 1b). Rates of change between survey periods have been established within ESRI ArcGIS, calculated by densifying the digitized shorelines to a measurement every 10 m before computing the mean Euclidean distance between positions (Figure 1c). The shoreline positions in the aerial images may also be influenced by small tidal variations or sea-level changes over the intervening period. Therefore, although beyond the error bands used (15 m), the coastline retreat rates should be interpreted with caution.
Historic aerial imagery (10 images with at least 60% overlap) of Peninsula Point from 1994 have been photogrammetrically analyzed in Agisoft Metashape © using Structure from Motion processing to establish the decadal scale volumetric changes occurring at the site. The 1994 photogrammetric DEM (11 points per m 2 ) has been compared to a 2004 lidar flight (collected at 1 point per m 2 ) and a 2018 UAV photogrammetric survey (400 points per m 2 ). The 2018 surface has been derived from 480 12 MP images (captured from a DJI Phantom 4 at a flying height of 70 m). Three-dimensional control points for the 1994 survey have been taken from common features with the 2004 lidar surface using areas away from any signs of activity, breaks of slope, or significant changes in topography. The large amounts of genuine change posed significant challenges to the number (6) and distribution of checkpoints achieved, selected in the same way but excluded from the registration process; in reality these points were often close to the control points used. A root mean square check point error of 0.7 m was recorded, but the 1994 photogrammetric surface was consistently higher on the landward (assumed stable) side than in 2004 ( Figure 2e). It cannot be determined whether this is a genuine surface lowering (approximately 0.5 m) from 1994 to 2004 or the result of vegetation changes for example, but the trend does not appear from 2004 to 2018. Therefore, error bounds based on the standard deviation of elevation differences in "nonchange areas" have been used in the interpretation of change. This resulted in error thresholds of ±0.7 m in the 1994-2004 changes and ±0.2 m in the 2004-2018 changes.

Passive Seismic Characterization of Massive Ice Spatial Variability
Passive seismic surveys have been conducted in order to map and characterize the spatial distribution of subsurface layers across the site at Peninsula Point. Tromino ® 3G ENGY instruments record three component accelerometer measurements of the ambient noise that is generated by microseisms from the Earth's interior and surface waves, caused principally by wave action at the shoreline (Kumar et al., 2018). The microseisms and surface wave signals propagate through the near-surface layers. Notable differences between amplitudes of horizontal and vertical motions are caused by acoustic impedance contrasts at layer boundaries and differences between shear-wave and compressional-wave velocities within each surface layer. The spectral ratio of horizontal over vertical responses exhibits resonance frequencies (fr) that can be directly linked to the layer thickness (h). The layer thicknesses are estimated by the quarter wavelength resonance approximation that was initially proposed by Joyner et al. (1981) and then later optimized for a wide variety of ground-motion metrics in diverse ground conditions (Boore, 2003) but primarily for generic rock profiles rather than frozen ground. Measurements of visible contact layer depths exposed in permafrost cliff sections adjacent to a set of control reading locations have been used to characterize relative acoustic impedance contrast signatures and subsequently to calculate shear-wave velocity (V s ) through each of the main material  Figures 4 and 5, and the black dot shows the approximate location of a borehole recording massive ice thickness, presented in Mackay and Dallimore (1992). A UTM reference scale is presented on the outer axes.
layers (see supporting information). This has enabled the determination of two distinct high frequency signal peaks associated with specific ground layer controls. A single peak in the resonance frequency (usually around 20 Hz) denoted active layer or slurry contact directly onto massive ice and a single deep later reflector, whereas traces with a high frequency peak around 30-50 Hz (usually with additional lower frequency peaks of 1-5 Hz) represented a permafrost to ice contact. There was no validation for the ice base contact, which was not exposed in the cliff, but borehole data from Mackay and Dallimore (1992) indicate that the ice in the areas is between 10 and 20 m thick. Once the impedance contrast signatures were classified, the following formula was applied to calculate and validate layer boundary depths (h): A shear-wave velocity of 920 m s −1 for the permafrost was determined from cliff outcrop data at the site and was used to derive depth to the massive ice contact (Carcione & Seriani, 1998). Validation sites have also been used to make sure that layer properties were sufficiently consistent over the survey area (see supporting information). These shear-wave velocities matched the control site surveys to known contacts but are at the lower end of what would be expected for permafrost shear-wave velocities, potentially reflecting areas of talik formation or degradation within the permafrost body, changing phase behavior, and slowing the transmittance speeds (Carcione & Seriani, 1998). Sites were selected to reflect changes in topography, usually at the top, middle, and base of a slope. Further research is required on whether transmittance speeds vary through the thaw season. The permafrost to ice contact depths, where present, have been combined with active layer probing and subtracted from surface elevations derived from the high-resolution UAV DEM. The massive ice surface has been modeled across the active areas of the site using smooth interpolation to reflect the continuous undulations of the surfaces exposed in the cliff sections.

Patterns of Change and the Role of Massive Ice
The mapped positions of the shoreline at Peninsula Point reveal a relatively continuous retreat rate of 3.47 m a −1 since 1935 (Figure 1c). This rate is within the long-term range (1 m-5 m a −1 ) for the area (Solomon, 2005) but lower than the 5 m a −1 rate for ice face retreat the site determined by French et al. (1983) and significantly higher than rates from 1952 (0.5-1.3 m a −1 ) established to the west along the Yukon coastline (Lantuit & Pollard, 2008;Radosavljevic et al., 2016). Air temperatures modeled using the Canadian Global Coupled Model ensembles were also projected to rise at a relatively constant rate over this period, amounting to an approximate 4°C increase from 1960 to 2020 (Manson & Solomon, 2007). Where excess ground ice melts the ground surface subsides and conditions change. Heat fluxes and hydrology are altered significantly (Aas et al., 2019), and ice-rich permafrost in contact with rivers and coasts undergoes thermal erosion and denudation processes . The shape of the coastline at Peninsula Point switched from convex in 1935 (i in Figure 1b) to concave in 1985. Early reports from the site refer to a flattish ice-cored hill (French et al., 1983), and supported by the morphological changes, it appears that the center of the massive ice was seaward of its current position (i in Figure 1b) as suggested by Mackay and Dallimore (1992). Superimposed on this wide-scale morphodynamic transition are localized responses. For example, an exposed slump backwall (iii in Figure 1b) supplied high amounts of saturated material that flowed out to form a fan (ii in Figure 1b) that resulted in a localized progradation of the coastline.
Photogrammetric analysis of historic aerial imagery has enabled characterization of long-term change at the study site (Figure 2). A notable feature has been the subsidence of a large connected dome ( Figure 2d) and channel (Figure 2c) feature. Localized thresholds and feedbacks such as snow and water flow accumulation are likely to play a key role in determining subsidence rates (Aas et al., 2019;Seppälä, 2011). The sections through the different ages of surfaces (1994; 2004 and 2018) also show a high degree of consistency and progressive phased development in the trends of change. Over approximately decadal time steps, the island has subsided at differential rates to reach a similar base elevation (Figures 2d and 2e), switching its overall topography from a relatively coherent dome profile to a retrogressive thaw slump (Burn & Lewkowicz, 1990;Lantuit et al., 2012). The landward profile sequence (Figure 2e) shows that coincident to the subsidence, a seaward encroachment resulted from fans of saturated thermokarst material released into the nearshore. The fan was removed within 14 years, by which time wave erosion had truncated the toe of the slope.
Subsidence of ice centered "earth hummocks" in the area has been reported at rates from 0.2 to 1 m a −1 (French et al., 1983). Since 1994 the central area of the site (covered by Figure 3) Mackay and Dallimore (1992) for validation, locations that were massive ice dominated have subsided up to 5 times faster than the island as a whole, but the 60% reduction in the difference between subsidence and coastal recession likely reflects a reduction in the supply of massive ice over time. It is noteworthy that the borehole location (Figure 1b) in 1988 was approximately 40 m from the exposed cliff face and 10 m above sea-level, but that site is now only 3.7 m above sea-level and seaward of the exposed cliff. Accounting for both the total measured subsidence and the landward retreat of the coastline, the site has experienced an average loss of 59 m 3 for each meter width strip of coastline annually. Within these landscape shifts, the mechanics of change reflect the rapid lateral expansion of slump backwalls that expose ground layers to thermoerosional processes such as insolation and precipitation, but once sufficiently insulated by a drape of failed material the vertical erosion (in particular) is greatly reduced. These switches from supply to transport limited system behavior operate over a subdecadal timeframe (Farquharson et al., 2019).

Mapping Subsurface Structure of Massive Ice
This work demonstrates that passive seismic surveys can determine the presence and depth of massive ice within challenging permafrost coast ground conditions. Stratigraphic borehole data have shown that the base of the massive ice in the area is relatively consistent but that there is significant variability in the surface elevation of the ice (Mackay & Dallimore, 1992;Moorman, 1998). The passive seismic surveys show concentrations of massive ice rise to high points in the landscape (Figure 4a; darker purple areas). Profiles through the ice and ground surfaces demonstrate that topography is still strongly associated with massive ice presence (Figures 4b and 4c). However, areas where the ice rises close to the surface, which have experienced high recent rates of subsidence, are beginning to invert this pattern with localized subsided lows where near-surface ice has melted (Figure 4b). The channel feature that has subsided significantly over the last few decades (Figure 2c) is still underlain by near-surface massive ice and thus remains susceptible to further subsidence through local processes (Aas et al., 2019). The undulating surface of the ice appears to overlie a larger north-south trend as the ice contact dips to the south of the island, reflecting an increasing dominance of permafrost in the ground layers and potentially a reduced sensitivity to rapid subsidence (Obu et al., 2017). The spatial variability of ground ice and its direct impact on thaw-induced surface deformation has recently been demonstrated in both small scale (an 11 m × 10 m plot) field experiments (Wagner et al., 2018) and investigation of localized thermokarst signals (Paquette et al., 2020).
The wider impacts of massive ice presence and sensitivity to change are related to both the direct release of carbon bearing material (Couture et al., 2018) and to the indirect changes to the exposure of highly labile, nutrient-rich, and bio-available material across retrogressive thaw slump (Abbott & Jones, 2015) and nearshore regions (Tanski et al., 2019). As an ice-rich area sequences through rapid subsidence and backwall recession it produces an increasing transition zone that remains largely unvegetated while the supply of meltwater continues (Kokelj & Burn, 2003). Gas emissions from these retrogressive slump areas and rapid thaw areas have the potential to cause a disproportionate release of carbon, particularly in methane resulting from ground saturation (Turetsky et al., 2020). It is therefore essential to better identify the spatial variation of massive ice and to constrain the temporal responses of these climate sensitive areas as they transition through phases of dynamic mass wasting to ice exhaustion, stabilization, and revegetation (Farquharson et al., 2019). By subtracting our modeled massive ice layer from the ground surface, a spatial map of the remaining organic-rich material present at the site (1.1 M m 3 of the total 1.4 M m 3 within the surveyed area) has been produced ( Figure 5). The overall massive ice content within the ground is comparable to values used elsewhere (22% compared with an average of 19% used in Couture et al., 2018, on the Yukon Coastal Plain), but the spatial variability ranges from comprising effectively all of the above sea-level ground profile to being absent from it. These approaches allow spatial estimates of contributions from ice melt relative to permafrost thaw to be partitioned.

Conclusions
Massive ice collapse remains a critical but poorly quantified driver of ice rich permafrost coast erosion processes and as such is a threat to Arctic ecosystems, infrastructure, and biogeochemical cycles (O'Neill et al., 2019). Taking an ice-rich permafrost type-site as an example, with clear and variable exposures of massive ground ice, we set the long-term erosion trends in context and use photogrammetric analyses of historic imagery to provide new quantitative data on decadal scale geomorphic processes. A relatively constant rate of change has been derived from shoreline mapping, which may reflect the influence of a sustained presence of massive ice at the site; however, significant doubt is cast on the accuracy and validity of this metric for understanding and predicting geomorphic responses under changing conditions. The volumetric change within the surveyed area associated with subsidence has been over 3 times greater than the volumes lost to coastal erosion since 1994. Passive seismic mapping of the site has identified the remaining ice surface contact zone and revealed several remaining near-surface ice masses that are likely to be highly vulnerable to future subsidence (Farquharson et al., 2019). By accounting for massive ice variability, both future sensitivity to change and potential constituent fluxes from terrestrial and nearshore environments can be better accounted for.