Ground surface elevation changes over permafrost areas revealed by multiple GNSS interferometric reflectometry

Ground subsidence and uplift caused by the annual thawing and freezing of the active layer are important variables in permafrost studies. Global positioning system interferometric reflectometry (GPS-IR) has been successfully applied to retrieve the continuous ground surface movements in permafrost areas. However, only GPS signals were used in previous studies. In this study, using multiple global navigation satellite system (GNSS) signal-to-noise ratio (SNR) observations recorded by a GNSS station SG27 in Utqiaġvik, Alaska during the period from 2018 to 2021, we applied multiple GNSS-IR (multi-GNSS-IR) technique to the SNR data and obtained the complete and continuous ground surface elevation changes over the permafrost area at a daily interval in snow-free seasons in 2018 and 2019. The GLONASS-IR and Galileo-IR measurements agreed with the GPS-IR measurements at L1 frequency, which are the most consistent measurements among all multi-GNSS measurements, in terms of the overall subsidence trend but clearly showed periodic noises. We proposed a method to reconstruct the GLONASS- and Galileo-IR elevation changes by specifically grouping and fitting them with a composite model. Compared with GPS L1 results, the unbiased root mean square error (RMSE) of the reconstructed Galileo measurements reduced by 50.0% and 42.2% in 2018 and 2019, respectively, while the unbiased RMSE of the reconstructed GLONASS measurements decreased by 41.8% and 25.8% in 2018 and 2019, respectively. Fitting the composite model to the combined multi-GNSS-IR, we obtained seasonal displacements of − 3.27 ± 0.13 cm (R2 = 0.763) and − 10.56 ± 0.10 cm (R2 = 0.912) in 2018 and 2019, respectively. Moreover, we found that the abnormal summer heave was strongly correlated with rain events, implying hydrological effects on the ground surface elevation changes. Our study shows the feasibility of multi-GNSS-IR in permafrost areas for the first time. Multi-GNSS-IR opens up a great opportunity for us to investigate ground surface movements over permafrost areas with multi-source observations, which are important for our robust analysis and quantitative understanding of frozen ground dynamics under climate change.


Introduction
Permafrost is the ground that remains at or below 0°C for at least two consecutive years and has undergone widespread degradation with global warming in past decades (Biskaborn B Zhenhong Li zhenhong.li@chd.edu.cn 1 et al. 2019). The active layer, which overlies the permafrost, causes seasonal subsidence and uplift of the ground surface in freeze/thaw cycles. The dynamic of the frozen ground serves as an important indicator of short-term climate change and exerts significant impacts on geomorphological, hydrological, biogeochemical processes, and infrastructures (Berger and Iams, 1996;Kokelj and Jorgenson, 2013). Continuous measurements of frozen ground movements are important for our quantitative understanding of the permafrost dynamics and can offer insights into the permafrost response to climate change.
Conventional field methods used to measure the vertical movements of frozen ground are often required to insert mechanical instruments such as telescoping tubes (Smith 1985), magnet probes (Mackay and Leslie 1987), and heavemeter (Mackay 1983) into the ground. Then, the thaw subsidence or frost heave can be retrieved by measuring the differential movements of the instrument benchmarks. Such conventional methods are limited to point measurement and are time-consuming. Traditional leveling can also be applied to measure elevation changes of the frozen ground surface (Burn 1990;Mackay and Burn 2002) with an accuracy of millimeters, but it needs a stable reference which is difficult to determine in permafrost areas. Targeting efficient surveys for frozen ground deformations over large areas, scholars introduced remote sensing techniques such as interferometric synthetic aperture radar (InSAR), light detection and ranging (LiDAR), and differential Global Positioning System (GPS) into permafrost monitoring. Among these remote sensing methods, InSAR measures the Earth's surface displacements with a wide coverage (e.g., 250 km × 250 km) and high precision (e.g., a few millimeters) (Bürgmann et al. 2000) and has been used to map seasonal movements of the ground surface (Liu et al. 2010;Short et al. 2011;Daout et al. 2017). Note that, in permafrost areas, InSAR often suffers from interferometric coherence loss and it is difficult (if not impossible) to find a stable reference point, limiting its use to obtain multi-year elevation changes. With the GPS receiver placed on the permafrost ground surface, the elevation changes can be quantified by the normal differential GPS method (Little et al. 2003;Shiklomanov et al. 2013). But the stability of the GPS receiver is largely affected by the ground, which is the major disadvantage of the differential GPS method. LiDAR has been used to investigate the ground deformation by imaging the topography with modest areal coverage (Jones et al. 2013). However, the costly instrument limits its extensive use, resulting in coarse temporal resolution of measurements (annually or monthly).
GPS Interferometric Reflectometry (GPS-IR) is a relatively new technique to measure the vertical movements of frozen ground. This technique uses the reflected GPS signals, which is quite different from the common positioning applications (e.g., differential GPS and precise point positioning) of GPS with the direct signals. GPS-IR is based on the fact that the reflected GPS signals contain the physical information of the reflector, in which the distance between the receiver phase center and reflector can be retrieved by quantifying the oscillations of the recorded SNR data. GPS-IR has been used to retrieve the snow depth variation (e.g., Larson et al. 2009;Larson and Nievinski 2013a), soil moisture (e.g., Larson et al. 2008;Chew et al. 2016), water level (e.g., Larson et al. 2013b), and coastal vertical land motion (e.g., Karegar et al. 2020) with an accuracy of centimeters. In previous permafrost studies, GPS-IR has been successfully applied to investigate the ground surface movements at different temporal scales from daily to inter-annual in the permafrost areas around the Arctic Hu et al. 2018;Zhang et al. 2020b) and in northeastern Qinghai-Tibet Plateau (Zhang et al. 2021). Compared with the aforementioned methods, GPS-IR is able to provide continuous and consistent measurements of frozen ground deformations at unprecedented daily intervals. Moreover, the GPS-IR elevation changes can be fixed to a geodetic reference frame by receiver point positioning. However, previous studies only used SNR data from GPS due to the limited signaltracking ability of receivers in permafrost areas. SNR data from other similar navigation satellite systems such as Russian GLONASS and European Galileo have also been used to study the ground surface variations (Qian and Jin 2016). All navigation satellite systems are collectively referred to as Global Navigation Satellite System (GNSS), so the GPS-IR can be extended to the generic GNSS-IR using GNSS SNR signals. In recent GNSS-IR studies, measurements of snow depth and sea levels derived from multiple GNSS-IR (multi-GNSS-IR) have been proved more robust and precise than that of a single system (Tabibi et al. 2017;Wang et al. 2019), while the performance of multi-GNSS-IR in permafrost areas has not been evaluated. The continuous development of multiple navigation systems and upgrade of GNSS receivers offer us a great opportunity to apply the multi-GNSS-IR technique to the abundant SNR data for studying frozen ground dynamics.
In this study, using the multiple SNR data collected by a permafrost GNSS site SG27 in Utqiaġvik, Alaska, we attempt to obtain the daily elevation changes by multi-GNSS-IR for the first time. The performances of multi-GNSS-IR elevation changes are assessed through a developed composite model driven by thermal forcing. We demonstrate the significant periodic noises in GLONASS-and Galileo-IR daily measurements, which are induced by topography differences associated with the revisit period of the GNSS satellites. Therefore, we then propose a model-based method to suppress the periodic noises and reconstruct the time series of multiple GNSS-IR elevation changes in the snow-free seasons of years 2018 and 2019. The combination of multiple GNSS-IR measurements for seasonal deformation retrieval and hydrological effects on elevation changes are also discussed in this paper.

Background of GNSS
GPS was the first GNSS consisting of 32 satellites at an altitude of 20,183 km. GPS satellites are unevenly distributed in six orbital planes inclined at an angle of 55°to the equator and orbit the Earth with a period of half a sidereal day, resulting in the daily repeat of their ground tracks. GPS satellites broadcast L-band signals at L1 (1575.42 MHz), L2 (1227.6 MHz), and L5 (1176.45 MHz) frequencies. GLONASS consists of 24 satellites homogeneously distributed in three orbital planes inclined at an angle of 64.8°to the equator. Each GLONASS satellite takes 11 h and 15.73 min to orbit the Earth at an altitude of 19,100 km and repeats its ground track every 8 sidereal days. GLONASS satellite signals are transmitted at two bands with initial frequencies of 1602 MHz and 1246 MHz, which vary a little satellite by satellite. The Galileo constellation consists of 20 + satellites at an altitude of 23,222 km. The satellites are distributed in three orbital planes inclined at an angle of 56°to the equator. Galileo satellites take about 14 h to orbit the Earth, resulting in a satellite ground-track repeat cycle every 10 sidereal days. Galileo satellite signals are transmitted in five frequencies including E1 (1575.42 MHz), E5a (1176.45 MHz), E5b (1207.15 MHz), E5 (1191.79 MHz), and E6 (1278.75 MHz). The total number of satellites from GPS, GLONASS, and Galileo exceeds 70 and is progressively increasing, providing users with abundant and diverse signals for multi-GNSS studies.

GNSS site and study area
SG27 is a continuous operating reference station (71.3229°N, 156.6103°W) using a geodetic GNSS receiver, located in Northern Utqiaġvik (formerly Barrow) Alaska. The Arctic Ocean is less than 3 km northwest of the site. SG27 is part of the EarthScope Plate Boundary Observatory (PBO) network recently integrated into a single pan-American network called Network of the Americas (NOTA). The primary scientific objectives of NOTA include the precise determination of plate motion and transient deformation related to earthquakes and volcanic eruptions. On June 18, 2018, the SG27 underwent a major instrument update with the receiver NetRS replaced by PolaRx, extending its ability to receive multi-GNSS (i.e., GPS/GLONASS/Galileo) signals. Before that, SG27 was able to record GPS signals only. The antenna of SG27 is~3.8 m above the ground surface and the bottom of the receiver monument is deeply anchored by~5 m beneath the ground. SG27 has open views toward east and north directions, and the surrounding ground is flat and homogenous polygon-free upland (Fig. 1), making it suitable for GNSS-IR experiments.
Utqiaġvik is located on the northern edge of Alaska's Arctic Coastal Plain underlain by continuous permafrost. The topographic relief of the Utqiaġvik area is low with elevations ranging from 0 to 20 m, and the graminoid-moss tundra prevails on organic-rich silty loam (Hinkel et al. 2001). Dominated by Arctic maritime climate, the annual precipitation in this area is 106 mm and 63% of the precipitation comes from rain events in the summer seasons (Shiklomanov et al. 2010). Mean monthly temperatures at Utqiaġvik range from − 24.6°C in February to 4.5°C in July with a mean annual temperature of − 9.7°C (Streletskiy et al. 2017). The active layer thickness (ALT) in Utqiaġvik is landscape-specific and varies from 29 cm in polygonized uplands to 83 cm in beach ridges (Shiklomanov et al. 2010). The thaw season in a typical calendar year lasts from early June to mid-August, but the end of the thaw season has tended to be later (September) in recent years with the global warming trend (Cox et al. 2017).

Meteorological data
To quantify the thermal forcing and interpret the frozen ground dynamics, we used near-surface air temperature by the co-located Barrow Observatory (Fig. 1). Barrow Observatory is maintained by National Oceanic and Atmospheric Administration (NOAA) and measures 2-m air temperature at an hourly interval. In this study, the hourly records of 2-m temperatures were simply averaged to construct daily observations. Daily precipitation records obtained by a weather station (USW00027502) near Utqiaġvik airport were used to interpret the summer heave during thaw seasons. This weather station is about 7 km southwest away from the SG27 station ( Fig. 1a).

GNSS-IR and seasonal deformation model
GNSS receiver receives two types of observations: (1) direct signals from GNSS satellites and (2) reflected signals from various reflectors such as ground, water, and buildings. The direct signals interfere with reflected signals and are recorded as SNR, which is the primary observables for GNSS-IR. For a planar ground, after removal of the trend contributed from direct signals using high-order polynomial, the detrended SNR (for simplicity, the detrended SNR is called SNR hereinafter) can be expressed as the function of the satellite elevation e (Larson 2016): where A and φ are the amplitude (unit: volts/volts) and phase offset (unit: rad), respectively. h (unit: meter) is the reflector height denoting the distance from the receiver phase center to the ground surface. λ (unit: meter) is the wavelength of the GNSS signal. Taking sin(e) as the independent variable, the SNR rate f (unitless, not exactly a frequency in hertz) equals 2h λ . Thus, the reflector height can be obtained with h f λ 2 . To avoid reflections from the buildings and other obstructions in the west and north directions of station SG27 (Fig. 1b), we applied an elevation angle mask (5-20°) and an azimuth mask (90-180°) to the GNSS SNR data as in previous studies Hu et al. 2018). The footprint of GNSS-IR depends on the first Fresnel zones ( Figure S1 in the supplemental material) of the reflected Fig. 1 Topography map of the study area. a The relief distribution around site SG27 derived from Wilson et al. (2014) and the footprint of GPS-IR at L1 frequency for example. Inset shows the locations of SG27 (red dot) and weather station coded USW00027502 (red triangle) in Utqiaġvik. b The views of SG27 toward four directions. The building of the Barrow Observatory is marked GNSS signal, which can be calculated with the signal wavelength, antenna height, elevation, and azimuth (see Text S1 in the supplemental material). For example, the footprint of GPS L1 is 6540.4 m 2 in this study (Fig. 1a). We conducted Lomb-Scargle periodogram (LSP) analysis to obtain the dominant frequency of the SNR data and scale it to the reflector height h. The retrieved reflector heights from all satellite tracks from each day are averaged to produce daily measurements with the standard deviation of the mean as their uncertainties. We did the same quality controls adopted by previous studies Hu et al. 2018) such as the median filter and the required minimum number of satellite tracks. In this study, we used all available GNSS SNR data to retrieve multi-GNSS-IR measurements from 2018 to 2021. Table 1 summarizes the GNSS signals used in our study.
The ground surface elevation changes d opposite to the reflector height changes: In this study, we fixed the reflector height to the a priori antenna height h 0 (i.e., 3.8 m) to obtain relative ground surface elevation changes. It should be noted that groundtrack-specific terrain differences exist in elevation changes obtained by Eq. (2). For repeat ground tracks, the terrain differences remain unchanged and thus do not affect the analysis of the temporal elevation changes. Affected by the solid earth dynamics, the vertical movements of the GNSS receiver of SG27 were up to a few centimeters during the study period ( Figure S2 in the supplemental material). It should be noted that the monument of SG27 was deeply installed into the  Figure S3 in the supplemental material for more details).
To characterize the deformation pattern of the frozen ground, we used the composite model developed by Hu et al. (2018). This composite model is derived from the Stefan function (Nelson et al. 1997) with some simplifications and assumptions. Here we only give a brief introduction to this composite model. The composite model characterizes the seasonal surface elevation changes as: where E is the variable relating to the thermal conductivity of the ground and I c is the composite index that can be expressed as: A T (t) and A F (t) are the accumulated degree days (units:°C days) calculated using the daily 2-m air temperature. The subscripts "T " and "F" represent the thaw season and freezing season, respectively. α is the scaling factor relating to the soil thermal conductivity and n-factor, which is 1.44 for Utqiaġvik (Hu et al. 2018). Chen et al. (2018) extended the composite model (i.e., Equation 3) to incorporate the long-term trend by adding a linear component which worked well in their data analysis with InSAR observations in a permafrost study. But this model was not adopted in this study since the dataset in this study only spanned two years. Specifically, we focused on the interpretation of the seasonal elevation changes in this study. In practice, we normalized the composite model: where d s is the maximum seasonal elevation changes; I c is the composite index normalized by its maximum value. d 0 is the reference elevation change relative to the a priori surface elevation (i.e., Equation 2), equivalent to the terrain reference in relative elevation changes measurements.

Preliminary elevation changes from GNSS-IR
First, we obtained the multi-GNSS reflector height measurements at a daily interval from June 18, 2018, to December 31, 2021. Figure 2a shows the representative results derived from GPS L1 SNR. The daily reflector heights derived from other GNSS signals are shown in Figure S4 in the supplemental material. From Fig. 2, we can see the times series of reflector heights in 2018 and 2019 are continuous and consistent while significant data gaps exist in snow-free days in 2020 and 2021. These gaps are also found in other time series of GNSS-IR reflector heights ( Figure S2 in the supplemental with the purple line and gray bars, respectively. Red and blue dashed lines mean the end of the snow season and snow onset, respectively, indicated by the snow depth measurements. Event A means the timing of the receiver change (18 June 2018). Events B, C, and D mean the data gaps (yellow shadow) in reflector heights. Event E means the data gap in snow depth measurements material), indicating poor quality issues of the SNR data on corresponding days. Therefore, we selected the continuous and complete records of reflector height during snow-free days in 2018 and 2019 for further processing and analysis in the main text. Then, we fixed the daily reflector heights with an a priori antenna height of 3.8 m (i.e., Eq. 2) to obtain the daily relative elevation changes of the ground surface during the snow-free seasons in years 2018 (Fig. 3) and 2019 (Fig. 4). The incomplete results during the snow-free seasons in years 2020 and 2021 are shown in Figures S5 and S6, respectively, in the supplemental material. For simplicity, GLONASS and Galileo are abbreviated as "GLO" and "GAL," respectively, in the figures and tables throughout this paper. Figures 3 and  4 show the time series of elevation changes derived from multi-GNSS, as well as the corresponding thermal indices. In 2018, the GPS L1 measurements show that the ground started to subside on 19 June and gradually subsided during the summer with a displacement range of~3 cm. The onset of the thaw season was 4 days later than the date (15 June) indicated by the thermal index √ A T (t). GPS measurements at both L2 and L5 frequencies show temporal changes similar to GPS L1, but with slightly larger noises and uncertainties. Despite the noises at centimeters level, the GLONASS-and Galileo-IR elevation changes also show the subsidence trend in July and August. The elevation changes time series are similar for GLONASS R1 and GLONASS R2, as well as the elevation changes of Galileo at five frequencies. Among all measurements, GPS L1 has the most consistent time series of measurements, clearly indicating the subsidence/uplift transition on 20 September 2018. The transition means the end of the thaw season and the onset of the freezing season, which was 3 days later than the date (17 September) determined by the thermal indices. However, it is hard to determine the thaw/freeze transitions from other GNSS-IR results due to the significant noises and greater uncertainties. In 2019, GPS L1 results are still the most consistent among the multiple elevation changes series. GPS L2 and L5 results are noisier than GPS L1, but both measurements are more consistent than Galileo and GLONASS results. GPS L1 elevation changes show that the ground surface underwent a strong settlement, with net subsidence exceeding 7 cm during the thaw season. The onset of subsidence on 26 June was 11 days later than the date (15 June) indicated by thermal indices. Other GNSS elevation changes also show an overall subsidence pattern resembling the GPS L1 measurements. The ground surface reached the lowest elevation in early September and kept nearly invariant till October. The GPS L1 elevation changes indicate that the freezing season started on 4 October 2019, which was 3 days ahead of the thermal-index-derived date (07 October 2019).
The subsidence and uplift of the ground surface are controlled by the phase changes of the water within the active layer, which is driven by the thermal exchanges between the ground and the atmosphere. Figures 3 and 4 show that the thermal indices basically characterize the pattern of elevation changes. To obtain the seasonal deformation amplitude of the ground surface, we also fitted the GNSS-IR measurements to the composite model (Figs. 5 and 6) and the statistics are summarized in Table 2. During the snow-free days in 2018, the fittings suggest maximum displacement ranging from − 4.11 ± 0.90 cm to − 2.85 ± 0.20 cm for different GNSS-IR elevation changes. The uncertainties in seasonal deformation of GLONASS and Galileo are substantially higher than that of GPS. A good agreement with an R 2 of 0.63 is obtained by GPS L1, followed by GPS L5 (R 2 0.52) and GPS L2 (R 2 0.28). Poor agreements with R 2 ranging from 0.12 to 0.26 are obtained for GLONASS and Galileo. Such poor agreements are mainly attributed to centimeter scale noise in GLONASS and Galileo measurements, which are comparable to the seasonal subsidence amplitude. In 2019, all the GNSS measurements suggest a similar displacement amplitude of about 10 cm and good agreements (R 2 exceeding 0.70). However, the fitting of GPS L1 to the model still has the smallest uncertainties and the best R 2 among all GNSS results as in 2018.

Time series reconstruction for GLONASSand Galileo-IR elevation changes
The elevation changes derived from GLONASS-and Galileo-IR are significantly noisier than that of GPS-IR, preventing the interpretation of the real changes of the frozen ground surface. To overcome the problem, we proposed a method to suppress the noises in the time series of measurements for GLONASS and Galileo. Unlike GPS with daily repeat ground tracks, the ground tracks of GLONASS and Galileo repeat every 8 days and 10 days, respectively. Figure 7 shows the daily azimuths of ground tracks of GPS, GLONASS, and Galileo during the study period. It should be noted that the ground tacks of GPS L2 and L5 overlap with a portion of GPS L1 ground tracks since L2 and L5 signals are not fully implemented in all GPS satellites. Moreover, ground tracks of GLONASS at different frequency bands are identical as well as the Galileo satellite ground tracks, so only GLONASS R1 and Galileo E1 are shown as representatives in Fig. 7. We can see that the GPS L1 ground tracks are very consistent at a daily interval, while the ground tracks of GLONASS R1 and Galileo E1 are scattered. The inconsistent daily tracks result in terrain differences within daily footprints and thus cause elevation biases in the day-to-day measurements which are fixed to the a priori reflector height (i.e., Eq. 2). These terrain effects have been recognized in previous GNSS-IR studies for snow depth monitoring (Tabibi et al. 2017;Zhang et al. 2020a). In  this study, we developed a method to remove the terrain differences over the permafrost area and thus reconstruct daily records of elevation changes for GLONASS and Galileo. First, we grouped the time series of elevation changes into subseries with the interval of the revisit period. In other words, we obtained 8 and 10 subseries with temporal consistent ground tracks for GLONASS and Galileo, respectively. Considering the terrain differences, we fitted the subseries to the composite model: where the subscript n denotes the subseries number.
is the coefficients matrix of reference elevation changes for subseries, in which m is the number of measurements in each subseries. As the extension of Eq. (5), Eq. (6) could be used to estimate terrain deference (i.e., d 0, n ) for each subseries. From the fitting, we obtained a unified seasonal deformation for the total subseries for GLONASS or Galileo. Table 2 shows the statistics of the fittings. Compared with the preliminary fitting results, the seasonal deformations in 2018 were closer to the GPS L1 and the uncertainties were significantly reduced by an average of 48.3% and 43.7% for GLONASS and Galileo results, respectively. The R 2 of GLONASS and Galileo also significantly increased to an average of 0.79 and 0.76, respectively. The fitting results in 2019 show that the deformation uncertainties of GLONASS and Galileo declined to an average of 0.31 cm and 0.36 cm, respectively, with R 2 exceeding 0.88.
We conducted the subseries fitting to obtain the deformation amplitude and sequenced reference elevation changes (i.e., topography differences) for the subseries. In turn, the reference elevation changes were used to remove the topographic effects and to reconstruct the time series of GLONASS-and Galileo-IR measurements. We set the reference elevation change of the first subseries (i.e., d 0, 1 ) as the reference terrain. Then, the rest of the subseries were fixed to the reference terrain by subtracting the differences between their reference elevation changes as: In Eq. (7), d r , n and d n are the reconstructed and original elevation changes, respectively, for subseries n. The terrain differences in subseries were removed through the fixing with Eq. (7).
The whole process of reconstruction is summarized in the flowchart as shown in Fig. 8.   Figures S7 and S8 in the supplemental material. We can see the reconstructed method significantly reduced the noises, and the temporal variations were closer to the GPS L1 results. Figure 9 also shows the power spectrum of the time series. The original elevation changes of GLONASS R1 showed an 8-day cycle which was removed after the reconstruction. For Galileo E1, the period noises related to the revisit cycle (10-day) were also significantly suppressed. Significant reductions in noise were also found in the rest of the reconstructed GLONASS and Galileo results ( Figures S3 and  S4 in the supplemental material). Taking the GPS L1 elevation changes as references, after reconstruction, the unbiased RMSE of GLONASS results reduced by an average of 41.8% (1.25-0.72 cm) and 25.8% (0.80-0.59 cm) in 2018 and 2019, respectively. For Galileo, the unbiased RMSE decreases were 50.0% (1.99-0.99 cm) and 42.2% (1.34-0.77 cm) in 2018 and 2019, respectively. We found that the results in 2019 are overall better than in 2018 by a few millimeters, possibly indicating that the ground surface was slightly rougher in 2018, which is associated with the annual variability of frozen ground movements. Another possible reason for the smaller RMSE in Galileo results in 2019 is the change in the number of satellites in orbit. The average number of available Galileo satellite tracks increased by 36% (11-15) from 2018 to 2019 (Fig. 12 in Sect. 5.2), resulting in a more uniform sampling of the ground in 2019.

Effects of slope and soil moisture on GNSS-IR-retrieved reflector heights
The retrieval of reflector heights is based on the assumption of flat ground. Ground tilting angles will cause a height bias in the retrieved reflector height (Larson and Nievinski 2013a). Soil moisture is another factor that significantly influences GNSS multipath modulation by changing soil permittivity. Variations in soil moisture will cause changes in retrieved reflector height, which are termed as "compositional height changes" (Nievinski 2013;Liu and Larson 2018). We used a GNSS multipath simulator (MPSImulator) developed by Nievinski and Larson (2014) to simulate the effects of slope and soil moisture on retrieved elevation changes. The settings are listed in Table 3. It should be noted that the simulator is only available for GPS and GLONASS. Figure 10a shows the slope distributions around SG27. The ground is quite flat within the footprint with an average slope of 2.16°. Figure 10b shows the slope simulation results. We can see that the height biases caused by the slope are up to a few centimeters with elevation changes range of [− 15

Soil moisture simulation
Soil moisture content Varying from 0 to 0.5 0] cm. However, the biases are relatively stable with amplitude smaller than 5 mm. Considering the centimeter-level uncertainties of GNSS elevation change measurements, these errors have a negligible impact on the temporal variation in relative elevation changes and thus the interpretation of the amplitude of ground surface movements. We did not correct the slope effects in this study. Figure 11 shows the soil moisture simulations. We can see that the compositional height changes are less than one centimeter with the soil moisture variations varying from 0 to 50%. Unfortunately, we do not know the soil moisture content at SG27 during the study period. However, a nearby (~60 m south-southeast of SG27) Circumpolar Active Layer Monitoring (CALM) site U1-1 measured the soil moisture from 1997 to 2011, which can provide a reference for us to infer the magnitude of soil moisture content. The records at depth of 5 cm show that the soil was saturated with water content at~35% and the soil moisture typically varied from 15 to 35% in thaw seasons ( Figure S5 in the supplemental material). According to the simulations, the maximum height change is smaller than 8 mm when soil moisture changes from the saturated state (i.e., 35%) to zero (i.e., soil is frozen), which has a limited impact on the measurements of relative elevation changes with typical uncertainties of 1.5-2.5 cm.

Combination of multiple GNSS-IR
Since the distributions of ground tracks of the three satellite constellations vary in azimuths (Fig. 7), GNSS SNR observations could be combined to provide dense measurements of the surface elevation changes. However, at the daily scale, elevation changes derived from other GPS (i.e., GPS L2 and L5), Galileo, and GLONASS are noisier than those from GPS L1 (Fig. 9). The greater uncertainties in measurements are partly due to the fewer number of valid satellite tracks. Specifically, the daily usable ground tracks of GPS L2 and GPS L5 were about half and one-third of the GPS L1 (3 0) during snow-free days in 2018 and 2019 (Fig. 12). As for GLONASS and Galileo, the daily usable ground tracks were only about two-thirds and one-third of the number of GPSL1 ground tracks. Thus, the simple averaging of the measurements derived from three constellations would introduce greater uncertainties in daily measurements than the single GPS L1 measurements. Consequently, we did not combine measurements to produce daily time series. Instead, we fitted all measurements to the model (i.e., Equation 6) to infer the maximum displacements at the seasonal scale. We obtained seasonal displacements of − 3.27 ± 0.13 cm (R 2 0.763) and − 10.56 ± 0.10 cm (R 2 0.912) in 2018 and 2019, respectively, from the model fitting with the combined multi-GNSS-IR, which were more robust than that of any single signal (Table 2).

Moving average of GLONASS-and Galileo-IR elevation changes
The inconsistencies in the original daily GLONASS-and Galileo-IR elevation changes are related to the 8-day and 10day cyclic ground tracks. A direct and simple way to solve this inconsistency is using the specific window averaging schemes according to the repeat cycle as 8-day averaging for GLONASS and 10-day averaging for Galileo. Figure 8 shows the moving average results for GLONASS R1 and Galileo E1 measurements as representatives. The window sizes for GLONASS and Galileo are 8 days and 10 days, respectively, with an identical step of 1 day. We can see that the averaged time series are smooth as this method successfully filters out the cyclic period noises. The averaged time series of the GLONASS and Galileo elevation changes are consistent and show a clear trend that resembles the pattern of the composite indices in each year (Figs. 3 and 4). Specifically, the subsidence in 2019 is significantly stronger than in 2018, which can be explained by the atmospheric forcing. The study area underwent an extremely warm summer in 2019 with a mean air temperature of 6.6°C, which is 2.2°C higher than that in 2018 (4.4°C). It indicates that the thaw front went deeper into the active layer under stronger atmospheric forcing in the 2019 summer, quantified by the large composite index (24.3 √ • C · day). The deeper thaw front caused more ice melt to water within the active layer, resulting in the strong subsidence of the ground surface in 2019. Compared with the GPS L1 results, the detailed temporal information at high frequency in the averaged time series is also removed. Moreover, the moving average method introduces high autocorrelations in the daily time series. The averaged time series can be used for studying long-term trends (seasonal to inter-annual trends) but is insensitive to short-term variations such as summer heave (Fig. 13).

Implication of hydrological effects on the elevation changes
The general trend of elevation changes driven by the thermal forces is well characterized by the composite model. However, the model could not explain the short-term summer heave revealed by the GNSS-IR or more specifically the GPS L1 measurements. We investigated the precipitation records obtained by the weather station USW00027502 and found that the significant summer heave coincided with the occurrences of rain events (Fig. 14), implying hydrological effects on the temporal uplifts in the thaw season (Gruber 2020). The precipitation could wet the soil and cause changes in soil permittivity. However, the increase in soil moisture will result in an increase in retrieved reflector height as discussed in Sect. 5.1. In other words, the wetting of soil would cause "subsidence" instead of "uplift" in the elevation measurements, but the magnitude is smaller than 4 mm given an extreme soil moisture variation of 20% ( Figure S9 in the supplemental material) during thaw seasons at SG27. We infer that the rain infiltrated the soil to the depth of the thawing front and was frozen to the ice with volume expansion, contributing to the temporal heave of the ground surface. The early uplift on October 4, 2019, was also possibly caused by the rain infiltration and refreezing, partly explaining the difference in freezing onset determination between the thermal index and GNSS-IR measurements (Fig. 4). Since the active layer is dominantly organic-rich soil in Utqiaġvik (Shiklomanov et al. 2010), the potential swelling of the soil under rain events may also contribute to the summer heave. We could not quantify these effects due to the lack of information on soil moisture within the active layer. However, the findings highlight the importance of considering the rain effects in the interpretation of the temporal deformation of the frozen ground. The rain effects are ignored in the simple composite model we used in this study, which could explain the non-linear trends (Fig. 6) in the fittings in 2019 with frequent rain events. Therefore, it is necessary to incorporate the hydrological effects to develop a complete permafrost model.

Conclusions
GPS-IR is a useful technique to retrieve continuous elevation changes of permafrost ground during the thaw and freezing seasons. In this study, we applied the multi-GNSS-IR (i.e., GPS-, GLONASS-and Galileo-IR) to the SNR observations collected by a permafrost GNSS site SG27 in Utqiaġvik, Alaska. The daily measurements of elevation changes in 2018 and 2019 were obtained for all GNSS signals. During the snow-free periods, multiple GNSS-IR measurements showed similar general movement trends of the permafrost ground surface, which was characterized by the composite model driven by thermal indices. Compared with the GPS L1, which has the most consistent measurements among all GNSS results, Galileo and GLONASS daily measurements showed obvious noises caused by the long repeat periods of the ground tracks. We proposed a method to remove the terrain differences and thus reconstruct the GLONASS and Galileo elevation changes time series, and the reconstructed time series agreed well with the GPS L1 results with the major noises suppressed. After reconstruction, the unbiased RMSE of GLONASS results reduced 41.8% (1.25-0.72 cm) and 25.8% (0.80-0.59 cm) in 2018 and 2019, respectively, while the unbiased RMSE of Galileo decreased 50.0% (1.99-0.99 cm) and 42.2% (1.34-0.77 cm) in 2018 and 2019, respectively. Combining the multi-GNSS-IR measurements and fitting them to the composite model, we obtained seasonal displacements of − 3.27 ± 0.13 cm (R 2 0.763) and − 10.56 ± 0.10 cm (R 2 0.912) in 2018 and 2019, respectively. The concurrences between the summer heave and rain events indicate the hydrological effects on the elevation changes.
Our study is the first attempt to use multi-GNSS SNR signals for permafrost dynamic study, which demonstrates the feasibility of GLONASS-and Galileo-IR for deformation monitoring over permafrost areas. We also found the evident hydrological effects on the elevation changes revealed by GPS-IR results. The periodic noise caused by the topographic variations is highlighted in this study, and the reconstruction method is recommended for multi-GNSS-IR measurements in other permafrost studies.
Author contributions The first author YH made the following contributions: conceptualization; methodology; formal analysis; investigation; validation; writing-original draft; writing-review and editing; project administration. JW as second author contributed to the paper as follows: formal analysis; validation; visualization. As third author, ZL made contributions to the following fields: conceptualization; funding acquisition; validation; writing-review and editing. JP as fourth author contributed to: resources; writing-review and editing; supervision.

Data availability
The raw GNSS data collected by SG27 and vertical position data are available from UNAVCO (http://pbo.unavco.org); the 2-m air temperature measurements are available from NOAA Global Monitoring Laboratory (https://gml.noaa.gov); the CALM soil moisture data are publicly available at https://www2.gwu.edu/~calm/data/ webforms/u1_f.htm; the dataset generated in this study are available upon request.

Conflicts of interest
The authors declare no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.