Assimilation of synthetic observations of radar backscatters at Ku-band improves SWE estimates

In cold regions, snow serves as the primary water source for downstream rivers and lakes. Accurate gridded snow water equivalent (SWE) estimation is hindered by the sparse ground observation network and the low resolution of satellite passive microwave products. To address this, Environment and Climate Change Canada (ECCC), the Canadian Space Agency (CSA), and Natural Resources Canada (NRCan) are developing the Terrestrial Snow Mass Mission (TSMM), a dual Ku-band satellite mission designed to measure backscatter at 13.5 and 17.25 GHz across the Northern Hemisphere at a 500 m spatial resolution with a weekly temporal resolution. This study assesses the feasibility of assimilating Ku-band backscatter to enhance SWE estimates in a synthetic experiment. We used the Soil-Vegetation-Snow version 2 (SVS2) land surface model, which incorporates the snowpack model Crocus, coupled with the Snow Microwave Radiative Transfer model (SMRT). Synthetic observations of SWE and backscatter extracted at weekly intervals from synthetic truths (model simulations) were assimilated with a particle filter at point-scale. This was done at three sites representing three different Canadian climates (Arctic, humid continental, Alpine) over three winter seasons. Meteorological forcing derived from the high-resolution Canadian meteorological model was perturbed to generate ensembles of snow simulations for assimilation. Results indicate that assimilating synthetic observations of backscatter improved SWE estimates at the Arctic and humid continental sites, reducing the mean continuous ranked probability score (CRPS) by up to 32 % compared to the open-loop ensemble. This performance was comparable to assimilating the SWE synthetic observations with observation errors larger than 20 %. Assimilating synthetic observations of backscatter at the Alpine site only improved the SWE estimates by 5 % as backscatter signals seemed to lose sensitivity to SWE values greater than ∼300 kg m⁻² in our experimental setup. Assimilating backscatter and SWE synthetic observations also improved the estimations of vertical profiles of snow density and specific surface area. These findings demonstrate the potential of direct assimilation of Ku-band backscatter to enhance both estimates of SWE and snowpack properties.

1 Introduction

Gridded snow water equivalent (SWE) estimates are an essential component for water and food security, ecosystem sustainability, and predicting flood and drought risks in cold regions (Pomeroy et al., 2016; Wagner et al., 2017; Sturm et al., 2017). As the climate is changing, snow is being altered both spatially and temporally, which is likely to increase avalanche risk and flooding in both the short and long term (Haeberli and Whiteman, 2015; Hamlet and Lettenmaier, 2007; Beniston, 2012; Wagner et al., 2017). Distributed SWE estimates from snowpack models help understand the spatial distribution of the snow properties over a domain compared to sparse observations, and can be used to initialize land surface and hydrological forecasts (Orio-Alonso et al., 2023; Garnaud et al., 2021). However, SWE estimates from snow models are subject to uncertainties linked to different parameterization of snow processes (Lafaysse et al., 2017) and uncertainties in the meteorological forcings, in particular in the precipitation amount and phase (Günther et al., 2019; Raleigh et al., 2015; Vionnet et al., 2022; Leroux et al., 2023). There is a critical need for improved estimates of SWE at high temporal and spatial resolutions to advance reanalysis products, operational models, and climate projections (Girotto et al., 2020).

Data assimilation is a commonly used method to decrease uncertainty in snow predictions by combining observations and model estimates. Observations used for snow data assimilation can originate from in-situ manual snow observations (Orio-Alonso et al., 2023; Magnusson et al., 2014; Oberrauch et al., 2024; Piazzi et al., 2018; Slater and Clark, 2006; Brasnett, 1999), remote sensing (Cluzet et al., 2020; De Lannoy et al., 2012; Durand et al., 2009; Larue et al., 2018a; Aalstad et al., 2018; Li et al., 2017; Su et al., 2008, 2010; Lahmers et al., 2022), or retrieved SWE products (Shrestha and Barros, 2025b; De Lannoy et al., 2012). Remote sensing observations for snow data assimilation predominantly utilize passive microwave sensors (Durand et al., 2009; Larue et al., 2018a; Lemmetyinen et al., 2018; Li et al., 2017; Dumont et al., 2012) or snow cover fraction and albedo products derived from optical sensors (De Lannoy et al., 2012; Garnaud et al., 2021; Aalstad et al., 2018; Su et al., 2008, 2010). Additional satellite-derived and airborn remote sensing datasets have been explored for snow data assimilation, including terrestrial water storage (TWS) information from the Gravity Recovery and Climate Experiment (GRACE) satellites (Su et al., 2010), snow depth from Pleiades imagery (Marti et al., 2016; Shaw et al., 2020), snow depth retrieved from ICESat-2 (Mazzolini et al., 2025), and SWE estimates from the NASA Airborne Snow Observatory (Lahmers et al., 2022).

Data assimilation of both in-situ and satellite observations can improve SWE estimates, though each observation type has different advantages and limitations. In-situ snow depth measurements provide better SWE estimates compared to satellite optical data when assimilated (e.g. Charrois et al., 2016), but the sparse distribution of in-situ measurements limits spatial representativeness (e.g. Alonso-González et al., 2023) and poses challenges when propagating information to unobserved regions (Pflug et al., 2024; Cluzet et al., 2022). Satellite observations offer broader spatial coverage but present other limitations. Passive microwave data provide global coverage with 40 year time series (e.g. Parkinson, 2022), but their coarse spatial resolution (e.g Larue et al., 2018b) makes derived snow products unreliable in mountainous terrain (Luojus et al., 2021). Active microwave sensors can retrieve SWE at higher spatial resolution (<500 m) and for deeper snowpacks (Tsang et al., 2022), providing an alternative to passive systems. Optical sensors face limitations from cloud cover, which affects data availability in the visible and near-infrared ranges (De Lannoy et al., 2012). Both optical and microwave remote sensing data have higher uncertainties for wet snow conditions and forested landscapes, which can limit accurate snow-covered area classification (Gascoin et al., 2024; Muhuri et al., 2021).

To address the need for more accurate, high spatial resolution SWE estimates from remote sensing data, Environment and Climate change Canada (ECCC), in partnership with the Canadian Space Agency (CSA) and Natural Resources Canada (NRCan), is developing the Terrestrial Snow Mass Mission (TSMM) (Derksen et al., 2021). This mission will be the only active synthetic aperture radar (SAR) mission dedicated to snow, aiming to launch a dual-frequency Ku-band satellite (13.5 and 17.25 GHz) with a weekly revisit frequency and a resolution of 500 m. TSMM is designed to complement existing passive microwave observations, thereby enhancing SWE retrieval capabilities (Lemmetyinen et al., 2018). Recent studies highlight the significance of Ku-band response to variations in SWE and snow microstructure, which are critical factors for SWE retrieval approaches utilizing Ku-band backscatter measurements (Lemmetyinen et al., 2018; Montpetit et al., 2025; Tsang et al., 2022). Previous research exploring the potential of assimilating retrieved SWE from Ku-band and/or X-band volume-scattering data has shown promise in improving snow profile properties (Phan et al., 2014) and SWE estimates (Cho et al., 2023; Pflug et al., 2024; Shrestha and Barros, 2025b; Garnaud et al., 2019). Specifically, Cho et al. (2023) demonstrated that improvements in SWE prediction root-mean-square error (RMSE) were more pronounced during the melting season than during the accumulation period in the mountainous regions of Colorado. Their work also indicated that lower assimilation RMSE was expected if retrieval algorithms are able to retrieve a wider range of SWE values, while higher tree cover fractions negatively impacted assimilation performance. To mitigate this, Pflug et al. (2024) proposed a method to estimate forest SWE from adjacent forest-free areas to improve SWE prediction in forested environments via data assimilation. More recently, Shrestha and Barros (2025b) showed that assimilating retrieved SWE from X- and Ku-band backscatter not only reduced snow depth and SWE model biases and errors in vertical snow density profiles, but also improved the forward simulation of volume backscatter. They also noted that larger observation errors in SWE retrievals led to increased RMSE in SWE prediction.

Durand et al. (2009) demonstrated that directly assimilating passive microwave radiance observations provides more accurate snow depth predictions compared to assimilating retrieved SWE. Similarly, in soil moisture research, assimilating brightness temperature measurements from passive microwave satellite observations is often preferred over retrieved soil moisture data to enhance estimation accuracy (e.g. De Lannoy and Reichle, 2016; Carrera et al., 2019). A significant limitation of retrieved products (Level 2 data) is their delayed availability, often hours after the original observations (Level 1 data), which poses challenges for near-real-time operational assimilation systems. To date, snow data assimilation studies have not yet explored the potential of assimilating Ku-band backscatter observations to improve SWE estimates. This study aims to address this gap in the literature.

Point-scale synthetic experiments are a common approach to evaluate the feasibility of assimilation schemes, utilizing synthetic observations derived from model runs (Charrois et al., 2016; Durand and Margulis, 2006; Larue et al., 2018a; Revuelto et al., 2021). The objective of this work is to demonstrate the potential of direct Ku-band backscatter assimilation to enhance SWE estimation through a synthetic point-scale experiment. This is made possible by the recent development of a new land surface scheme (Soil-Vegetation-Snow version 2, SVS2) at ECCC that includes the multi-layered snow model Crocus (Vionnet et al., 2025a). SVS2 can provide the necessary snow inputs to a forward radiative transfer model, such as the Snow Microwave Radiative Transfer model (SMRT) (Picard et al., 2018) to estimate backscatter at Ku-band. SVS2/Crocus coupled with SMRT were included in the Multiple Snow Data Assimilation System (MuSA) platform (Alonso-González et al., 2022) to develop the synthetic data assimilation experiment detailed in this paper. Due to the strong non-linearity of multi-layered snow models, the particle filter is used for the data assimilation (Cluzet et al., 2021; Charrois et al., 2016; Larue et al., 2018a; Revuelto et al., 2021). In this study, SWE estimates derived from assimilating synthetic observations of backscatter are compared against SWE estimates obtained by assimilating synthetic SWE observations with varying levels of uncertainty. Section 2 presents the study sites and data used, the models applied, and the synthetic experiment design. The synthetic experiment results on bulk snow properties prediction (SWE and snow depth) and vertical snow properties (density and snow specific area, SSA) are presented in Sect. 3 and discussed in Sect. 4.

2 Design of the synthetic experiments

2.1 Study Sites

The synthetic experiments are run at three study sites that span different Canadian climates: Trail Valley Creek (TVC) is an Arctic site in the Northwest Territories (68.74, −133.5°, elevation of 91 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$, above sea level), Rogers Pass is an Alpine site in British Columbia (51.23, −117.71°, elevation of 1905 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$), and Powassan is an agricultural site, in a humid continental climate in Ontario (46.08, −79.36°, elevation of 256 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$) (Fig. 1). These sites were chosen because field experiments for science readiness activities of TSMM were conducted there, gathering in-situ snow observations (snow pits and snow courses), airborne and ground-based Ku-band backscatter, and meteorological observations between 2018–2024 (Montpetit et al., 2024; Kelly et al., 2024; Madore et al., 2023). The study period includes three consecutive winters, from September 2020 to August 2023. Table 1 summarizes the total snowfall amounts and mean air temperature between 1 September–30 June at each site for each winter season extracted from the High Resolution Deterministic Prediction System (HRDPS) (Milbrandt et al., 2016) and the Canadian Precipitation Analysis System (CaPA) (Khedhaouiria et al., 2020; Lespinas et al., 2015) using Wang et al. (2019a) for the precipitation phase partitioning (Sect. 2.3.1). TVC receives less snow than the other sites, while Rogers Pass received significantly more snow. Powassan is the warmest site with mean air temperatures above freezing and TVC the coldest site. Meteorological data, including air temperature, relative humidity, and wind speed, were measured at Powassan during the 2022–2023 winter season (Kelly et al., 2024), and were used to inform the forcing perturbations needed for the ensemble generation (Sect. 2.3.1).

Figure 1 Map showing the locations of the three sites: Trail Valley Creek (TVC), Rogers Pass, and Powassan. Basemap: ESRI world imagery | Powered by Esri.

Table 1 Total snowfall (water equivalent, w.e.) and mean air temperature (1 September–30 June) at each study site.

2.2 SVS2/Crocus and SMRT

The land surface scheme SVS2 developed at ECCC (Vionnet et al., 2025a) contains the snowpack model Crocus (Vionnet et al., 2012; Lafaysse et al., 2017, 2025). Crocus is a one-dimensional snowpack model that simulates the seasonal evolution of the physical properties of the snowpack and its vertical layering. For each snow layer, Crocus simulates the evolution of the thickness, density, liquid water content, temperature, age, and snow microstructure represented by the snow specific surface area (SSA) and the snow grain sphericity (Brun et al., 1992; Vionnet et al., 2012; Carmagnola et al., 2014). Crocus was initially developed to simulate the properties of alpine snow in the context of avalanche hazard forecasting (e.g. Durand et al., 1993). Recently, Woolley et al. (2024) proposed an Arctic configuration of SVS2/Crocus that improves the simulations of Arctic snowpack properties through a better representation of wind-packing and inclusion of the effect of basal vegetation on snow compaction. This Arctic version of Crocus was used at TVC and the default version of Crocus was used at Powassan and Rogers Pass (see Table A1 for the Crocus parameterizations). Montpetit et al. (2025) used the Arctic version of Crocus to successfully retrieve SWE at TVC from Ku-band Synthetic Aperture Radar (SAR) measurement. In our study, a maximum of 20 snow layers was specified to simulate the evolution of the snowpack properties with SVS2/Crocus.

SMRT is a snow radiative transfer model (Picard et al., 2018) that was coupled with SVS2 (Meloche et al., 2025). SMRT was used to compute backscatter signal (σ) at the two TSMM frequencies, 13.5 and 17.25 GHz in VV polarization, from simulated snow layered properties (density, thickness, SSA, and temperature) from SVS2/Crocus. The DORT solver was used with the Improved Born Approximation and an exponential microstructure. An incidence angle of 35° was assumed. The simulated soil temperature and water content from SVS2 in the upper 5 cm of the soil were used to calculate the soil permittivity using the Mironov model (Mironov et al., 2019).

2.3 Assimilation design

SVS2, coupled with SMRT, was included in MuSA version 2.0 (Alonso-González et al., 2022). MuSA is a Python toolbox that can generate ensembles of snow simulations and offers a wide selection of assimilation methods. Given its compatibility with non-linear models, such as multilayered snow models, the particle filter method was adopted in this study. This sequential data assimilation approach updates the snow state estimates at each observation time by weighing and resampling particles according to their likelihood given the observations. The particle filter has been used extensively with complex snow models (Cluzet et al., 2021; Magnusson et al., 2017; Revuelto et al., 2021; Larue et al., 2018a; Charrois et al., 2016).

2.3.1 Generation of the Ensemble

Because uncertainties in the snow simulations were assumed to originate from uncertainties in the meteorological forcing (Raleigh et al., 2015; Günther et al., 2019), the snow ensemble were generated by only perturbing the meteorological forcing. Meteorological forcing for SVS2/Crocus was obtained at each site from the HRDPS at 2.5 km grid spacing. These forcing included air temperature, specific humidity, wind speed, surface pressure, and incoming longwave and shortwave radiation. Successive short-term HRDPS forecasts were combined to generate continuous hourly meteorological forcing. HRDPS forecasts initialized every 6 h (00:00, 06:00, 12:00, 18:00 Z) were used, with forecast lead times of 7–12 h extracted from each run and concatenated to create a continuous hourly time series. Air temperature and humidity were at 2 $\mathrm{m}\mathrm{a}.\mathrm{g}.\mathrm{l}.$ (above ground level) and wind speed was at 10 $\mathrm{m}\mathrm{a}.\mathrm{g}.\mathrm{l}.$ The precipitation amount was taken from the CaPA 2.5 km (Khedhaouiria et al., 2020; Lespinas et al., 2015). Once the meteorological forcings were perturbed, the phase of the precipitation was determined using the approach of Wang et al. (2019b) that relies on the near-surface wet-bulb temperature.

Separate perturbations are applied to different meteorological forcing in order to represent errors in the HRDPS forcing data compared to local observations. The time evolution of each perturbation follows a first-order auto-regressive model describing the time evolution of an error as in Magnusson et al. (2017):

$$\begin{array}{}\text{(1)}& {\displaystyle }{\displaystyle }{q}_k=\mathit{\alpha }{q}_{k-\mathrm{1}}+\sqrt{\mathrm{1}-{\mathit{\alpha }}^{\mathrm{2}}}{\mathit{\omega }}_k\end{array}$$

where q_k is the error at time k, ω_k is a white noise with mean of 0 and standard deviation of 1 that changes with k, and α is a function of the decorrelation time length τ:

$$\begin{array}{}\text{(2)}& {\displaystyle }{\displaystyle }\mathit{\alpha }=\mathrm{1}-{\displaystyle \frac{\mathrm{\Delta }t}{\mathit{\tau }}}\end{array}$$

with Δt being the model time step. α was determined by calculating the autocorrelation of the error between the HRDPS and the observations with a lag of 1 and τ was determined from the calculated α values.

Random perturbations were applied to the meteorological inputs either as additive or multiplicative perturbations. Following Charrois et al. (2016), Larue et al. (2018b), the additive perturbations were drawn from a normal distribution and were applied to the air temperature and incoming longwave radiation forcing while the multiplicative perturbations were drawn from a log-normal distribution and were used for the precipitation, wind speed, and shortwave radiation forcing. The perturbations applied to the meteorological forcing (n_k) at time k were as follow:

$$\begin{array}{}\text{(3)}& {\displaystyle }& {\displaystyle }{n}_k={\mathit{\mu }}_f+q_k{\text{SD}}_f& {\displaystyle }& {\displaystyle }\text{if additive perturbation}\text{(4)}& {\displaystyle }& {\displaystyle }{n}_k=\text{exp}({\mathit{\mu }}_f+q_k{\text{SD}}_f)& {\displaystyle }& {\displaystyle }\text{if multiplicative perturbation}\end{array}$$

For additive perturbations, μ_f and SD_f are the mean and standard deviation of the normal distribution, respectively. For multiplicative perturbations, they are the parameters of the underlying normal distribution in the log-normal formulation. For the normal distributions, we assumed μ_f equal to 0 and SD_f is taken as the standard deviation of the differences between HRDPS model predictions and actual meteorological observations at Powassan during 2022–2023 (Sect. 2.1). For the log-normal distributions (Eq. 4), we calculated SD_f as the standard deviation of the logarithm of the ratio between HRDPS predictions and observations, and we adjusted μ_f so that the mean of n_k equaled 1. This approach made it possible to base the uncertainty of our perturbations on discrepancies in the HRDPS forcing data. For simplicity, these discrepancies between observations and HRDPS were assumed to hold true at the other two sites. Table 2 summarizes the parameters used to generate the perturbations. As direct longwave radiation measurements were not available to determine its perturbation, a linear regression was determined using the HRDPS forcing data, modelling changes in longwave radiation as a function of changes in air temperature (Table 2). This method propagates the air temperature uncertainty to longwave radiation based on their HRDPS-derived correlation. Based on this perturbation strategy, a total of 100 members were generated for each assimilation experiment as it was found suitable for snow assimilation with the particle filter (Piazzi et al., 2018).

Table 2 Parameters of the perturbation applied to the meteorological forcing to generate the ensemble.

^* The linear regression between the longwave radiation perturbation (ϵ_k,LW) and air temperature perturbation (${\mathit{ϵ}}_{k,T_{\mathrm{a}}}$) is ${\mathit{ϵ}}_{k,\text{LW}}=\mathrm{3.7}{\mathit{ϵ}}_{k,T_{\mathrm{a}}}$ with ${\mathit{ϵ}}_{k,T_{\mathrm{a}}}$ in K.

2.3.2 Data Assimilation Experiments

This study focuses on idealized experiments, in which synthetic observations were generated from reference model runs and subsequently assimilated (e.g. Revuelto et al., 2021; Durand and Margulis, 2006; Larue et al., 2018b). For each winter season and site, we generated 10 reference runs of SVS2/SMRT using perturbed meteorological inputs (Sect. 2.3.1). These runs were used as synthetic true states providing reference snow states and backscatter synthetic values. We used multiple reference runs to evaluate data assimilation performance across different snowpack conditions within each winter season, following the approach of Revuelto et al. (2021). Because the main source of uncertainties is assumed to originate from the meteorological forcing (Sect. 2.3.1), the model parameterization between these reference runs and the ensemble members was identical. SWE and backscatter values were extracted from the reference runs at weekly intervals on Mondays at 12:00 UTC, which is the expected revisit period of TSMM (but not particularly the expected date of observation at the sites and the time was chosen to avoid wet snowpack observations potentially occurring later in the day) so long as the following criteria are met: (1) bulk liquid water content within the snowpack below 1 % of mass as TSMM measurements for wet snowpacks would be discarded, and (2) backscatter values at 13.5 and 17.25 GHz above −24 and −26 dB, respectively, which correspond to the lowest snow backscatter observations at these frequencies from a tower-based radar system deployed in Sodankylä, Finland used in Lemmetyinen et al. (2018), Pan et al. (2024). Table 3 presents the number of synthetic observations for each winter season across the 10 reference runs. TVC and Rogers Pass have the highest number of synthetic observations. Rogers Pass has the longest snow seasons but has more occurrences of wet snow than at TVC.

Table 3 Sample sizes by site and winter season across the 10 reference runs.

A series of assimilation experiments were conducted to evaluate the impact of different synthetic observation types and their associated uncertainties. The synthetic observations were generated by adding random noise to the values extracted from the reference runs, which are considered to be the true snowpack states. The noise was drawn from a normal distribution with a standard deviation equal to the observation uncertainty. The following synthetic observations were assimilated: (1) SWE with uncertainty ranging from 5 % – the best expected accuracy of manual measurements Beaudoin-Galaise and Jutras, 2022) – to 30 % – corresponding to typical uncertainties in radar-based SWE retrievals (e.g. Montpetit et al., 2025; Pomerleau et al., 2020; Cho et al., 2023); (2) backscatter at 13.5 GHz with a 1 dB uncertainty; (3) backscatter at 17.25 GHz with a 1 dB uncertainty; (4) both 13.5 and 17.25 GHz backscatter simultaneously, with a 1 dB uncertainty for each frequency; and (5) the backscatter difference between 13.5–17.25 GHz, with a 1.4 dB uncertainty, as sometimes employed in brightness temperature assimilation to enhance snow information extraction (e.g Larue et al., 2018b). The assumed 1 dB uncertainty reflects the target measurement accuracy of individual TSMM observations. However, this does not include additional errors that would arise in practice, such as spatial mismatch between observation and model grid (representativeness error), inaccuracies in converting model states to backscatter (observation operator errors), or other systematic biases.

Following each assimilation step (i.e. observation time step), a new 100-member ensemble was generated from the assimilated particles using the particle filter's resampling process (Sect. 2.3.3). The forcing perturbations for the new 100 members were generated by drawing new samples from the noise distributions (Sect. 2.3.1).

Figure 2 shows an example of one randomly chosen reference run used to generate synthetic observations for assimilation. In total, 10 reference runs were generated at each site for each winter season (Figs. S1–S3 in the Supplement). By construction of the synthetic experiment, the reference runs were mostly within the generated 100-member ensembles. The snowpack was shallowest in TVC, deeper in Powassan, while the Alpine snowpack in Rogers Pass reached depths exceeding 3 m. The backscatter at 17.25 GHz was slightly higher than at 13.5 GHz. While the backscatter at Powassan and TVC increased with SWE, it saturated early in the season at Rogers Pass, reaching a consistent value that barely changed with SWE.

Figure 2 Spread of the open loop (OL) ensemble composed of a 100 members (between 5th–95th percentiles) and one randomly chosen reference run at (a, d, g, j) Powassan, (b, e, h, k) TVC, and (c, f, i, l) Rogers Pass for (a–c) backscatter at 13.5 GHz, (d–f) backscatter at 17.25 GHz, (g–i) snow depth, and (j–l) SWE for the winter 2020–2021. Similar figures but with all 10 reference runs and for the three winter seasons can be found in the Supplement (Figs. S1–S3).

The spread-skill (ensemble spread divided by RMSE) of the ensembles and the climatological variance condition (mean of the member variances divided by the variance of the synthetic truths at the observation times) were calculated to assess ensemble reliability (Fortin et al., 2014; Dirkson and Buehner, 2025a, b; Johnson and Bowler, 2009). The ensemble size was accounted for in the calculation of the RMSE to eliminate its effect on the spread-skill scores (Dirkson and Buehner, 2025b). The reference runs (the truths) were used in the calculation of these two scores. Table 4 presents the spread-skills over the three winter seasons and the reference runs for backscatter (13.5 and 17.25 GHz), snow depth, and SWE. Most values are close to 1, indicating that the ensembles are reliable and that the ensemble spread accurately captures the predicted errors in ensemble mean.

Table 4 Spread skills and climatological variance conditions over the 10 reference runs and the three winter seasons at each site for backscatter (σ) at 13.5 and 17.5 GHz, SWE, and snow depth.

2.3.3 The Particle Filter

The particle filter with Sequential Importance Resampling (PF-SIR) (Gordon et al., 1993) was used for the assimilation of different snow variables described in Sect. 2.3.2. In the particle filter, the prior distribution (also referred to as background) of the model states x_t at time t, $p\left({\mathit{x}}_t\right|{\mathit{y}}_{\mathrm{1}:t-\mathrm{1}})$, is compared to the observations at time t, y_t, to estimate a posterior distribution, p(x_t|y_t), from calculated weights, w, between the prior distribution and the observations:

$$\begin{array}{}\text{(5)}& {\displaystyle }{\displaystyle }p\left({\mathit{x}}_t\right|{\mathit{y}}_t)\propto \mathit{w}p({\mathit{x}}_t\left|{\mathit{y}}_{\mathrm{1}:t-\mathrm{1}}\right)\end{array}$$

The prior distribution is composed of an ensemble of model states (particles). During the first step of the PF-SIR, the prior particles are weighted based on their distance to the observation. The weight for each ensemble member i is calculated as:

$$\begin{array}{}\text{(6)}& {\displaystyle }{\displaystyle }{w}_i={\displaystyle \frac{\text{exp}\left(-\frac{\mathrm{1}}{\mathrm{2}}[\mathit{y}-\widehat{{\mathit{y}}_{\mathit{i}}}]{\mathbf{R}}^{-\mathrm{1}}[\mathit{y}-\widehat{{\mathit{y}}_{\mathit{i}}}]\right)}{{\sum }_{j=\mathrm{1}}^{N_{\mathrm{e}}}\text{exp}\left(-\frac{\mathrm{1}}{\mathrm{2}}[\mathit{y}-\widehat{{\mathit{y}}_j}]{\mathbf{R}}^{-\mathrm{1}}[\mathit{y}-\widehat{{\mathit{y}}_j}]\right)}}\end{array}$$

where R is the error covariance matrix that accounts for the uncertainty in the observations, N_e is the total number of members, and ${\widehat{y}}_i=\mathit{\text{H}}\left(x_i\right)$ with H(⋅) is the observation operator which maps from the model state to the observation space. During the second step of the PF-SIR, a resampling of the particles is conducted to select the particles with the highest weights. In our study, a systematic resampling is used, selecting particles at evenly spaced intervals along the cumulative weight distribution and is a commonly used approach in particle filter applications (e.g. Kuptametee and Aunsri, 2022). This generates a new ensemble of N_e members (100 members in this study), with selection proportional to particle weights.

The particle filter is known to suffer from degeneracy when all ensemble members collapse to a few particles (e.g. Moradkhani et al., 2012). To monitor degeneracy, the effective sample size (N_eff) was calculated after the calculation of the weights as:

$$\begin{array}{}\text{(7)}& {\displaystyle }{\displaystyle }{N}_{\text{eff}}={\displaystyle \frac{\mathrm{1}}{{\sum }_{i=\mathrm{1}}^{N_{\mathrm{e}}}{w}_i^{\mathrm{2}}}}\end{array}$$

In our study, we considered that degeneracy happened when $N_{\text{eff}}/N_{\mathrm{e}}$ is below 20 %, similarly to Larue et al. (2018a) and Cluzet et al. (2021), which considered thresholds between 14 %–20 %.

2.3.4 Metrics for Evaluation

The RMSE and the continuous ranked probability score (CRPS) were computed to evaluate the performance of the different assimilation experiments. The RMSE is calculated as:

$$\begin{array}{}\text{(8)}& {\displaystyle }{\displaystyle }\mathrm{RMSE}=\sqrt{{\displaystyle \frac{{\sum }_{i=\mathrm{1}}^M(\overline{x}-y_i{)}^{\mathrm{2}}}{M}}}\end{array}$$

where $\overline{\mathit{x}}$ is the mean of the ensemble members at all the observation times, M is the number of observations, and y_i are the values of the reference runs at the observations times.

The CRPS is a verification metric that measures the difference between a predicted probability distribution and an observed value, with lower values indicating better forecast skills (Hersbach, 2000). The CRPS is calculated as:

$$\begin{array}{}\text{(9)}& {\displaystyle }{\displaystyle }{\text{CRPS}}_t=\int \left(F_t\right(x)-O_t(x){)}^{\mathrm{2}}\mathrm{d}x\end{array}$$

where CRPS_t is the CRPS at time t, F_t(x) is the cumulative distribution function (CDF, i.e. the probability that a variable is less than or equal to a given value) of the ensemble forecast at time t and O_t(x) is the CDF of the reference run (the truth) at the time of the observation (a step function that equals 0 for x<observed value and 1 for x≥observed value). The CRPS is the integrated squared difference between the predicted and observed CDFs across all possible values. CRPS_t values were then averaged over time for each assimilation experiment to obtain a mean CRPS, which represents the overall forecast accuracy.

These scores were calculated for the open loop ensemble, the background particles, and the assimilated particles (the analysis). The scores obtained from the open loop served as a baseline for assessing improvements in snow predictions (SWE and snow depth) achieved through assimilation. Comparing the scores between the analysis and the background particles highlighted the performance of the assimilation at each assimilation step. The RMSE was calculated for the mean of the ensemble and the observations to evaluate the performance of the ensemble means, while the CRPS was applied to all ensemble members. A single RMSE and CRPS value were computed for each distinct experimental configuration, which comprised a unique reference run and year.

The RMSE and CRPS of the assimilation results were normalized against those of the open-loop or the background particles. The normalized CRPS against the open loop (CRPS_{norm, OL}) was calculated as:

$$\begin{array}{}\text{(10)}& {\mathrm{\text{CRPS}}}_{\text{norm, OL}}=\mathrm{100}\times ({\text{CRPS}}_{\text{OL}}-{\text{CRPS}}_{\text{exp}})/{\text{CRPS}}_{\text{OL}}\end{array}$$

where CRPS_OL is the CRPS of the open loop and CRPS_exp is the CRPS of the assimilation for a same reference run and winter. Negative values of CRPS_{norm, OL} meant that the assimilation performed worse than the open-loop and positive values showed an improvement of the assimilation over the open-loop. Similarly, we normalized CRPS scores against the background particles (CRPS_norm,back).

We evaluated the simulated vertical profiles of snow density and SSA, both with and without assimilation, against those from the reference runs. First, the height of all vertical profiles was normalized between 0–1 to allow a direct comparison as in Woolley et al. (2024). We then divided this normalized height into equally spaced layers, each with a thickness of 0.005 (Viallon-Galinier et al., 2020). Subsequently, the vertical density and SSA values were compared at each of these 0.005 layers, and the RMSE and CRPS for each layer were averaged per profile. This analysis aimed at quantifying improvements in estimated snowpack vertical properties via the assimilation method compared to the open-loop ensemble members.

3 Results

In this section, we first present the results of the assimilation experiments on bulk SWE and snow depth estimates. These results are compared across the different types of synthetic observations being assimilated. Finally, we examine how well each synthetic observation type improves vertical snow profile estimation (density and SSA) beyond the open loop baseline.

3.1 Estimation of Snow Bulk Properties

Figure 3 presents seasonal evolution of SWE at the three sites in 2022–2023 obtained from three assimilation experiments, using a single reference run across all three sites. These experiments include: assimilation of 13.5 GHz backscatter synthetic observations, assimilation of 17.25 GHz backscatter synthetic observations, and assimilation of SWE synthetic observations with a 10 % uncertainty. Generally, the assimilation experiments showed improved SWE predictions by reducing the spread of the ensembles compared to the open loops. At all sites, backscatter assimilation revealed a spread increasing during the accumulation period and diminishing during the melt season. At Powassan and Rogers Pass, SWE predictions after assimilating backscatter synthetic observations at 17.25 GHz presented a slightly smaller spread than those from backscatter at 13.5 GHz, although this narrower spread did not always encompass the reference run. The assimilation of backscatter synthetic observations at Rogers Pass barely improved SWE estimates compared to the open loop, despite this site having the highest number of weekly observations throughout the winter that passed the selection criteria specified in Sect. 2.3.2. On the other hand, assimilating synthetic SWE observations with an error of 10 % provided an analysis more centered around the reference run with a much reduced spread than when assimilating either backscatter. Figure S4 presents the assimilation results for the same winter season and reference run, but considers additional assimilation configurations: the simultaneous assimilation of both backscatter frequencies, the assimilation of the backscatter frequency difference, and the assimilation of SWE with an observation error of 20 %. For this case, the assimilation results obtained using the various backscatter frequency combinations were largely consistent with those derived from single-frequency assimilation, except for improved SWE estimates at Powassan when both frequencies were assimilated simultaneously compared to assimilating the 13.5 GHz backscatter alone. Furthermore, the assimilation of SWE with a 20 % observation error resulted in a larger ensemble spread in the SWE estimates relative to the case in which a 10 % error was prescribed. Overall, assimilating SWE synthetic observations provided SWE estimates with a narrower spread around the reference run compared to assimilating backscatter synthetic observations.

Figure 3 Results of the assimilation of the reference run #1 for the 2022–2023 snow season at (a, d, g) Powassan, (b, e, h) TVC, and (c, f, i) Rogers Pass. (a–c) are the results of assimilating backscatter synthetic observations at 13.5 GHz, (d–g) for the assimilation of backscatter synthetic observations at 17.25 GHz, and (g–i) for the assimilation of SWE synthetic observations with an error of 10 %. The gray envelope shows the spread between the 5th–95th percentiles of the open loop (OL) without assimilation, while the blue envelope shows the results of the assimilation (5th to 95th percentiles) and the corresponding reference run is in red.

Figure 4 shows the variability in seasonal CRPS scores, where each value represents the mean CRPS for a single reference run over one winter season. In most cases, assimilating the individual backscatter synthetic observations improved SWE estimates compared to the open loop (Fig. 4a). The TVC site consistently shows larger improvements with mean normalized CRPS values of SWE estimates exceeding 30 % compared to Powassan (∼23 %) and Rogers Pass (between 0 %–3 %). The assimilation of backscatter at either frequency provided similar results in terms of mean normalized CRPS, but the assimilation of the backscatter at 17.25 GHz slightly improved the estimates at Powassan (+1.5 % compared to 13.5 GHz for SWE estimates) while the assimilation of backscatter at 13.5 GHz worked better at Rogers Pass (+3 % compared to 17.25 GHz for SWE estimates) and at TVC (+1.5 % compared to 17.25 GHz for SWE estimates). Very little improvements in SWE estimates were found over the open loop at Rogers Pass when assimilating the backscatter at 17.25 GHz. The assimilation of individual backscatter showed greater improvements for SWE predictions at Powassan and TVC than for snow depth predictions (Fig. 4b), with an increase of ∼10 % and ∼6 % of mean normalized CRPS, respectively. The normalized RMSE were qualitatively similar to the normalized CRPS (Fig. S5), showing the positive impact of the assimilation of backscatter to improve the mean prediction accuracy of SWE and snow depth estimates.

Figure 4 Normalized seasonal CRPS relative to the open loop over all the different runs for the three winter seasons (sample size of 30 for each box plot) for (a) SWE prediction and (b) snow depth prediction at the three sites based on different synthetic observations being assimilated. Box plots show median (center line), interquartile range (box), 10th–90th percentiles (whiskers), and mean (+). No outliers are shown for clarity.

Assimilating both frequencies of backscatter synthetic observations demonstrated the highest gains in prediction accuracy at Powassan and TVC, outperforming the assimilation of individual frequencies with an increase in mean normalized CRPS of 10 % for SWE and snow depth estimates at Powassan and between 3 %–5 % at TVC (Fig. 4). This improvement was not systematic at Rogers Pass. In addition, the combined assimilation approach led to a reduction in the spread of normalized scores at both Powassan and TVC (as shown by the size of the boxplots), indicating a more robust and less variable assimilation outcome for these locations. In contrast, Rogers Pass exhibited a wider spread of normalized scores under the dual-frequency assimilation, suggesting less consistent performance. Lastly, assimilating the difference between the frequencies showed slight improvements in normalized scores at Powassan and Rogers Pass for both SWE and snow depth estimates compared to assimilating individual frequencies (+2 %–4 % in mean normalized CRPS).

Assimilating SWE synthetic observations with an uncertainty of 5 % showed the best estimates of SWE and snow depth, with mean normalized CRPS values mostly above 75 % and 60 % for SWE and snow depth estimates, respectively. Improvements at TVC were greater than those of the other sites. As the uncertainty in SWE synthetic observation increased, the improvements of the assimilation compared to the open loop were reduced. Assimilating SWE with uncertainties equal to or below 20 % performed better than assimilating the backscatter synthetic observations, with improvements in mean normalized CRPS commonly above 45 % for SWE estimates and 35 % for snow depth estimates. However, when assimilating SWE with an uncertainty of 30 %, the results were similar to assimilating the backscatter, except at Rogers Pass where it still performed better.

To understand the effectiveness of the assimilation method on reducing forecast uncertainty at each assimilation time step, the CRPS of the analyses were compared to those of the background particles. The performance of the particle filter in narrowing the spread of the background ensemble is evaluated for different periods of the winter season (Fig. 5). At all sites, the assimilation algorithm did not improve on the background particles in the middle of the winter when the SWE at each site was the highest, i.e. between January–February for Powassan and between January–April for TVC and Rogers Pass, which have longer winter periods than Powassan. The assimilation method performed best during the accumulation and melt periods when SWE values were the lowest. It is important to note that the number of observations during the melt period was low as only dry snow conditions were assimilated. In the middle of the winter season, the particle filter failed to reduce the spread of background particles when assimilating the different backscatter observations. Assimilating both backscatter synthetic observations performed best at Powassan across the winter and at TVC early in the winter season. Assimilating the difference of backscatter outperformed the other assimilation of backscatter at Rogers Pass during the melting period, with mean normalized CRPS against the background particles of ∼40 % for both SWE and snow depth estimates. The particle filter performed significantly better across the winter when assimilating the SWE synthetic observations, with greater results for lower SWE uncertainties, but once again, the improvements were lower in the middle of the season. It can be noted that the spread of the SWE estimates from assimilating the synthetic SWE observations sometimes decreased with increasing SWE uncertainty. This highlights a trade-off in assimilation experiments: low uncertainties and therefore high observational constraints maximize potential improvements but increase sensitivity to observation errors (deviation of the observation from the truth), whereas higher uncertainties in the observations sacrifice peak performance for greater reliability across diverse conditions. The normalized CRPS against the open loop (Fig. S6) did not show this dependency on seasonality, meaning that the SWE or snow depth analysis consistently improved upon the open loop across the season. These results suggest that assimilating observations early in the winter had the greatest impact, enhancing SWE and snow depth estimates for the rest of the season.

Figure 5 Normalized CRPS values at each synthetic observation time against the background particles over all the runs and the three winter seasons for (a, c, e) SWE prediction and (b, d, f) snow depth prediction for (a, b) Powassan, (c, d) TVC, and (e, f) Rogers Pass based on the month of the observations. Box plots show median (center line), interquartile range (box), 10th–90th percentiles (whiskers), and mean (+). No outliers are shown for clarity. The number of observations (n) within each boxplot is shown in the graphs.

3.2 Impact of Assimilation on the Vertical Snowpack Properties

The vertical profiles of density and SSA were compared for the different assimilation experiments with the corresponding reference run profiles. Figure 6 shows examples of density and SSA profiles at TVC on 27 December 2021, obtained with the experiments assimilating backscatter synthetic observations at 13.5 and 17.25 GHz and SWE synthetic observations with an uncertainty of 10 %. In this example, assimilating SWE observations provided an estimate of the density and SSA profiles closer to the reference run than assimilating the individual backscatter observations, which only moderately improved the profiles of SSA and density. Figure 7 shows the normalized seasonal CRPS values (relative to the open loop) for density and SSA profiles. For density profiles, improvements relative to the open loop were greatest at TVC compared to the other two sites, with the best estimates obtained when assimilating both frequencies of backscatter synthetic observations (mean seasonal normalized CRPS up to 9 %). Some improvements in density were observed at Rogers Pass when assimilating backscatter synthetic observations at 17.25 GHz and the backscatter difference (mean values of 5 %), while at Powassan, only modest improvements were noted when assimilating the backscatter difference (mean values of 3.5 %). When assimilating SWE synthetic observations, density profiles were generally improved at all three sites, with the largest improvements occurring at Rogers Pass. The range of normalized seasonal CRPS values for density profile estimates at TVC exhibited large variability, indicating inconsistent assimilation performance at this site. For SSA profiles, the distributions of normalized CRPS showed improvements when assimilation backscatter synthetic observations at TVC and Rogers Pass, with positive mean values up to 5 % and 17 %, respectively. At TVC, the improvements were greatest when assimilating backscatter synthetic observations at 17.25 GHz and the backscatter difference. No improvements over the open loop in SSA profiles were observed at Powassan when assimilating backscatter synthetic observations, with degradation occurring when assimilating backscatter synthetic observations at either 13.5 or 17.25 GHz. When assimilating SWE synthetic observations, SSA profiles were generally improved at all three sites, with the largest improvements overall occurring at Rogers Pass and for lower observation uncertainties.

Figure 6 Results of the assimilation of the reference run #1 at TVC on (a–c) the vertical profile of density and for (d–f) SSA on 27 December 2021 12:00 UTC. (a, d) are the results of assimilating backscatter (σ) at 13.5 GHz, (b, e) of assimilating σ at 17.25 GHz, and (c, f) of assimilation SWE with an error of 10 %.

Figure 7 Normalized seasonal CRPS (relative to open loop) for (a) vertical density profiles and (b) vertical SSA profiles, at each of the three sites (Powassan, TVC, and Rogers Pass) aggregated over the three winter seasons and the 10 reference runs (sample size of 30 for each boxplot). Box plots show median (center line), interquartile range (box), 10th–90th percentiles (whiskers), and mean (+). No outliers are shown for clarity.

4 Discussion

4.1 Backscatter Assimilation Improves Snow Estimates

The SWE estimates showed positive improvements after data assimilation of backscatter synthetic observations, with site-specific variations across sites representing different climate zones in Canada. At TVC, located in the Arctic, SWE estimates from individual backscatter frequency assimilation demonstrated most promising results, with mean CRPS improvements ranging between 30 %–35 % over the open loop ensemble. In contrast, Powassan, situated in a humid continental climate, showed slightly lower improvements, with improvements of ∼23 % over the open loop of mean CRPS. These differences can be attributed to several key factors. First, the SVS2-Crocus configuration at TVC used in this study came from Woolley et al. (2024) as its best represented the vertical profiles of density and SSA at the Arctic site, while the default SVS2-Crocus configurations were specified at Powassan. In addition, TVC experiences colder winter conditions, providing between 20–30 valid observations during the snow season. On the other hand, the snowpack at Powassan experiences more frequent melt and rain-on-snow events, limiting synthetic observation availability to between 9–16 data points per winter. Consequently, the higher number of assimilated synthetic observations at TVC compared to Powassan may have contributed to the improved SWE and snow depth estimates. The winter conditions at Powassan resulted in frequent melting and precipitation events producing complex vertical snow profiles, particularly with the formation of ice lenses. While these features may not have impacted the assimilation results in this synthetic experiment, it is expected that in real-world data assimilation, complex vertical snow layering due to melt events significantly impact microwave signals and pose considerable challenges for radiative transfer model simulations (Bartsch et al., 2007; Dolant et al., 2018; Picard et al., 2018; Montpetit et al., 2013). Also, backscatter assimilation results may have been affected when assimilation occurred during or shortly after snowfall, since backscatter has little sensitivity to fresh snow particles. Indeed, fresh snow consists of particles with small optical diameters (Mätzler, 2002), poorly scattering microwaves (Chang et al., 1976).

Assimilation strategies combining the individual backscatter frequencies were tested as often done in snow assimilation studies (e.g. Revuelto et al., 2021; Durand and Margulis, 2006) but yielded mixed results. At Powassan, combining both backscatter frequencies greatly enhanced SWE and snow depth estimates, with mean normalized CRPS improved by ∼10 %. However, this approach did not consistently improve results across all sites. Assimilating the frequency difference only slightly improved snow depth estimates at TVC. These findings reflect the complex and sometimes contradictory results observed in previous studies. Larue et al. (2018a), Larue et al. (2018b) encountered mixed results when assimilating brightness temperatures at different frequencies, with improvements in SWE RMSE when assimilating the difference of the frequencies in their synthetic experiment but the opposite when using real data.

Rogers Pass, located in an Alpine climate, presented unique challenges for data assimilation. Despite having the highest number of synthetic backscatter observations (between 24–29), SWE estimates showed little improvement over the open loop ensemble when assimilating backscatter. This seemed to be caused by backscatter signal saturation over a SWE threshold of ∼300 kg m⁻² (Fig. 3). Backscatter simulated by SMRT fed by SVS2 snow outputs does not depend only on SWE, but also on the vertical snow profile properties, particularly the SSA. This low threshold of SWE above which backscatter saturates can be caused by several factors: (1) the SSA in SVS2/Crocus might be under-estimated in the simulations, which can then led to overestimated modeled volume scattering in SMRT (Woolley et al., 2025; Vionnet et al., 2025a) and (2) SMRT needs further testing and development to estimate backscatter from deep alpine snowpacks at the considered frequencies. Future studies will investigate the SWE threshold at which backscatter saturates based on SSA profiles. In addition, this experiment should be conducted with cross-polarization backscatter that could have a stronger response to SWE than co-polarization (Borah et al., 2022). Ongoing radar tower-based field experiments tend to indicate that the radar signal, even at 17.25 GHz, can penetrate the maximum SWE observed at Rogers Pass as shown by Madore et al. (2023) who were able to retrieve SWE values up to a 1000 kg m⁻² with radar measurements at 24 GHz.

The assimilation results from assimilating backscatter synthetic observations were compared to assimilating SWE synthetic observations. Assimilating SWE synthetic observations with uncertainties less than 10 % provided the best estimates of SWE and snow depth at all sites; however, such small uncertainties are only realistic when observations are collected in-situ (e.g. Beaudoin-Galaise and Jutras, 2022). From SWE retrieval algorithms, larger uncertainties can be expected, and can be particularly high with large forest cover fractions (Foster et al., 2005; Cho et al., 2020; Pflug et al., 2024; Bonnell et al., 2024). When assimilating SWE with uncertainties below 20 %, the results in terms of bulk SWE and snow depth estimates were better than when assimilating the backscatter synthetic observations. However, the results were when assimilating SWE synthetic observations with an uncertainty of 30 %, except at Rogers Pass where it still performed better. This decrease in the quality of SWE estimates as SWE observation uncertainty increases is consistent with Shrestha and Barros (2025b), who also looked at SWE retrieval uncertainties between 5 %–30 %.

The assimilation of backscatter synthetic observations demonstrated site-dependent effectiveness in improving vertical snow property estimates. For density profiles, TVC showed the most promising results, with assimilation of both backscatter frequencies simultaneously outperforming single-frequency assimilation. Assimilating the difference of backscatter frequency improved the density profiles estimates at both Powassan and Rogers Pass. For SSA profiles, Rogers Pass exhibited the strongest performance, with backscatter assimilation achieving improvements comparable to direct SWE assimilation. The improvements in density and SSA profiles at TVC are particularly encouraging, as current multi-layered snowpack models struggle with representing Arctic snowpack stratigraphy (e.g. Woolley et al., 2025; Vionnet et al., 2025a). The limited improvements at Powassan, even when assimilating SWE synthetic observations, suggest challenges in constraining SSA at this site. Overall, the improvement of estimated vertical snow properties after assimilating either backscatter or SWE synthetic observations is consistent with findings from Shrestha and Barros (2025a).

4.2 Performance of the Particle Filter Algorithm

Despite the limited performance of the particle filter for high SWE values, the particle filter algorithm proved to greatly reduce the ensemble spread composed of the background particles when SWE values at the observation times were below 200 kg m⁻² (Fig. S7). Noticeable improvements of mean CRPS often above 40 % over the background ensemble were found for SWE values below 50 kg m⁻² at Powassan and Rogers Pass. This translated in the assimilation algorithm performing best during the early accumulation period and the melting period (Fig. 5). Similarly, Revuelto et al. (2021) showed that assimilating MODIS-reflectance with the particle filter performed better for shallower snowpacks. This highlights the value of assimilating observations early in the winter season, during the beginning of the snowpack accumulation phase. Early assimilation helps establish more accurate initial conditions, which in turn improves the effectiveness of assimilation throughout the season, especially during mid-winter.

Assimilation with the particle filter is inherently prone to degeneracy. Figure S8 illustrates this, presenting the effective sample sizes (N_eff, Eq. 7) averaged over the 10 reference runs across different assimilation time steps, experiments, and sites. Overall, the particle filter was not prone to degeneracy during the middle of the snow season, with N_eff values usually above 70 %, signifying good particle diversity and greater consistency between observations and ensemble predictions. In contrast, we observed that the algorithm typically degenerated either early in the snow season or towards the end of the melt season, with N_eff values dropping below the 20 % threshold (Sect. 2.3.3). These periods of low N_eff largely coincided with ephemeral snow events occurring before or after the main seasonal snowpack (Figs. 2, S1–S3). This is caused by two main factors: (i) precipitation events leading to ensemble members with different precipitation amounts and phases, widening the ensemble spread and reducing the number of particles close to the observations to be resampled, and (ii) melt events generate wet snowpacks, reducing the number of synthetic observations to be assimilated and modifying soil moisture content, which strongly affects backscatter values and introduces additional variability in the ensemble, further reducing the number of particles close to the observations. These factors explain why effective sample sizes for backscatter assimilation are sometimes lower than those for SWE assimilation. To mitigate degeneracy, systematic resampling was employed within the PF-SIR, and new perturbations were applied after each assimilation time step. This strategy allowed the algorithm to successfully recover by propagating new particles with increased ensemble spread. On average, assimilating the individual backscatter observations yielded slightly higher N_eff values (+3 to +6) than assimilating SWE with 10 % uncertainty, indicating that SWE observations were more discriminating and informative for constraining the particle ensemble.

4.3 Current Limitations

The assimilation method used in this study presents some limitations. The parameterizations of the different snow processes simulated by SVS2/Crocus were kept the same when generating the different members of the open loop ensembles, the reference runs, and the ensembles in the assimilation experiments. This was based on the assumption that meteorological forcings are the main source of uncertainty in snowpack modelling (Günther et al., 2019; Raleigh et al., 2015). Future studies will consider an ensemble that accounts simultaneously for the uncertainties in the meteorological forcing as well as the snowpack model as done in Cluzet et al. (2021), Deschamps-Berger et al. (2022). This could improve the assimilation of backscatter observations at the Alpine site and reduce the backscatter saturation with observed SWE by having potentially more realistic SSA vertical profiles. The assimilation experiments were limited to 100 members as it has been found to be suitable for snow assimilation experiments with the particle filter (e.g. Piazzi et al., 2018). Although some studies have shown the possibility to apply other assimilation methods with multi-layered snowpack models, such as the Ensemble Kalman Filter and the 1D-Var (Shrestha and Barros, 2025b; Phan et al., 2014), they were not considered in our study. Finally, a systematic bias in the observations was not considered, but it may affect the assimilation outcomes (e.g. Revuelto et al., 2021). When assimilating actual data, it may be necessary to implement a bias correction procedure for the observations.

5 Conclusions

This study investigates the potential of assimilating backscatter observations at two frequencies (13.5 and 17.25 GHz) in preparation for the TSMM satellite mission (Derksen et al., 2021). The synthetic assimilation experiments were conducted at three different sites spanning different Canadian climates, including an Arctic site (TVC), a humid continental site (Powassan), and an Alpine site (Rogers Pass). To test the assimilation for different snowpack conditions within each climate regime, three winter seasons at each site were considered, each with 10 random reference runs, which are assumed to represent the true snowpack state, from which the synthetic observations were extracted at a weekly interval (proposed measurement frequency of TSMM) when the snowpack was dry. The results of the synthetic experiments are as follows:

• Assimilating individual backscatter improves SWE and snow depth estimates at all sites, with mean normalized CRPS against the open loop of SWE estimates up to 32 % at TVC, up to 23 % at Powassan, and up to 3 % at Rogers Pass.

• Assimilating both frequencies at the same time performed better at Powassan and TVC with a mean normalized CRPS of SWE estimates improved by ∼10 % compared to the assimilation of individual frequencies.

• Assimilating the difference of the frequencies at the same time performed better at Rogers Pass with a mean normalized CRPS of SWE estimates over the open loop of 5 %.

• The SWE estimates obtained with backscatter assimilation were comparable with estimates from assimilating SWE observations with high uncertainties (>20 %), which can be expected from radar-based SWE retrievals (e.g. Pomerleau et al., 2020; Shrestha and Barros, 2025a).

• This study also showed that assimilating backscatter or combination of backscatter frequencies can improve the estimates of snow vertical profile properties, such as density and SSA, in some cases by up to ∼12 % in terms of improved mean CRPS over the open loop.

This study provides initial insights into the assimilation of backscatter data directly within a snowpack model to improve predictions of SWE, snow depth, and vertical snow properties. The coupling of SVS2/Crocus and SMRT within the MuSA framework facilitates experimentation and advancements in snow data assimilation. Consequently, our work establishes the foundation for an assimilation scheme tailored to future Ku-band SAR missions. A major limitation found in this study is the saturation of backscatter for SWE values above ∼300 kg m⁻², resulting in quasi-constant backscatter values for SWE values above this threshold. This resulted in the limited improvements in SWE and snow depth estimates shown at the Alpine site. Work is currently done to overcome this limitation, particularly to improve the parameterization of SSA in SVS2 as well as to further test SMRT against field observations.

In the last few years, field studies have been conducted to gather data in preparation of TSMM. Kelly et al. (2024) obtained measurements of backscatter at 13.5 GHz over the Powassan site during the 2022–2023 winter with the CryoSAR in conjunction with in-situ snow measurements of bulk snow properties and vertical snow properties. Similar data are available at TVC (Montpetit et al., 2024, 2025). The next steps of this study will include using these field data to test the methodology with real data at point-scale, comparing the performance of assimilating direct backscatter against retrieved SWE following the method of Montpetit et al. (2025), running an OSSE over a 2D domain using a simulator created to generate synthetic data using the TSMM satellite orbit, and developing the 2D assimilation of backscatter within the Canadian Land Data Assimilation Scheme (CaLDAS) (Carrera et al., 2015).

Appendix A: Parameterization Used for the Crocus Simulations

Table A1 Parameterization used for the SVS2/Crocus simulations. See Lafaysse et al. (2017), Vionnet et al. (2012), and Woolley et al. (2024) for a description of the parameterizations.