Modelling the L-Band Snow-Covered Surface Emission in a Winter Canadian Prairie Environment

: Detailed angular ground-based L-band brightness temperature ( T B ) measurements over snow covered frozen soil in a prairie environment were used to parameterize and evaluate an electromagnetic model, the Wave Approach for LOw-frequency MIcrowave emission in Snow (WALOMIS), for seasonal snow. WALOMIS, initially developed for Antarctic applications, was extended with a soil interface model. A Gaussian noise on snow layer thickness was implemented to account for natural variability and thus improve the T B simulations compared to observations. The model performance was compared with two radiative transfer models, the Dense Media Radiative Transfer-Multi Layer incoherent model (DMRT-ML) and a version of the Microwave Emission Model for Layered Snowpacks (MEMLS) adapted speciﬁcally for use at L-band in the original one-layer conﬁguration (LS-MEMLS-1L). Angular radiometer measurements (30 ◦ , 40 ◦ , 50 ◦ , and 60 ◦ ) were acquired at six snow pits. The root-mean-square error (RMSE) between simulated and measured T B at vertical and horizontal polarizations were similar for the three models, with overall RMSE between 7.2 and 10.5 K. However, WALOMIS and DMRT-ML were able to better reproduce the observed T B at higher incidence angles (50 ◦ and 60 ◦ ) and at horizontal polarization. The similar results obtained between WALOMIS and DMRT-ML suggests that the interference phenomena are weak in the case of shallow seasonal snow despite the presence of visible layers with thicknesses smaller than the wavelength, and the radiative transfer model can thus be used to compute L-band brightness temperature. measurements as well as the three snow microwave emission models and the soil emission model. The model parameterizations are presented, after which we present results and a comparison of the model performance.


Site and Data
During the 2014-2015 winter, a ground-based L-band radiometer measurement campaign was conducted at the Kernen Crop Research Farm (KCRF; 52.149 • N; 106.545 • W), a 380 ha property within the city of Saskatoon owned and operated by the University of Saskatchewan, Canada. L-band radiometer measurements and coincident snow pit and meteorological observations were performed. The study area and the campaign and datasets are described in detail in Reference [17].
At KCRF, tree scenes were located adjacent to each other within the same field to ensure similarity in background emission, with only the overlying snow conditions altered. The three scenes include: Scene 1-Undisturbed snow: A scene of naturally accumulating snow-covered ground; Scene 2-Snow free: snow was removed on a weekly basis to maintain bare ground; Scene 3-Artificially compacted snow: a scene with deep and dense snow. Additional snow was manually added to Scene 3 and compacted on December 10th, 19th, 2014 and January 11th 2015 and then left to evolve naturally for the rest of the season. As this study focuses on snow emission modelling, Scene 2 was not used in this study. The scenes were characterized by silt-loam bare soil conditions. Wheat residue was noted and not disturbed during the study. Surface roughness was derived using a terrestrial Light Detection and Ranging (LiDAR) system using a surface roughness tool called Roughness from Point Cloud Profiles (RPCP) [26] implemented in Whitebox Geospatial Analysis Tools (GAT) software [27]. Surface roughness had a root-mean-square height (RMSH) of 1.78 cm and 1.64 cm within Scene 1 and 3, respectively, at the beginning of the study and remained almost unchanged throughout the study (RMSH = 1.79 cm and 1.75 cm).
L-band measurements were acquired by a surface-based hyperspectral dual polarization L-band Fourier transform radio-frequency interference (RFI) detecting radiometer with 385 channels designed for a frequency range from 1400 MHz to ≈1550 MHz. The radiometer antenna is a 19-element air loaded conformal muffin tin design that has a 30 • half-power (−3 dB) beamwidth. A method was developed for separating out the thermal spectrum from RFI-contaminated channels to get unique RFI-free T B from the measured spectrum [28]. Only the protected radio-astronomy frequency spectrum of 1400-1427 MHz was used to calculate the T B . The radiometer was set 2.75 m above the surface, and measurements at the angles 30 • , 40 • , 50 • , and 60 • relative to nadir were taken of the three scenes on a weekly basis. On 9 November 2014, radiometric measurements were taken while the soil was frozen and snow-free. From December 2014 to March 2015, six radiometric measurements were taken, coincident with manual snow pit measurements in the vicinity of Scene 1 and 3. The snow pits included documenting the snow stratigraphy, including the presence of ice lenses. Profiles of snow temperature and snow density were taken for the observed snow layers. Mass density was measured using a 100 cm 3 density cutter, and samples were weighed with a digital scale with an accuracy of ±0.1 g. The snow and soil temperature at 2.5 cm intervals were measured with a digital temperature probe (±0.1 • C). Soil was frozen at each visit. Figure 1 shows snow pits performed close to Scene 1 and 3 during each visit. Note that on 7 December 2014, only a single snow pit is available and refers to both scenes. Snow pits in the vicinity of Scene 1 were generally shallow ( Table 1) and composed of a depth hoar layer at the bottom and a high-density rounded grain winds slab snow layer at the surface. One or two high-density melt/ice crust layers and/or ice lenses were present within the snowpack, resulting from mid-winter melt events (see in Reference [17] the Figure 4 and details). Note that this strong stratification between the top and bottom of the snowpack made the snow density measurements a challenge because of the hardness (surface wind slab and melt/ice crust) and the instability (depth hoar) of the snow layers. Because of the artificial compaction of snow, the snow density of the bottom layer is higher in Scene 3. There was still a high snow density observed in February showing that the artificial high-density snow/ice crust made up a large proportion of the lower 10 cm of the snowpack in Scene 3. However, as the season progressed and metamorphism continued within the snowpack, there was a decrease in the density of the snow/ice crust layers found within the bottom layers. Note that all air temperature measurements below −6 • C ensure that the snow was dry during each visit (Table 1).

Emission Models
All three snow microwave emission models and the soil emission model used in this study are already well described in detail (see previously provided references). Accordingly, we only recall here the principal components of each model, the model inputs and adjustments made for this study.

WALOMIS
The WALOMIS [18] coherent snow emission model is based on a wave approach, i.e., solving Maxwell's equation for a multi-layered medium [29,30]. Each layer is characterized by thickness, temperature, and density. The most important simplification in this model is to neglect scattering by snow grains. This assumption is invalid for high microwave frequencies, however, in the case of Lband, scattering by grains is insignificant in comparison with absorption and reflection at the interfaces between layers due to the L-band wavelength being several orders of magnitude larger than snow grain size. Under these assumptions, the vertically and horizontally polarized TB of a given snowpack is calculated with the propagation-matrix derived from Reference [31].
WALOMIS was initially implemented to investigate the microwave emission at L-band for semiinfinite snow-firn in Antarctica [18,25]. In the case of the Antarctic ice-sheet, the soil emission can be ignored because of the high ice thickness (>1000 m). Thus, the lowest layer of the model is considered a semi-infinite ice layer. In contrast, in the case of seasonal snowpack, the soil emission is not negligible for the total emission of the snow-covered ground. Therefore, for the present study WALOMIS was adapted to take into account the soil emission from below the snowpack replacing the semi-infinite bottom ice layer by a soil layer characterized by the observed temperature and permittivity.
Because of the high sensitivity of the interference phenomena to the layer thickness with the wave approach (which is not the case with the non-coherent radiative transfer approach), the result obtained for a specific snowpack configuration (i.e., a given set of inputs) may differ considerably

Emission Models
All three snow microwave emission models and the soil emission model used in this study are already well described in detail (see previously provided references). Accordingly, we only recall here the principal components of each model, the model inputs and adjustments made for this study.

WALOMIS
The WALOMIS [18] coherent snow emission model is based on a wave approach, i.e., solving Maxwell's equation for a multi-layered medium [29,30]. Each layer is characterized by thickness, temperature, and density. The most important simplification in this model is to neglect scattering by snow grains. This assumption is invalid for high microwave frequencies, however, in the case of L-band, scattering by grains is insignificant in comparison with absorption and reflection at the interfaces between layers due to the L-band wavelength being several orders of magnitude larger than snow grain size. Under these assumptions, the vertically and horizontally polarized T B of a given snowpack is calculated with the propagation-matrix derived from Reference [31].
WALOMIS was initially implemented to investigate the microwave emission at L-band for semi-infinite snow-firn in Antarctica [18,25]. In the case of the Antarctic ice-sheet, the soil emission can be ignored because of the high ice thickness (>1000 m). Thus, the lowest layer of the model is considered a semi-infinite ice layer. In contrast, in the case of seasonal snowpack, the soil emission is not negligible for the total emission of the snow-covered ground. Therefore, for the present study WALOMIS was adapted to take into account the soil emission from below the snowpack replacing the semi-infinite bottom ice layer by a soil layer characterized by the observed temperature and permittivity.
Because of the high sensitivity of the interference phenomena to the layer thickness with the wave approach (which is not the case with the non-coherent radiative transfer approach), the result obtained for a specific snowpack configuration (i.e., a given set of inputs) may differ considerably from those Remote Sens. 2018, 10, 1451 5 of 15 obtained with a slightly different snowpack. To account for the variable nature and the imperfect layering of the snowpack within the footprint of the radiometer, it is essential to average a large number of simulations using inputs that represent natural snow variability. As thousands of simulations are required, it would be impossible to obtain the input profiles from direct measurements. Several studies suggested stochastic methods to generate such profiles from measurements in Antarctica (e.g., Reference [18,29,32]). A similar procedure used here is described in Section 4.2. The output of the model is the average T B from all the generated profiles.

DMRT-ML
The DMRT-MultiLayer (DMRT-ML) is an incoherent model that describes the snowpack as a multilayer medium, where each snow layer is characterized by its thickness, temperature, density, grain optical radius, stickiness parameter, and liquid water content. The model is available from http://gp.snow-physics.science/dmrtml. It is based on the DMRT theory [30]. In this study, stickiness is not investigated because scattering by grains is negligible, and this parameter has no effect at L-band (typically less than 0.1 K). Because all the measurements were made in cold conditions with dry snow, the liquid water content was considered to be zero. For each layer, the effective dielectric constant is represented using the first order quasi-crystalline approximation and the Percus-Yevik approximation for spherical grains. The absorption and scattering coefficients are calculated assuming a medium of "ice spheres in air background" and the emission and propagation of radiation through the snowpack are computed using the Discrete Ordinate Method (DISORT: [33]) with 64 streams, which takes multiple scattering between the layers into account, but not the interferences.

LS-MEMLS-1L
The LS-MEMLS [19] model estimates L-band microwave emission from a ground surface covered by a layer of dry snow. This emission model is based on parts of MEMLS, [12] with the assumptions of no absorption and no volume scattering in dry snow, which are applicable to the L-band frequencies in dry snow. Once the interface reflectivities are known, the Kirchhoff coefficients associated with a single (snow) layer above an infinite half-space (ground) are computed to derive T B . Snow is characterized only by its permittivity, controlled by the dry snow density. Schwank et al. [8] assumed a single snow layer with a homogeneous density distribution, which allowed a numerical inversion of the model with minimal a priori information, for purposes of retrieval of snow and ground parameters. Although e.g., Reference [34] applied the model also in a configuration exhibiting a vertical distribution of snow densities, in this study LS-MEMLS is applied in the original one-layer configuration (LS-MEMLS-1L) to evaluate its applicability for snow density retrievals [8,9].

Soil Emission Model
At L-band, soil emission has a significant contribution to the signal emerging from the surface in environments with seasonal snow [35]. Hence, a soil reflectivity model is a prominent component of seasonal snow microwave-emission models. In this study, the soil reflectivity is calculated from the Fresnel equations and the roughness is considered as negligible. A specular soil reflectivity model is used in this study because WALOMIS needs electric field reflectivity between layers, while known rough soil emission models (i.e., Reference [36]) provide only the power reflectivity without phase information. Because the main purpose of this study is to evaluate the performance of snow emission models, it is important that the same soil emission model is used for the three snow emission models in order to avoid any bias in simulations that come from soil emission modelling. The same specular soil reflectivity model is thus used with each of the three snow emission models. The hypothesis of a specular soil is plausible in our case because the root-mean-square height (RMSH: 1.79 cm and 1.75 cm; see Section 2) of the soil measured with the LiDAR is much smaller than the L-band wavelength (RMSH < λ/12). Fresnel equations calculate the soil reflectivity from the permittivity of the frozen soil and the permittivity of the layer on top (air or snow). Snow permittivity is calculated from snow emission models, but frozen soil permittivity (ε g ) was not measured and remains an unknown. Hence, ε g was inferred from frozen ground snow-free radiometric measurements taken on 9 November (see Table 1). An iterative process with an increment of 0.1 was used to calculate the frozen soil permittivity that minimized the root-mean-square error (RMSE) between the measured T B (T B mes ), and simulated T B (T B sim ) at vertical (V-pol) and horizontal polarization (H-pol) at the four measured incidence angles such that: The optimization of ε g was done on Scene 1 and 3 separately. In this case, ε g should be considered as an effective parameter that allows representation of frozen soil emission for the three models, but can also partially compensate for the specular assumption.

Frozen Soil Permittivity
The optimized ε g was calculated from T B simulations for ε g ranging from 1 to 15 and using November snow-free frozen soil parameters. RMSE εg was computed from angular radiometer measurements performed on 1st November 2015 on Scene 1 and 3 at vertical and horizontal polarization ( Figure 2). The optimal value of ε g was 4.6 and 4.9 for Scene 1 and 3, respectively. RMSE εg of 10.7 K and 8.2 K, respectively, was observed for Scene 1 and 3, but an important component of this error was the poor simulation performance at 60 • V-pol (discussed in more detail later). RMSE εg computed without the T B at 60 • are 6.9 K and 5.7 K for an optimized ε g of 4.3 and 4.8 for Scene 1 and 3, respectively. The small differences between both sites are not significant, and could be related to differences in soil RMSH, which is not considered in the soil model. These optimized ε g were used in the following simulations. Fresnel equations calculate the soil reflectivity from the permittivity of the frozen soil and the permittivity of the layer on top (air or snow). Snow permittivity is calculated from snow emission models, but frozen soil permittivity (εg) was not measured and remains an unknown. Hence, εg was inferred from frozen ground snow-free radiometric measurements taken on November 9th (see Table  1). An iterative process with an increment of 0.1 was used to calculate the frozen soil permittivity that minimized the root-mean-square error (RMSE) between the measured TB (TB mes), and simulated TB (TB sim) at vertical (V-pol) and horizontal polarization (H-pol) at the four measured incidence angles such that: The optimization of εg was done on Scene 1 and 3 separately. In this case, εg should be considered as an effective parameter that allows representation of frozen soil emission for the three models, but can also partially compensate for the specular assumption.

Frozen Soil Permittivity
The optimized εg was calculated from TB simulations for εg ranging from 1 to 15 and using November snow-free frozen soil parameters. RMSEεg was computed from angular radiometer measurements performed on 1st November 2015 on Scene 1 and 3 at vertical and horizontal polarization ( Figure 2). The optimal value of εg was 4.6 and 4.9 for Scene 1 and 3, respectively. RMSEεg of 10.7 K and 8.2 K, respectively, was observed for Scene 1 and 3, but an important component of this error was the poor simulation performance at 60° V-pol (discussed in more detail later). RMSEεg computed without the TB at 60° are 6.9 K and 5.7 K for an optimized εg of 4.3 and 4.8 for Scene 1 and 3, respectively. The small differences between both sites are not significant, and could be related to differences in soil RMSH, which is not considered in the soil model. These optimized εg were used in the following simulations.

WALOMIS Gaussian Noise Parameterization
As is true of any model based on the wave approach, the result obtained for a specific snowpack configuration (i.e., a given set of inputs) may differ considerably from those obtained with a slightly different snowpack, which is not the case with the radiative transfer approach. This is due to the high sensitivity of interference phenomena to layer characteristics. Hence, on the basis of Reference [18], Remote Sens. 2018, 10, 1451 7 of 15 10,000 snow density profiles were generated by adding a Gaussian noise (σ d ) to the measured density profile and T B was obtained from the average of 10,000 WALOMIS simulations performed with these profiles. Figure 3 shows the effect of the Gaussian noise for an example of simulations for Scene 1 on 4 March 2015 with fixed layer thickness. Even with a very high Gaussian noise of σ d = 90 kg m −3 , the simulations show a wavy angular pattern owing to the high sensitivity of interference phenomena, very different to the measured angular spectrum. These results suggest that for shallow seasonal snow, adding Gaussian noise to density profiles does not reproduce the variability of snow cover within the radiometer field of view.

WALOMIS Gaussian Noise Parameterization
As is true of any model based on the wave approach, the result obtained for a specific snowpack configuration (i.e., a given set of inputs) may differ considerably from those obtained with a slightly different snowpack, which is not the case with the radiative transfer approach. This is due to the high sensitivity of interference phenomena to layer characteristics. Hence, on the basis of Reference [18], 10,000 snow density profiles were generated by adding a Gaussian noise (σd) to the measured density profile and TB was obtained from the average of 10,000 WALOMIS simulations performed with these profiles. Figure 3 shows the effect of the Gaussian noise for an example of simulations for Scene 1 on 4 March 2015 with fixed layer thickness. Even with a very high Gaussian noise of σd = 90 kg m −3 , the simulations show a wavy angular pattern owing to the high sensitivity of interference phenomena, very different to the measured angular spectrum. These results suggest that for shallow seasonal snow, adding Gaussian noise to density profiles does not reproduce the variability of snow cover within the radiometer field of view. In situ measurements performed during the campaign revealed strong variability in the thickness of internal layers within the snowpack (see Reference [17]), which could better represent interference phenomena in the model. WALOMIS simulations were also performed from 10,000 profiles of layers thickness by adding a Gaussian noise (σh) to the measured layer thickness, keeping the measured snow density. Figure 4 shows that with increasing σh, the angular pattern of the simulation gets closer to the TB measurements. Increasing the variability of layer thickness to 2-4 cm results in agreement with the TB observations (Figure 4). With a σh = 2 cm and 4 cm, the simulated Hpol also capture the TB decrease with incidence angle. However, contrary to observations which slowly decrease with incidence angle, simulations at V-pol tend to increase with incidence angle, before decreasing at 60°. As expected, because interference phenomena are highly sensitive to optical path-length across layers, a Gaussian noise of σh = 2 cm, was found to give the best agreement with measurements and was used for the subsequent simulations. In situ measurements performed during the campaign revealed strong variability in the thickness of internal layers within the snowpack (see Reference [17]), which could better represent interference phenomena in the model. WALOMIS simulations were also performed from 10,000 profiles of layers thickness by adding a Gaussian noise (σ h ) to the measured layer thickness, keeping the measured snow density. Figure 4 shows that with increasing σ h , the angular pattern of the simulation gets closer to the T B measurements. Increasing the variability of layer thickness to 2-4 cm results in agreement with the T B observations (Figure 4). With a σ h = 2 cm and 4 cm, the simulated H-pol also capture the T B decrease with incidence angle. However, contrary to observations which slowly decrease with incidence angle, simulations at V-pol tend to increase with incidence angle, before decreasing at 60 • . As expected, because interference phenomena are highly sensitive to optical path-length across layers, a Gaussian noise of σ h = 2 cm, was found to give the best agreement with measurements and was used for the subsequent simulations.

Footprint Integration
Because of the large antenna full beamwidth (30°) and the footprint geometry of the surfacebased radiometer, the simulated TB of a single directional incidence angle (θ) might not be representative of the measured TB over a large footprint, especially at higher incidence angles [37]. Hence, in this study, a weighting function was computed to estimate the integrated TB within the footprint for a range of incidence angles. For this estimation, the area included in θ ± 15° was considered. A Gaussian weighting was applied to θ with a standard deviation of 7.5° in order to represent the antenna directional power sensitivity pattern. Then, a factor 1/r 2 , with r the distance from the radiometer, was applied to the obtained coefficients to attenuate the contribution as a function of the location within the footprint. The coefficients used for the weighting are illustrated in Figure 5, normalized by the maximum. When the weighting function was applied to the simulations of the three snow emission models at Scene 3 on 4 March 2015, there is a decrease of TB at both H-pol and V-pol mostly at incidence angles higher than 55° ( Figure 6). This decrease in TB slightly improved the results for all three models at high incidence angles, thus it was used in the following simulations.

Footprint Integration
Because of the large antenna full beamwidth (30 • ) and the footprint geometry of the surface-based radiometer, the simulated T B of a single directional incidence angle (θ) might not be representative of the measured T B over a large footprint, especially at higher incidence angles [37]. Hence, in this study, a weighting function was computed to estimate the integrated T B within the footprint for a range of incidence angles. For this estimation, the area included in θ ± 15 • was considered. A Gaussian weighting was applied to θ with a standard deviation of 7.5 • in order to represent the antenna directional power sensitivity pattern. Then, a factor 1/r 2 , with r the distance from the radiometer, was applied to the obtained coefficients to attenuate the contribution as a function of the location within the footprint. The coefficients used for the weighting are illustrated in Figure 5, normalized by the maximum.

Footprint Integration
Because of the large antenna full beamwidth (30°) and the footprint geometry of the surfacebased radiometer, the simulated TB of a single directional incidence angle (θ) might not be representative of the measured TB over a large footprint, especially at higher incidence angles [37]. Hence, in this study, a weighting function was computed to estimate the integrated TB within the footprint for a range of incidence angles. For this estimation, the area included in θ ± 15° was considered. A Gaussian weighting was applied to θ with a standard deviation of 7.5° in order to represent the antenna directional power sensitivity pattern. Then, a factor 1/r 2 , with r the distance from the radiometer, was applied to the obtained coefficients to attenuate the contribution as a function of the location within the footprint. The coefficients used for the weighting are illustrated in Figure 5, normalized by the maximum. When the weighting function was applied to the simulations of the three snow emission models at Scene 3 on 4 March 2015, there is a decrease of TB at both H-pol and V-pol mostly at incidence angles higher than 55° ( Figure 6). This decrease in TB slightly improved the results for all three models at high incidence angles, thus it was used in the following simulations. When the weighting function was applied to the simulations of the three snow emission models at Scene 3 on 4 March 2015, there is a decrease of T B at both H-pol and V-pol mostly at incidence angles higher than 55 • (Figure 6). This decrease in T B slightly improved the results for all three models at high incidence angles, thus it was used in the following simulations.

Snow Emission Model Intercomparison
The effect of dry snow at L-band is mostly related to refraction and impedance matching [16]. Impedance matching by dry snow reduces dielectric gradients and consequently increases thermal emission of the scene, while refraction caused by the snow layer in contact with the ground surface leads to a steeper incidence angle at the ground surface in comparison with the observation angle [16]. The impact is in general higher at H-pol, since at V-pol these effects are partly compensatory, and even fully compensatory around an incidence angle of 51°. At H-pol, the effect of snow typically increases with incidence angle. At the plot scale, the presence of snow can change TB at H-pol from 5 K at 30° up to 20 K at 60° [17]. The specular soil reflectivity model with optimized permittivity and the optimized weighting function were used to simulate the TB at the Scene 1 and 3 and for the six sampling periods with the three snow emission models (Figure 7). The three models show very similar overall RMSE at V-pol and H-pol with values ranging between 7.2 K and 10.5 K. LS-MEMLS-1L gives lower RMSE at H-pol while DMRT-ML gives slightly better results at V-pol. Note that for the three models and both polarizations, the RMSE increases with increasing incidence angle ( Table  2). The worst results are obtained at 60°, which has a strong impact on the overall RMSE (noted "All" in Table 2). It is thus not possible to state that a specific snow emission model gives better results overall. Nevertheless, at 50° and 60°, LS-MEMLS-1L clearly underestimates the variability of TB Hpol, with a standard deviation of the simulations much lower than the standard deviation of observations. Simulations at 60° ranged between 203.9 K and 213.9 K, while the measured TB H-pol ranged between 192.8 K and 225.8 K. On the other hand, the standard deviations of WALOMIS and DMRT-ML simulations at TB H-pol are much closer to the measurements (Table 3). These results suggest that snow layers significantly impact the TB H-pol. The multi-layer model configurations are able to better capture this effect, but the complex interaction of the radiation within the layers and the difficulty to precisely measure the snow layer characteristics in the field at the meter scale (footprint) leads to an overall RMSE comparable to LS-MEMLS-1L applied in a 1-layer configuration. We thus face the problem where more complex and more sensitive radiometric models require precise in situ information for comprehensive evaluation, a condition that will provide limitations for more general use, especially at the satellite scale.

Snow Emission Model Intercomparison
The effect of dry snow at L-band is mostly related to refraction and impedance matching [16]. Impedance matching by dry snow reduces dielectric gradients and consequently increases thermal emission of the scene, while refraction caused by the snow layer in contact with the ground surface leads to a steeper incidence angle at the ground surface in comparison with the observation angle [16]. The impact is in general higher at H-pol, since at V-pol these effects are partly compensatory, and even fully compensatory around an incidence angle of 51 • . At H-pol, the effect of snow typically increases with incidence angle. At the plot scale, the presence of snow can change T B at H-pol from 5 K at 30 • up to 20 K at 60 • [17]. The specular soil reflectivity model with optimized permittivity and the optimized weighting function were used to simulate the T B at the Scene 1 and 3 and for the six sampling periods with the three snow emission models (Figure 7). The three models show very similar overall RMSE at V-pol and H-pol with values ranging between 7.2 K and 10.5 K. LS-MEMLS-1L gives lower RMSE at H-pol while DMRT-ML gives slightly better results at V-pol. Note that for the three models and both polarizations, the RMSE increases with increasing incidence angle ( Table 2). The worst results are obtained at 60 • , which has a strong impact on the overall RMSE (noted "All" in Table 2). It is thus not possible to state that a specific snow emission model gives better results overall. Nevertheless, at 50 • and 60 • , LS-MEMLS-1L clearly underestimates the variability of T B H-pol, with a standard deviation of the simulations much lower than the standard deviation of observations. Simulations at 60 • ranged between 203.9 K and 213.9 K, while the measured T B H-pol ranged between 192.8 K and 225.8 K. On the other hand, the standard deviations of WALOMIS and DMRT-ML simulations at T B H-pol are much closer to the measurements (Table 3). These results suggest that snow layers significantly impact the T B H-pol. The multi-layer model configurations are able to better capture this effect, but the complex interaction of the radiation within the layers and the difficulty to precisely measure the snow layer characteristics in the field at the meter scale (footprint) leads to an overall RMSE comparable to LS-MEMLS-1L applied in a 1-layer configuration. We thus face the problem where more complex and more sensitive radiometric models require precise in situ information for comprehensive evaluation, a condition that will provide limitations for more general use, especially at the satellite scale.  Table 2. Root-mean-square height (RMSE) of all the dates and Scene 1 and 3 together. "All" is the RMSE calculated with incidence angles of 30°, 40°, 50° and 60°.

Soil Permittivity Parameterization
The optimized values of frozen soil permittivity are close to other studies [16,17,35,38,39]. Rautiainen et al. [35] obtained L-band real part soil permittivity values between 3.3 to 3.8 in boreal forest frozen soils; Schwank et al. [16] showed that for L-band, the real part of the frozen soil permittivity varies from 3.5 to 4.5, while Hallikainen et al. [38] showed that it varies from 5 to 8 in the 10 to 18 GHz frequency range. At the same frequency range, Mironov et al. [39] developed a temperature dependent permittivity model and showed that the permittivity could vary from 3 to 4.5 for a frozen soil at −25 °C. The optimized values are similar to the ones obtain at the same site ( [17]; εg = 5.1), but using a two parameter (εg and snow density) retrieval method [8]. The plausible optimized εg values thus show that the assumption of a flat soil at L-band was reasonable for our site and the impact of soil roughness is reasonably included in the εg. However, in this study, soil emission was empirically parameterized to minimize its impact on the simulation of snow covered ground, but it remains that the snow-free frozen soil RMSE are similar to the snow simulations. In particular, Figure 6 shows that at V-pol, there is a clear underestimation of the TB at an incidence angle of 30°, while there is an overestimation of TB at higher incidence angle (50° and 60°). There is thus an opposite trend in the angular spectrum between the measurements and the observations, an

Soil Permittivity Parameterization
The optimized values of frozen soil permittivity are close to other studies [16,17,35,38,39]. Rautiainen et al. [35] obtained L-band real part soil permittivity values between 3.3 to 3.8 in boreal forest frozen soils; Schwank et al. [16] showed that for L-band, the real part of the frozen soil permittivity varies from 3.5 to 4.5, while Hallikainen et al. [38] showed that it varies from 5 to 8 in the 10 to 18 GHz frequency range. At the same frequency range, Mironov et al. [39] developed a temperature dependent permittivity model and showed that the permittivity could vary from 3 to 4.5 for a frozen soil at −25 • C. The optimized values are similar to the ones obtain at the same site ( [17]; ε g = 5.1), but using a two parameter (ε g and snow density) retrieval method [8]. The plausible optimized ε g values thus show that the assumption of a flat soil at L-band was reasonable for our site and the impact of soil roughness is reasonably included in the ε g . However, in this study, soil emission was empirically parameterized to minimize its impact on the simulation of snow covered ground, but it remains that the snow-free frozen soil RMSE are similar to the snow simulations. In particular, Figure 6 shows that at V-pol, there is a clear underestimation of the T B at an incidence angle of 30 • , while there is an overestimation of T B at higher incidence angle (50 • and 60 • ). There is thus an opposite trend in the angular spectrum between the measurements and the observations, an effect most likely related to the soil model. It should be note that those error impact substantially the results with snow-covered surface. Simulations and observations at V-pol will have to be investigated further.

WALOMIS Gaussian Noise Parameterization
This study shows that an approach developed by Reference [18], which applies Gaussian noise to in situ snow density measurements to estimate T B variability is insufficient to smooth the snow layer interference phenomena in WALOMIS. However, Gaussian noise applied to snow depth layer thickness of σ h = 2 cm leads to an angular spectrum comparable to the observations. A value of σ h = 2 cm is plausible considering the high standard deviation of snow depth in the area and the variable internal layering of the snowpack. Snow depth surveys around the area suggest that despite the surface appearing relatively homogeneous, the standard deviation of the snow depth is in the range of 28 to 56% ( Table 1). The high variability in snow depth is related to wind redistribution, a key process in this prairie environment [40]. Strong winds also lead to a high snow density layer observed at top of the snowpack due to compaction processes [41], thus resulting in the optimized snow depth layer thickness variability σ h of 2 cm at the plot scale.

Footprint Integration
Roy et al. [37] showed that spatially distributed surface-based radiometer observations within a SMOS pixel gives constantly lower T B (between 8.5 and 22.9 K) than SMOS observations at 60 • V-pol (not apparent at lower incidence angles). It was hypothesized that these discrepancies were due to the large beamwidth (30 • ) of the surface-based radiometer, which produced an incidence angle in the far range of up to 75 • . This may exaggerate the influence of high incidence angles on reducing the magnitude of V-pol T B . The calculation of the weighted footprint integration in this paper shows that the effective incidence angle can explain a small part of the ground-based radiometer bias to lower T B at 60 • V-pol ( Figure 6). For the three models, RMSE was reduced 2.4 K at 60 • V-pol because of the footprint integration, but no clear improvement was observed at H-pol. The footprint integration cannot fully account for the entire difference between simulations and measurements at higher incidence angles. Thus, part of the difference may be the result of the sky contribution (low T B ) captured by antenna side lobes.

Snow Emission Model Intercomparison
Results show that all models give similar overall RMSE values between 7.2 K and 10.5 K. It is thus difficult to identify a preferred model for this environment. Those snow-covered RMSE are very similar to snow-free RMSE, which suggest that a large part of the errors come from the soil parameterization (see Section 5.1). However, our results suggest that the multi-layer models, DMRT-ML and WALOMIS, are able to better simulate the range of observed T B compared to the one-layer model. DMRT-ML and WALOMIS simulations are similar, and suggest that the interference phenomena, simulated by WALOMIS, have a marginal impact on the simulations. Leduc-Leballeur et al. [14] showed that the coherence effect considered in WALOMIS significantly improve the simulation for the semi-infinite snow layered medium in Antarctica, which is not the case for shallow seasonal snowpack. It seems that DMRT-ML, by considering the refraction between layers incorporate the main features needed to reproduce variations in the T B similar to the observations. However, the snowpack during the winter was shallow with very strong meter scale spatial variability and a stratigraphy that made the snow density measurements challenging to document within the field of view of the radiometer (Table 1). Furthermore, because of the several melt/refreeze events that occurred during the winter several ice layers were present in the snowpack resulting in substantial spatial variability. Because of this high spatial variability in the presence of ice layers and snow density within the snowpack, it is unlikely that individual snow pit measurements are representative of the radiometer footprint, and thus a further source of model uncertainty. Therefore, it is difficult to precisely quantify the effect of coherence in an ice layer, even in a relatively controlled experiment such as that described. A similar experiment but with deeper and more homogeneous snow could help to improve our understanding of snow layer effects on L-band emission. The higher complexity of simulations with WALOMIS may not be necessary, unless coherence and layering effects are very strong.
Although the one-layer model does not capture stratigraphy and coherence effects captured in the more detailed multi-layered models, the overall performance was similar. Hence, in a snow density inversion scheme such as that proposed by Reference [8], the one-layer model is a practical solution because it includes only a small number of free parameters, which is almost a prerequisite to achieve unique values of respective retrievals.

Summary and Conclusions
Recent studies show the non-negligible impact of snow on L-band emission [8,17,19] even when the snowpack is shallow. In this study, the Wave Approach for LOw-frequency MIcrowave emission in Snow was adapted and parameterized for seasonal snow using ground-based L-band radiometer observations in a prairie environment. A specular soil reflectivity model was added to WALOMIS, and frozen soil permittivity (ε g ) of 4.3 and 4.8 used for Scene 1 and 3, respectively, was estimated using an optimisation scheme that compared, observed and simulated T B (DMRT-ML) from snow free frozen ground measurements.
Gaussian noise of snow depth layers of σ h = 2 cm obtained a comparable multi-angular T B response but still underestimated the observations at H-pol. The same simulations also underestimated at lower incidence angles and overestimated at higher incidence angles compared to V-pol observations. The calculation of a weighted footprint integration shows that the effective incidence angle can explain part, but not all of the ground-based radiometer bias to lower T B at 60 • V-pol, but no clear improvement was observed at H-pol. The WALOMIS simulations were then compared to two other radiative transfer models, the DMRT-ML and a simplified adaptation of MEMLS for L-band (LS-MEMLS-1L), used in a single-layer configuration. RMSE between simulated and measured T B were similar for the three models with overall RMSE between 7.2 and 10.5 K. However, WALOMIS and DMRT-ML were able to better reproduce the range of T B of the observations at higher V-pol incidence angles (50 • and 60 • ) and across all incidence angles at H-pol.
The impact of snow on surface emission at L-band is relatively small compared to higher frequencies, where the wavelength is similar to snow grain size, inducing scattering. However, there is still valuable and complementary snow information in L-band observations, such as sensitivity to snow density [9] and ice layers [6]. Hence, it is important to develop and assess existing emission models that will help to better quantify these different effects for applications such as the observation of the soil freeze/thaw status. In this study, an effort was made to develop and parameterize WALOMIS for seasonal snow. While the results of all three models are similar, the wave approach of WALOMIS suggests that it is an appropriate new tool to better understand and model L-band emission when interference and snow layering is present.