Global Positioning System interferometric reflectometry (GPS-IR) measurements of ground surface elevation changes in permafrost areas in northern Canada

Global Positioning System interferometric reflectometry (GPS-IR) is a relatively new technique which uses reflected GPS signals to measure surface elevation changes to study frozen-ground dynamics. At present, more than 200 GPS stations are operating continuously in the Northern Hemisphere permafrost areas, which were originally designed and maintained for tectonic and ionospheric studies. However, only one site in Utqiaġvik, Alaska (formerly Barrow), was assessed to be usable for studying permafrost by GPS-IR. Moreover, GPS-IR has high requirements on the ground surface condition, which needs to be open, flat, and homogeneous. In this study, we screen three major GPS networks in Canada and identify 12 out of 38 stations located in permafrost areas as useful ones where reliable GPS-IR measurements can be obtained. We focus on the five Canadian Active Control System stations and obtain their daily GPS-IR surface elevation changes. We find that the ground surface subsided in Alert, Resolute Bay, and Repulse Bay respectively by 0.61±0.04 cm yr⁻¹ (2012–2018), 0.70±0.02 cm yr⁻¹ (2003–2014), and 0.26±0.05 cm yr⁻¹ (2014–2019). At the other two sites of Baker Lake and Iqaluit, the trends are not statistically significant. The linear trends of deformation were negatively correlated with those of thaw indices in Alert, Resolute Bay, and Repulse Bay. Furthermore, in Resolute Bay, we also find that the end-of-thaw elevations during 2003–2012 were highly negatively correlated with the square root of thaw indices. This study is the first one using multiple GPS stations to study permafrost by GPS-IR. It highlights the multiple useful GPS stations in northern Canada, offering multi-year, continuous, and daily GPS-IR surface deformation, which provides new insights into frozen-ground dynamics at various temporal scales and across a broad region.

1 Introduction

Since the International Polar Year (2007–2009), permafrost has undergone a warming trend globally, with an average increase in ground temperature at or near the depth of zero annual amplitude by 0.29±0.12 ^∘C during 2007–2016 (Biskaborn et al., 2019; Romanovsky et al., 2010; Smith et al., 2010). Warming permafrost causes ground ice melting, active-layer thickening, and the release of previously sequestered carbon (Brown et al., 2000; Trucco et al., 2012). It affects hydrological, geomorphological, and biogeochemical processes (Mackay, 1966; Shur and Jorgenson, 2007; Lantuit and Pollard, 2008; Kokelj and Jorgenson, 2013). Measuring and quantifying permafrost changes are crucial for understanding the dynamics of the active-layer and near-surface permafrost (collectively called as frozen ground in this paper), studying the response of permafrost environments to climate change, and assessing the risk of permafrost changes to infrastructures.

Surface elevation changes can serve as an indicator of frozen-ground changes. The freeze–thaw of frozen ground is associated with the phase transition of soil moisture, leading to a ∼9 % change of ice volume due to the density difference between water and ice. Such volume change in the freeze–thaw cycle causes the ground surface to uplift or subside seasonally. Surface deformation can be measured by either traditional benchmark-based methods or modern geodetic and remote-sensing ones. The traditional methods use vertical tubes or pipes, anchored deep into the permafrost, as datum references of ground surface for repeat leveling surveys (Mackay, 1983). Modern methods include interferometric synthetic aperture radar (InSAR), light detection and ranging (LiDAR), and global navigation satellite system (GNSS) positioning. InSAR has been used to measure and quantify surface subsidence in various permafrost landforms (Liu et al., 2010, 2014, 2015; Chen et al., 2018). However, InSAR suffers from coarse temporal resolutions and interferometric coherence loss. Furthermore, InSAR measurements need reference points where the surface deformation is known or assumed to be zero. LiDAR has been used to construct differential elevation models to investigate surface deformation (Jones et al., 2015). However, LiDAR surveys are usually conducted at annual or multi-annual intervals. GNSS positioning has also been used to measure and quantify surface subsidence and uplift (Little et al., 2003; Shiklomanov et al., 2013; Streletskiy et al., 2017). However, those GNSS surveys are usually conducted at the beginning or end of thaw seasons.

Global Positioning System interferometric reflectometry (GPS-IR) is a technique which uses reflected GPS signals to measure ground surface changes, such as elevation, soil moisture content, and vegetation growth condition (Larson, 2016, 2019). GPS-IR has been successfully used to study frozen-ground dynamics by measuring surface deformation at one station in Utqiaġvik, Alaska (Liu and Larson, 2018; Hu et al., 2018). Compared with the aforementioned modern methods, GPS-IR measurements of surface elevation changes have higher temporal resolutions, usually at daily intervals. Their accuracies are on the order of a few centimeters (typically ∼2 cm). Their spatial coverage is antenna-height dependent, e.g., 1000 m² for a 2 m high antenna. Such spatial coverage fills a gap between regional-scale satellite observations and in situ point measurements. Furthermore, GPS-IR measurements are free of solid-Earth movement, such as glacier isostatic adjustment and plate movement (Liu and Larson, 2018). GPS-IR measurements are converted from the vertical distances between the antenna and the reflecting surface. As both the antenna and the surface are equally affected by solid-Earth movement, GPS-IR measurements can directly reflect frozen-ground dynamics.

However, some limitations exist when using GPS-IR. This technique can only be usable in certain surface conditions. It needs the reflecting surface to be open, nearly flat, and relatively homogeneous (Larson, 2016). More than 200 GPS stations are continuously operating in permafrost areas in Northern Hemisphere (Fig. 1). However, only one station, located in Utqiaġvik, Alaska, has been examined and proved for the use of GPS-IR (Liu and Larson, 2018). The underutilization of the GPS-IR technique to a large number of existing stations motivates us to assess the usability of all of these sites. We chose Canada as our study area due to the public accessibility of GPS data, abundant weather records, and detailed geological surveys. We first design a three-step framework to identify useful stations under the same protocols and to ensure the reliability of measurements. We then screen all of the major public GPS networks in Canada to identify the usable stations. We then estimate surface deformation by GPS-IR and use them in turn to study frozen-ground dynamics.

Figure 1 Locations of continuously operating and open-data GPS stations in the permafrost areas north of 50^∘ N. The permafrost zonation, represented by various colors, is based on Brown et al. (1997).

The significance of this study relies on that it first can provide usable GPS sites in permafrost areas for GPS-IR studies, which are complementary to the existing permafrost monitoring programs such as the Circumpolar Active Layer Monitoring (CALM) program and Global Terrestrial Network for Permafrost (GTN-P). Permafrost changes have large spatial heterogeneity due to location, topography, precipitation, and vegetation. Despite the significant increase in the number of in situ sites in the past decades, the CALM and GTN-P sites are still sparse and unevenly distributed. The identified suitable GPS stations can fill in the spatial gaps of the CALM and GTN-P sites. Moreover, GPS-IR measurements are typically continuous and span multiple years. They can be used to study permafrost in a detailed manner and provide new insights into frozen-ground dynamics.

In Sect. 2, we describe the mechanism of GPS-IR to measure surface elevation changes, our proposed framework to identify useful GPS stations, and the datasets we used in this study. In Sect. 3, we present basic information of the identified useful GPS stations (e.g., monument material, foundation depth, and antenna height above ground) and environment conditions of the study sites, such as biome and surficial material. We then show the results of surface elevation changes in thaw seasons at the study sites in Sect. 4. In Sect. 5, we interpret the GPS-IR results in revealing frozen-ground dynamics on various temporal scales, discuss the possible error sources of the results and their limitations in permafrost studies, and present their potential for validating and calibrating space-borne InSAR measurements. We conclude by summarizing the results and findings.

2 Methodology

2.1 GPS interferometric reflectometry (GPS-IR)

Larson (2016, 2019) presented the principle of GPS-IR and its applications in measuring snow depth, surface soil moisture, vegetation growth condition, and water level. However, the use of GPS-IR for studying permafrost was not explicitly presented. So, here, we describe how GPS-IR retrieves surface elevation changes and their link to ground deformation in permafrost areas.

The input of GPS-IR is signal-to-noise ratio (SNR) data of GPS signals, one of the observables recorded by GPS receivers. It represents the strength of the received signal. SNR series at low satellite elevation angles (e.g., 5–20^∘ used in this study) oscillate with respect to the elevation angle due to the interference between direct and reflected signals. The oscillating frequency mainly depends on the vertical distance between the antenna and the reflecting surface (called reflector height and denoted as H). If a GPS station is located above a smooth horizontal surface (e.g., Fig. 2), SNR can be expressed by a sinusoidal function of elevation angle e (Larson, 2016):

$$\begin{array}{}\text{(1)}& \mathrm{SNR}=A\left(e\right)\sin \left({\displaystyle \frac{\mathrm{4}\mathit{\pi }H}{\mathit{\lambda }}}\sin e+\mathit{\varphi }\right),\end{array}$$

where A(e) is the oscillation amplitude, also varying with e, λ is the wavelength of the carrier wave of GPS signal, and ϕ is the phase. When taking sin e as an independent variable, the oscillating frequency is

$$\begin{array}{}\text{(2)}& f={\displaystyle \frac{\mathrm{2}H}{\mathit{\lambda }}}.\end{array}$$

If f is determined, H can be obtained as

$$\begin{array}{}\text{(3)}& H={\displaystyle \frac{f\mathit{\lambda }}{\mathrm{2}}}.\end{array}$$

In practice, we divide the SNR series into individual segments corresponding to rising and setting satellite tracks. Then we remove their second-order polynomial fits and use the residual ones, which are mainly contributed from the reflected signals. For simplicity, we use SNR series hereafter to denote the residual SNR series. We conduct the Lomb–Scargle periodogram (LSP) analysis on any given SNR series to obtain its frequency spectrum. Then we use the peak value of the spectrum to represent the frequency f and obtain H by using Eq. (3). The oversampling parameter of LSP can be determined based on the expected resolution (e.g., 1 mm in this study) of the estimated reflector height. The programs for data processing are available in the software tools of GNSS interferometric reflectometry (Roesler and Larson, 2018).

Figure 2 Schematic diagram showing the geometry of the GPS antenna, GPS signals, and the ground surface in permafrost area. H, or reflector height, is the vertical distance between the GPS antenna and the surface, and e is the satellite elevation angle.

If the monument is deeply anchored (e.g., Fig. 2), the GPS antenna is stable with respect to the permafrost. The variation in the distance H only depends on the change in surface elevation. The change of H is opposite to that of surface elevation; i.e., surface uplift (subsidence) leads to decreasing (increasing) H (Liu and Larson, 2018). For the daily measurements of reflector height, we first assign minus signs to them and then remove the average to represent surface elevation changes.

2.2 A framework for identifying useful GPS stations for studying permafrost by GPS-IR

GPS-IR requires the ground surface to be open and relatively flat and smooth. To identify suitable ones from the existing GPS stations under the same protocols and to ensure the reliability of GPS-IR measurements, we design a three-step framework, which is described in detail as follows.

• Step 1: selecting GPS stations in permafrost areas.

  We first check whether permafrost is present where the GPS station is located. This step aims to identify the GPS stations in permafrost areas. We use the International Permafrost Association map compiled by Brown et al. (1997), which shows the spatial distribution of permafrost in the Northern Hemisphere.

• Step 2: estimating an azimuth range with an open, flat, and homogeneous ground surface.

  In this step, we aim to estimate an azimuth range, where the surface is open, nearly flat, and relatively homogeneous at each selected station in step 1. Normally, ground photos of a GPS station are taken as a part of metadata. In practice, we use ground photos and Google Earth images of a GPS site to check its surrounding environment and then estimate an azimuth range free of obstructions. We also use these images to recognize whether the selected surface is nearly planar and smooth. In Fig. S1 in the Supplement, we present ground photos of two GPS stations as typical positive and negative examples.

• Step 3: ensuring high reliability of GPS-IR measurements.

  At present, 32 operational GPS satellites orbit around the Earth twice daily. Therefore, multiple SNR series are available within a day. In practice, for any given day, we first process all SNR series within the determined azimuth range and elevation angle range to obtain their H using the standard method summarized in Sect. 2.1. Then, we calculate their median and discard the ones deviating from the median by 0.25 m or more. Then we compute the mean value and the standard deviation (σ) of the remaining H and remove those H values deviating from the mean value by larger than 3σ as outliers. The final retained H values and their corresponding SNR series are regarded as reliable. We average the finally retained H (denote the average as $\overline{H})$ to represent the vertical distance between the antenna and the reflecting surface on that day. The uncertainty of $\overline{H}$ is represented by its standard deviation, i.e., the standard deviation of H divided by the square root of the sample size. To further ensure the reliability of $\overline{H}$, a minimal number of 10 pieces of reliable SNR series are required.

2.3 Dataset and information

We use SNR data of L1 C/A signals of the identified GPS sites. L1 C/A is the legacy civilian code broadcasted by all the GPS satellites. By using them, we can obtain GPS-IR measurements spanning over several years to more than a decade, which enables us to study permafrost on various temporal scales and reveal its response to the changing climate.

To understand and interpret GPS-IR results, we use air temperature and snow depth at the study sites. These measurements are recorded by the nearest weather stations of GPS sites and can be downloaded from Environment Canada (http://climate.weather.gc.ca/historical_data/search_historic_data_e.html, last access: June 2020). We also use borehole ground temperatures at the study sites, which are provided by GTN-P (http://gtnpdatabase.org/boreholes, last access: June 2020).

We also summarize the climate and environment information of the study sites, such as mean annual air and ground temperature (Ednie and Smith, 2015; Environment Canada, http://climate.weather.gc.ca/climate_normals/index_e.html, last access: June 2020), surficial material (Cruishank, 1971; Dredge, 1994; Taylor et al., 1982; Throop et al., 2010), and ground ice content of near-surface permafrost (O'Neill et al., 2019), to provide background information for reference.

3 Identified GPS stations and study sites

We have screened all of the three major GPS networks in Canada, namely the Canadian Active Control System (CACS), Canadian High Arctic Ionospheric Network (CHAIN), and Portable Observatories for Lithospheric Analysis and Research Investigating Seismicity (POLARIS) (Fig. 3). CACS is a nationwide network and is maintained by the Geodetic Survey Division in conjunction with the Geological Survey of Canada (Lahaye et al., 2001). It serves to build and maintain the Canadian Spatial Reference System, which is fundamental for mapping, navigation, and studying crustal deformation. CHAIN was designed to investigate the impact of solar output on planetary environment (Jayachandran et al., 2009). The network is operated by the University of New Brunswick. It consists of 25 GPS stations, of which three (KUGC, REPC, and QIKC) are shared with CACS. It is important to note that most of the receiver antennas of CHAIN stations are anchored onto the roofs of buildings. Consequently, the monuments may move due to the foundation instability and thermal expansion or contraction of buildings. When using these stations for GPS-IR studies, corrections for such instability should be conducted. POLARIS, operated by the University of Western Ontario, was initiated for mapping solid Earth's structure and assessing earthquake hazards (Eaton et al., 2005). It includes seven geodetic-quality GPS stations.

Figure 3 Locations of GPS stations in the Canadian permafrost areas. The identified stations are labeled by their four-character IDs.

Following the framework in Sect. 2.2, we identified 12 GPS stations out of 38 ones located in permafrost areas as suitable ones for GPS-IR studies. Table 1 gives their basic information, including locations, monument types, foundation conditions, data time spans, and spatial coverage of GPS-IR measurements. The receiver and antenna types are listed in Table S1 in the Supplement. Five of them are from CACS, and the rest are from CHAIN. None of the POLARIS stations were identified as suitable.

Table 1 Basic information of the identified GPS stations. NA – not available

Given that the GPS-IR measurements of the CHAIN stations might be affected by the unstable buildings, in this study we present and interpret the measured elevation changes at the five identified CACS stations. Their monuments are all anchored into bedrock (Table 1). Figure 4 shows their ground photos, and Fig. S2 shows examples of their SNR series and corresponding LSP spectrum analysis.

Figure 4 Ground photos of the identified CACS GPS stations. Source: https://webapp.geod.nrcan.gc.ca/geod/data-donnees/cacs-scca.php?locale=_en.

These five sites are all located in the Canadian Arctic. The climate in this region is dominantly polar climate due to high latitude. The biomes are mainly tundra and Arctic desert. Permafrost is continuous in this area, and normally its thickness increases with latitude. In the far northern latitude of 75^∘, permafrost can be thicker than 500 m (Sladen, 2011). Ground temperatures at or near the depth of zero annual amplitude ranged from colder than −15 ^∘C to warmer than −2 ^∘C, and they decreased northward in concert with climate (Smith et al., 2013). During 2008–2014, ground temperatures at the depth of 15 m increased at an average rate of ∼0.17 ^∘C yr⁻¹ at 10 extensively distributed sites in the Canadian Arctic (Ednie and Smith, 2015). Thawing ice-rich permafrost has initiated wide-spread development of thermokarst landforms in this region, such as retrogressive thaw slump (Lantuit and Pollard, 2008; Kokelj et al., 2015) and active-layer detachment (Lewkowicz and Harris, 2005; Lewkowicz and Way, 2019).

We summarize the basic regional information of the five sites respectively in Table 2, including biome, land cover, ground ice content of near-surface permafrost, mean annual air temperature (MAAT), and mean annual ground temperature (MAGT). In Alert and Resolute Bay, the biomes are both arctic desert due to the high latitude, and the land surfaces are dominantly bare soil. The biomes at the other three sites are all tundra. But, due to their specific locations, the ground surface is mainly bare soil in Repulse Bay but is covered by a peat layer in Baker Lake and is sparsely vegetated in Iqaluit.

Table 2 Regional background of the study sites.

^a The land-cover information is for the areas around the boreholes, which are close to the identified GPS stations. ^b The ground ice contents are surficial material unit-based,
which are simulated with surficial geology, deglaciation, paleo-vegetation, glacial lake and marine limits, and modern permafrost distribution (O'Neill et al., 2019). ^c MAAT
refers to mean annual air temperature during 1981–2010 (Environment Canada, http://climate.weather.gc.ca/climate_normals/index_e.html). ^d MAGT refers to mean annual
ground temperature at or near the depth of zero annual amplitude, except Repulse Bay (Smith et al., 2013). MAGT at Repulse Bay was at the depth of 15 m (Ednie and
Smith, 2015). ^e Note that it is contrary to the field observations (Taylor et al., 1982) that found ground ice exists in the active-layer and near-surface permafrost in Alert.

4 Results: surface elevation changes measured by GPS-IR

We obtain multi-year and seasonal time series of surface elevation changes at the five CACS sites. In this study, we only present the measurements in thaw seasons, when air temperature is above 0 ^∘C and ground is not covered by snow.

We build the best linear fit to the thaw-season measurements and obtain the trends at the five sites (Fig. 5). We find that in Alert, Resolute Bay, and Repulse Bay, the ground surface subsided at a rate of 0.61±0.04 cm yr⁻¹ (2012–2018), 0.70±0.02 cm yr⁻¹ (2003–2014), and 0.26±0.05 cm yr⁻¹ (2014–2019), respectively. However, at the other two sites, the displacements of ground surface were 0.04±0.02 cm yr⁻¹ in Baker Lake during 2010–2017 and $-\mathrm{0.05}\pm \mathrm{0.02}$ cm yr⁻¹ in Iqaluit during 2010–2019. These last two trends were not statistically significant (t test, α=0.05).

Figure 5 Time series and the best-linear-fit dashed lines of surface elevation changes in thaw seasons at the five CACS sites. For clarity, we do not show the error bars. For the y axis, “relative” means that the presented elevation changes are referenced to the mean value of the entire records at each site.

In Fig. 6, we present the seasonal surface elevation changes in Resolute Bay. The seasonal results of the other sites can be found in Zhang et al. (2019). During a thaw season, the ground surface typically subsides progressively and reaches its lowest position at the end of season. However, at the Resolute Bay site, the surface elevation changes in the thaw seasons were irregular. The surface uplifted abnormally and significantly within the thaw seasons, for instance in 2003 and 2007. A similar phenomenon was also observed at a site near Yellowknife in Canada by Gruber (2020) using an inclinometer. This phenomenon could be due to the refreezing of soil moisture which migrated from the thawed active layer or the swelling of soil when it became wet. However, we lack measurements of soil moisture and ice content to investigate the cause of the observed uplift. In addition to such abnormal changes, the elevation changes among the thaw seasons were inconsistent. Given the complexity of these seasonal elevation changes, we turn to investigating the interannual variability and linear trends of surface deformation in Resolute Bay.

Figure 6 Surface elevation changes in each thaw season in Resolute Bay during 2003–2014. Red dots denote the measurements in the thaw seasons. Grey error bars denote the uncertainties. The mean value of the measurements has been removed. The shorter thaw season in 2013 was due to the late thawing onset on DOY 181 and early freezing onset on DOY 227, estimated from air temperature and snow depth records.

5 Discussion

In this section, we first interpret GPS-IR-measured surface elevation changes in Resolute Bay, as they are the longest among the five sites. We then qualitatively study the linear trends of surface deformation at the five sites. We also discuss the possible error sources of these GPS-IR measurements, the limitations of using GPS-IR measurements in permafrost studies, and their capability in validating and calibrating space-borne InSAR observations.

5.1 Interannual variability in end-of-thaw elevations in Resolute Bay

Net seasonal subsidence is an effective indicator of the response of frozen ground to the atmosphere, as it mainly depends on the soil moisture content within the active layer and the heat from the atmosphere. But, as shown in Sect. 4 and Fig. 6, it is challenging to reliably obtain seasonal subsidence in Resolute Bay due to the irregularity and inconsistency of surface elevation changes in thaw seasons. As an alternative, we use the end-of-thaw-season surface elevations to investigate the frozen-ground dynamics.

The end-of-thaw elevation is determined as the mean value of the elevations at the last 7 d of a thaw season, since the thawing front moves slowly at the end of thaw and the further surface deformation is limited. According to the Stefan equation, active-layer thickness is approximately proportional to the square root of the thaw index (Brown et al., 2000; Smith et al., 2009; French, 2017). The thaw index is represented by the degree days of thawing (DDT) derived by the accumulation of daily air temperatures above 0 ^∘C until the end of thaw season. As surface subsidence is mainly caused by ice melting within the active layer, we compare the end-of-thaw elevations to the square root of the annual thaw indices (Fig. 7).

Figure 7 (a) Time series of the end-of-thaw-season elevations and $\sqrt{\mathrm{DDT}}$ during 2003–2014. The right vertical axis for $\sqrt{\mathrm{DDT}}$ has been reversed to show the correlation between $\sqrt{\mathrm{DDT}}$ and the end-of-thaw elevations. (b) Scatter plot of the end-of-thaw elevations versus $\sqrt{\mathrm{DDT}}$. The red dashed line is the best-linear-fit line. (c) Scatter plot and the best-linear-fit line of the end-of-thaw-season elevations vs $\sqrt{\mathrm{DDT}}$ after removing the measurements of 2013 and 2014.

In Fig. 7a, the end-of-thaw-season elevations and $\sqrt{\mathrm{DDT}}$ were highly negatively correlated between 2003 and 2012, whereas the end-of-thaw elevations were low in cool summers in 2013 and 2014. To further investigate their correlation, we draw a scatter plot of end-of-thaw-season elevations versus $\sqrt{\mathrm{DDT}}$ (Fig. 7b) but find that the linear line fitted poorly. After removing the measurements in 2013 and 2014, the R² and root-mean-square error (RMSE) of the best linear fit improve significantly, from 0.24 to 0.83 and 2.57 to 1.19 cm, respectively (Fig. 7c).

We postulate that the highly negative correlation between the end-of-thaw elevations and $\sqrt{\mathrm{DDT}}$ during 2003–2012 was due to thickening active layer. A larger DDT indicates that more heat is available to penetrate into the deeper part of the frozen ground, leading to active-layer thickening, more ice melting within the frozen ground, and thus larger subsidence and lower surface elevation. This assumption of a thickening active layer during 2003–2012 is consistent with the borehole ground temperatures during 2008–2012 (Fig. 8). The ground temperatures showed that the thawing front (i.e., the 0 ^∘C isotherm) deepened and exceeded 1 m depth in 2011.

Figure 8 Time series of monthly ground temperatures at depths of 0.5, 1.0, and 1.5 m from June 2008 to September 2012 (Ednie and Smith, 2015). In August 2011, the ground temperature at 1.0 m depth was above 0 ^∘C.

However, in 2013 and 2014, the end-of-thaw elevations were low, even in the relatively cool summers (corresponding to low DDT). This is possibly due to the Markovian behavior of the active layer. Markovian behavior describes the reset of the active layer's response to air temperature after an extremely warm or cold summer, and this new response regime will last until the next extreme thaw season (Nelson et al., 1998). In Resolute Bay, the year 2011 had the warmest summer, with the DDT of 542.9 ^∘C d, more than 4 times larger than that in 2004 (132.3 ^∘C d). After 2011, the response of the active layer to the atmospheric forcing may have changed due to the changes in thermal properties of the active layer and ice content at the permafrost table. So, even with low DDT, the maximal thaw depths were still larger than expected, resulting in low end-of-thaw-season surface elevations. Yet, ancillary data such as thermal properties, ice content, soil moisture, and thaw depth are needed to test these postulated changes in the active layer.

5.2 Linear trends of surface deformation at the CACS sites

The ground surface deformed differently among the five sites. In this subsection, we study the possible links between linear trends of surface deformation and air temperature, land cover, and ground ice near the permafrost table.

We make basic statistics of the annual thaw indices during the study periods at the sites of Alert, Resolute Bay, and Repulse Bay (Table 3). All of these sites had warming thaw seasons, with trends of 9.35 ^∘C d yr⁻¹ in Alert during 2012–2018, 8.17 ^∘C d yr⁻¹ in Resolute Bay during 2003–2014, and 66.41 ^∘C d yr⁻¹ in Repulse Bay during 2014–2019, respectively.

Table 3 Basic statistics of annual DDT at Alert, Resolute Bay, and Repulse Bay.

^* Negative means subsidence, and vice versa.

The ground surface underwent subsidence with increasing DDT in Alert, Resolute Bay, and Repulse Bay. At these three sites, the surficial materials are sandy soil and barely vegetated (Table 2). Due to the lack of an insulating organic layer, bare soil facilitates the heat transfer between the atmosphere and the ground. When the climate was warming, the transient layer (i.e., the layer between the active-layer and long-term permafrost table and subjected to freeze and thaw seasonally to centennially; Shur et al., 2005) started to thaw with ground ice melting and surface subsidence, such as that seen in Alert, Resolute Bay, and Repulse Bay, even though they have low ice content in near-surface permafrost (Table 2).

Liu and Larson (2018) obtained surface elevation changes during 2004–2015 at Utqiaġvik, Alaska, by using GPS-IR and found a subsidence trend of 0.26±0.02 cm yr⁻¹. During the same time span, the thaw season in Utqiaġvik also had a warming trend, with 4.79 ^∘C d yr⁻¹. The results of Liu and Larson (2018) are consistent with ours: warming thaw seasons lead to surface subsidence. These findings in Utqiaġvik and our sites indicate that permafrost in high latitudes was degrading and air temperature is the dominant driver.

5.3 Possible error sources of GPS-IR-measured surface deformation

GPS-IR measurements of surface elevation changes might be affected by the surrounding environment (e.g., troposphere, vegetation, and soil moisture) and instruments (including antenna and monument). In this section, we discuss the impact of these variables on GPS-IR measurements and their magnitudes.

GPS signals refract when they propagate through the troposphere, leading to changes of propagating velocity and direction. Such refraction effects change the geometry among the direct and reflected signals and the receiver antenna and then introduce bias into reflector height retrievals. Tropospheric bias mainly depends on the antenna height and atmospheric conditions at a given elevation angle (Williams and Nievinski, 2017). In our study, because (1) all sites are located in the Canadian Arctic, characterized by a dry and cold climate, and (2) their antenna heights are ∼2 m (Table 1), the tropospheric biases at these sites are expected to be limited. More quantitatively, we calculate the tropospheric biases at RESO in the thaw season in 2014 by using the astronomical refraction model of Bennett (1982) and present them in Fig. S3. The magnitudes of tropospheric biases are ∼1.6 cm and stay relatively steady during the thaw season. As the magnitudes of biases are comparable to the uncertainties of GPS-IR measurements and we focus on the temporal variations in reflector heights, it is not necessary to correct for them.

Soil moisture also affects GPS-IR measurements of surface elevation changes through impacting phases of SNR series. For any given SNR series, soil moisture has a slightly different influence on the phase of each point; i.e., ϕ in Eq. (1) is also a function of elevation angle. Taking ϕ as a constant in practice introduces bias. Such bias is called compositional reflector height, as it manifests itself as a part of reflector height (Nievinski, 2013). Liu and Larson (2018) simulated the compositional height and found that they are less than 2 cm and their variation range is less than 1 cm, given a variation range of soil moisture between 15 % and 40 %. In this study, the compositional heights and their variation range are expected to be limited, as the precipitation is light and limited due to the cold and dry polar climate. Moreover, as we focus on the temporal variations in reflector heights at interannual and multi-annual timescales, we expect negligible impact of compositional heights on our results and interpretation.

Regarding the vegetation, at the study sites, the biomes are Arctic desert or tundra. The ground is bare or sparsely vegetated, and the vegetation is short enough, i.e., less than the wavelength of L-band GPS signals. The vegetation is approximately transparent for GPS signals. The impact of vegetation on GPS-IR measurements is therefore negligible.

The antenna gain pattern also impacts GPS-IR measurements. As the GPS stations used in this study were installed originally for geodetic or ionospheric studies, the receiver antennas were designed to favor direct signals with high elevation angles and suppressing signals with low and negative elevation angles by using asymmetric antenna gain patterns. During the data time span, the antennas are not replaced. The impacts of antenna gain patterns can be regarded as system biases and barely impact the GPS-IR results.

As for the monuments of the identified CACS stations, their material is galvanized or stainless steel and aluminum. The coefficients of linear thermal expansion measured at 20 ^∘C of steel and aluminum are 11–$\mathrm{13}\times {\mathrm{10}}^{-\mathrm{6}}$ and $\mathrm{23.1}\times {\mathrm{10}}^{-\mathrm{6}}$ m (m ^∘C)⁻¹, respectively. Given a temperature variation range of 20 ^∘C in the thaw season, for a 2 m high aluminum or steel monument, the magnitude of thermal expansion is less than 1 mm, at least 1 order of magnitude smaller than the elevation changes. The thermal expansion or contraction impact is ignorable for GPS-IR measurements.

5.4 Limitations of GPS-IR measurements of surface deformation in permafrost studies

GPS-IR-measured surface deformation has relatively large uncertainties, whose magnitudes are on the order of a few centimeters (i.e., ∼2 cm in Resolute Bay). The uncertainties are mainly caused by the rugged surface, presence of vegetation, and other unexpected disturbances. Such uncertainties make it difficult to study the daily changes of surface elevation based on GPS-IR measurements and even the seasonal changes if their magnitudes are comparable to those of seasonal subsidence. Resolute Bay is such a case, where daily and seasonal elevation changes cannot be obtained reliably. However, 12-year-long measurements enable the interannual variability in end-of-thaw elevations and decadal linear trend to be obtained with high confidence.

Data gaps exist in GPS observations due to instrumental problems. GPS-IR measurements before and after the gaps are contaminated by the bias introduced by the replacement of instruments. The data gaps and bias hinder the study of permafrost with long-term, continuous, and consistent GPS-IR measurements.

The interpretation of GPS-IR measurements in permafrost areas needs ground observations, such as soil temperatures and moisture. However, these data are usually not available at GPS sites, as they were installed initially for tectonic and ionospheric research. Moreover, surface condition records are often brief or absent. This being the case, we usually only have GPS-IR measurements and lack ancillary data such as ground temperature or soil moisture to help interpret the GPS-IR results.

These limitations indicate that, in the future, better location choices and maintenance of GPS stations and other ground measurement sensors are needed to exploit the full potential of GPS-IR observations in permafrost studies.

5.5 Potential of linking GPS-IR measurements to large-scale mapping from InSAR

Both GPS-IR and InSAR can measure surface elevation changes. In Table 4, we summarize their typical temporal and spatial samplings rates, advantages, and limitations. As we mentioned in Sect. 1, GPS-IR measurements are at daily intervals and local scales. In contrast, space-borne InSAR observations have much coarser temporal resolutions (the shortest to date being 6 d) and larger spatial scales (covering tens of kilometers) and also require a reference point with known surface deformation or that is assumed to be stable. These characteristics make GPS-IR and InSAR measurements complementary to each other. GPS-IR measurements could be used to overcome the limitations of InSAR observations. In particular, as GPS-IR measurements are continuous and at daily intervals over a few years to decades, they can provide baseline information for reference and can validate InSAR observations.

Table 4 Comparing GPS-IR and space-borne InSAR for measuring surface elevation changes in permafrost areas

Several major research programs such as the Arctic Boreal Vulnerability Experiment (ABoVE), Next-Generation Ecosystem Experiments (NGEE), and European Space Agency Permafrost Climate Change Initiative (CCI) use remote-sensing elevation changes (e.g., InSAR) to investigate permafrost dynamics. GPS-IR measurements can be used to calibrate and validate them and provide baseline information for historical, current, and future remote-sensing measurements from air and space.

6 Conclusions

In this study, for the first time, we implement a framework for assessing useful GPS stations for GPS-IR studies in permafrost areas and identify 12 useful GPS stations extensively distributed across the Canadian permafrost areas. Our framework can be applied to GPS networks in other regions and nations to identify more usable GPS stations. Our identified useful stations and the potential ones are also complementary to existing monitoring networks such as the CALM and GTN-P programs.

This study is also the first one using multiple GPS stations to study permafrost by GPS-IR. At the five identified CACS sites, we obtain their time series of elevation changes. The ground surface subsided in Alert by 0.61±0.04 cm yr⁻¹ during 2012–2018, in Resolute Bay by 0.70±0.02 cm yr⁻¹ during 2003–2014, and in Repulse Bay by 0.26±0.05 cm yr⁻¹ during 2014–2019. At the other two sites of Baker Lake and Iqaluit, the linear trends are not statistically significant. The trends at Alert, Resolute Bay, and Repulse Bay are negatively correlated to those of annual thaw indices; i.e., warming thaw seasons lead to surface subsidence. These findings collectively indicate that frozen ground at the study sites is sensitive to air temperature changes.

In Resolute Bay, we also find a highly negative correlation between the end-of-thaw elevations and the square root of thaw indices during 2003–2012 and suspect that it was possibly due to active-layer thickening under the warming thaw seasons. And we also find that the end-of-thaw elevations were low even with cool summers in 2013 and 2014. Continuous and daily measurements reveal the complexity of frozen-ground dynamics, i.e., the irregularity and inconsistency of seasonal surface elevation changes and the summer heave in Resolute Bay. To further investigate the dynamics and mechanisms of frozen-ground changes, it is important to measure other variables such as ground temperature, soil moisture, and ground ice content.

Our discussion on error sources and limits of GPS-IR measurements recommends that better location choice and maintenance of GPS stations be conducted to fully use the potential of those stations in frozen ground. The multi-year, continuous, daily GPS-IR measurements with intermediate spatial coverage can validate or calibrate remote-sensing observations of elevation changes in permafrost areas.