Snow Depth Variations in Svalbard Derived from GNSS Interferometric Reﬂectometry

: Snow plays a critical role in hydrological monitoring and global climate change, especially in the Arctic region. As a novel remote sensing technique, global navigation satellite system interferometric reﬂectometry (GNSS-IR) has shown great potential for detecting reﬂector characteristics. In this study, a ﬁeld experiment of snow depth sensing with GNSS-IR was conducted in Ny-Alesund, Svalbard, and snow depth variations over the 2014–2018 period were retrieved. First, an improved approach was proposed to estimate snow depth with GNSS observations by introducing wavelet decomposition before spectral analysis, and this approach was validated by in situ snow depths obtained from a meteorological station. The proposed approach can e ﬀ ectively separate the noise power from the signal power without changing the frequency composition of the original signal, particularly when the snow depth changes sharply. Second, snow depth variations were analyzed at three stages including snow accumulation, snow ablation and snow stabilization, which correspond to di ﬀ erent snow-surface-reﬂection characteristics. For these three stages of snow depth variations, the mean absolute errors (MAE) were 4.77, 5.11 and 3.51 cm, respectively, and the root mean square errors (RMSE) were 6.00, 6.34 and 3.78 cm, respectively, which means that GNSS-IR can be a ﬀ ected by di ﬀ erent snow surface characteristics. Finally, the impact of rainfall on snow depth estimation was analyzed for the ﬁrst time. The results show that the MAE and RMSE were 2.19 and 2.08 cm, respectively, when there was no rainfall but 5.63 and 5.46 cm, respectively, when it was rainy, which indicates that rainfall reduces the accuracy of snow depth estimation by GNSS-IR.


Introduction
The Arctic region, one of the regions with the most significant temperature increases, has exhibited a significant response to climate change [1], and the snow and glaciers in Svalbard have changed sharply in recent years [2,3]. Snow is a critically important and rapidly changing feature of the Arctic [4], but snow depth observations in the Arctic are still insufficient in terms of quantity and quality. Ultrasonic sensors can only obtain the snow depth at certain points, and the discreteness and uneven distribution of sites severely restrict the real-time acquisition and application of snow information [5]. Remote sensing methods can obtain large-scale snow information, but the accuracy Remote Sens. 2020, 12, 3352 3 of 16

Data
A field experiment with GNSS-IR to sense snow depth variations was carried out in Ny-Alesund (78 • 55 N, 11 • 55 E), Svalbard, as shown in Figure 1. The local environment is ideal for snow depth estimation due to the radio silence in Ny-Alesund and flat ground around the GNSS tracking station. Considering a 1.67 m-high antenna, the specular reflection footprint is within a circle of~20 m in a planar area when the satellite elevation angle is 5 • -25 • . Therefore, the shape of the terrain can be taken as a small plane and has almost no influence on snow depth estimation. The GNSS station is equipped with a Unicore UR4B0 receiver and a Novatel GNSS-750 antenna, and GPS L1 SNR observations from 2014 to 2018 were collected at a sampling rate of 30 s. The meteorological station at Ny-Alesund, approximately 230 m from the GNSS station, provided in situ snow depth measurements and rainfall observations. The in-situ snow depth measurements from the meteorological station were collected at UTC 6 h every day. All the meteorological observations can be found at http://eklima.met.no, which contains historical data and real-time meteorological observations from the Norwegian Meteorological Institute (NMI).

Data
A field experiment with GNSS-IR to sense snow depth variations was carried out in Ny-Alesund (78°55′N, 11°55′E), Svalbard, as shown in Figure 1. The local environment is ideal for snow depth estimation due to the radio silence in Ny-Alesund and flat ground around the GNSS tracking station. Considering a 1.67 m-high antenna, the specular reflection footprint is within a circle of ~20 m in a planar area when the satellite elevation angle is 5°-25°. Therefore, the shape of the terrain can be taken as a small plane and has almost no influence on snow depth estimation. The GNSS station is equipped with a Unicore UR4B0 receiver and a Novatel GNSS-750 antenna, and GPS L1 SNR observations from 2014 to 2018 were collected at a sampling rate of 30 s. The meteorological station at Ny-Alesund, approximately 230 m from the GNSS station, provided in situ snow depth measurements and rainfall observations. The in-situ snow depth measurements from the meteorological station were collected at UTC 6 h every day. All the meteorological observations can be found at http://eklima.met.no, which contains historical data and real-time meteorological observations from the Norwegian Meteorological Institute (NMI).

Approach Based on LSP Spectral Analysis
Note that the signal-to-noise ratio, after removing the trend term, is a sine-like function that oscillates with the change in satellite elevation. The oscillation frequency can be obtained by differentiating the change in satellite elevation with time. Because the frequency characteristics contain height information, we can derive the snow depth estimation from the height of the antenna phase center relative to the reflector. In this paper, SNR observation data with a satellite elevation of 5°-25° were selected for further processing and analysis. Since the SNR sequence is not sampled at equal intervals, LSP spectral analysis is often used to process the non-uniformly sampled signal-tonoise ratio sequences after trend removal to obtain the main frequency [26,27]. Due to the influence of thermal and environmental noise, on some occasions, it is difficult to obtain the main frequency of the signal-to-noise ratio sequences via LSP spectral analysis. Figure 2 shows the spectral results of the de-trended SNR sequence for GPS pseudo-random number (PRN) 19 on February 28, 2018. This is a case that illustrates the failure of LSP in obtaining the main frequency, because two peaks can be found, which correspond to two main frequencies and can cause an error in snow depth estimation. Since it is particularly important to control the quality of the GNSS observations in snow depth estimation, an improved approach is proposed in the next section.

Approach Based on LSP Spectral Analysis
Note that the signal-to-noise ratio, after removing the trend term, is a sine-like function that oscillates with the change in satellite elevation. The oscillation frequency can be obtained by differentiating the change in satellite elevation with time. Because the frequency characteristics contain height information, we can derive the snow depth estimation from the height of the antenna phase center relative to the reflector. In this paper, SNR observation data with a satellite elevation of 5 • -25 • were selected for further processing and analysis. Since the SNR sequence is not sampled at equal intervals, LSP spectral analysis is often used to process the non-uniformly sampled signal-to-noise ratio sequences after trend removal to obtain the main frequency [26,27]. Due to the influence of thermal and environmental noise, on some occasions, it is difficult to obtain the main frequency of the signal-to-noise ratio sequences via LSP spectral analysis. Figure 2 shows the spectral results of the de-trended SNR sequence for GPS pseudo-random number (PRN) 19 on February 28, 2018. This is a case that illustrates the failure of LSP in obtaining the main frequency, because two peaks can be found, which correspond to two main frequencies and can cause an error in snow depth estimation. Since it is particularly important to control the quality of the GNSS observations in snow depth estimation, an improved approach is proposed in the next section.

Improved Approach Based on Wavelet Analysis
Commonly used methods for improving the signal-to-noise ratio include the empirical modal solution, independent component analysis, wavelet decomposition and so on. However, under low signal-to-noise ratio conditions, there may be aliasing effects with the empirical modal solution, the autocorrelation matrix of the signals may have zero eigenvalues, and the whitening process cannot be completed for independent component analysis. By contrast, the multi-resolution characteristics of wavelet decomposition can accurately represent the local information of the signal in the time and frequency domains and overcome the shortcomings of the short-time Fourier transform at the singleresolution scale. Wavelet decomposition is now widely used in time-frequency analysis, signal-tonoise separation and weak-signal extraction [28,29]. In this study, before performing LSP spectral analysis on the SNR sequence, wavelet decomposition was carried out to reduce the influence of the noise power. Grossmann and Morlet proposed the use of wavelet decomposition to express the signals that meet the conditions as a series of successive approximation expressions for the first time, and each of these expressions is the smoothed form of the original signal, corresponding to different resolutions [30]. Wavelet analysis has the advantage of multiple resolutions, and both the time and frequency windows can be dynamically adjusted according to the specific form of the signal. Multiresolution analysis is a set of functions that construct closure approximation L 2 (R), and the set of functions constitutes a normalized orthogonal basis. For any given signal, its expansion in wavelet space can be described as , , where Wj is the orthogonal complement of the closed subspace Vj on L 2 (R); Ψ(t) is a function in L 2 (R). The corresponding multi-resolution analysis is A peak corresponds to a main frequency, but two similar peaks mean the failure of LSP in obtaining the main frequency.

Improved Approach Based on Wavelet Analysis
Commonly used methods for improving the signal-to-noise ratio include the empirical modal solution, independent component analysis, wavelet decomposition and so on. However, under low signal-to-noise ratio conditions, there may be aliasing effects with the empirical modal solution, the autocorrelation matrix of the signals may have zero eigenvalues, and the whitening process cannot be completed for independent component analysis. By contrast, the multi-resolution characteristics of wavelet decomposition can accurately represent the local information of the signal in the time and frequency domains and overcome the shortcomings of the short-time Fourier transform at the single-resolution scale. Wavelet decomposition is now widely used in time-frequency analysis, signal-to-noise separation and weak-signal extraction [28,29]. In this study, before performing LSP spectral analysis on the SNR sequence, wavelet decomposition was carried out to reduce the influence of the noise power. Grossmann and Morlet proposed the use of wavelet decomposition to express the signals that meet the conditions as a series of successive approximation expressions for the first time, and each of these expressions is the smoothed form of the original signal, corresponding to different resolutions [30]. Wavelet analysis has the advantage of multiple resolutions, and both the time and frequency windows can be dynamically adjusted according to the specific form of the signal. Multi-resolution analysis is a set of functions that construct closure approximation L 2 (R), and the set of functions constitutes a normalized orthogonal basis. For any given signal, its expansion in wavelet space can be described as where Wj is the orthogonal complement of the closed subspace Vj on L 2 (R); Ψ(t) is a function in L 2 (R). The corresponding multi-resolution analysis is Remote Sens. 2020, 12, 3352

of 16
where ϕ is the scale function, j and j' are the resolution scale, J is the maximum resolution scale, and k is the translation scale. For the signal-to-noise ratio sequences after trend removal, the Mallat algorithm was used for multi-resolution analysis, and the features of interest are extracted by reconstructing the original signal structure [31]. The corresponding frequencies of different resolution scales can be calculated as where wj is the center frequency and ∆wj is the frequency radius. Since wj is an orthogonal complement of Vj, w is also an orthogonal sequence set. After a comparative experiment, a compactly supported orthogonal wavelet decomposition was performed with a wavelet order of 6 [32]. The results of the discrete wavelet decomposition of the SNR sequences of the PRN17 satellite on February 28 of 2018 are shown in Figure 3.
Remote Sens. 2020, 12, x FOR PEER REVIEW 5 of 16 where φ is the scale function, j and j' are the resolution scale, J is the maximum resolution scale, and k is the translation scale. For the signal-to-noise ratio sequences after trend removal, the Mallat algorithm was used for multi-resolution analysis, and the features of interest are extracted by reconstructing the original signal structure [31]. The corresponding frequencies of different resolution scales can be calculated as where wj is the center frequency and ∆wj is the frequency radius. Since wj is an orthogonal complement of Vj, w is also an orthogonal sequence set. After a comparative experiment, a compactly supported orthogonal wavelet decomposition was performed with a wavelet order of 6 [32]. The results of the discrete wavelet decomposition of the SNR sequences of the PRN17 satellite on February 28 of 2018 are shown in Figure 3. The wavelet decomposition of the signal has the following structure: [cD1, cD2, cD3, cD4, cD5, cD6, cA6], where cD is the detailed coefficient of the sequence, and cA is the approximation coefficient of the sequence. The wavelet decomposition of the signal has the following structure: [cD1, cD2, cD3, cD4, cD5, cD6, cA6], where cD is the detailed coefficient of the sequence, and cA is the approximation coefficient of the sequence. Different resolution scales correspond to completely different frequencies, and the higher the resolution scale, the narrower the frequency window. The coefficients obtained by wavelet decomposition correspond to different signal frequencies, such as the receiver thermal-noise frequency, environmental-noise frequency and signal frequency related to snow depth. To further determine the depth-related frequency in the original signal, LSP spectral analysis was performed on the coefficients Remote Sens. 2020, 12, 3352 6 of 16 at each resolution scale; the results corresponding to the different resolution scales are shown in Figure 4.
Remote Sens. 2020, 12, x FOR PEER REVIEW 6 of 16 Different resolution scales correspond to completely different frequencies, and the higher the resolution scale, the narrower the frequency window. The coefficients obtained by wavelet decomposition correspond to different signal frequencies, such as the receiver thermal-noise frequency, environmental-noise frequency and signal frequency related to snow depth. To further determine the depth-related frequency in the original signal, LSP spectral analysis was performed on the coefficients at each resolution scale; the results corresponding to the different resolution scales are shown in Figure 4. .43 cm, respectively. According to the measurements provided by the Ny-Alesund weather station, the snow depth measured at 6:00 UTC on February 28, 2018, was 27 cm. Therefore, cD4+cD5 is more effective for denoising than the other coefficients. For cD4, although also close to the true values, the results do not show good denoising effects. After multi-resolution analysis, the reconstructed spectral-analysis results are in good agreement with the measured snow depth.
After estimating the snow depth from individual SNR sequences, a filtering procedure was applied to obtain the daily average snow depth by using observations of all the satellites over the whole azimuth in one day. Results with random errors over three standard deviations were eliminated, as were negative estimations of daily snow depth.

Daily Averaged Snow Depths Over 5 Years
Since a filtering procedure was applied to the individual snow depth retrievals in this study, we first checked the accuracy of individual satellite-derived and measured snow depths. Figure 5 presents two cases of the comparison between the individual satellite-derived and measured snow .43 cm, respectively. According to the measurements provided by the Ny-Alesund weather station, the snow depth measured at 6:00 UTC on 28 February 2018, was 27 cm. Therefore, cD4+cD5 is more effective for denoising than the other coefficients. For cD4, although also close to the true values, the results do not show good denoising effects. After multi-resolution analysis, the reconstructed spectral-analysis results are in good agreement with the measured snow depth.
After estimating the snow depth from individual SNR sequences, a filtering procedure was applied to obtain the daily average snow depth by using observations of all the satellites over the whole azimuth in one day. Results with random errors over three standard deviations were eliminated, as were negative estimations of daily snow depth.

Daily Averaged Snow Depths over 5 Years
Since a filtering procedure was applied to the individual snow depth retrievals in this study, we first checked the accuracy of individual satellite-derived and measured snow depths. Figure 5 presents two cases of the comparison between the individual satellite-derived and measured snow depths in one day. There were more than 30 operational satellites in the GPS constellation during the field experiment, and GPS satellites orbit the earth every 12 hours, with one rising arc and one setting arc, which means more than 120 dSNR sequences can be used to estimate the snow depth in one day. As shown in Figure 5, although more than 100 individual satellite-derived snow depths were acquired in one day, they were remarkably scattered due to the effect of noise.
Remote Sens. 2020, 12, x FOR PEER REVIEW 7 of 16 depths in one day. There were more than 30 operational satellites in the GPS constellation during the field experiment, and GPS satellites orbit the earth every 12 hours, with one rising arc and one setting arc, which means more than 120 dSNR sequences can be used to estimate the snow depth in one day. As shown in Figure 5, although more than 100 individual satellite-derived snow depths were acquired in one day, they were remarkably scattered due to the effect of noise.
(a) (b) Figure 5. Scatterplot of individual satellite-derived snow depths versus in-situ snow depths on Jan. 30, 2015 (a) and Apr. 30, 2015 (b). Red squares represent the measured snow depths from field work, collected at UTC 6 h every day. Blue circles correspond to snow depths derived from different satellites. Figure 6 and Figure 7 shows the comparison between individual satellite-derived and measured snow depths in 2014 and 2015, which can fully represent the effects of the filtering procedure under different snow packs (shallow, medium and thick). Similar to those in Figure 5, the snow depths obtained by using individual satellites only have an accuracy of a few tens of centimeters, consistent with the results of [15,23]. Compared to the soil surface, the snow surface is significantly more spatially heterogeneous [33]. Therefore, a filtering procedure for the daily average is necessary because of the low accuracy of the snow depth estimated from individual satellites.   Figure 5, the snow depths obtained by using individual satellites only have an accuracy of a few tens of centimeters, consistent with the results of [15,23]. Compared to the soil surface, the snow surface is significantly more spatially heterogeneous [33]. Therefore, a filtering procedure for the daily average is necessary because of the low accuracy of the snow depth estimated from individual satellites.
Remote Sens. 2020, 12, x FOR PEER REVIEW 7 of 16 depths in one day. There were more than 30 operational satellites in the GPS constellation during the field experiment, and GPS satellites orbit the earth every 12 hours, with one rising arc and one setting arc, which means more than 120 dSNR sequences can be used to estimate the snow depth in one day. As shown in Figure 5, although more than 100 individual satellite-derived snow depths were acquired in one day, they were remarkably scattered due to the effect of noise.
(a) (b) Figure 5. Scatterplot of individual satellite-derived snow depths versus in-situ snow depths on Jan. 30, 2015 (a) and Apr. 30, 2015 (b). Red squares represent the measured snow depths from field work, collected at UTC 6 h every day. Blue circles correspond to snow depths derived from different satellites. Figure 6 and Figure 7 shows the comparison between individual satellite-derived and measured snow depths in 2014 and 2015, which can fully represent the effects of the filtering procedure under different snow packs (shallow, medium and thick). Similar to those in Figure 5, the snow depths obtained by using individual satellites only have an accuracy of a few tens of centimeters, consistent with the results of [15,23]. Compared to the soil surface, the snow surface is significantly more spatially heterogeneous [33]. Therefore, a filtering procedure for the daily average is necessary because of the low accuracy of the snow depth estimated from individual satellites.   The estimated daily snow depths were obtained by averaging individual satellite-derived snow depths for one day. The comparison between the estimated daily snow depth and measured daily snow depth in Ny-Alesund during the 2014-2018 period is presented in Figure 8, and the performance of LSP and the proposed approach is also compared in Figure 8. The data from January to June of each year were analyzed in detail because the snow depth changes remarkably during this period, while the snow depth in other periods is very small or even 0.  The estimated daily snow depths were obtained by averaging individual satellite-derived snow depths for one day. The comparison between the estimated daily snow depth and measured daily snow depth in Ny-Alesund during the 2014-2018 period is presented in Figure 8, and the performance of LSP and the proposed approach is also compared in Figure 8. The data from January to June of each year were analyzed in detail because the snow depth changes remarkably during this period, while the snow depth in other periods is very small or even 0.  Figure 8 shows that the snow depth in Ny-Alesund reaches a maximum in spring, with slow accumulation in winter and fast melting in summer. Compared to that in 2014, the amount of snow in Ny-Alesund has declined remarkably in recent years, and the snow depth has generally been lower than 40 cm. Moreover, frequent fluctuation of the snow depth can be observed in Figure 5. The temperature rise in the Arctic region, one of the most significant regions of temperature increase in the 21st century, is aggravating the loss of snow cover year by year. Although there is a certain systematic error between the estimated and measured snow depth, the two methods generally remain consistent with the field measurements. The estimated daily snow depths were obtained by averaging individual satellite-derived snow depths for one day. The comparison between the estimated daily snow depth and measured daily snow depth in Ny-Alesund during the 2014-2018 period is presented in Figure 8, and the performance of LSP and the proposed approach is also compared in Figure 8. The data from January to June of each year were analyzed in detail because the snow depth changes remarkably during this period, while the snow depth in other periods is very small or even 0.  Table 1. Figure 8 shows that the snow depth in Ny-Alesund reaches a maximum in spring, with slow accumulation in winter and fast melting in summer. Compared to that in 2014, the amount of snow in Ny-Alesund has declined remarkably in recent years, and the snow depth has generally been lower than 40 cm. Moreover, frequent fluctuation of the snow depth can be observed in Figure 5. The temperature rise in the Arctic region, one of the most significant regions of temperature increase in the 21st century, is aggravating the loss of snow cover year by year. Although there is a certain systematic error between the estimated and measured snow depth, the two methods generally remain consistent with the field measurements.

Performance of the Improved Approach
To verify the reliability and applicability of the proposed approach, the standard deviation (STD), root mean square error (RMSE) and mean absolute error (MAE) of the two approaches were Red line represents the measured snow depth from field work. Green and blue lines denote the snow depths estimated with LSP and LSP after wavelet decomposition, respectively. The error bars show the daily standard deviations of the GNSS-based estimations, and the average standard deviations were ∼7.13 and 6.80 cm for LSP and LSP after wavelet decomposition, respectively, as also shown in Table 1.

Performance of the Improved Approach
To verify the reliability and applicability of the proposed approach, the standard deviation (STD), root mean square error (RMSE) and mean absolute error (MAE) of the two approaches were analyzed in detail. The RMSE and MAE can be used to evaluate the difference between the estimated value and referenced value, while the STD is used to quantify the variation or dispersion of a set of data [34][35][36]. Table 1 compares the two methods for snow depth estimation. The STD is a good indicator of the stability of the results. It can be observed that the proposed approach has higher precision than the LSP method in terms of the STD, indicating that its results are more robust. Furthermore, the MAEs of the estimated and measured values are 5.72 and 5.57 cm, and the RMSEs are 7.52 and 7.37 cm, respectively, for the LSP approach and the proposed approach, which means that the proposed method is more consistent with the field measurements than the original approach. By comparing the accuracy in different years, we can see that the improvement is limited because the wavelet decomposition transform is only a denoising method and does not change the frequency characteristics contained in the original signal. When the noise is large, as shown in Figure 2, the improvement effect is obvious; however, when the noise is small, the accuracy of the two methods is roughly equivalent.
As discussed in the method section, after estimating the snow depth from individual SNR sequences, the daily average snow depth was obtained using observations of all the satellites over the whole azimuth in one day. For the daily snow depth estimation, snow depths derived from individual satellites with random errors over three standard deviations were eliminated during data processing, so data utilization is an important indicator for evaluating the reliability of the improved method. Figure 9 represents the data utilization rates of the two methods over 5 years, and Table 2 shows the annual mean data utilization rates of the two methods. Compared with the LSP method, the improved wavelet-based method has significantly higher data utilization, especially in 2014, 2016 and 2018, and the mean data utilization rate over 5 years increased from 72.18% to 80.94% for the improved method, which means that the improved method is more robust and reliable.
Remote Sens. 2020, 12, x FOR PEER REVIEW 10 of 16 As discussed in the method section, after estimating the snow depth from individual SNR sequences, the daily average snow depth was obtained using observations of all the satellites over the whole azimuth in one day. For the daily snow depth estimation, snow depths derived from individual satellites with random errors over three standard deviations were eliminated during data processing, so data utilization is an important indicator for evaluating the reliability of the improved method. Figure 9 represents the data utilization rates of the two methods over 5 years, and Table 2 shows the annual mean data utilization rates of the two methods. Compared with the LSP method, the improved wavelet-based method has significantly higher data utilization, especially in 2014, 2016 and 2018, and the mean data utilization rate over 5 years increased from 72.18% to 80.94% for the improved method, which means that the improved method is more robust and reliable. 2018 (e). Blue and red lines denote the data utilization rates for snow depth estimation with LSP and LSP after wavelet decomposition, respectively. The data utilization rate refers to the ratio of individual-satellite-derived snow depths actually used to calculate the daily averaged snow depth.

Effects of Snow-Surface Characteristics on Estimated Snow Depth
The estimated snow depth is strongly affected by the reflected snow surface characteristics. The influence of dry snow conditions were discussed by Jacobson [10], Larson [14] and Henkel [37]. Najibi further investigated the interactions of snow accumulation, snow melting, bare soil, and fixed snow depth [38]. Koch derived snow height under three periods, including the beginning of the snowcovered season, snow accumulation period, and snow melting period [39]. In this study, given the complex snow conditions in Ny-Alesund, the changes in snow depth were divided into the following three stages: the snow-ablation, snow-accumulation and snow-stabilization stages. The differences between the estimated and measured values at different stages were studied in detail. The snowaccumulation stage is defined as an increase in snow by 2 cm for three consecutive days, while the snow-ablation stage is defined as a decrease in snow by 2 cm for three consecutive days. By contrast, the snow-stabilization stage refers to a change in snow depth within 1 cm for three consecutive days. The results are shown in Table 3. 2018 (e). Blue and red lines denote the data utilization rates for snow depth estimation with LSP and LSP after wavelet decomposition, respectively. The data utilization rate refers to the ratio of individual-satellite-derived snow depths actually used to calculate the daily averaged snow depth.

Effects of Snow-Surface Characteristics on Estimated Snow Depth
The estimated snow depth is strongly affected by the reflected snow surface characteristics. The influence of dry snow conditions were discussed by Jacobson [10], Larson [14] and Henkel [37]. Najibi further investigated the interactions of snow accumulation, snow melting, bare soil, and fixed snow depth [38]. Koch derived snow height under three periods, including the beginning of the snow-covered season, snow accumulation period, and snow melting period [39]. In this study, given the complex snow conditions in Ny-Alesund, the changes in snow depth were divided into the following three stages: the snow-ablation, snow-accumulation and snow-stabilization stages. The differences between the estimated and measured values at different stages were studied in detail. The snow-accumulation stage is defined as an increase in snow by 2 cm for three consecutive days, while the snow-ablation stage is defined as a decrease in snow by 2 cm for three consecutive days. By contrast, the snow-stabilization stage refers to a change in snow depth within 1 cm for three consecutive days. The results are shown in Table 3. Table 3 shows that the proposed approach performs better than the traditional approach under most snow states. At the snow-accumulation stage, the MAE and RMSE of the estimated and measured results decreased by 16.83% and 26.90%, respectively, while at the snow-accumulation stage, they decreased by 12.18% and 12.48%, respectively. At the snow-stabilization stage, they decreased by 9.01% and 12.53%, respectively. The accuracy of the two approaches in the stable-snow period is significantly higher than that in the changing-snow period, and the improved accuracy of the proposed approach is more obvious in the changing-snow period. These results are possibly because the physical parameters of the surface-such as the dielectric constant, snow-particle size and surface roughness-are different during different snow periods. The periods will have a certain effect on the received signal power, and the penetration of the L-band signal will also be different. According to a previous study [38], although there is a difference in the amplitude of the accumulating and melting trends, the effect of snow accumulation and snow melting on reflected GPS signals has similar pattern when compared to each other. Koch also pointed out that comparing the GPS-derived snow height with the validation data, the agreement for dry snow in accumulation period (RMSE: 0.13 m) and wet snow in melting period (RMSE: 0.14 m) is very good [39]. By contrast, under stable-snow condition, the snow depth derived from GNSS-IR can be more accurate. Although the observation noise is large in the snow-accumulation and ablation periods, wavelet decomposition shows good denoising effects and can improve the estimation in this period. During the stable-snow period, the physical characteristics of the snow surface are stable and minimally affected by the observation noise. The wavelet decomposition does not change the frequency composition in the original signal, so there is no obvious difference between the two approaches at this point. Upon comparing the estimation results of the two approaches for different years, it can be found that, except for 2015, the proposed approach significantly improves the accuracy when the snow changes. A possible reason is that the snow changes in 2015 were relatively slow, while the snow depths in other years experienced relatively large fluctuations. In summary, the proposed approach can significantly improve the accuracy of snow depth estimation when the snow depth changes drastically.

Effects of Rainfall on Estimated Snow Depth
Rainfall accelerates the melting or freezing of snow by changing its density and dielectric constant. In this study, the local rainfall data were used to analyze the influence of rainfall on GNSS-IR-based snow depth estimation for the first time. The meteorological station in Ny-Alesund, Svalbard, also provides rainfall observations, which makes the research on the effects of rainfall feasible. Figure 10 shows several typical examples of the impact of rain, where Figure 10a,b are cases without rainfall, and Figure 10c,d are cases with rainfall. A statistical comparison of the snow depth estimation with and without rainfall is shown in Table 4. feasible. Figure 10 shows several typical examples of the impact of rain, where Figure 10a,b are cases without rainfall, and Figure 10c,d are cases with rainfall. A statistical comparison of the snow depth estimation with and without rainfall is shown in Table 4.  From Figure 10 and Table 4, it can be observed that when the snow depth variation is stable, the estimations agree well with the field measurements. When there is no rainfall, the differences between the estimations and field measurements are smaller; however, when the rainfall changes greatly, the estimations fluctuate greatly and decrease in accuracy. Based on all the results over 5 years, it can be found that the MAEs of the estimation results without and with rainfall are 2.19 and 5.63 cm, and the RMSEs are 2.08 and 5.46 cm, respectively, which means that the estimation accuracy is significantly higher when there is no rainfall. A possible reason is that rainfall affects the snow's dielectric constant and physical characteristics, such as size and density. From Figure 10 and Table 4, it can be observed that when the snow depth variation is stable, the estimations agree well with the field measurements. When there is no rainfall, the differences between the estimations and field measurements are smaller; however, when the rainfall changes greatly, the estimations fluctuate greatly and decrease in accuracy. Based on all the results over 5 years, it can be found that the MAEs of the estimation results without and with rainfall are 2.19 and 5.63 cm, and the RMSEs are 2.08 and 5.46 cm, respectively, which means that the estimation accuracy is significantly higher when there is no rainfall. A possible reason is that rainfall affects the snow's dielectric constant and physical characteristics, such as size and density.

Conclusions
GNSS-IR is a new remote sensing technique for snow depth estimation. In this study, long-term GNSS observations from 2014 to 2018 in Ny-Alesund were collected and used to retrieve snow depths with GNSS-IR. First, an improved approach based on wavelet decomposition was proposed to improve the accuracy of snow depth estimation with GNSS-IR. Due to the good denoising effect of wavelet decomposition, the noise power and signal power can be accurately separated without changing the frequency composition in the original signal, which can significantly improve the accuracy and reliability of the estimations. The errors with the proposed method were lower, and the data utilization was better. The validity of the proposed approach is particularly obvious when the snow depth changes sharply. Second, the accuracy of snow depth estimation during the three periods of the accumulation, ablation and stabilization of snow cover is discussed in this paper. The results show that snow depth estimation is more accurate during the snow-stabilization period, and more complicated snow surfaces are present during the snow-accumulation period. The difference in estimation accuracy during different snow-cover periods is quite obvious. Finally, combined with the local meteorological observations at Ny-Alesund, the influence of rainfall on the accuracy of GNSS-IR estimation was studied for the first time. The results show that the accuracy of snow depth estimation without rainfall is significantly higher than that with rainfall.
With the development of global navigation satellite systems, future research can focus on combinations of observations at different frequencies and from different navigation satellite constellations. In addition, the snow water equivalent can be further calculated from snow depth variation once the snow density is determined. Future work could also focus on the influences of snow type, snow wetness, grain size and snow density on GNSS-IR snow depth estimations to improve the retrieval precision further.