Snow Density Retrieval in Quebec Using Space-Borne SMOS Observations

Abstract

Snow density varies spatially, temporally, and vertically within the snowpack and is the key to converting snow depth to snow water equivalent. While previous studies have demonstrated the feasibility of retrieving snow density using a multiple-angle L-band radiometer in theory and in ground-based radiometer experiments, this technique has not yet been applied to satellites. In this study, the snow density was retrieved using the Soil Moisture Ocean Salinity (SMOS) satellite radiometer observations at 43 stations in Quebec, Canada. We used a one-layer snow radiative transfer model and added a $\tau -\omega$ vegetation model over the snow to consider the forest influence. We developed an objective method to estimate the forest parameters ($\tau$, $\omega$) and soil roughness ($S_D$) from SMOS measurements during the snow-free period and applied them to estimate snow density. Prior knowledge of soil permittivity was used in the entire process, which was calculated from the Global Land Data Assimilation System (GLDAS) soil simulations using a frozen soil dielectric model. Results showed that the retrieved snow density had an overall root-mean-squared error (RMSE) of 83 kg/m³ for all stations, with a mean bias of 9.4 kg/m³. The RMSE can be further reduced if an artificial tuning of three predetermined parameters ($\tau$, $\omega$, and $S_D$) is allowed to reduce systematic biases at some stations. The remote sensing retrieved snow density outperforms the reanalysis snow density from GLDAS in terms of bias and temporal variation characteristics.

1. Introduction

Snow cover plays a critical role in terrestrial hydrological, climatological, and ecological processes. It influences the energy balance on the land surface, based on its high albedo and low thermal conductivity [1]. The measurement of snow water equivalent (SWE) is important to understand the timing and magnitude of snowmelt runoff [2]. However, when there have been multiple continental or large-scale snow depth remote sensing products, the corresponding snow density product required to convert snow depth to SWE is still limited. Snow density, SWE, and snow pressure are not as widely measured as snow depth at stations. When the passive microwave remote sensing technique at the Ku- to Ka-band was used to retrieve the snow depth for more than four decades [3], observations at these frequencies were rarely used to estimate snow density. Instead, after the snow depth is retrieved, the error in snow density auxiliary information becomes an important source of SWE uncertainty [4]. At a high frequency, snow volume scattering cannot be neglected. The influences of the snow depth and snow grain size on the brightness temperature (T_B) are much higher than that of snow density. The coherency between scattered waves from snow particles [5] results in a sophisticated relationship between T_B and snow density, which is differently understood by different snow radiative transfer models [6]. Therefore, to estimate snow density, the lower-frequency channel is better suited for its retrieval as it is more sensitive to snow refraction rather than snow volume scattering [7,8].

Snow density varies spatially, temporally, and vertically, influenced by the snow compaction rate and snow compaction time [9,10,11]. The use of a fixed snow density (for example, 240 kg/m³) will result in an overestimated SWE in the early snow season and an underestimated SWE in the late snow season [12,13,14]. The CCI+ Version2 (CCIv2) SWE product [13] utilized time-varying snow density functions according to snow classes [15]. The empirical method that calculates the snow density as a regression function of the snow depth, day of year, and snow classes [16] is unable to capture the spatial and inter-annual variability [17]. In a warming climate, the direct observation of snow density can be of great value to detect the increased occurrence of rain-on-snow events.

In the field, the most commonly used method to measure snow density is to weigh the snow in a container of a fixed volume [18]. Dielectric permittivity measurement [19], micro-computed tomography [20,21], and neutron-scattering probes [22] can also be used to measure snow density. After being inter-calibrated against the gravimetric measurements, these techniques are either fast in the field or of fine vertical resolution [23,24,25]. Snow density can also be estimated by L-, C-, and X-band synthetic aperture radar (SAR) data. The study in [26] utilized the change in surface scattering caused by snow, and estimated the snow density using a parameterized backscattering model simplified from the Integral Equation Model (IEM) [27]. In recent years, with the development of full-polarized SAR, surface [28,29] and volume scattering [30,31] components were extracted from different signal decomposition methods to establish empirical or physical relationships with snow density.

After the launch of the Soil Moisture Ocean Salinity (SMOS) [32] satellite for soil moisture observation purposes in November 2009, the multiple-angle L-band radiometer became an additional type of sensor applicable for snow density estimation. The revisit frequency of SMOS is two to three days, which makes it suitable for studying temporal snow density variation and supporting daily SWE products. Preliminary studies on the physical [8] and the practical feasibility [7] of snow density retrieval using ground-based radiometer observations have been conducted. It was found that the existence of snow on the bare soil changes both the propagation angle of microwaves inside the snow medium and the refraction at the air–snow and snow–soil boundaries. With the increase in snow density, the propagation angle becomes steeper, the emissivity at the horizontal polarization increases, and the emissivity at the vertical polarization increases at small incidence angles and decreases at large incidence angles. Because of the large contrast between the L-band wavelength and snow grain size, volume scattering generated by snow particles can be neglected [7,8,26]. This means that the computational cost to calculate snow T_B and retrieve the snow density can be largely reduced. Using the simplified snow radiative transfer model, the studies in [7,33,34] successfully estimated the snow density in tower-based radiometric experiments. However, if the snow density retrieval algorithm is moved from ground-based to space-borne instruments, it will be challenged by stronger observation noise, emission from other land surface types, and the Radio Frequency Interference (RFI) from human activities. Whether the snow density is retrievable becomes uncertain.

Therefore, the scientific goal of this article is to retrieve the snow density from SMOS satellite observations, according to the theory in [7,8]. The study area is Quebec, Canada, determined by the availability of ground measurements and the appropriate terrain condition. Forests in Quebec should be considered using space-borne observations in theory. To minimize the impact of forests during snow density retrieval, the τ-ω vegetation model was utilized. The prior knowledge of soil permittivity was introduced and calculated from the Global Land Data Assimilation System (GLDAS) simulated soil temperature and total soil water content [35,36] using a frozen soil dielectric model [37]. We also developed a method based on objective optimization to determine $\omega$, $\tau$, and soil roughness ($S_D$) using the SMOS observations in the snow-free period and applied them subsequently to estimate the snow density in the snow season. With the predetermined parameters and the prior knowledge, the ill-posed problem of snow density retrieval was greatly relieved.

This paper is organized as follows. Section 2 introduces the study area and datasets utilized. Section 3 describes the forward emission model and the retrieval algorithm. Section 4 shows the snow density retrieval results. Section 5 and Section 6 are the Discussion and Conclusions, respectively.

2. Study Area and Datasets

2.1. Study Area

The study area is located in Quebec, Eastern Canada (45°N–63°N, 57°W–80°W), with elevation above sea level from −3 to 1081 m (Figure 1). Characterized by cool temperatures in summer and abundant snowfall in winter, this area includes eight different terrestrial ecozones [38]. Snow in this region has large spatial variability, with a snow cover duration ranging from 120 days in Southern Quebec to 240 days in Northern Quebec, and an annual maximum SWE from less than 100 mm at low altitudes to more than 300 mm at high altitudes [39,40].

2.2. Datasets

2.2.1. The SMOS Brightness Temperature Dataset

The SMOS mission launched by the European Space Agency [32] provides L-band T_B at multiple observation angles and dual polarizations. The SMOS Level 1C (L1C) product is provided in the ISEA-4H9 grid (icosahedral Snyder equal area earth fixed) with a spatial resolution of 43 km on average [41]. However, the observations in some areas are affected by the RFI. We employed a two-step regression approach, developed in [42], to smooth the angular dependence of the SMOS L1C T_B product. The method can largely remove outliers and reduce the observation uncertainty, which was validated in [43]. In the post-processed product, the SMOS T_B was resampled to fixed incidence angles at 40° and from 2.5° to 62.5° with an interval of 5°. Observations at all angles from ascending orbits with an overpass time of 6:00 AM were used for snow density retrieval. The study period was from October 2019 to June 2020.

2.2.2. In-Situ Snow Measurements

The in-situ snow measurements come from the Canadian historical Snow Water Equivalent dataset (CanSWEv3.0) [44,45], which contains both the snow depth and snow water equivalent at the same station. The entire dataset has measurements from 2832 stations in the 1928–2021 period. In our study area and period, there are 43 stations with elevation < 1000 m. According to the International Geosphere-Biosphere Programme (IGBP) classification acquired from the MODIS MCD12Q1 product [46], stations in the south are located in highly forest-covered areas, including evergreen needleleaf forests, mixed forests, and woody savannas, whereas the remaining are mostly located in the savanna region. Stations falling within grids with a dominant IGBP type of water bodies were removed.

2.2.3. Other Auxiliary Datasets

To calculate the ground permittivity, we introduced the simulated total soil water content and soil physical temperature at 0–10 cm from the Global Land Data Assimilation System (GLDAS) Noah model product [35,36]. A soil dielectric model was used to calculate the frozen soil permittivity [37], and the required soil texture (sand, silt, and clay contents) and soil bulk density were extracted from the Harmonized World Soil Database [47]. In the T_B calculation, the vegetation temperature was considered the same as the air temperature from GLDAS, whereas the snow temperature was calculated as the mean of the GLDAS soil temperature and air temperature. The snow temperature has a very small influence on the imaginary part of snow permittivity. To calculate the influence of the forest canopy, we introduced the forest cover fraction from the MODIS MOD44B product [48].

3. Methods

3.1. Forward Emission Model

The radiative transfer model used to describe the emission of the soil–snow–vegetation system is an empirical rough soil reflectivity model [49], coupled with a simplified snow emission model neglecting absorption and scattering coefficients [7] and a $\tau -\omega$ vegetation model [50]. For an SMOS grid with a forest cover fraction of $F_{C }$, the measured T_B can be written as in Equation (1), which includes an open snow component (of areal fraction of $1-F_{C }$) and a forest-covered snow component (of areal fraction of $F_{C }$):

$$T_B^p{=T}_{B,f}^p F_C{+T}_{B,S}^p (1-F_{C })$$

(1)

where $T_B^p$ is the satellite-observed T_B at polarization $p$ ($p$= H or V), $T_{B,S}^p$ is the open snow T_B, and $T_{B,f}^p$ is the snow-emitted T_B observed above the forest canopy. The SMOS grid is considered fully snow-covered.

$T_{B,f}^p$ is expressed as a three-component model as [50]:

$$T_{B,f}^p{=T}_{B,S}^p{ \gamma +T}_{C }\left(1-\omega \right) \left(1-\gamma \right){+T}_{\mathrm{C}} \left(1-\omega \right) \left(1-\gamma \right) r^p \gamma$$

(2)

where $T_{C }$ is the vegetation physical temperature. $\omega$ is the effective scattering albedo accounting for the volume scattering of the canopy [51], and $\gamma$ is the vegetation attenuation factor. The first component is the snow emission attenuated by the vegetation canopy, the second component is the thermal emission from the forest propagated upward, and the third is the thermal emission from the forest propagated downward, reflected by the snow and attenuated by the vegetation canopy. $\omega$ and $\gamma$ are assumed independent of polarization in this study. $\gamma$ is expressed as $\gamma =\mathit{exp} (-\tau \mathit{sec}\theta )$, where $\tau$ is the vegetation optical depth (VOD) and $\theta$ is the observation angle. $r^p$ is the effective surface reflectivity of the snow–soil system.

Due to the low-loss characteristic of snow cover in the L-band [52], the simulation of ${ T}_{B,S}^p$ in this study utilized the simplified snow emission model neglecting snow absorption and volume scattering, as [7]

$$T_{B,S}^p{=a}_G^p T_G{+a}_S^p T_S{+a}_{sky}^p T_{sky}$$

(3)

where $T_G$, $T_S$, and $T_{sky}$ are the soil physical temperature, snow physical temperature, and downwelling sky T_B, respectively. If the snow is homogeneous, the Kirchhoff coefficients $a_G^p$, $a_S^p$, and $a_{sky}^p$ can be expressed as

$$\begin{gathered} a_G^p=\frac{(1-s_G^p) (1-s_S^p{) t}_S}{{1+r}_S^2{ s}_G^p s_S^p-{ r}_{S }{(s}_G^p{+s}_S^p)-t_S^2{ s}_G^p{ s}_S^p} \\ {a}_{S }^p=\frac{(1-s_S^p) (1-r_S-t_S) (1-r_S{s}_G^p{+t}_S{s}_G^p)}{(1-r_S{s}_G^p) (1-r_S{s}_S^p)-t_{S }^2{s}_G^p{ s}_S^p} \\ {a}_{sky}^p=\frac{(1-r_S{s}_G^p{) (s}_S^p{+r}_S (1-{2s}_S^p{\left)\right)+s}_G^p (1-{2s}_S^p{) t}_S^2}{(1-r_S{s}_G^p) (1-r_S{s}_S^p)-{ t}_S^2{ s}_G^p s_S^p} \end{gathered}$$

(4)

where $s_G^p$ is the snow–soil interface reflectivity, $s_S^p$ is the specular air–snow interface reflectivity, and $r_{S }$ and $t_S$ are the reflectivity and transmissivity of the snow layer, respectively.

When there is no absorption, no volume scattering assumptions are made for snow, $t_{S }$= 1 and $r_{S }$= 0. The Kirchhoff coefficients can be simplified to

$$\begin{gathered} a_G^p=\frac{(1-s_G^p) (1-s_S^p)}{1-s_G^p s_S^p} \\ {a}_{S }^p=0 \\ {a}_{sky }^p=1-a_G^p \end{gathered}$$

(5)

where $a_G^p$ can be used as $1-r^p$ in Equation (2).

We evaluated the reasonability of neglecting absorption and volume scattering at the L-band using the Microwave Emission Model of Layered Snowpacks (MEMLS) based on the Improved Born Approximation (IBA) [53,54], and found that the sensitivity of the snow-emitted T_B to snow depth is only 0.018–0.072 K per m, using an exponential correlation length of 0.18 mm.

To calculate $s_G^p$, a semi-empirical QHN model [55,56] with frequent-independent coefficients [49] was used (Equations (6) and (7)). It allows the use of one parameter, the standard deviation of height ($S_D$), to model the rough soil surface reflectivity, which is easier to fit than models using two parameters:

$$s_G^p=\left(\left(1-Q_R\right) s_G^{p\ast }{+Q}_R s_G^{q\ast }\right) \mathit{exp} (-H_R{\mathit{cos}}^{N_{\mathrm{Rp}}}\theta )$$

(6)

$$H_R=(\frac{a_1{S}_D}{a_2{S}_D{+a}_3}{)}^6$$

(7)

where $s_G^{p\ast }$ and $s_G^{q\ast }$ are calculated from the soil permittivity using Fresnel equations [57] at polarization $p$ and $q$ (e.g., $q$ = V if $p$ = H and vice versa). The soil dielectric model developed in [37] was used to calculate the soil permittivity. The remaining parameters were optimized as $Q_R$ = 0.075, $N_{RV}$ = 1.503, $N_{RH}$ = 0.131, and $H_R$ is a function of $S_D$, with coefficients as $a_1$ = 0.887, $a_2$ = 0.796, $a_3$ = 3.517 [49].

3.2. Retrieval of Predetermined Parameters ($\tau$, $\omega$, $S_D$) in Snow-Free Period

To reduce the difficulty in snow density estimation, three unknown parameters ($\tau$, $\omega$, $S_D$) were predetermined based on SMOS observations in the snow-free period. We assumed that the values of the three parameters remain constant over time and are independent of polarization. Using the soil permittivity from GLDAS, we calculated the error between the simulated and the observed SMOS T_B during the snow-free period, using different combinations of $\tau$, $\omega$, and $S_D$.

$$T_{B, \mathit{error}}\left({\tau , \omega , S}_D\right)={{\displaystyle \sum }}_{t=1}^{t_n}{{\displaystyle \sum }}_{{\theta }_k=\theta }^{{\mathsf{\theta }}_n}{{\displaystyle \sum }}_{p=H,V}{\left(T_{B, \mathit{obs}}^p\left({\theta }_k,t\right)-T_{B, \mathit{mod}}^p\left({\theta }_k{, t, \tau , \omega , S}_D\right)\right)}^2$$

(8)

$T_{B, \mathit{error}}$ was calculated from days within 2 weeks before the onset or after the end of the snow season. The searching ranges of $\tau$, $\omega$, and $S_D$ were $\tau$ from 0 to 0.5 with a step of 0.01, $\omega$ from 0 to 0.4 with a step of 0.01, and $S_D$ from 0 to 100 mm with a step of 1 mm. It was found that a few combinations of $\tau$, $\omega$, and $S_D$ gave similarly small $T_{B, \mathit{error}}$. We extracted the first 0.1% smallest $T_{B, \mathit{error}}$ and found the corresponding $\tau$, $\omega$, and $S_D$ combinations (see Figure 2 as an example for the HQ-CM4L station). Later, the values of these $\tau$, $\omega$, and $S_D$ combinations were averaged, and the single combination closest to the averaged value was chosen as the final estimate.

3.3. Retrieval of Snow Density

After determining the three unknown parameters, it is possible to calculate soil and vegetation emissions. Subsequently, the snow density can be retrieved using a cost function as

$${ \mathit{CF}(\rho }_S)={{\displaystyle \sum }}_{{\theta }_k=\theta }^{{\theta }_n}{{\displaystyle \sum }}_{p=H,V}{\left(T_{B, \mathit{obs}}^p\left({\theta }_k\right)-T_{B, \mathit{mod}}^p\left({\theta }_k{, \rho }_S \right)\right)}^2$$

(9)

The searching range for snow density (${\rho }_S$) was 50 to 500 kg/m³ with a step of 1 kg/m³. ${\rho }_S$ with the smallest ${\mathit{CF}(\rho }_S)$ will be considered as the retrieved snow density. Figure 3 is the flow chart for the entire snow density retrieval process.

3.4. Objective Postprocessing Method to Reduce Retrieval Uncertainty

Snow density retrieval is sensitive to the predetermined parameters $\tau$, $\omega$, and $S_D$. In our study, we found that using SMOS measurements from different periods to fit $\tau$, $\omega$, and $S_D$ can lead to different results. Better performance can be achieved by using SMOS measurements before or after the snow season, which varies by station. This is due to changes in bias in the simulated snow-free T_B using the same GLDAS dataset and predetermined parameters, with the impacting factors unknown. To reduce uncertainty and ensure objectivity, we generated three sets of $\tau$, $\omega$, and $S_D$ fitted from before, after, and before and after the snow season, and we generated three sets of snow density time series as candidates. Abnormal snow density time series can be recognized by counting the occurrence of abnormal daily snow density values reaching the 50 or 500 kg/m³ searching boundary limit. The remaining snow density time series will be averaged to obtain the final retrieval result.

4. Results

4.1. Performance of Multiple-Angle Brightness Temperature Simulation

To examine the behavior of the forward model in matching the SMOS observations, Figure 4 shows examples of simulated and observed T_B for one station during the snow density retrieval. It shows that for the HQ-CM4E station, the simulated T_B matched well with the SMOS-observed T_B, except for large incidence angles at horizontal polarization and unstable snow conditions (Figure 4a). On 29 November 2019 (Figure 4a), the relatively low-biased T_B at small angles and the complex snow condition during this season resulted in an underestimation of snow density. However, in Figure 4b–f, the errors between the retrieved snow density and in-situ measurements were within 45 kg/m³. For example, on 6 December 2019, the estimated snow density was 155 kg/m³ with an error of 3.5 kg/m³, and on 16 February 2020, the estimated snow density was 206 kg/m³ with an error of 5.9 kg/m³. On 16 February 2020 (Figure 4d), the increased snow density compared to previous months was retrieved from an increased T_B in small incidence angles, a decreased T_B in large angles at vertical polarization, and an increased T_B at horizontal polarization. On 8 April 2020 (Figure 4f), the model-simulated soil permittivity and the retrieved snow density significantly increased compared to previous examples, which may have been caused by melt–refreeze events in this season.

4.2. Performance of Snow Density Retrieval at Example Stations

Figure 5 and Figure 6 show the time series of retrieved snow density using the retrieval method described in Section 3.2 and Section 3.3 at six different stations, and we compared them with the CanSWE ground measurements and reanalysis snow density from GLDAS. The GLDAS product provides the snow depth and snow water equivalent, which allows us to calculate the average snow density along the snow profile. Among 43 stations, some have better retrieval performance than others. Therefore, we have selected and presented three examples with good performance in Figure 5 and three examples with poor performance in Figure 6.

As shown in Figure 5, the retrieved snow density captured the increasing trend of the observed snow density from 100~200 kg/m³ in the early and mid- snow season to 300~400 kg/m³ in the late snow season. The root-mean-squared error (RMSE) varies from 52.7 kg/m³ to 67.9 kg/m³ for different stations. The mean bias is smaller than 25 kg/m³, which is approximately 10% of the mean snow density. If we only consider the retrievals during relatively stable snow conditions (December to March), the RMSE can be reduced to 36.7~42.1 kg/m³, but the changes in R and bias are complex, because the snow density variation range is reduced. In the late snow season, the remote sensing retrieved snow density is noisier than the in-situ measurements. However, in the early and mid-snow season, the GLDAS reanalysis dataset underestimates the snow density compared to the measurements. Therefore, the SMOS-retrieved snow density outperforms the GLDAS reanalysis in the time series variation characteristic. In the late snow season after April, the snow density increases due to melt–refreeze events and snow compaction, and it decreases due to new snowfalls. The observed, SMOS-retrieved, and GLDAS snow density fluctuate at the same pace but with different magnitudes. The shape of the SMOS-retrieved snow density is closer to that of GLDAS, whereas the variation in the observed snow density is smaller in general.

As shown in Figure 6, the snow density retrieval results with poor performance show systematic biases, with a mean bias ranging from −77.5 kg/m³ to 108.7 kg/m³. However, although the RMSE increases to 90.9~131.6 kg/m³, the unbiased RMSE (ubRMSE) is within 47.5~74.2 kg/m³, which is comparable to the results in Figure 5. The ubRMSE can be used to remove the contribution of systematic bias to RMSE. The small ubRMSEs imply that the variation trend of the observed and the retrieved snow density is still consistent. We found that all three stations in Figure 6 are located close to water bodies, which may be frozen or unfrozen in winter depending on the air temperature and the size of the lake. Additionally, it is possible that the lake was unfrozen when the vegetation and soil roughness parameters were fitted in the snow-free season, but frozen when we were estimating the snow density. In this case, the influence of the lake on snow density retrieval is quite complex and will need to be explored in future research. We also found that the SMOS T_B at the stations in Figure 6a,b is similar to that of the stations in Figure 5. However, the SMOS T_B at the station in Figure 6c is much lower than that of the stations in Figure 5, for approximately 20–30 K; therefore, the retrievals were underestimated instead of overestimated. When the retrieved snow density shows positive or negative biases at the three stations near the lakes, the GLDAS snow density consistently biases low in the early and mid-snow season.

4.3. Validation of Retrieved Snow Density at All Stations

In Figure 7, the snow density retrieval performance is summarized from 43 stations to test the robustness of our retrieval algorithm. Figure 7a,b show the scatterplots of the SMOS-retrieved and GLDAS snow density against the measured snow density, respectively. It shows, from all validation stations, that the retrieved snow density has an RMSE of 83 kg/m³ and a Pearson correlation coefficient of 0.5 with measurements, whereas the GLDAS snow density has a lower RMSE of 76 kg/m³ and a higher correlation of 0.7. However, the GLDAS snow density has an overall underestimation of 48 kg/m³, and the absolute bias is significantly higher than that of the retrieved snow density (9.4 kg/m³). If only the December–March period is summarized, the RMSE of the retrieved snow density improves to 72.3 kg/m³, which is comparable to that of GLDAS at 71.9 kg/m³. Near 200 kg/m³, GLDAS underestimated the snow density, whereas our algorithm overestimated the snow density. The underestimation of GLDAS is from the early and mid-snow season, as shown in Figure 5 and Figure 6. It implies that although the retrieved snow density contains larger noise than the reanalysis, it still has a more trustworthy temporal variation trend and a generally unbiased characteristic. Figure 7c,d show that both the retrieved and the GLDAS snow density biases are insensitive to the snow depth.

Figure 8 shows the distribution of validation metrics across each station on a geographic map. It shows that stations located in the savannas in the north of Quebec tend to have higher correlations with measurements and lower ubRMSEs. Conversely, stations in forests and woody savannas have lower correlation coefficients and higher ubRMSEs. However, in savannas, some stations exhibit strong overestimation and large RMSEs, with seven out of nine located near large water bodies. Therefore, as summarized in Figure 9 and

Table 1. The validation metrics of retrieved snow density compared with in-situ measurements summarized based on the MCD12Q1 IGBP classification.

MCD12Q1 IGBP Classification	R		Bias (kg/m3)		RMSE (kg/m3)		ubRMSE (kg/m3)		Measurement Number		Station Number
	October – May	December – March	October – May	December – March	October – May	December – March	October – May	December – March	October – May	December – March
evergreen needleleaf forest	0.35	−0.19	−33.44	−12.95	78.55	52.97	71.08	51.36	93	62	1
woody savannas	0.55	0.39	−16.81	22.27	85.81	77.07	84.15	73.78	1011	680	4
mixed forest	0.47	0.48	12.5	3.67	76.59	49.53	75.56	49.39	393	258	10
savannas	0.5	0.25	−11.37	20.01	82.8	73.76	82.01	70.99	2664	1816	28
ALL SITES	0.5	0.29	9.44	18.33	82.89	72.31	82.35	69.95	4161	2816	43

, the average values of the validation metrics for different land cover types are similar, because both forested and non-forested regions have their own challenges. More statistics on snow density retrieval accuracy at each station are provided in Appendix A.

5. Discussion

In this paper, the snow density retrieval accuracy is influenced by the following factors. Firstly, soil permittivity is a key parameter that affects the L-band multiple-angle T_B. However, our retrieval algorithm uses the soil permittivity calculated from the reanalysis dataset as a known variable to support snow density estimation. Secondly, the vegetation and soil roughness parameters ($\tau$, $\omega$, $S_D$) were assumed the same in snow-free and snow-covered periods. However, at least $\tau$ and $\omega$ were found to vary with the air temperature in subzero conditions, especially at low frequencies [58]. Thirdly, natural snow is a layered medium. However, we neglected the possible vertical variation in snow density in the snow profile, although the study in [7] mentioned that the retrieved snow density may be closer to the bottom snow density instead of the average. Fourthly, factors such as wet snow can make snow density retrieval impossible or unreliable. The ice lens and snow crust formed during the freeze–thaw process can complicate the snow stratigraphy, increase the refraction inside the snow, and make the current retrieval algorithm less accurate. These factors were not considered in this paper, as simplifications were made to make snow density retrievable from satellites. Additionally, the snow density retrieval may be influenced by SMOS observation errors, footprint differences in different incidence angles, lakes (frozen or unfrozen), ice, and complex landscapes, etc.

In-situ measurements used to calculate the error of our retrieval algorithm were from the point scale, with a single site inside each SMOS grid. Therefore, differences in scale between the retrieved and the measured snow density may result in an overestimation of the uncertainty of our algorithm.

In the previous section, we noted that the retrieved snow density showed stronger biases at some stations. Figure 10 shows that the systematic bias can be corrected if the three predetermined parameters ($\tau$, $\omega$, $S_D$) are artificially tuned to match the measured snow density. However, this result assumes that the target snow density is provided. As shown in Figure 10, the bias can be reduced to below 15 kg/m³ at all three stations.

To analyze the GLDAS snow density bias compared to the measurements, we examined the relationships between GLDAS and the in-situ snow depth (SD) and snow water equivalent (SWE), as seen in Figure 11a,b, respectively. The figures reveal that many GLDAS SD values between 50 and 100 cm are unbiased but the corresponding SWEs are underestimated. Additionally, some SDs in the deep snow range are seriously underestimated. Figure 11c provides an example of one station to explain this phenomenon. GLDAS SD coincides well with in-situ measurements in the early snow season and begins to underestimate after mid-January due to the underestimated snowfall amount in GLDAS before the snowmelt. After the snowmelt onset (which occurs almost simultaneously for GLDAS and the station in mid-March), the GLDAS snowpack melts faster than the station, because it has less SWE to melt and the single-layer snow scheme utilized in the Noah model overestimates snowmelt speed [59]. Melt–refreeze events can dramatically increase the snow density. Therefore, during the snowmelt period, GLDAS overestimates the snow density compared to in-situ measurements.

6. Conclusions

This study conducted snow density retrieval experiments based on L-band multiple-angle SMOS satellite observations and compared the results with the in-situ measurements from 43 CanSWE stations in Quebec, Canada. A forward model was used to describe the emission of the soil–snow–vegetation system. The vegetation and soil roughness parameters were objectively determined using SMOS T_B in the snow-free period and applied to estimate the snow density. The new retrieval method achieved bias of 9.4 kg/m³ and an RMSE of 83 kg/m³ for snow density at all stations. Currently, some stations show large systematic biases, but these biases can be reduced.

The need for prior information about soil permittivity could be a challenge for satellite-based snow density retrieval. Similarly, to estimate the unfrozen soil water content under the snowpack, a previous study [60] introduced the snow density from the Snow Data Assimilation System (SNODAS) reanalysis dataset [61]. Whether we can develop a method to simultaneously estimate the snow density and soil permittivity from satellites, as with ground-based experiments, remains to be seen.

This paper extended the snow density retrieval from a ground-based radiometer to satellite application and established a method to solve the forest influence problem at coarse resolutions. Although the satellite-retrieved snow density is noisier than that of GLDAS, it is generally less biased and shows better time series variation characteristics. Thus, estimating the snow density from satellites using SMOS has scientific value. We anticipate a lower RMSE for SMOS-retrieved snow density if the vegetation and soil roughness parameters can be better estimated and the water body influence can be properly treated.