Seasonal Surface Subsidence and Frost Heave Detected by C-Band DInSAR in a High Arctic Environment, Cape Bounty, Melville Island, Nunavut, Canada

Abstract

Differential interferometry of synthetic aperture radar (DInSAR) can be used to generate high-precision surface displacement maps in continuous permafrost environments, capturing isotropic surface subsidence and uplift associated with the seasonal freeze and thaw cycle. We generated seasonal displacement maps using DInSAR with ultrafine-beam Radarsat-2 data for the summers of 2013, 2015, and 2019 at Cape Bounty, Melville Island, and examined them in combination with a land-cover classification, meteorological data, topographic data, optical satellite imagery, and in situ measures of soil moisture, soil temperature, and depth to the frost table. Over the three years studied, displacement magnitudes (estimated uncertainty ± 1 cm) of up to 10 cm per 48-day DInSAR stack were detected. However, generally, the displacement was far smaller (up to 4 cm). Surface displacement was found to be most extensive and of the greatest magnitude in low-lying, wet, and steeply sloping areas. The few areas where large vertical displacements (>2.5 cm) were detected in multiple years were clustered in wet, low lying areas, on steep slopes or ridges, or close to the coast. DInSAR also captured the expansion of two medium-sized retrogressive thaw slumps (RTS), exhibiting widespread negative surface change in the slump floor.

1. Background

Approximately one-quarter of the land area in the northern hemisphere is underlain with permafrost [1], defined as a layer of ground which lies below a seasonally freezing and thawing active layer and remains frozen continuously for at least two years [2]. The Canadian Arctic Archipelago (CAA) is one such environment. Repeated freezing and thawing of the active layer, along with subsurface ice loss, can lead to distinctive geomorphological processes such as large-magnitude localized surface subsidence, landslide-like active layer detachments (ALDs) [3], and long-lived mass-movement events known as retrogressive thaw slumps (RTSs) [4]. Collectively, processes associated with the melting of massive ice are described as thermokarst [5]. In addition to this, more widespread, ‘isotropic’ subsidence occurs seasonally due to thaw consolidation in areas of continuous permafrost, typically with smaller magnitudes [6]. Correspondingly, freeze-up in the autumn can cause surface uplift or frost heave due to the expansion of water undergoing a phase transition [7] and through the upward expansion of excess ice. Specifically, if there is high water availability at depth, ice lenses can form by meteoric water infiltrating soil and becoming trapped near the active layer (permafrost layer transition zone) [8]. Ground ice content [9,10], soil moisture [11], and soil composition [7] are the primary predictors of an area’s susceptibility to these disturbances. The local surface energy flux, controlled by air temperature [12], snow cover [13], and microtopography [14] is also important as it controls active layer thaw penetration in any given season [15]. The timing of such surface displacements has also been shown to be closely related to ground temperature dynamics at the surface and at depth [16]. Precipitation is also a component of the direct heat exchange with frozen ground and has been incorporated as a climatic factor in models monitoring surface deformation in permafrost environments [17]. Vegetation cover can also affect surface deformation, as dense vegetation and the related soil litter will have an insulating effect, reducing the impact of spring thaw and fall freeze-up on active layer thickness (ALT) [18]. Recent research has shown that rapid heave events are synchronous with precipitation events, suggesting that water balance and possibly lateral water movement in the subsurface are factors involved in these processes and pointing to a ‘thermo-hydro-mechanical’ system far more complex than previously thought [19]. Furthermore, processes on several different time frames are at play and interact with each other. While the main freeze-up and thaw sequence occurs annually, freeze–thaw transitions at the transient layer (between the active layer and permafrost table) may occur with a periodicity of multiple years or multiple decades [20].

These processes not only have implications for the carbon and hydrological balance in this region [21] but make engineering works, particularly the construction and maintenance of roads and air strips (infrastructure essential to remote Northern communities), challenging to design and expensive to implement [22] as surficial disturbances can cause built materials to buckle and crack. Furthermore, an estimated 2000 Pg (20 billion tonnes) of carbon is presently stored in the ground at high latitudes [23] where permafrost currently impedes its decomposition in the subsurface and into the atmosphere. Hence, carbon released from thawing permafrost as the mean annual air temperature (MAAT) rises is likely to be substantial and contribute to global warming [20,23], which exacerbates the problem of thawing permafrost. This intensification of permafrost thaw during the past two decades has been observed in the Canadian Arctic [10]. It is therefore of great utility and importance to be able to effectively monitor and predict the surface displacement resulting from permafrost degradation in order to assist with the planning of engineering work and with carbon budgeting in Arctic ecosystems.

Vertical deformation in permafrost environments has previously been measured at a fine scale using instruments such as heave sleeves [24] and tilt arms [19]. These point measurements can be difficult to interpret in highly heterogeneous areas and suffer from potential measurement errors associated with thermal expansion and contraction of the instruments themselves. Remote sensing and surveying methods, which allow for repeated monitoring of larger areas often at a significantly lower cost, have also been used, namely differential global positioning systems (DGPS) [6,25], light detection and ranging (LiDAR) [26], and structure from motion (SfM) photogrammetric techniques [27].

Another technique is differential interferometry of synthetic aperture radar (DInSAR), which can generate surface displacement maps by converting phase differences in SAR backscatter signals collected on two different dates to line-of-sight (LOS) displacement along the vector connecting the ground to the sensor [28]. A weather-independent, high-spatial resolution SAR backscatter signal gives DInSAR data analysis the potential to monitor surface change over large areas with high spatial detail and at repeated intervals. Snowmelt, alongside large precipitation events, can cause challenges when examining DInSAR in High Arctic environments, causing incoherence in InSAR pairs due to increased soil moisture and the large difference in backscatter characteristics between snow and the soil surface [29,30]. Low levels of vegetation make tundra environments fertile ground for DInSAR, reducing the potential for loss of coherence due to complex radar-vegetation interactions [31]. DInSAR has previously been used to map permafrost-related surface displacement for: (i) high latitudes [15,32,33,34], including specifically the mapping of retrogressive thaw slumps [35]; (ii) high altitude environments [36,37,38]; and (iii) built environments [39,40]. Attempts are also being made to develop fine-scale, pan-Arctic ground ice maps using the distribution and magnitude of subsidence detected by InSAR [41].

This research examines repeat-pass DInSAR with ultrafine beam C-band (λ = 5.5 cm) Radarsat-2 data with a spatial resolution of approximately 3 m to generate surface displacement maps with high precision for the late summer seasons of 2013, 2015, and 2019 at the Cape Bounty Arctic Watershed Observatory (CBAWO), Melville Island, Nunavut. The displacement maps are interpreted alongside optical satellite imagery, a systematic mapping of land cover [42], and field soil moisture and temperature measurements in order to address the principal objective of describing the land cover and topographic characteristics of those areas that display surface displacement in one or multiple years. This is done with the aim of improving the understanding of both seasonal and secular land movement processes occurring at the CBAWO and informing the input variables to future permafrost disturbance susceptibility models, which may be scaled up to the region and the wider North American Arctic. This paper extends previous research in this field by: (i) utilizing the high spatial resolution and high coherence potential of C-band ultrafine beam SAR in a heterogeneous tundra environment, together with a detailed land cover map, to study how displacement patterns vary by land cover type, topography, and in proximity to ALDs; and (ii) testing, as a proof of concept, the utility of two simple indices derived from DInSAR maps to characterize longer term displacement patterns over multiple summers in the case where multi-year, coherent DInSAR stacks are unavailable.

2. Methods

2.1. Overview

The principal objective of this work is to relate surface displacement as measured by DInSAR to landscape characteristics and climate conditions during each of the study years. Therefore, the subsequent methods consist of three main parts: (i) collection and analysis of field data; (ii) acquisition and processing of optical remote sensing data and the derivation of wetness and topographic indices; and (iii) acquisition of SAR data and calculation of surface displacement maps via DInSAR.

2.2. Study Area

Fieldwork was completed in July and August 2020 at the CBAWO located on the southern coast of Melville Island, Nunavut, Canada (74°55′N, 109°35′W), 425 km west of Resolute Bay (Qausuittuq) (Figure 1). The site features gently sloping topography and a paired watershed system with two large lakes (West Lake and East Lake, unofficial names). The area is underlain with continuous permafrost, with a typical ALT of 0.5–1 m and measured temperatures at the top of the permafrost (TTOP) between −13 and −15 °C [43]. Vegetation rarely grows to more than a few centimeters and is primarily composed of prostrate dwarf shrub tundra [44]. Melting temperatures typically persist from June through August [45], with an average July air temperature of 4.0 °C, and 24-hour daylight persists from mid-April to mid-August. The area is geomorphologically active, featuring more than 100 large ALDs, the majority of which occurred during an unusually warm summer in 2007 [21], and several active RTSs. This site has been intensively studied since 2003, ensuring long-established weather, hydrological, soil and geomorphological observations and records are available. Ground penetrating radar (GPR) has revealed highly variable ground ice composition at the CBAWO, even within localized areas [46]. A high resolution (1.6 m) land-cover classification has been produced for the CBAWO, based on topography-derived indices, multi-spectral satellite imagery, and multi-year field observations [42]. In addition to identifying areas of water and snow, the classification divides the landscape further into three broad vegetation classes: polar (semi) desert (35.1%), mesic tundra (34.7%), and wet sedge (16.6%) (Figure 1). In areas of polar (semi) desert -which are found mostly in upland areas—vegetation can only be found in the moist, nutrient-rich depressions around the edges of polygons, which increases the ease at which polygon features can be delineated on the landscape. In contrast, areas classified as mesic tundra and wet sedge are characterized by a high abundance of graminoids and bryophytes, with wet sedge areas generally having higher phytomass and soil moisture content [47]. The bare ground class is generally characterized by areas devoid of vegetation or by areas of exposed bedrock, gravel, or frost shattered rock (felsenmeer).

2.3. Field Measurements: Frost Depth, Soil Moisture, and Meteorological Data

Soil moisture data aids our interpretation of the DInSAR-generated surface deformation maps because radar backscatter is highly dependent on soil moisture [49,50]. Soil moisture also contributes to surface energy flux and active layer development [51]. Thus, in conjunction with soil temperature at depth [16], it is anticipated that areas with high soil moisture content will be more susceptible to displacement. Soil moisture and temperature were recorded at a depth of 5 cm at 3-hour intervals via Decagon em50 loggers equipped with 5TM moisture/temperature/conductivity sensors available for polar (semi) desert (Plateau, 2015-19), mesic tundra (Upper goose and Upper ptarmigan, 2013-19), and wet sedge (Muskox, 2015) sites (Figure 1). Additionally, depth to the frost table measurements taken at different times and in different conditions reveal the extent of thaw development in any given season. Frost depth was measured on 9 July 2019 by pushing a metal probe to resistance at 33 sites around the CBAWO and repeated on 2 August 2019, to align with SAR acquisition dates (Figure 1). Five measurements (four measurements each 5 m from a central point measurement) were taken at each site and averaged.

Shielded air temperature and precipitation were recorded hourly at the WestMet station, located in an upland area of the CBAWO at approximately 90 m a.s.l (Figure 1). Daily temperature and precipitation averages were calculated for the summers of 2013, 2015 and 2019. Meteorological data from the CBAWO was unavailable after 11 August 2019. As a replacement, temperature readings from the Environment and Climate Change Canada weather station at Mould Bay, Prince Patrick Island (76°14′N, 119°39′W), which is highly correlated to conditions at Cape Bounty [52], were used. To better characterize temperature conditions for each season, temperature time series were converted to thawing degree day (TDD) totals, calculated annually as the sum of daily average temperatures for all days with a daily average above 0 °C [53]. Similarly, freezing degree days (FDD) were calculated for each season. However, since the displacement maps will only capture heave associated within the stack date range, FDD was only calculated during the interval of the stack (i.e., stack FDD or SFDD).

2.4. Optical Satellite Data

To help identify landscape features, interpret the generated displacement maps, and calculate a spectral wetness index (Section 3.3), multispectral WorldView-2 images of the study site with eight bands ranging from 400–1040 nm were acquired on 12 August 2015 and 4 August 2019. The 2015 image is cloud-free and the 2019 image has minimal cloud cover (6%) mostly constrained to the far north-east corner of the study area. Orthorectification of the 2019 image was completed using the geometric function suite within ArcMap [54] using the rational polynomial coefficient .rpb file provided with the image product. The quality of the orthorectification procedure was assessed using six control points with known positions; an average positional error of 1.61 m was determined. The 2015 image (1.84 m resolution) was co-registered to the 2019 orthorectified image (1.24 m resolution) in ArcMap using the nearest neighbour method. An absolute radiometric correction was carried out for each band to convert pixel values (DN) to top of atmosphere reflectance via the method presented in Kuester (2016) [55].

2.5. Wetness Indices

High soil moisture conditions can increase the susceptibility of an area to disturbance, especially via terrain destabilization on slopes due to higher soil pore water pressure [56] and can lead to ALD formation [11] (Lewkowicz and Harris, 2005). High water availability, delivered to the base of the active layer via runoff and infiltration, may increase the magnitude of frost heave through added ice volume [24], which will be captured as uplift in the DInSAR. To estimate soil moisture across the study area, a linear combination of the TOA reflectance values of all eight bands of WorldView imagery with band coefficients (see supplementary Table S1) for wetness/shadows from Yarbrough et al. (2014) [57] was constructed via the Kauth–Thomas method [58] for the 2015 image. Additionally, a topographic wetness index (TWI) based on the effects of local topography on flow direction and runoff accumulation, was calculated via Equation (1):

$$TWI=\ln \left(\frac{\alpha }{\tan \beta }\right),$$

(1)

where $\alpha$ is a scaled flow accumulation, found by multiplying the flow accumulation by the digital elevation model (DEM) pixel size, and $\beta$ is the local slope in radians [59]. Flow direction, flow accumulation, slope, and TWI were calculated in Whitebox GAT [60] based on a 1 m-resolution DEM of the CBAWO derived from a summer 2012 WorldView-2 stereopair [47]. The study area was segmented into zones of low, medium, and high TWI based on 33rd and 66th percentiles of TWI.

2.6. SAR Imagery

Six Radarsat-2 ultrafine beam C-HH (λ = 5.5 cm) SAR images were acquired for the 2019 season, in addition to three previously acquired images for each of 2013 and 2015. All were captured at approximately 19:14 local time from a right-looking, ascending orbit, with an incidence angle of 33°, the same as the 2013 and 2015 images. The footprint of each SAR image was large enough so that only one image was required to cover the entire study area and no image mosaicking was required. Radar wavelengths within the C-band have been shown to work well for differential interferometry at high latitudes [61], and experience significantly less ionospheric refraction than the L-band signal [62], which would result in a higher signal-to-noise ratio. The temporal coverage was designed to align with the summer field season and was chosen to minimize loss of coherence from snowmelt in the spring (May-June) and fresh snowfall in the fall (late September). The images were delivered with an approximate spatial resolution of 4 m and in single-look complex (SLC) format, where each pixel contains information about both the magnitude and the phase of the radar return signal.

2.7. DInSAR Processing

Differential interferometry with Radarsat-2 images processed in GAMMA v.20191209 software [63] was used to generate summer season displacement maps for 2013, 2015, and 2019. ArcticDEM tiles (2 m resolution) [64] were downloaded for use in the topographic correction during DInSAR processing as they provided full coverage across the SAR swath, while GDAL Warp was used to mosaic these tiles together. Next, 2 × 3 multilook intensity images (MLIs) were generated for each SAR image to reduce speckle [65]. Perpendicular baselines were calculated using the orbital data delivered with the imagery and examined, and the image of 16 July 2013 was chosen as the reference image. A simulated SAR intensity image was created based on the DEM and the reference image co-registered to it. Offsets were calculated for the MLI images and combined with the DEM to resample the SLCs using look-up tables to complete co-registration. Differential interferogram pairs for every available pair of dates were generated using baseline refinement and adaptive filtering, before being unwrapped using a minimum cost flow (MCF) algorithm. Given the high spatial resolution of the displacement map, and the objective to identify more significant spatial patterns, an adaptive and weighted interpolation procedure using a search radius of 6 pixels was applied to the interferograms to remove any small patches of incoherence. These areas of incoherence can be useful for informing the interpretation of displacement maps (especially when identifying areas of snow [66]), so the interpolation parameters were set such that all larger areas of incoherence were preserved. Based on the coherence patterns, as well as the presence or absence of fringe discontinuities within each interferogram pair, a qualitative assessment was made to select the high-quality pairs (

Table 1. Interferometric pairs selected for stacking.

Year	Stack Dates	Pairs Used	Temporal Baseline (days)
2013	16 July–2 September	16 July–9 August 9 August–2 September	24 24
2015	30 July–16 September	30 July–23 August 23 August–16 September	24 24
2019	9 July–26 August	9 July–2 August 2 –26 August	24 24
2019	2 August–19 September	2 –26 August 26 August–19 September	24 24

) for use in the creation of the 48-day displacement maps by way of a stacking procedure in GAMMA. This process assumes that atmospheric statistics remain stationary through the period covering all interferograms, and so stacked phase rates (radians/year) may be calculated using N interferograms using Equation (2):

$${\phi }_{rate}=\frac{{{\displaystyle \sum }}_{j=1}^N\mathsf{\Delta }{t}_j{\phi }_j}{{{\displaystyle \sum }}_{j=1}^N\mathsf{\Delta }{t}_j^2},$$

(2)

where j is each interferometric pair, $\mathsf{\Delta }t$ is the temporal baseline of the pair, and $\phi$ is the phase difference within the pair. Generally, all phases should be unwrapped with reference to a point with known displacement or to a point which is not expected to be displaced. The selection of this point is important, as all displacements on the map are not absolute but rather relative to the motion of this ‘stable’ point [7]. A known area of exposed bedrock (WGS84, UTM Zone 12X: 541577, 8313180), previously used in DInSAR processing for this area [34], was chosen as the reference point. Finally, displacement maps were generated which convert the stacked phase changes into a vertical deformation rate via Equation (3):

$$\mathsf{\Delta }{z}_{rate}=\frac{\lambda }{4\pi \cos \theta }{\mathsf{\phi }}_{rate},$$

(3)

where $\lambda$ is the radar wavelength (here 5.5 cm), $\theta$ the incidence angle (here 33°) of the incoming radar beam, and $\mathsf{\phi }$ the calculated phase difference. The resulting map of deformation rate has units (m/year) so a scale factor is applied equal to the fraction of a year that the stack dates represent in order to yield an average displacement value in metres. The vertical displacement maps are calculated from LOS maps with the assumption that all motion happens normal to the ground surface. The displacement maps were then geocoded from Doppler radar coordinates to UTM with a spatial resolution of approximately 5.5 m.

Noise and atmospheric effects introduce error to the calculated surface deformations. This error was estimated by extracting the average absolute measured vertical change within six control polygons at known stable points around the CBAWO (Table S2), namely areas of bare rock or felsenmeer (rubble) which should in principle show zero vertical change [7]. The largest average absolute displacement found across all six polygons was 0.97 cm in the 2019 stack. Therefore, throughout interpretation of the displacement maps, any net vertical change of <±1 cm will be considered to be within the measurement error and described as ‘stable’. This value is consistent with DInSAR error in high latitude environments [30,34].

2.8. Derivation of Metrics to Characterize Long-Term Spatial Patterns of Displacement (2013-19)

Since high quality interferometric pairs spanning multiple years were unavailable due to temporal decorrelation, two DInSAR-derived rasters were generated to characterize long-term spatial patterns of surface displacement for 2013-19, which are designed to offer an approach for characterizing long-term landscape change when only non-consecutive and disjointed SAR datasets are available. First, the aggregate (net) displacement from the 2013, 2015, and 2019 (late season) stacks was calculated by simply summing the values for each pixel using the raster calculator in ArcMap; and second, the 2013, 2015, and 2019 DInSAR displacement maps were used to generate zonal maps showing pixels which displayed significant net displacement (>2.5 cm) in zero, one, two, or all three of the DInSAR maps using reclassification in combination with the raster calculator tool in ArcMap. This 2.5 cm threshold was selected as it is approximately two standard deviations from zero of the measured displacement across the entire scene for all years and is designed to highlight the most extreme surface change in a given stack. Within each zone, 1000 random points were generated and the values of slope, elevation, vegetation classification, and wetness indices were extracted for their use in ANOVA and chi-squared tests to provide a statistical characterization of the landscape features within each zone.

3. Analysis

3.1. Meteorological and Soil Conditions

Based on in situ measurements collected on 2 August 2019, sites situated in wet sedge meadows had the shallowest average frost depth (56.7 cm), while polar (semi) desert sites had the deepest (73.2 cm). Across all sites, the mean frost depth increased by 23% from 51.7 cm to 63.8 cm in the 24-day period between 9 July and 2 August 2019. Based on long-term data from soil stations, 2013 saw an early onset of frozen soil conditions (Figure 2) and an extended ‘zero-curtain’ period (where soil temperatures fluctuate closely around freezing for several days or longer while latent heat transfer occurs) [67], which would imply higher availability of unfrozen water in the active layer at freeze-up creating a large distributed latent heat effect [68]. The onset of soil thawing conditions also came later in 2013 as compared to 2015 (June 5 and May 12 respectively). Soil temperature and moisture data for 2019 were unavailable for dates after 10 August, so did not capture the freeze-up period. Differences in measured soil moisture at the stations between years was generally small, although at the sites located in mesic tundra areas, the 7 day mean soil water content was 13% higher in early August 2019 (0.25 m³/m³), compared to 2013 (0.22 m³/m³), and 20% higher when compared to a similar date in 2015 (0.21 m³/m³).

3.2. Early versus Late-Season Displacement (2019)

The dates of the 2019 SAR acquisitions allowed for the creation of two 48-day displacement maps, an early stack (9 July–26 August) and an overlapping later stack (2 August–19 September) (Figure 3). Radarsat-2 acquisitions were also made on 22 May and 15 June, but interferograms involving these early season SLCs had insufficient coherence due to spring snowmelt. In these two stacks, vertical displacement greater than 5 cm was detected, with a maximum value of −11.3 cm, detected in the late stack. Mean displacement across all pixels was −1.04 cm in the early stack and −0.31 cm in the late stack, with a standard deviation of 1.3 cm in both cases. As expected, the late season map captured freeze-up into September, displaying higher levels of surface uplift (> 1 cm) than the earlier stack (11.4% of total area versus 5.5%, excluding bodies of water) (Figure 3). The majority of the uplift detected in the early stack was located around streams and in streambeds, the areas which are most susceptible to phase unwrapping errors. However, some of the positive surface displacement captured here may be due to the aggregation of sedimentary materials in streambeds. We consider it possible that a large portion of this uplift is due to phase unwrapping errors emerging from the MCF algorithm. Conversely, a larger percentage of the total land area displayed subsidence (> 1 cm) in the early stack (43%) as compared to the later stack (27.4%). Areas of net subsidence captured in the late season stack may have an excess contribution from August thaw settlement, and reduced by September frost heave, or may be due to thaw settlement occurring after the initial freeze-up associated with hydraulic processes in the subsurface [19]. Generally, the magnitude of thaw settlement is larger than that of surface uplift, indicating that between early July and mid-September, the average overall net vertical displacement across the scene is likely negative. There is also a higher percentage of the scene below the coherence threshold in the early season stack as compared to the later stack (7.3% versus 3.4% of the study area, respectively). This change may relate to changes in snow cover between dates, especially in incised stream channels, but likely underestimates the change in snow cover away from streams as additional incoherence will be introduced in the late season due to increased soil moisture and standing water as a result of this snowmelt and high precipitation levels during late summer 2019 [30].

3.3. Comparison of Single-Season Displacement Maps (2013, 2015, 2019)

Displacement maps for the late summers of 2013 (16 July–2 September), 2015 (30 July–16 September), and 2019 (2 August–19 September) (Figure 4) were compared and examined alongside meteorological and soil conditions. The 2013 and 2015 displacement maps for the same dates were analyzed by Rudy et al. (2018) [34]. Across the CBAWO, vertical ground displacement (> 1 cm) was most widespread in 2019 and 2013 (42.4% and 36.8% of the land area showing either subsidence or uplift, respectively), and less evident in 2015 (21.4%) (Figure 4). Summer air temperature is likely a principal control of the magnitude and extent of subsidence and heave in any given summer [15]. Summers 2013 and 2019 were similar in terms of both TDD (228°C and 255.9°C, respectively) and SFDD (46.6°C and 55.9°C, respectively), whereas summer 2015 was much warmer (TDD = 352.9 °C), with a late onset of freezing conditions in mid-September leading to a low SFDD of 20.1 °C. It can be assumed that the majority of large surface displacement associated with the relatively warm conditions in 2015 would have occurred after initial thaw in mid-June (Figure 2) and in July where most TDDs are accumulated, prior to the first acquisition date, and consequently were not captured, while the late onset of frost lowered the potential for widespread detectable uplift in the late season.

Air temperature conditions were also reflected in the soil temperature at a depth of 5 cm, with 2013 being characterized by an early onset of frozen ground conditions (7 August) and a long zero-curtain period as compared to 2015 (Figure 2). In addition, comparatively low amounts of precipitation fell throughout summer 2015 (a total of 28 mm from 1 May–7 August). In the same interval the total precipitation was 42 mm in 2013 and 82 mm in 2019 (Figure 5a). The extremely wet conditions in 2019 were a result of several large and persistent rainfall events throughout the summer, causing the ground to become saturated. Soil moisture was an average of 13% and 20% higher in 2019 compared to 2013 and 2015, respectively, across all four logger stations. The zero-curtain period of 2013 was also long, suggesting a high water content in the active layer [34] which would increase the potential for ice formation and heave in late 2013. Taken together, this would imply that both air temperature throughout the summer and summer precipitation are important predictors of the levels of displacement (Figure 5a).

This examination only describes the seasonal picture. However, there are also likely longer-term processes affecting changes year-to-year. As an example, examination of optical satellite imagery of the CBAWO from 2003 to 2019 (unpublished) shows a notable decline in late July snow cover during this period. This change will likely impact the water availability and thermal regime of some areas normally fed by late-season snowmelt [68]. The reduced insulation effect by summer snowpack can also lead to both a longer thaw duration and a greater potential for subsidence and cooling along with potential aggregation of permafrost during winter. What is worth noting is that even though 2015 was the warmest of the three summers studied based on total thawing degree days (TDD = 353 °C), it was only the fifth warmest summer on record since 2003 (with 2011 being the warmest; TDD = 465 °C). Assuming these processes are controlled primarily by air temperature [12], years predating 2013 could have seen much larger levels of surface change. In fact, this was observed in the warm summer of 2007, during which deep ALTs were measured and large ALDs occurred at the CBAWO [21].

3.4. Displacement Maps in Relation to Environmental and Topographic Variables

Next, three variables (land-cover type, topographic wetness index, and slope) were examined in relation to the three seasonal displacement maps in order to discern the general characteristics of areas susceptible to large vertical displacements. Generally, the largest displacements were found in areas of wet sedge, in wet, low lying areas, and on steep slopes.

3.4.1. Land-Cover Type

Vegetation cover has been used as an indicator of geomorphic stability in Arctic study sites [69] since vegetation may be removed in the development of ALDs and RTSs, and recover in the period after stabilization. The percentage of pixels which showed larger vertical displacements (>2.5 cm) was calculated for each land-cover type for each year with substantial differences observed across land-cover types (

Table 2. Percentage of total area displaying surface displacement greater than 2.5 cm for each summer, within each land-cover type, topographic wetness index level, and slope.

	Land-Cover Type				TWI			Slope
Summer Stack	Bare Ground	Mesic Tundra	Polar (semi) Desert	Wet Sedge	Low	Med	High	<2°	2–5°	5–8°	>8°
2013	0.9	1.3	0.5	3.1	3.2	2.2	2.6	2.4	2.8	3.0	3.1
2015	0.2	0.2	0.1	2.0	1.1	1.0	1.2	0.8	1.0	1.3	2.3
2019	3.9	2.1	2.4	7.1	7.2	6.8	6.2	4.1	5.3	7.6	19.8

). Across all years, wet sedge exhibited the most displacement, with 3.1%, 2.0%, and 7.1% of this class area displaying large displacements in the 2013, 2015, and 2019 stacks, respectively. In contrast, only 2.4% of polar (semi) desert land cover was displaced (>±2.5 cm) in 2019. Additionally, a one-way analysis of variance (ANOVA) using a sample of 10,000 random pixels showed significant differences in absolute displacement between land-cover types for all three stacks (p < 0.05). However, the significance of each pair-wise comparison between land-cover types varied from year to year. Differences between vegetation types may be even greater in reality, as DInSAR is likely to underestimate displacement in saturated soil conditions since radar penetration depth is strongly dependent on soil water content [70]. Thus, DInSAR-measured displacement may be capturing the movement of the wet/dry interface within the soil between acquisition dates rather than true ground surface movement [33,34]. Moreover, McFadden (2019) [71] detected millimetre-scale displacement with inclinometer instruments associated with water table changes in the active layer. Further, the smoothing effect of the adaptive filter used in processing may ‘share’ phase change in one pixel with neighbouring areas of incoherence, thereby diminishing the overall displacement signal in those areas [72]. Vegetation type is broadly constrained by soil water content at the CBAWO [43] and average volume water content at wet sedge sites exceeded 0.6 m³/m³ (essentially saturated) around the acquisition dates in 2019. This means that wet sedge sites are likely to be highly susceptible to these particular underestimations and that the differences between displacement magnitudes in wet sedge compared to drier land-cover classes is likely even higher than presented here. Although measured frost depth did vary significantly by cover type, across the 30 in situ measurement locations, no strong relationship was found (r = 0.20) between the change in measured frost depth during the season and the DInSAR-measured vertical displacement

3.4.2. Topographic Wetness Index

The percentage of the area within areas of low, medium, and high TWI displaying moderate and large amounts of vertical displacement was calculated (Table 2). A one-way ANOVA of 10,000 randomly sampled pixels indicated that absolute displacement did not vary significantly by TWI zone in 2013 (F = 2.16, p = 0.115) but did in 2015 (F = 5.07, p = 0.006) and 2019 (F = 22.12, p < 0.001). The lowest levels of absolute displacement were in areas of moderate TWI, suggesting displacement is more likely to occur both in well-drained upland areas and low-lying (wet) areas around streams.

3.4.3. Slope Angle

A similar analysis was conducted for slope, using the raster generated as an intermediate step in the TWI generation. In 2013, no significant differences were found in the extent (Table 2) or the magnitude of displacements (one-way ANOVA: F = 0.428, p = 0.733) between zones representing low, medium, high, and very high slope gradients. In contrast, in 2019, mean absolute displacement was 37% higher in steep areas than in gently sloping areas. These particularly large differences associated with slope in 2019 may be due to widespread ground saturation compared to the other years, leading to high susceptibility for slope failure in steeply sloping areas [11].

In the absence of more extensive field data and DInSAR stacks that go beyond these three seasons, it is difficult to make conclusions that go beyond these general characterizations. However, it would seem that areas at the CBAWO which are most susceptible to late-season surface displacement are steeply sloping areas with high moisture content, specifically in areas classified as wet sedge. Positive displacement is enhanced during wet seasons, and during seasons with early onset of thaw and associated high FDD throughout late August and into September.

3.5. Surface Displacement Related to Landscape Features

On average, absolute surface displacement was found to be greater inside previously mapped ALDs than outside. Additionally, C-band DInSAR captured the development of a medium sized RTS feature within the study area.

Previous DInSAR displacement maps generated at the CBAWO (2013 and 2015) detected some signs of uplift in the compression zone at the base of ALDs which formed in summer 2007 [34]. This may be due to the DInSAR capturing the growth of syngenetic permafrost associated with the material deposited in the initial ALD slide and associated solifluction [73]. There was some evidence of uplift in some of the detachment compression zones in the 2019 displacement map, but not in all. However, via a two-sample t-test comparing a random 2500-pixel sample inside previously mapped disturbances [48] to the same number of pixels outside (but within 200 m) of the disturbances, in all three late summer stacks the average absolute displacement values were significantly higher inside the disturbed areas (

Table 3. Mean vertical displacement and mean absolute vertical displacement inside and outside former landscape disturbances, and accompanying t-test results.

	Mean Vertical Surface Displacement inside and outside Mapped Disturbances (cm)			Mean Absolute Vertical Surface Displacement inside and outside Mapped Disturbances (cm)
Period	Outside	Inside	t-Test	Outside	Inside	t-Test
16 July–2 September 2013	−0.06	−0.35	t = −12.79, p < 0.001	0.79	0.85	t = 3.14, p = 0.002
30 July–16 September 2015	−0.12	−0.29	t = −7.27, p < 0.001	0.57	0.63	t = 4.14, p < 0.001
2 August–19 September 2019	−0.42	−0.85	t= −14.95, p < 0.001	0.83	1.02	t = 9.27, p < 0.001

). Moreover, the magnitude of net displacement when both negative and positive displacement values were considered was also significantly greater inside disturbed areas. Since many of these areas are now visibly stable, the higher levels of subsidence may be due to bare earth where vegetation was disturbed, leading to different active layer and thawing conditions as compared to the areas immediately surrounding the ALDs. Additional field observations would be required to identify slope failures which are still active, and to test the effectiveness of DInSAR in monitoring their expansion and recovery. Processes associated with permafrost growth due to the deposits in the ALD toe may be slow, thereby failing to give a clear displacement signal in a single 48-day stack.

RTS development has been monitored previously using InSAR techniques [74]. These slumps form when the ablation of massive ground ice causes the overlying soil material to fall and be removed as a liquefied slurry in the slump floor, forming a characteristic steep headwall scarp [75]. Once formed, thawing of the newly exposed ice can cause the headwall to retreat upslope on the order of meters per year [4,76]. In one area just north of the CBAWO (WGS84, UTM Zone 12X: 544000, 8318245), the late season 2019 DInSAR stack captured unstable ground behind the slump headwall of two active medium-sized RTSs, which may be indicative of thawing ground ice (Figure 6). No other major expression of displacement was observed in the surrounding area, and all displacement was coincident and within the RTS limits. This demonstrates that DInSAR may have the potential to assist other techniques such as photogrammetry and light detection and ranging (LiDAR) [77] to estimate material loss in RTSs via an area-to-volume calculation. However, field measurements would be required for verification of any slumping captured by DInSAR.

3.6. Metrics to Characterize Spatial Patterns of Displacement over Multiple Summers (2013-19)

In the absence of continuous SAR data covering multiple years, by using two simple indices of multiple-year change we identified several patterns of displacement unfolding over these three years. Statistical tests are used to identify which landscape parameters (land cover type, topography, wetness) may best characterize areas susceptible to large surface displacement in multiple summer seasons.

InSAR pairs with temporal baselines of sufficient length to span from one summer to another (e.g., 2015-19) were generated but had insufficient coherence to provide meaningful insight into long-term, multi-year change. Consequently, two simple metrics were developed to provide insight into spatial patterns of displacement and, potentially, the longer-term processes at the CBAWO: (i) frequency of displacement; and (ii) aggregate seasonal displacement. Although it does not provide information concerning absolute surface change from year to year (as in [7,30]), this analysis helps to identify the characteristics of those areas which are likely to change significantly over multiple years.

3.6.1. Assessing Frequency of Displacement

The resultant frequency raster (Figure 7) indicates that the vast majority of pixels which showed significant net vertical displacement at any point only did so in a single year (

Table 4. Percentage of pixels where significant uplift, significant subsidence or both were captured by DInSAR at any time during the stack periods of 2013, 2015, and the later 2019 stack.

Net Vertical Subsidence (>2.5 cm) (%)			Net Vertical Uplift (>2.5 cm) (%)			Subsidence One Year, Uplift in Another (%)	No change >±2.5 cm Observed (%)
1 Season	2 Seasons	3 Seasons	1 Season	2 Seasons	3 Seasons
6.18	0.89	0.06	1.93	0.02	<0.01	0.06	90.87

), while areas that indicated change during two or three seasons are sparse (~1% of CBAWO) but are clustered in a number of highly localized locations. The locations were mostly in low-lying areas, either close to streams (Figure 7 regions of interest (ROI) 3, 4, 8, 10, 11, 12 and 15), in coastal areas (ROIs 6, 7, 10, 14, 16, 17), or on steep slopes or ridges (ROIs 5 and 13). In ROI 2, high displacement frequencies are coincident with patterned ground indicative of thermokarst development [2]. Another area, ROI 9, was the only major displacement cluster at the CBAWO to show significant uplift in one year and significant subsidence in another. It is characterized by saturated ground to the point of accumulating standing water, and wet sedge cover. It is possible that the wet conditions and water availability from an upslope snowbank at ROI 9 has led to thaw consolidation and ice aggregation at the base of the active layer during the study period.

Areas which were displaced (>2.5 cm) in multiple years were found to be characterized by different topographic, vegetation, and wetness conditions than those which were not (

Table 5. Statistical tests identifying characteristics of areas of frequent change.

			Mean Value in Each Frequency Zone (Years Where >2.5 cm Displacement Occurred)
Variable	Stat. Test	Test Result	0	1	2	3
Wetness (Spectral)	ANOVA	F = 7.98, p < 0.001 *	−1.23	−1.28	−1.34	−1.29
Wetness (TWI)	ANOVA	F = 0.94, p = 0.423	6.10	6.02	6.21	6.08
Slope	ANOVA	F = 39.5, p < 0.001 *	5.80	8.78	8.10	9.72
Elevation	ANOVA	F = 525.0, p < 0.001 *	64.15	48.81	23.76	18.42
Land-cover class	Chi-Squared	χ2 = 570.8, p < 0.001 *	N/A	N/A	N/A	N/A

* Significant test statistic.

). Areas with high-frequency displacement are far more likely to be found in wet, low-lying areas (<50 m a.s.l.) and marginally more likely on steep slopes (Figure S1). Additionally, there were differences in the land-cover make-up of areas between zones (Figure S2). As an example, despite constituting only 14% of the area of the CBAWO, wet sedge cover made up 36% and 57% of the areas which changed in any two and any three seasons, respectively. The distribution of pixels which were displaced in one or fewer seasons across land-cover types was approximately representative of the CBAWO as a whole. This, again, indicates that wet sedge is the most dynamic land-cover type at the study site. No significant difference in TWI (F = 0.94, p = 0.423) was found between zones and, assuming that wetness is a significant predictor of displacement, this suggests that a spectral measure of wetness (such as the Kauth–Thomas method presented here) may be a better approach when using remote sensing imagery to assess susceptibility to landscape disturbance in relation to soil moisture conditions.

These analyses portray a landscape where persistently active areas are highly localized and share similar characteristics (i.e., wet, low-lying, and sloping). Ground ice content at depth has been found to be highly variable around the CBAWO [46], which may explain the highly fragmented spatial distribution of significant displacement throughout the CBAWO. The sparsity of areas changing in multiple seasons and the apparent unpredictability of displacement may indicate that reliable hazard susceptibility models may be challenging to develop for similar environments, although more extensive field data, and more than three seasonal DInSAR stacks, would likely be required to test the efficacy of such models.

3.6.2. Assessing Aggregate Displacement of Multiple Individual Seasonal Maps

Due to the periodic nature of ground displacement expected due to the seasonal thaw and freeze cycle [2], the aggregate values of three stacks do not represent the absolute change in elevation from 2013 to 2019, but do allow for a better characterization of displacement patterns across the entire period. The results indicate broad partitioning of the CBAWO landscape with large aggregate uplifts constrained to upland areas and on slopes, and large aggregate subsidence in low-lying areas (Figure 8). As in the assessment of displacement frequency, significant uplift and subsidence are almost mutually exclusive (only 0.06% of the CBAWO displayed both in the study years). This upland/lowland divide could reflect the early Holocene marine limit, which has been estimated at 35–90 m a.s.l. [78] and 45–90 m a.s.l. [79] for southern Melville Island. Higher levels of subsidence in low-lying areas may be associated with thermokarst development due to more abundant segregation of ice in fine grained marine sediments [10]. The upland areas with a large aggregate uplift map closely to an abundance of mud ejections systematically mapped at the CBAWO in 2012-13 [80]. These ejections of slurry from the subsurface are hypothesized to be closely associated with high pore water pressure, which may imply poor drainage and, as a consequence, cause ice growth during freeze-up, which may increase the potential for detecting uplift in these upland areas. Additionally, large amounts of aggregate uplift were detected on north-northwest slopes. These areas, associated with a low potential incoming solar radiation, have previously been speculated to accumulate shallow ice [81] in several locations in the High Arctic.

5. Conclusions

We have demonstrated that DInSAR stacks using C-band Radarsat-2 data can effectively map vertical surface change in a High Arctic environment, capturing surface subsidence associated with thaw in the early summer, and uplift, associated with frost heave in the late season. Displacement maps across a 48-day period were generated for late summer 2013, 2015, and 2019 and an additional early season stack was generated for July–August 2019. Displacements of up to 10 cm were detected, but the vast majority of pixels displayed far smaller displacements (±2 cm) per 48-day stack. Summer 2019, characterized by moderate TDD but a high cumulative precipitation and correspondingly wet soil moisture conditions, showed both the most extensive and the highest magnitude average displacement. Correspondingly, 2015, characterized by low cumulative precipitation, but an early onset of thaw conditions and high TDD, was relatively inactive. This points to the importance of precipitation and soil condition in processes which result in surface displacement in any given summer, although the late onset of freezing conditions in 2015 at the end of the stack interval meant that only a small fraction of potential frost heave in this year would have been captured by DInSAR. Hence, additional stacks would be required to examine this relationship more thoroughly.

As expected, DInSAR captured less uplift (associated with frost heave) in the early (9 July–26 August) stack compared to the later stack (2 August–19 September) (5.5% of total area showing >1 cm uplift versus 11.4%). Areas which displayed large absolute displacement (>2.5 cm) in all three seasons were sparse, but highly localized around streams, on steep slopes and ridges, and in coastal areas. This indicates a landscape with highly variable isotropic displacement, which may make it challenging to predict via susceptibility modelling, but with a number of highly active hotspots. Downslope transport likely also contributes to the regular nature of displacement activities in many of these areas. When all three years are taken together, large aggregate subsidence was exclusive to low-lying areas, and large aggregate uplift was found exclusively in upland areas. This upland/lowland transition occurred through 50–80 m a.s.l., which is consistent with previous estimates of the early Holocene marine limit in southern Melville Island, and thus may be associated with differing sedimentary depositions (each with different associated frost susceptibilities and depths of ice formation). Mean surface displacement magnitude was found to be significantly greater within previously mapped landscape disturbances such as ALDs across all years. DInSAR also captured cohesive and large subsidence within the slump floors of two medium-sized retrogressive thaw slumps.

Overall, this research has taken a multi-faceted approach to demonstrate the utility of DInSAR in monitoring landscape change in a High Arctic environment. Although the two simple metrics we developed to characterize change across all three years provided some insight, more frequent SAR acquisitions with fewer data gaps between seasons will be required to: (i) more extensively examine surface displacement in relation to environmental and topographic variables; and (ii) to understand long-term landscape change and to the timing of displacement each summer with relation to the onset and intensity of thawing and freezing conditions.