Vulnerable top-of-permafrost ground ice indicated by remotely sensed late-season subsidence

Ground ice is foundational to the integrity of Arctic ecosystems and infrastructure. However, we lack fine-scale ground ice maps across almost the entire Arctic, chiefly because ground ice cannot be observed directly from space. Focusing on northwestern Alaska, we assess the suitability of late-season subsidence from Sentinel-1 satellite observations as a direct indicator of vulnerable excess ground ice at the top of permafrost. The idea is that, towards the end of an exceptionally warm summer, the thaw front can penetrate materials that were previously perennially frozen, triggering increased subsidence if 5 they are ice rich. For locations independently determined to be ice rich, the late-season subsidence in an exceptionally warm summer was 4–8 cm (5th–95th percentile), while it was lower for ice-poor areas (-1–2 cm). The distributions overlapped by 2%, demonstrating high sensitivity and specificity for identifying top-of-permafrost excess ground ice. The strengths of late-season subsidence include the ease of automation and its applicability to areas that lack conspicuous manifestations of ground ice, as often occurs on hillslopes. One limitation is that it is not sensitive to excess ground ice below the thaw front 10 and thus the total ice content. Late-season subsidence can enhance the automated mapping of vulnerable permafrost ground ice, complementing existing (predominantly non-automated) approaches based on largely indirect associations with vegetation cover and periglacial landforms. Improved ground ice maps will prove indispensable for anticipating terrain instability in the Arctic and sustainably stewarding its ecosystems.

To estimate surface displacements at a resolution of 60 m, we used Sentinel-1 observations between early June and mid-September 2017-2019. The Sentinel-1 observations (Torres et al., 2012), acquired at a frequency of 5 GHz with VV polarization 90 in the Interferometric Wide mode, were available at 12-day intervals from an easterly look direction (37 degrees incidence angle; path 15, frame 367). There was one exception in 2017, during which a gap of 18 days occurred.

Estimating subseasonal subsidence time series
To estimate a subseasonal displacement d time series from the Sentinel-1 observations Copernicus Sentinel (2020), we applied Short BAseline Subset (SBAS) processing (Berardino et al., 2002). The rationale of this Interferometric SAR (InSAR) 95 technique is to derive displacement time series from redundant interferograms. We first formed interferograms with temporal separations of up to 24 days, using spectral diversity techniques for coregistration (Scheiber and Moreira, 2000) and removing the topographic phase contribution using the TanDEM-X DEM (Rizzoli et al., 2017;TanDEM-X DLR, 2020). After multilooking, we unwrapped the interferograms using SNAPHU (Statistical-Cost, Network-Flow Algorithm for Phase Unwrapping; Chen and Zebker (2001)). Ionospheric phase corrections were deemed to be unnecessary. We then estimated the displacement 100 time series, with a temporal sampling of 12 days, from the interferogram stack. We used a weighted least squares approach based on the singular value decomposition (Berardino et al., 2002), with the weights determined by the Cram'er-Rao phase variance estimate (Tough et al., 1995).
The time series are reported as displacements d along the line-of-sight direction, with positive values corresponding to increasing distance. Owing to the hilly terrain, we chose not to convert them to vertical displacements. However, we did 105 not discern aspect-dependent trends that are associated with downslope movements. For simplicity's sake, we will refer to displacements with increasing distance as subsidence. To emphasize the late-season subsidence, we set the zero point of the relative subsidence time series to be at the beginning of the late season t 1 .

Referencing and assessing its quality
To reduce long-wavelength atmospheric errors in the subsidence observations, we spatially referenced the raw time series at 110 multiple locations with outcropping bedrock or a thin rubble veneer. These locations were assumed to be stable (Reger and Solie, 2008;Antonova et al., 2018).
To assess the quality of the referencing, we chose eight bedrock validation points distributed across the study region. Assuming these points to remain stable, the subsidence observations at these locations are a measure of the observational uncertainty . This uncertainty estimate subsumes decorrelation errors and uncompensated atmospheric contributions 115 under a single number.

Late-season subsidence from spline fitting
We estimated the late-season subsidence d l by fitting spline basis functions to the referenced subsidence time series. The advantage of fitting a flexible and yet simple spline function is that measurement noise, such as residual atmospheric errors, can be reduced (Berardino et al., 2002). To capture a range of subseasonal subsidence patterns, we used three cardinal quadratic 120 B-splines for the subsidence rate (first derivative), corresponding to the cubic spline basis functions for the subsidence time series shown in Fig. 3a).  The late-season subsidence d l was defined to be the cumulative subsidence between 10 August (t 1 ) and 10 September (t 2 ) for all years, the end point having been chosen to minimize the confounding impact of diurnal frost heave (Chen et al., 2020).
We derived d l from the spline fit. 125 We also estimated total subsidence during the early and mid thaw season (10 June to 10 August), d e . While d e is not expected to contain direct information about permafrost ground ice, being sensitive to soil texture and interannual variations in active layer moisture (Lewkowicz, 1992;Harris et al., 2011), it provides a simple reference for assessing the late-season speed-up characteristic of melting top-of-permafrost excess ground ice.
3.2 Assessment with independent ground ice data 130

Manual mapping
We assessed the suitability of late-season subsidence as an indicator of permafrost ground ice content by comparing the satellite observations to an independently derived ground ice map. The map comprised two primary categories: ice rich and ice poor.
The sensitivity and specificity of late-season subsidence were assessed from the observed distributions conditional on the location being ice rich and ice poor, respectively.

135
The independently derived ground ice map was obtained using manual interpretation and expert knowledge, drawing on field observations and high-resolution (∼1 m) satellite imagery. The mapped focus area, 8 km by 8 km in size, was chosen because of the wide range of ecotypes and the availability of field observations. A drawback of the map are its gaps, as ice-rich permafrost cannot be identified reliably in areas that lack unambiguous indicators. Furthermore, six percent of the area were discarded in the comparison to avoid unrepresentative values over lakes and infrastructure.

140
The key considerations in the manual mapping were: 1. Ice-rich permafrost was assigned to ecotypes where high-resolution imagery and field observations revealed direct indicators of excess ice at the top of permafrost. Ice wedge polygons -some only barely visible or inconspicuous except along lakeshores and beaded streams, others in an advanced state of degradation -are widespread (Fig. 4). They are abundant in old alluvial and lacustrine deposits (DOWL Engineers, 1994;Shannon & Wilson, Inc., 2006), as well as in 145 colluvial sediments that are rich in retransported silt (DOWL Engineers, 1994). Pingos in drained lake basins are also direct indicators of excess ground ice (Mackay, 1973). 2. Exposed bedrock and rubble-covered surfaces (Fig. 4b) were classified as ice poor. This classification is generally supported by geotechnical investigations and field evidence from the study area (Pewe et al., 1958;DOWL Engineers, 1994).
Deviations from this general pattern cannot be excluded (Robinson and Pollard, 1998;French and Shur, 2010). The mapping process is inherently subjective. The discretization of top-of-permafrost excess ground ice into a limited number of categories present challenges (Tryck Nyman Hayes, 2006;Paul et al., 2020). The most difficult decisions were about where to draw the line between the ice-rich/ice-poor and the indeterminate category.

165
We also compared the late-season subsidence observations to geotechnical assessments of permafrost ground ice content from Although coring is the most reliable method for determining ground ice content, the spatial and temporal representativeness of the cores need to be considered. The point observation a single core represents could paint a misleading picture of ground ice content (Morse et al., 2009;Kanevskiy et al., 2017) when compared to a ∼60 m resolution cell. The cores were further taken a decade before the remote sensing observations. However, none of the locations showed signs of severe disturbance, and We classified the cores as ice poor and ice rich based on the descriptions in the geotechnical report (Shannon & Wilson, Inc., 2006). The ground ice content of each core was summarized verbally as ice rich or ice poor, complemented by pictures and estimates of the visible ice content. All cores but one were ice rich with visible ice contents >30% in the form of segregated 180 and massive ice. The only ice-poor core was a gravelly soil grading into weathered bedrock (5% visible ice content in joints).
Two cores taken from within the same ice wedge polygon were combined for the comparison exercise. Tab. S1 summarizes the locations and ground ice properties of all 13 coring locations.

Results
4.1 Spatial and temporal variability of late-season subsidence 185 4.1.1 Regional variability The observed late-season subsidence in our study area showed two distinct peaks in 2019 (Fig. 5a). One corresponded to no or very small subsidence (-1 -1 cm). The other corresponded to regions with elevated subsidence (∼5 cm) in 2019. In 2017, the subsidence in these regions was lower (1-4 cm). There were no notable instances of pronounced negative late-season subsidence. Late-season subsidence was consistently low in elevated uplands (e.g. bedrock or talus without vegetation cover) and along rivers (recent floodplains without or with dense vegetation cover; Fig. 5a-b). This association supports the assumption that the reference and validation points on bedrock (circles in Fig. 5) by circles, could be assumed to be stable. The early-mid and late-season subsidence observations at the validation points, shown in Fig. 3b) had a root-mean-square deviation of 0.6 cm.
Other low-lying areas and many hillslopes exhibited elevated late-season subsidence, especially in the exceptionally warm

Temporal variability
The inter-annual variability of the observed subsidence time series is illustrated by three examples in Fig. 6. The late-season subsidence in 2019 exceeded 5 cm in the ice-rich (ice wedges) deposits shown in panels a and b, respectively. According to the 200 radar observations, there was an acceleration in the rate of subsidence in the late thaw season, as the rate more than doubled.
The rate and the total subsidence during the late season is approximately a factor of three smaller in the years 2017 and 2018.
For the floodplain in Fig. 6c, the interannual variability was less than 1 cm. The late-season subsidence was of a magnitude comparable to the observational uncertainty in all years at this ice-poor location.

Suitability in an exceptionally warm summer
Late-season subsidence in the hot summer of 2019 was markedly different for ice-rich and ice-poor areas (Fig. 7a). It exceeded 4 cm at all cored locations rich in top-of-permafrost ground ice. For the ice-rich areas, as determined by independent manual mapping (Fig. 7c), d l varied between 4 and 8 cm (5th and 95th percentile; Fig. 7b. For for ice-poor areas the range was -1 to 2 cm. Only in the extreme tails (2%) do the distributions overlap, ensuring a robust separability of the two classes.  The triangles indicate the location of the ice cores (cf. Fig. 5): the one on the bench was ice poor, the others ice rich. Ice wedge polygons were observed in the field ∼300 m upslope from the bench.
The separability based on the 2019 late-season subsidence d l is better than that based on the early-mid-season subsidence d e . The d e distributions of ice-rich and ice-poor areas overlap (Fig. 8c), whereas the d l distributions are concentrated around two separate peaks (see also Fig. 7a).
The candidate relocation site Tatchim Isua was characterized by a narrow zone with low late-season subsidence (Fig. 9).
This ∼100 m wide zone largely coincides with the gravel-covered bench that a single core from 2005 indicates to be ice Further examples from a range of geologic settings serve to illustrate the suitability of d l for identifying ice-rich permafrost. 220 Fig. 10a-d shows instances of ice-rich permafrost with ice wedge polygons. They all exhibited elevated d l of 4-8 cm, corresponding to an increased subsidence rate during the late season. Conversely, the observed subsidence was below 1 cm for the points shown in Fig. 10i-l, which were independently determined to be poor in ground ice because of their young age, such as active and inactive floodplains, or because the formation of excess ice is impeded by their composition.
The most interesting cases are those areas where the manual ground ice mapping was indeterminate because there was 225 no strong evidence for either category. They exhibited a bimodal distribution of d l (Fig. 7a). The larger mode d l ∼ 5 cm, comparable to those of ice-rich areas. Examples include colluvial hillslope deposits without conspicuous ground ice indicators or ecosystem-driven perennial ground ice near the surface (Jorgenson et al., 2010;Kanevskiy et al., 2012Kanevskiy et al., , 2014Paul et al., 2020), but see Burn (1997). Ice-rich permafrost may occur at depth, perhaps under a thick talik. The thermal regime of such 290 permafrost is more complex (Jorgenson et al., 2010;Connon et al., 2018), and melt-induced subsidence does not necessarily occur primarily at the end of the thaw period. The sensitivity may further be diminished by sinkhole formation and piping underneath cohesive materials such as peat (Osterkamp et al., 2000); or by the retardation of thaw consolidation due to inefficient drainage of the excess melt water (Morgenstern and Nixon, 1971).
The specificity is impaired when unrelated processes induce late-season subsidence. Gradual subsidence due to processes 295 such as organic layer degradation (Stephens et al., 1984) may be distinguished from the ground ice signals shown in Fig. 6a-b by the late-season acceleration typical for ice-rich permafrost. Elevated late-season subsidence may also reflect ice content at the bottom of the active layer. Cold permafrost and ample moisture supply promote ice segregation at the base of the active layer, which can be continuous with the intermediate layer (Mackay, 1981). Deformation related to precipitation events and lateral flow may pose additional challenges, in particular in peatlands and on slopes (Roulet, 1991;Matsuoka, 2001;Gruber, 300 2020; Zhang et al., 2020). In hilly terrain it will be advantageous to resolve the downslope and surface-normal movement components, but we currently lack adequate satellite observations to do so routinely. The limited downslope movements in our study area likely contributed to the high specificity of late-season subsidence for mapping ice-rich permafrost.
A final limitation of this study and the preceding discussion is that the complexity of ground ice was simplified to just two categories: ice rich and ice poor. In reality, however, excess ground ice content is a continuous parameter (Morse et al., 2009;305 Kanevskiy et al., 2012; Paul et al., 2020) whose magnitude could be constrained using late-season subsidence observations.
A quantitative assessment at the appropriate observational scale will require densely sampled ground ice cores (Morse et al., 2009;Jorgenson et al., 2010).

Enhancing automated ground ice mapping
Late-season subsidence can enhance the automated mapping of vulnerable permafrost ground ice. Remotely sensed late-season 310 subsidence can be mapped on pan-Arctic scales, thanks to the global availability of Sentinel-1 data. The mapping can be automated, as no manual interpretation an no calibration using in-situ cores are required.
Late-season subsidence is complementary to state-of-the-art mapping approaches based on visible manifestations of ground ice and indirect associations with vegetation cover or topographic variables. The greatest contribution of late-season subsidence will likely be where field observations are sparse or where ground ice is indistinct and poorly correlated with surface charac-315 teristics. Examples includes uplands and hillslopes (Fig. 9), as well as areas underlain by ice wedges from the early Holocene or the Pleistocene (Pewe et al., 1958;Dredge et al., 1999;Reger and Solie, 2008;Morse et al., 2009;Jorgenson et al., 2015;Farquharson et al., 2016).
Incorporating geological constraints can counteract weaknesses of remotely sensed late-season subsidence. These include larger errors in vegetated areas (Wang et al., 2020) or false negatives when the thaw front fails to penetrate deep into the ice-rich 320 materials (Shur et al., 2005). Most importantly, geological and other observational constraints will be needed to estimate total excess ice contents.