Seasonal and Interannual Ground-Surface Displacement in Intact and Disturbed Tundra along the Dalton Highway on the North Slope, Alaska

Spatiotemporal variation in ground-surface displacement caused by ground freeze–thaw and thermokarst is critical information to understand changes in the permafrost ecosystem. Measurement of ground displacement, especially in the disturbed ground underlain by ice-rich permafrost, is important to estimate the rate of permafrost and carbon loss. We conducted high-precision global navigation satellite system (GNSS) positioning surveys to measure the surface displacements of tundra in northern Alaska, together with maximum thaw depth (TD) and surface moisture measurements from 2017 to 2019. The measurements were performed along two to three 60–200 m transects per site with 1–5 m intervals at the three areas. The average seasonal thaw settlement (STS) at intact tundra sites ranged 5.8–14.3 cm with a standard deviation range of 2.1–3.3 cm. At the disturbed locations, averages and variations in STS and the maximum thaw depth were largest in all observed years and among all sites. The largest seasonal and interannual subsidence (44 and 56 cm/year, respectively) were recorded at points near troughs of degraded ice-wedge polygons or thermokarst lakes. Weak or moderate correlation between STS and TD found at the intact sites became obscure as the thermokarst disturbance progressed, leading to higher uncertainty in the prediction of TD from STS.


Introduction
Frost heave is the upward or outward movement of the ground surface caused by the formation of ice lenses in porous materials (e.g., [1]). This is a ubiquitous phenomenon in cold regions, where the ground surface experiences freezing. The fundamental mechanism of frost heave has been recognized since the 1920s as the formation of segregated ice lenses by water migration to the freezing front in the frost-susceptible soil (e.g., [2,3]). The heaved ground surface, in turn, settles during thawing seasons, and the seasonal thaw settlement (STS) depends on the amount of ice lenses formed in the previous years and the maximum thaw depth (TD) of the observation year. Frost heave redistributes soil water and displaces components of the soil layer (e.g., [4]). Researchers and practitioners of civil engineering have studied the mechanism of frost heave and modeled its behavior (e.g., [5][6][7][8]) because the surface displacement damages infrastructure, such as roads, railroads, and buildings, through repeated seasonal heave and settlement. Ice segregation is a fundamental driving force to develop geomorphological features such as palsas [9] and lithalsas (e.g., [10]) in permafrost regions. The physical processes underlying the phenomenon of frost heave have been studied by a number of researchers, as summarized by Dash et al. [11].
Ice segregation in the active layer and the resulting differential seasonal ground displacement through frost heave and thaw settlement produce highly heterogeneous soil structures and microrelief in Arctic ground surfaces underlain by permafrost. At the Land 2021, 10, 22 2 of 18 same time, ice segregation damages fine plant roots and seedlings (e.g., [12,13]), alters vertical moisture distribution [14,15] in freezing soil layers, and restricts water accessibility for aboveground vegetation and belowground microbial activity at different soil depths. Quantifying these key processes is critical to understanding cold-region ecosystems that are experiencing rapid climate change.
The development of automatic measurement systems with sensors (e.g., strain gauge and potentiometer) and loggers has provided continuous and high-resolution measurements of ground-surface displacement (e.g., [16,17]). Gruber [18] introduced a low-cost and simple-structured "tilt arm" measurement system. Although these automatic systems can provide continuous measurements, there still remains a problem of limited spatial extent of the measurements.
Little et al. [19] demonstrated the ability of a differential global positioning system (DGPS) to measure ground-surface displacement caused by frost heave and thaw settlement in wide areas more rapidly than classical methods. A similar DGPS displacement measurement in Alaskan tundra was continued by Shiklomanov et al. [20] and Streletskiy et al. [21]. The recent technological development of RTK (real-time kinematic)/PPK (postprocessed kinematic) GNSS (global navigation satellite system) positioning enabled much more rapid measurements with a vertical precision of 10-20 mm. However, these measurements are also episodic and require the physical presence of surveyors.
As a solution to the limited spatial extent of in situ displacement measurement, a radar remote-sensing technique has been deployed to estimate the spatial variation of active layer thickness in broad areas. For example, Liu et al. [22], Schaefer et al. [23], and Michaelides et al. [24] estimated spatial variation in the active layer in Alaskan permafrost regions using the SAR (synthetic aperture radar) interferometry technique (InSAR). In developing the algorithm to determine active layer thickness from remotely sensed variables, they presumed a relationship between TD and STS because greater TD results in a larger water capacity that potentiates larger frost heave and STS. However, supporting direct field validation for the remote-sensing algorithm is still missing, due to the scarcity of studies on the TD-STS relationship.
Furthermore, recent InSAR studies applied to permafrost regions have demonstrated the ability to capture high-resolution spatial variation in seasonal and interannual ground displacement caused by ground freezing and thawing, for example, in Alaska (e.g., [25,26]), Siberia [27][28][29], and Himalayan highland [30]. Few studies conducted comparisons between InSAR measurements and in situ ground survey; however, the numbers of measuring points and observation periods were limited by various logistical reasons.
Our study focuses on the knowledge gap between the high-resolution spatial distribution of ground-surface displacement and remote sensing in permafrost regions. We conducted a three year in situ measurement of ground-surface displacement on the North Slope, Alaska, deploying a GNSS positioning system. The objectives of our field campaign were (1) to capture spatial and temporal variations in ground-surface displacement in an ice-rich permafrost region in a broader and denser manner than previously reported, (2) to better understand the magnitude of disturbance impact on the freeze-thaw-related ground-surface displacement, and (3) to examine the relationship between TD and STS that was presumed in past remote-sensing works without field measurements.

Study Sites
The study sites are located about 300-500 m west of the Dalton Highway, midway between the Brooks Range and Deadhorse ( Figure 1). All study sites are located in Arctic bioclimate subzone E covered with dominant vegetation unit 8 (tussock sedge, dwarf shrub, and moss tundra) according to the CAVM (Circumpolar Arctic Vegetation Map) Team [31] and Walker et al. [32]. No tall shrubs (>0.5 m) were found, and no patterned ground was discernable in the areas of our survey transects, except for faint ice-wedge polygons at site IC. Land cover consists of acidic dry tussock-sedge tundra developed on silty loam to silty clay loam soils (e.g., [33,34]) with occasional gravels or cobbles. According to Alaska Paleo-Glacier Atlas Version 2 [35], site SG (Sagwon) is situated a few kilometers out of the maximum glaciation range, while sites HV (Happy Valley), and IC (Ice Cut) are located in the area between the maximum glaciation and early Wisconsin glaciation boundaries. The landscape of these study sites is gently rolling hills of highly ice-rich permafrost sediment (Yedoma Ice complex) with an elevation range of 300-400 m.
Site SG (Figure 2a,b, Figures A1 and A2) is located on a relatively flat summit area of a Yedoma remnant hill. A portion of site SG experienced surface disturbance between 1995 and 2007, according to available high-resolution images, and thermokarst development is ongoing. Transects SG_T1 and SG_T3 intersect the disturbed and undisturbed tundra, while the entire transect SG_T2 is situated in intact tundra. We subdivided site SG into SG_disturbed and SG_intact for the later descriptions and discussion. Site HV (Figure 2c Small tussocks with a height up to 20 cm are sporadically distributed at these sites, and 12%, 28%, and 16% of the total survey stakes were installed on the tussocks at SG, HV, and IC, respectively. Mean annual air temperatures at the NOAA (National Oceanic and Atmospheric Administration) Sagwon station (GHCND: USS0048U01S) ranged from −6 to −9 • C and the mean annual precipitation was 238 mm in the period 2002-2018 [36]. Daily mean air temperature ranged from −40 to 23 • C ( Figure 3). Other meteorological statistics for our study period using daily data at the Sagwon station [37] are summarized in Table 1. For the calculation of TDD Menne thawing degree days) and FDD (freezing degree days), we used periods of June-September and the previous October-May for the respective year.   We conducted an RTK (real-time kinematic)/PPK (postprocessed kinematic) GNSS (global navigation satellite system) positioning survey to measure seasonal and interannual ground-surface displacement at targeted tundra sites (Figure 2). At each site, we installed transect lines permanently marked by surveyor stakes at intervals of 1 or 5 m. The measurements were repeatedly performed at the marked ground surfaces at the beginning and end of thawing season from 2017 to 2019. We used a Trimble R9s GNSS receiver as a base station and the R2 receiver placed on the top of a carbon-fiber survey rod with a flat topo shoe (Seco, 5192-02) as a rover. The nominal vertical accuracies of static and RTK/PPK GNSS positioning were 5 mm + 0.5 ppm RMS (root mean square) and 20 mm + 1 ppm RMS, respectively (Trimble Inc., USA).
The height of the measuring ground surface was defined as the height at which the rover rod's bottom settled by the entire rover weight. To avoid ambiguity in selecting the ground surface around a survey stake by different surveyors at timings, we installed the survey stakes at the center of flat tundra areas with an approximate diameter larger than 10 cm. The GNSS positioning was performed for more than 10 s at each measurement point. The measured points with vertical precision larger than 15 mm were discarded and remeasured when the RTK survey was performed. Postprocessing of PPK or static measurements was carried out using the software Trimble Business Center. Postprocessed points with low vertical precision (larger than 15 mm) were removed from further analysis during postprocessing.
A permanent stake was installed in the permafrost layer at least 1 m deep at each site and repeatedly used as a reference point. While we can safely assume that the reference height was stable during a thawing season, frost jacking of the stakes was sometimes unavoidable. To minimize the error in vertical height of the reference rods raised by frost jacking, we used individual postprocessed base locations measured every year. The base locations were measured for at least 24 h to obtain a vertical measurement precision of static GNSS positioning less than 20 mm. Combining the measurement errors of the rover and base, we estimated the overall GNSS positioning accuracies to be less than 21 and 29 mm for STS and interannual displacement, respectively. STS was determined by the difference in the ground-surface heights measured at the beginning and end of thaw season each year. Interannual surface displacement for the individual interval was determined according to the height difference between measurements of the ground-surface heights for the corresponding time interval.

Thaw Depth and Surface Moisture
We measured TD and surface moisture at the end of thawing seasons nearby the survey points within a few days of the GNSS survey. TD was determined by pounding a steel rod into the ground until frost was encountered. As an indicator of volumetric water content in the ground surface, the travel time of electromagnetic waves (TDR period) at depths of 6, 12, and 20 cm was measured by HydroSence II (Campbell Scientific Logan).

Statistical Analysis
Statistical testing and calculations were done using R version 3.6.1 software [38] with the "corrplot" package [39].

Variation in the Maximum Thaw Depth
For each study site, the maximum thaw depth (TD) measurements with confidential intervals are summarized in Table 2 and Figure 3. The mean TDs at SG_disturbed were the largest (62-73 cm) throughout the study period, and were about 1.5 times larger than the values at other intact sites. The TDs at SG_intact and HV were in a similar range (41-48 cm), while those at IC had an intermediate value ( Figure 3, Table 2). Our observed TD ranges at the intact locations were comparable to the TD (about 50 cm) measured at a site a few hundred meters from our sites by Raynolds et al. [33] and Walker et al. [40]. Variations in TD within the areas were largest at SG_disturbed (Figure 4). The standard deviations of TD at SG_disturbed were markedly higher (20-22 cm) than those at IC and HV (9 cm), reflecting the influence of past surface disturbance (Table 2). Although mean TD at SG_intact was as shallow as that at HV, the larger variations in TD compared to other intact sites indicate the nearby influence of surface disturbance at SG_disturbed, and a potential increase of TD in the near future.  The overall trend of interannual change in TD was an increase at all sites from 2017 to 2019. There were monotonic increases in mean TD at SG_disturbed and IC, while the TD at SG_intact and HV showed stagnation in 2018. The largest TD in 2019 was mainly attributed to the warmest summer conditions and the thicker snow depth in the previous winter (Table 1). The development of thermokarst may also cause an increase in TD. However, analysis of the thermokarst processes at the individual site was beyond the scope of this study.
Our TD measurements were conducted at the end of the thawing season each year. TDD 1/2 values on the days of TD measurements were 29, 26, and 31, in 2017, 2018, and 2019, respectively. TDD 1/2 is often used as an analytical estimate of the progress of seasonal thaw depth in permafrost regions, according to the simplified Stefan equation (e.g., [41][42][43]) as well as for the interannual comparison of active layer thickness (e.g., [20,21]). TDD 1/2 values on our measurement dates were 94%, 96%, and 97% of the annual maximum values, respectively ( Figure 3). Therefore, we regarded our measured TD as maximum seasonal values for an interannual comparison.

Variation in Surface Displacement
Our GNSS surveys captured the variation in seasonal thaw settlement (STS) during thawing seasons and frost heave during freezing seasons, together with interannual ground-surface displacement at each site (  (Table 2). On the other hand, the ground surface showed monotonic interannual subsidence at SG. We observed that fall season measurements had a much larger average total subsidence at SG_disturbed (18.2 cm) than at SG_intact (11.6 cm). It is highly probable that the thermokarst process is ongoing at SG due to the past surface disturbance, resulting in monotonic interannual subsidence.
Surface disturbance and the thermokarst that followed enhanced the magnitude and spatial variation of ground-surface displacement due to the freeze-thaw cycle in the active layer. Large seasonal and interannual displacement (up to 44 and 56 cm for seasonal and interannual, respectively) was observed in concave areas (arrows in Figure 5) due to thermokarst subsidence in the disturbed zones along transects SG_T1 and SG_T3 (Figure 2a). Standard deviations of STS at SG_disturbed (6.5-9.3 cm) were about three times larger than values at SG_intact, HV, and IC (2.1-3.3 cm) ( Table 2).
Characteristics of the fine-scale spatial variation in ground-surface displacement tended to be preserved over the years (Figures 5-7). For example, the points with larger STS and interannual subsidence in a certain year tended to have relatively larger values regardless of the observed year. The dependency of the amount of ground-surface displacement on the fine-scale locations obviously indicates the importance of microtopography and the local difference in soil properties, as demonstrated by Walker et al. [45], Romanovsky et al. [46], and Walker et al. [40].

Correlation between Seasonal Thaw Settlement and Maximum Thaw Depth
Pearson's positive correlations between STS and TD were found to be statistically significant in all groups (Figure 8). Moderate correlations (R = 0.41-0.45, p < 0.001) were found at IC, HV, and SG_intact, and low correlation (R = 0.37, p < 0.001) was found at SG_disturbed. An analysis of covariance (ANCOVA) showed insignificant differences for slopes at the location, but significant (p < 0.001) differences in intercepts of the regression curves; i.e., the differences in average TD among the sites. ANCOVA also showed an insignificant interaction between the presence of tussocks and STS. Our observations showed that the hypothesis of a positive relationship between STS and TD generally holds, but features high uncertainty in the use of point prediction. Despite the significant moderate correlation between STS and TD at intact sites, the uncertainty in prediction (width of the 95% confidence interval) using linear regression curves was high (at least ±17 cm at HV). As expected from the low correlation at SG_disturbed, the uncertainty for this site was the largest (±42 cm).
The correlation between STS and TD may be useful for ensembles of measurements to predict from one parameter to the other. At the same time, the uncertainty range over 34 cm for the TD prediction from a single STS observation demonstrates the difficulty in the use of regression curves for individual point measurements. A possible quantitative explanation for the differences in the prediction uncertainty ranges among the intact locations is the degree of thermokarst or proximity to disturbed tundra. HV was located farthest from the Dalton Highway, and there were no visible thermokarst depressions at or near the site. SG_intact was next to the SG_disturbed locations, and IC could be treated as intermediately affected by thermokarst because slight thermokarst depressions were recognizable near the site, and the site was close to the highway.  Overall, as the thermokarst disturbance progressed, the correlation between STS and TD became obscure, and the uncertainty in the prediction of TD from STS became larger. The lower correlation between STS and TD in thermokarst-affected locations was reasonable. Assuming similar weather conditions and frost susceptibility of the soils in the study area, differences in soil moisture and freezing rate play important roles in the amount of seasonal ground-surface displacement (frost heave and thaw settlement) (e.g., [11]). While an increase in TD in thermokarst-affected locations expands the potential STS range, larger spatial (horizontal and vertical) variations in soil moisture and altered thermal regime in the active layer increase the range in frost heave amount compared to the intact tundra with its relatively uniform soil properties. Our observations suggest that thermokarst, by melting segregation ice and massive ground ice in the transient (shielding) layer of permafrost, shifted hydrology in the active layer (e.g., lower water table and enhanced drainage) and altered the correlation between STS and TD in the intact tundra.

Surface Moisture Distribution and Ground-Surface Displacement
Low to moderate correlations (significance level of 99.9%) between STS and surface moisture at a depth of 12 cm were found at intact tundra sites ( Figure 9). However, moisture at the shallower (6 cm) and deeper (20 cm) levels was uncorrelated with STS in our study. In the active layer of the tundra with shallow TD (about 50 cm), water drainage from the soil layer close to the frost table is usually hindered by the relatively flat topography and underlying frozen ground. In comparison, the top few centimeters of the active layer typically consist of highly permeable live plants, moss, or peat, which rapidly drains precipitation water downward. As a result, the surface moisture in the upper and lower active layer tends to be less variable in space. The negative correlation between STS and surface moisture at IC, HV, and SG_intact was counterintuitive because we expected larger frost heave where soil moisture is higher, which is essential for the formation of ice lenses upon freeze-up. To explain this negative correlation, we may need to compare soil texture at a fine spatial resolution; however, these ancillary spatial data were not available for this study. On the other hand, no significant correlation was shown at SG_disturbed between STS and surface moisture (only measurements at a depth of 12 cm were available). There was also no significant correlation between TD and surface moisture at all sites.

Conclusions
The spatial variation in ground-surface displacement at intact tundra sites without marked patterned ground features and a disturbed site was characterized using repeated precise GNSS positioning surveys. The average seasonal thaw settlement (STS) at intact tundra sites (HV and IC) ranged 5.8-7.2 cm and 10.2-14.3 cm in average climate summers (2017 and 2018) and a warmer summer (2019), respectively. Standard deviations were relatively stable (2.1-3.3 cm) over the observation years. At the disturbed site (SG_disturbed), STS was highest in all observed years and among all sites (12.6, 8.6, and 19.9 cm in 2017, 2018, and 2019, respectively). The spatial variation in ground-surface displacement among the observed points roughly persisted during our observation period. There were few variations in meteorological and ground-surface conditions, including vegetation cover, among the study sites. Hence, the observed variations in ground-surface displacement were presumably governed by in situ physical properties and conditions of the active layer.
Low to moderate positive Pearson's correlations between maximum thaw depth (TD) and STS were found in intact areas. On the other hand, a lower and nonsignificant correlation was found at SG_disturbed, where ongoing thermokarst was prominent. The lower correlation suggests shifts of hydrology in the active layer (e.g., lower water table and enhanced drainage) and increased contribution of thermokarst subsidence by melting segregation ice and massive ground ice in the transient (shielding) layer of permafrost. However, soil moisture at a depth of 12 cm negatively correlated with STS at all sites, and no correlations were found at other depths.
The TD at SG_disturbed was greatest (62-73 cm) throughout the study period and about 1.5 times larger than that at other intact sites (41-48 cm). The interannual changes in TD and STS could be explained by the weather conditions (snow depth, freeze/thaw indices) at the intact locations. On the other hand, the interannual changes at SG_disturbed were less dependent on atmospheric forcing, indicating significant influence by thermokarst.
Overall, as the thermokarst disturbance progressed, the correlation between STS and TD became obscure, and the uncertainty in the prediction of TD from STS became larger. Our field surveys captured the spatial and temporal variations in freeze-thaw-related ground-surface displacement at intact tundra sites and a disturbed area in a typical tundra landscape underlain by ice-rich permafrost.