Subgridding high-resolution numerical weather forecast in the Canadian Selkirk mountain range for local snow modeling in a remote sensing perspective

Snow water equivalent (SWE) is a key variable in climate and hydrology studies. Yet, current SWE products mask out high-topography areas due to the coarse resolution of the satellite sensors used. The snow remote sensing community is hence pushing towards active-microwave approaches for global SWE monitoring. Designing a SWE retrieval algorithm is not trivial, as multiple combinations of snow microstructure representations and SWE can yield the same radar signal. Retrieval algorithm designs are converging towards forward modeling approaches using an educated first guess on the snowpack structure. Snow highly varies in space and time, especially in mountain environments where the complex topography affects atmospheric and snowpack state variables in numerous ways. In Canada, automatic weather stations are too sparse, and high-resolution numerical weather prediction systems have a maximal resolution of 2.5 km × 2.5 km, which is too coarse to capture snow spatial variability in a complex topography. In this study, we designed a subgridding framework for the Canadian High Resolution Deterministic Prediction System (HRDPS). The native 2.5 km × 2.5 km resolution forecast was subgridded to a 100 m × 100 m resolution and used as the input for snow modeling over two winters in Glacier National Park, British Columbia, Canada. Air temperature, relative humidity, precipitation, and wind speed were first parameterized regarding elevation using six automatic weather stations. We then used Alpine3D to spatialize atmospheric parameters and radiation input accounting for terrain reflections, and we performed the snow simulations. We evaluated modeled snowpack state variables relevant for microwave remote sensing against simulated profiles generated with automatic weather station data and compared them to simulated profiles driven by raw HRDPS data. The subgridding framework improves the optical grain size bias by 18 % on average and the modeled SWE by 16 % compared to simulations driven with raw HRDPS forecasts. This work could improve the snowpack radar backscattering modeling by up to 7 dB and serves as a basis for SWE retrieval algorithms using forward modeling in a Bayesian framework.

1 Introduction

Seasonal snow governs several feedback loops that directly affect our planet's climate and plays a major role in its hydrological dynamics. With its high albedo, snow reflects a large proportion of the incoming solar radiation, which in turn helps to mitigate global warming (IPCC, 2019). Furthermore, snow insulates the underlying soil, affecting the microbial activity, carbon fluxes, and permafrost freeze–thaw cycles (Natali et al., 2019; Biskaborn et al., 2019). Seasonal snowmelt provides connected watersheds with freshwater, sustaining natural ecosystems and human infrastructure. Finally, extreme precipitation events and resulting snowmelt can cause devastating floods (Pomeroy et al., 2016; Vionnet et al., 2020), so managing runoff would highly benefit both society and the economy (Sturm et al., 2017). Yet, snow mass (or snow water equivalent, SWE) remains poorly characterized, especially in mountainous regions where a significant amount of SWE is stored at the continent scale (Wrzesien et al., 2018). Global SWE products inferred from passive microwave observations are available at a 25 km resolution (Luojus et al., 2021), which is too coarse to capture the SWE spatial variability (Derksen et al., 2021), so mountains are simply omitted or masked out. Both observations from passive microwaves and modeling efforts yield negative biases when estimating mountain or deep-snow SWE on the global scale (Vuyovich et al., 2014; Wrzesien et al., 2018; Pulliainen et al., 2020). Hence, the snow remote sensing community is promoting active remote sensing, which provides higher-spatial-resolution information compared to passive microwave products (Tsang et al., 2022; Rott et al., 2010; Derksen et al., 2021). The sensitivity of the synthetic aperture radar (SAR) signal to SWE has been proven at the Ku-band (King et al., 2015; Lemmetyinen et al., 2016). Recent studies suggest that the C-band could also be used for snow depth retrieval (Lievens et al., 2019, 2022), a key parameter for SWE retrieval, although it is contrasting with previous research (Dozier and Shi, 2000). Linking SWE to SAR backscattering is complicated as it depends on more than solely SWE (which is a function of snow height and density) and also on the snow microstructure. Several combinations of SWE and snowpack microstructures can yield similar backscattering values, creating multiple inversion solutions (Tsang et al., 2022). As a result, recent inversion algorithms tend towards a Bayesian framework where a forward scattering model is used to generate possible backscattering values, and a weighted cost function allows finding the model that fits best (Lemmetyinen et al., 2018; King et al., 2018; Zhu et al., 2021). So far, these studies only paired airborne radar observations with fields measurements, but coupling a radiative transfer model with a snow physics model still has to be explored in the active-microwave domain.

Advanced thermodynamic multilayered snow models such as Crocus or SNOWPACK produce SWE and microstructure parameter estimates (Brun et al., 1992; Vionnet et al., 2012; Lehning et al., 2002). Such models can be driven by either automatic weather station (AWS) measurements, atmospheric models, or reanalysis products. On the one hand, weather stations provide very accurate measurements of the atmospheric conditions at the local scale. However, they need human maintenance, are subject to outages and local biases, and usually undersample the spatial heterogeneity of the processes at stake, especially in complex terrain. As a result, AWS spatial interpolation in mountainous areas is not always accurate (Lundquist et al., 2019). On the other hand, the Canadian High Resolution Deterministic Prediction System (HRDPS) (Milbrandt et al., 2016) produces numerical forecasts at 2.5 km resolution and is known for its negative bias in precipitation, yielding a negative bias for snow depth and SWE (Bellaire et al., 2011, 2013; Côté et al., 2017).

Several numerical weather prediction (NWP) downscaling schemes have already been proposed. Liston and Elder (2006) introduced the MicroMet model, which is now widely used and is part of several more recent models. In MicroMet, a high-resolution DEM (30 m to 1 km) is used to generate the overlying atmospheric forcing from a coarser grid or a sparse network of automatic weather stations. This allows a physically sound downscaling to be produced when compared to naive interpolation methods but without the need to run a computationally intensive fully dynamic atmospheric model at the local scale. In MicroMet, lapse rates are used for air temperature, dew point temperature (for relative humidity), and precipitation. The algorithm for wind speed takes terrain slope and curvature into account. Incoming solar radiation is split between direct and diffuse radiation and adjusted with cloud cover and terrain shading. Fiddes and Gruber (2014) developed the TopoSCALE model, which can be seen as an iteration over the MicroMet model. The main difference with MicroMet lies in the precipitation subgridding that considers wind redistribution by altering the precipitation field with climatology data after applying the lapse-rate correction from Liston and Elder (2006). In the Canadian Hydrology Model (CHM), Marsh et al. (2020) take this idea one step further by adding snow modeling to the atmospheric model subgridding. In their study, the high-resolution DEM used for subgridding is first transformed into an unstructured triangular mesh (or “triangulated irregular network”, TIN). The input meteorology can be either real AWSs or an array of “virtual stations” extracted from any atmospheric model and defined by latitude, longitude, and elevation. The provided atmospheric forcing is spatialized over the study area and then used to run snow simulations using either iSnobal, SNOWPACK, or Crocus (Marks et al., 1999; Lehning et al., 2002; Brun et al., 1992). Vionnet et al. (2021) used the CHM with a novel wind-downscaling strategy to subgrid forecasts from the HRDPS and simulate snow conditions at 50 m during one snow season using two-layer iSnobal within CHM as the snow model. With snow hydrology as a main application, the evaluation for CHM in Marsh et al. (2020) and Vionnet et al. (2021) is naturally focused on SWE and snow depth. However, for remote sensing applications, specifically for the SAR signal inversion, snow microstructure and layering truly matter (King et al., 2018; Zhu et al., 2021; Tsang et al., 2022). Remote sensing products are written in a gridded raster format; the TIN mesh used in CHM, although very efficient, becomes problematic when pairing the model's output with satellite imagery. The Alpine3D model is a spatially distributed 3D model, which allows the vertical 1D multi-layer snow model SNOWPACK to be run over a gridded DEM, considering the spatial processes affecting atmospheric variables (Bartelt and Lehning, 2002; Lehning et al., 2002, 2006). Weather data are spatialized using the MeteoIO library (Bavay and Egger, 2014). However, MeteoIO is geared towards AWS spatialization, and it excludes an atmospheric model subgridding scheme. This highlights the current community needs for both the design and the evaluation of an atmospheric model subgridding framework to perform snow modeling in a SAR remote sensing coupling context. Our research should answer the following questions:

1. How do subgridded HRDPS forecasts compare to reference automatic weather stations in the simulation domain?

2. Do the resulting atmospheric forcings lead to an improvement in snowpack modeling, especially for critical snow parameters in remote sensing applications?

3. Which degree of spatial variability with regards to snow parameters can be reached by such a subgridding framework?

We thus first built a subgridding module to downscale HRDPS grids as a virtual weather station array. Second, we spatialized atmospheric parameters and performed snow simulations on the study area using the Alpine3D model over two consecutive winters (2018–2019 and 2019–2020). Weather parameter subgridding and snowpack state parameters were assessed at three reference weather stations using an array of statistical criteria and a dynamic time warping algorithm (Hagenmuller and Pilloix, 2016; Hagenmuller et al., 2018; Herla et al., 2021). Finally, we assessed the spatial variability capacity of the proposed subgridding framework over the whole simulation domain and within one HRDPS grid cell.

2 Study area

This study was conducted in the Rogers Pass area of Glacier National Park (GNP), British Columbia, Canada (Fig. 1), which is part of the Selkirk range in the Columbia Mountains. The pass is used as a transportation corridor by the Trans-Canada Highway and the Pacific Railway, making it the busiest transport corridor in western Canada (Bellaire et al., 2016). The pass is exposed to 144 avalanche paths, and as a result, Rogers Pass hosts the largest avalanche control operation in Canada (Delparte et al., 2008). The operation has been ongoing since 1965, and the site has been used for snow research ever since, making this area the longest record of mountain snow in western Canada (Fitzharris, 1987; Bellaire et al., 2016; Madore et al., 2022). The study area is 18 km by 16 km, covering 288 km² of complex topography, with elevations ranging from 840 m a.s.l at the valley bottom to 3284 m a.s.l. In winter, the Columbia Mountains snowpack is characterized as a transitional snowpack with a maritime influence. Hence, westerly fluxes coming from the Pacific mainly govern the precipitation pattern. Occasionally, dryer and colder systems from the northeast hit the range, bringing some continental influence to the east of the study area. On average, the snowpack reaches 3.2 m at its peak, usually around the end of March and early April.

Figure 1 Glacier National Park, British Columbia, Canada.

The Park has eight automatic weather stations (AWSs) at different elevations around the Highway corridor. However, Heather Hill and Rockfall stations are outside of the main study area and are more subject to long outages. As a result, we removed them from the AWS set used in this study. The measured variables are air temperature (TA, °C), relative humidity (RH, %), wind speed (VW, m s⁻¹), wind direction (DW, degrees), precipitation (PSUM, mm), incoming longwave radiation (ILWR, W m⁻²), and incoming shortwave radiation (ISWR, W m⁻²). Some stations include a snow height (HS, cm) sensor. Table 1 summarizes the set of meteorological variables available for each AWS. This work mainly aims at providing a realistic first guess of the snowpack structure in the context of SAR remote sensing signal inversion algorithm development. At relevant frequencies (Ku-band, X-band, C-band), the snowpack becomes opaque to microwaves when wet. As a result, this study focuses on the accumulation period, and we used two winter time series: from September to April, 2018–2019 and 2019–2020. The 2018–2019 season had overall colder temperatures and was relatively dry. The 2019–2020 season had milder temperatures and abundant precipitation. As a result, in 2019–2020, the snowpack was deeper and mostly composed of rounded grains where, in 2018–2019, the shallower snowpack and colder temperatures led to a mostly faceted snowpack. Both seasons had rain-on-snow episodes in the early season, which created melt–freeze crusts at the bottom of the snowpack.

Table 1 Inventory of instruments on the weather stations of the study site. TA stands for air temperature, RH for relative humidity, VW for velocity of wind, DW for direction of wind, PSUM for precipitation water equivalent, ILWR for incoming longwave radiation, and ISWR for incoming shortwave radiation.

3 The numerical weather prediction downscaling processing chain design

Figure 2 summarizes the numerical weather prediction (NWP) downscaling processing chain.

Figure 2 Numerical weather prediction downscaling scheme flowchart.

3.1 HRDPS subgridding and Alpine3D simulations

The High Resolution Deterministic Prediction System (HRDPS) produced by the Meteorological Service of Canada provides a 2.5 km gridded hourly forecast of atmospheric variables for most of Canada (Milbrandt et al., 2016). Atmospheric variables are computed for each pixel at a reference elevation provided by the underlying 2.5 km resolution digital elevation model. This study is based on a grid of 70 HRDPS prediction cells overlying the study area and including the following variables: TA, RH, VW, DW, PSUM, ISWR, and ILWR. First, using the 20 m Canadian Digital Elevation Model (CDEM), we transformed each HRDPS cell data into a virtual weather station. To do so, each cell's centroid coordinates were recomputed, minimizing the hypotenuse distance between each underlying CDEM pixel centroid and the HRDPS centroid, in an 800 m radius around the original HRDPS centroid. Then, TA, RH, PSUM, and ILWR were parameterized to correct for the model's biases and account for the elevation discrepancy between the HRDPS cell elevation and the new centroid CDEM elevation. VW was parameterized to account for the topography underlying the 2.5 km resolution grid unresolved by HRDPS.

Based on all weather stations in the park, the bias in air temperature has a nonlinear relationship with the elevation difference between the station elevation and the original HRDPS cell elevation over the 2018–2020 period. We generated a training set by randomly selecting 75 % of this dataset uniformly across elevations, and the remaining 25 % served as validation set. The data were transposed into logarithmic space to perform a linear regression. The resulting logarithmic fit was then applied over the TA dataset when the elevation difference between the virtual weather station and its overlying HRDPS cell was over 100 m.

$$\begin{array}{}\text{(1)}& {\text{TA}}_{\mathrm{p}}=\left\{\begin{array}{ll}{\text{TA}}_{\mathrm{hrdps}}+\ln (-\mathrm{0.012}{\mathrm{\Delta }}_{\mathrm{E}}+\mathrm{6.89})& \text{if }\left|{\mathrm{\Delta }}_{\mathrm{E}}\right|>\mathrm{100}\\ {\text{TA}}_{\mathrm{hrdps}}& \text{otherwise }\end{array}\right.\end{array}$$

Here Δ_E corresponds to the elevation difference between the CDEM cell and the HRDPS cell, TA_p is the parameterized air temperature, and TA_hrdps is the raw HRDPS air temperature.

RH was corrected by first converting relative humidity to dew point temperature, as described in Liston and Elder (2006). This dew point temperature was then adjusted using the logarithmic fit presented above and converted back to relative humidity.

Snowfall was first parameterized using an elevation lapse-rate correction. This lapse rate was computed by performing a simple linear regression of precipitation as a function of elevation. We used a dataset of 4 weeks of manual SWE measurements on four conventional HN24 precipitation boards placed between 1330 and 1920 m at Mount Fidelity, all placed in flat and open areas sheltered from the wind.

$$\begin{array}{}\text{(2)}& {\text{PSUM}}_{\mathrm{p}}={\text{PSUM}}_{\mathrm{hrdps}}+\mathrm{0.0011}\times {\mathrm{\Delta }}_{\mathrm{E}}\times {\text{PSUM}}_{\mathrm{hrdps}}\end{array}$$

Finally, the HRDPS ILWR was downscaled using the lapse-rate correction as highlighted by Marty et al. (2002), and VW was downscaled to the 20 m CDEM resolution at each new centroid position using the sky view factor approach (Helbig and Löwe, 2014; Helbig et al., 2017).

These virtual weather stations were then spatially interpolated on a 100 m grid via MeteoIO (Bavay and Egger, 2014) using the CDEM grid resampled to 100 m. TA was spatialized using a simple lapse rate computed from the AWS data and inverse distance weighting (IDW). RH, VW, and DW were spatialized using the MicroMet algorithms described in Liston and Elder (2006). To spatialize precipitation, we used topographic parameters and prevailing winds to alter the precipitation field to account for wind snow redistribution (Winstral et al., 2002). Incoming shortwaves were spatialized considering terrain shading, slopes, and reflections from neighboring cells. Finally, ILWR was spatialized using IDW. All the spatial interpolation algorithms mentioned above are a part of the MeteoIO library, which is integrated into the Alpine3D model.

Alpine3D is a spatially distributed 3D model, which allows running the vertical 1D snow model SNOWPACK over an area while taking into account the spatial processes affecting the input atmospheric variables, such as terrain shadowing (Lehning et al., 2006). SNOWPACK is a detailed multi-layer thermodynamic finite-element model of snow microstructure and metamorphism. In this model, the snow microstructure is represented by four main variables: grain size, bond size, dendricity, and sphericity for each snow layer. In addition, the model simulates several metrics of interest when monitoring the evolution of the snowpack, such as height of snow, SWE, density, optical grain size, or snow temperature (Bartelt and Lehning, 2002; Lehning et al., 2002). To do so, the model is fed with three text files describing the weather parameters on the time domain of the simulation, the initial state of the soil layers on which the snow is going to develop (and initial snow layers if relevant), and finally the configuration of the simulation. Alpine3D uses a DEM and a land-use layer to properly initiate each SNOWPACK cell. Depending on the land-use category each cell falls into, canopy information is provided for forested cells to represent snow interception and forest snow processes. The model was run at Rogers Pass on the same 100 m grid described above, over an area of 18 km × 16 km (288 km²) centered on the Highway 1 corridor for winters 2018–2019 and 2019–2020. The snowdrift scheme was turned off, and we generated outputs for three reference stations: Fidelity, Hermit, and Abbott. To assess the spatial variability capacity of the subgridding framework, the model was run on the whole simulation domain, and we also generated outputs at six points within the same cell for intra-cell variability assessment. The specific cell was chosen because it is the only cell in the simulation domain that features a north and south slope with elevations ranging from below the treeline to the alpine on both aspects. No glacier is present in the area. Table 2 summarizes the topographic characteristics for the chosen intra-cell spatial variability points.

Table 2 Summary of the topographic characteristics for the chosen intra-cell spatial variability points. BTL stands for below treeline, TL for treeline, and ALP for alpine.

3.2 Validation data and atmospheric parameter subgridding evaluation

To compare with snow simulations driven by the raw HRDPS and the subgridding framework, we performed SNOWPACK simulations driven with AWS data at Fidelity, Hermit, and Abbott stations. We filtered weather station data to remove outliers; data gaps smaller than 6 h were linearly interpolated, while larger gaps were filled using parameterized HRDPS data. Note that Abbott station is located in a thunderstorm-prone area. Hence, the station is shut down all summer and is only turned back on mid-October. Parameterized HRDPS data were again used to fill in this gap. Moreover, SR50 snow depth measurements at Fidelity and Abbott stations were compared against each snow simulation approach.

The numerical weather forecast subgridding quality was statistically assessed using three well-known criteria: bias, mean absolute error (MAE), and Spearman R correlation coefficient. These indicators allowed respectively quantifying the systematic difference between the models and the ground-truth measurements at the AWSs, the prediction accuracy, and the strength of the association between modeled variables and the ground truth. To smooth small time lags between modeled meteorological events and measurements at the AWSs, we averaged the meteorological time series over 2 h time steps, and reaccumulated precipitation over the same period.

3.3 Snow modeling evaluation

Dynamic time warping (DTW) is an algorithm developed to measure the similarity between two sequences. In a nutshell, DTW computes the optimal match between two signals while allowing for an elasticity in time (or space, in the case of snow profiles). It first resamples the two sequences on 1D grids of the same elemental size and length. Then, a local cost matrix D is built, summarizing the distance between every elemental pair. From there, an accumulated cost matrix G is built by computing the accumulated cost to iterate from one element of D to the next one, respecting a predefined constraint set. The optimal alignment is found by minimizing the alignment accumulated cost.

Although originally designed for speech recognition (Sakoe and Chiba, 1978), DTW is extensively used in time series analysis, and it has recently received increased interest in the snow community (Hagenmuller and Pilloix, 2016; Hagenmuller et al., 2018; Viallon-Galinier et al., 2020; Herla et al., 2021). In the snow science community, DTW has only been used so far in an avalanche forecasting perspective, focusing on aligning standard snow parameters (e.g., grain type, hardness, liquid water content). In this study, we present a new development to the open-source DTW snow profile alignment package written by Herla et al. (2021), allowing aligning snow profiles on remote-sensing-oriented snow parameters, namely layer density and optical grain size (OGS), two key parameters in snow radiative transfer modeling. To do so, an alternative cost function was added to compute the local cost matrix D.

$$\begin{array}{}\text{(3)}& D_{i,j}=w_{\mathrm{d}}{d}_{\mathrm{d}}(q_i^{\mathrm{d}},r_j^{\mathrm{d}})+w_{\mathrm{ogs}}{d}_{\mathrm{ogs}}(q_i^{\mathrm{ogs}},r_j^{\mathrm{ogs}}),\end{array}$$

where w_d and w_ogs are averaging weights respectively applied to density and OGS ($w_{\mathrm{d}}+w_{\mathrm{ogs}}=\mathrm{1}$), $r_n^k$ denotes the nth element of the reference profile R, and $q_n^k$ denotes the nth element of the query profile Q. Finally, d_d() and d_ogs() correspond to the distance function for density and OGS respectively, which is simply the absolute difference between the two elements, normalized over the entire vertical profile.

The space elasticity in the alignment algorithm allows us to find the best match for a layer of the query profile in the same depth range in the reference profile. The algorithm constraints define the amount of elasticity allowed, i.e., the warping window definition and the local slope constraint. These constraints are essential to keep the algorithm from degenerating and from generating irrelevant alignments. However, HRDPS tends to underestimate precipitation in mountain environments (Bellaire et al., 2011, 2013; Côté et al., 2017). Hence, profiles generated from atmospheric models and ground truth profiles can significantly diverge in HS. Physically matching layers are then too far apart, according to the algorithm's constraints, preventing the algorithm from generating relevant matches. Therefore, we artificially inflated profiles modeled using NWP, each layer being multiplied by the height ratio with the station profile. This allowed us to rely solely on snow microstructure parameters for the alignment, assessing only the microstructure representation. The generated aligned (or warped) profile was then used to compute a mean bias for density and OGS with respect to the ground truth. As precipitation is usually underestimated by HRDPS, HS should be underestimated as well, which should impact the overburden pressure on basal layers. A small negative bias on density might result with regards to AWS-driven SNOWPACK runs, depending on the amount of missing snow. For OGS, the temperature gradient in this region is low and metamorphism mainly happens through gravitational settling, leading to little variability in OGS in the snowpack (Madore et al., 2018). As a result, we do not expect much impact of the inflation approach on this microstructure parameter, as the main discrepancies should come from offsets in rain-on-snow modeling and melt/percolation events. Height of snow was compared to the station's SR50 measurements when available, and the Nash–Sutcliffe model efficiency coefficient (Nash and Sutcliffe, 1970) allowed us to assess the SWE modeling quality using HRDPS and subgridded HRDPS data versus the station runs over the two seasons. For more details on the DTW implementation used in this paper, see Herla et al. (2021).

4 Results

4.1 Numerical weather forecast subgridding performance

Figure 3 summarizes the performances of the subgridding framework (denoted as SGF in the figures) applied to the 2018–2019 and 2019–2020 HRDPS time series. The SGF delivers a mixed performance for TA in 2018–2019. The HRDPS model shows a 1.8 °C negative bias at Abbott, which is reduced by 1.2 °C when using the subgridding framework. However, the negative bias increases by 0.6 °C at the Fidelity station and increases very slightly at the Hermit station (< 0.1 °C). On the other hand, the mean absolute error (MAE) and the Spearman R coefficient slightly increase for most of the validation stations. TA shows a very strong correlation with station measures (R>0.8).

Relative humidity subgridding yields good performances as the HRDPS model bias and MAE are reduced by 1 % to 6 % at all sites. The bias is constant at Fidelity station, as is the MAE at Abbott. The Spearman correlation coefficient also slightly increases at Hermit and Fidelity stations and only slightly decreases at Abbott. Finally, modeled RH shows a strong correlation with station values ($\mathrm{0.6}<R<\mathrm{0.8}$).

Overall, the SGF performs best at subgridding VW. HRDPS is constantly overestimating wind speed by 1.5 cm s⁻¹. This bias is considerably reduced to 0.25 m s⁻¹ on average, and MAE is reduced by around 1 m s⁻¹. However, the correlation with station values is overall weak ($\mathrm{0.2}<R<\mathrm{0.4}$), and the subgridding workflow seems to have even weakened this correlation, except for Hermit station, which shows a negligible correlation (R<0.2).

Finally, subgridding shows a good performance for PSUM as well. In agreement with the literature, the bias in modeled precipitation shows a general lack of precipitation in the HRDPS model and ranged from 0.05 to 0.2 mm. For 2018–2019, MAE values range between 0.25 and 0.35 mm. With the subgridding, bias and MAE decrease at Abbott and slightly increase at Hermit, and bias decreases at Fidelity, while the MAE remains constant. Correlation with station values is moderate at Abbott and Hermit ($\mathrm{0.4}<R<\mathrm{0.6}$) and strong at Fidelity ($\mathrm{0.6}<R<\mathrm{0.8}$). Finally, subgridding slightly improves the correlation with AWS measurements at the former sites and stayed constant at the latter.

Figure 3 Atmospheric parameter evaluation for the seasons 2018–2019 and 2019–2020. The first row shows biases with respect to AWS measurements for each atmospheric parameter, the second row shows mean absolute error, and the third shows Spearman R coefficients. Blue and orange bars refer to values for 2018–2019 raw HRDPS and subgridding framework respectively. Green and red bars refer to values for the 2019–2020 season, in the same order.

The 2019–2020 season yields very similar results to those for 2018–2019; the bias and MAE are corrected on the same scale, and the Spearman R coefficient is in the same range for each variable. The only notable difference is that the PSUM bias at Abbott is positive for this season, meaning that HRDPS overestimated precipitation, which is highly unusual. As a result, the SGF introduces even more precipitation bias (+0.07 mm).

4.2 Subgridding performance for snow modeling

Figures 4 and 5 summarize the subgridding framework performances for the seasons 2018–2019 and 2019–2020. From here, snow simulations are denoted as SGF-SNOWPACK, HRDPS-SNOWPACK, and AWS-SNOWPACK when driven respectively by subgridded atmospheric parameters from the SGF, raw HRDPS forecasts, and AWS measurements. For 2018–2019, the snowpack similarity behaves identically at every site and for both HRDPS-SNOWPACK and SGF-SNOWPACK. The season begins with average similarity values (around 0.5), and then it plummets to low values in mid-October (<0.5) before improving to higher levels of similarity in November and for the rest of the season (0.6 < sim < 0.8). In general, HRDPS-SNOWPACK tends to have a closer similarity to AWS-SNOWPACK early in the season. SGF-SNOWPACK tends to score higher similarities than HRDPS-SNOWPACK in the mid-season before converging back with HRDPS-SNOWPACK in the spring. For the 2019–2020 season, similarity is again highly variable around 0.5 at Abbott and Hermit for both HRDPS-SNOWPACK and SGF-SNOWPACK. The early season at Fidelity shows very low similarities for SGF-SNOWPACK and HRDPS-SNOWPACK. Then, starting in November, the similarities stabilize and slowly rise throughout the season to 0.8 at every site. Again, SGF-SNOWPACK shows a higher similarity than HRDPS-SNOWPACK during the mid-winter period at Abbott and Hermit. HRDPS-SNOWPACK reaches the same level of similarity by early spring or the end of the winter. At Abbott, HRDPS-SNOWPACK and SGF-SNOWPACK similarities are very close, though the former shows more fluctuations during most of the winter and early spring. The average similarity for all sites and seasons is 0.8 for SGF-SNOWPACK and 0.75 for the HRDPS-SNOWPACK. Overall, SGF-SNOWPACK increased snowpack similarity by 7 % at Abbott, by 2 % at Hermit, and by 6 % at Fidelity with respect to HRDPS-SNOWPACK and when compared against AWS-SNOWPACK. The mean error in OGS shows a similar behavior at every site in 2018–2019 for both approaches. The error peaks and fluctuates strongly in the early season for both approaches and then stabilizes around mid-November. Considering the high variability in the fall for both seasons (September to November included), the first 3 months were considered a spin-up phase for the model to initiate a proper snowpack. Hence, the numerical analysis of the results was carried out starting on the first of December. Both HRDPS-SNOWPACK and SGF-SNOWPACK are slightly overestimating OGS with respect to AWS-SNOWPACK, but overall, the proposed method allows the bias to be decreased by 0.04 mm at Fidelity. However, OGS bias remains constant at Abbott and Hermit in the same period. The same pattern repeats for the 2019–2020 season at all sites, where SGF-SNOWPACK decreases on average the OGS bias by 0.07 mm at Fidelity and 0.09 mm at Hermit. The bias in OGS remains unchanged at Abbott once again that season.

Figure 4 Snowpack similarity assessment for the season 2018–2019. The first row shows DTW similarity for each station, the second shows bias in OGS with respect to AWS-SNOWPACK, the third shows density bias with respect to AWS-SNOWPACK, and the last row shows simulated HS (and measured when available). The blue curves refer to SGF-SNOWPACK, the orange curves to HRDPS-SNOWPACK, and the green curves to AWS-SNOWPACK. In the HS plots for Abbott and Fidelity, the red curves refer to SR50 measurements at the station plot.

Figure 5 Snowpack similarity assessment for the season 2019–2020. The first row shows DTW similarity for each station, the second shows bias in OGS with respect to AWS-SNOWPACK, the third shows density bias with respect to AWS-SNOWPACK, and the last row shows simulated HS (and measured when available). The blue curves refer to SGF-SNOWPACK, the orange curves to HRDPS-SNOWPACK, and the green curves to AWS-SNOWPACK. In the HS plots for Abbott and Fidelity, the red curves refer to SR50 measurements at the station plot.

Similarly to OGS, in 2018–2019, the mean error in density shows higher values in the early season and stronger fluctuations for both HRDPS-SNOWPACK and SGF-SNOWPACK. Variations tend to stabilize by mid-November, and the error increases again towards the end of the season. Again, the numerical analysis was performed from the first of December until the end of the simulation for both seasons. Generally, HRDPS-SNOWPACK seems to underestimate snow density. SGF-SNOWPACK brought density 3.17 kg m⁻³ closer to AWS-SNOWPACK at Abbot and 6.79 kg m⁻³ closer at Fidelity for the 2018–2019 season. However, in the same period at Hermit, the density bias increases by 1.71 kg m⁻³ on average. In 2019–2020, the density bias decreases on average by 4.65 kg m⁻³ at Fidelity and by 3.62 kg m⁻³ at Hermit. However, the density bias increases by 0.9 kg m⁻³ over the same period at Abbott.

Finally, SGF-SNOWPACK mean error in modeled HS is 35 cm in 2018–2019 (20 cm improvement when compared to HRDPS-SNOWPACK) and 29 cm in 2019–2020 (29 cm improvement) when compared with the SR50 measurement at Fidelity. For reference, the modeled HS with AWS-SNOWPACK shows a mean error of 8 cm in 2018–2019 and of 14 cm in 2019–2020. However, SGF-SNOWPACK seems to degrade the quality of the HS modeling regarding HRDPS-SNOWPACK when compared to the SR50 at Abbot station, overestimating HS for each season. Yet, the AWS-SNOWPACK run at Abbott shows a high discrepancy with the SR50 measurements as well, overestimating HS (especially in 2018–2019). Finally, HS remains relatively unchanged at Hermit for 2018–2019 as the framework does not bring a substantial improvement when compared to the AWS-SNOWPACK-modeled HS.

SWE modeling is considerably improved at all stations except at Abbott in 2019–2020 (Fig. 6) when compared to AWS-SNOWPACK. Table 3 summarizes the Nash–Sutcliffe efficiency coefficient (NSE) values at each site and for each season. On average, SGF-SNOWPACK improves the SWE NSE by 13 %, up to 57 % at Hermit in 2019–2020.

Figure 6 SWE modeling for the seasons 2018–2019 and 2019–2020. The first row represents SWE simulations at each site for the 2018–2019 season, and the second row represents SWE simulations for the 2019–2020 season. The blue curves refer to SGF-SNOWPACK, the orange curves refer to HRDPS-SNOWPACK, and the green curves refer to AWS-SNOWPACK.

Table 3 Nash–Sutcliffe model efficiency coefficient for SWE at each site for each season. Values in bold denote the highest scoring method at each suite and for each modeled season.

4.3 Intra-cell spatial variability provided by the subgridding framework: a case study

Figure 7 summarizes the main atmospheric parameter values for the six spatial variability sites, averaged per month (and accumulated for precipitation). First, as a direct consequence of the applied lapse rate for TA downscaling and spatialization, the general rule of thumb that TA should get colder with elevation is respected. However, the selected point on the north aspect is 100 m lower than the south aspect point, and as a result, the figure shows slightly warmer temperatures consistently throughout the season on the north aspect. Second, the aspect gradient is respected with lower incoming shortwave radiation and slightly lower temperatures in the north aspects. Moreover, wind direction is mostly coming from the south and southwest. South slopes are thus more exposed to the wind and north aspects are more sheltered. This is reflected by wind speed values being higher in the south aspects, especially in the alpine region. Wind-generated snow redistribution is accounted for, leeward slopes getting more snow than windward slopes. Finally, the altitudinal precipitation rate gradient is also respected by the subgridding framework, with precipitation rates getting higher with elevation. This atmospheric parameter variability is then propagated to the modeled snowpacks (Fig. 8). Again, the altitudinal gradient in HS is present, with a deeper snowpack from below the treeline to the alpine region. Moreover, the melt onset date is a few days earlier on the south aspect than on the north aspect, and water is percolating faster and deeper in the snowpack on the south aspects.

Figure 7 Intra-cell spatial variability in subgridded atmospheric parameters for the season 2018–2019. Each row represents monthly averages for each atmospheric variable and cumulative precipitation for each elevation band. Blue bars correspond to the north aspect, and orange bars correspond to the south aspect.

Season 2019–2020 results for spatial variability appear in Appendix A, Figs. A1 and A2. The subgridding framework creates the same altitudinal and aspect gradients, with more precipitation overall, milder temperatures, and stronger wind speeds on average.

Figure 8 Snowpack spatial variability assessment for the season 2018–2019. The left column represents snow profiles on the south aspect in the alpine, treeline, and below-treeline elevation bands. The right column represents snow profiles at the same elevation bands but on the north aspect. Colors represent grain type, which are denominated according to the international classification of Fierz et al. (2009).

4.4 Spatial variability in snow water equivalent provided by the framework over the simulation domain

Figure 9 shows the evolution of simulated SWE averaged by elevation band and aspect. The elevation gradient is well represented over all four quadrants. In the east and west aspects, which face dominant winds in the area, the high alpine elevation is showing SWE equal to or lower than that in the alpine elevation band. This reflects the effect of wind on ridges, modeled in the SGF as an alteration of the PSUM field with wind speed and terrain features. Figure 10 shows the variability in the modeled SWE within each HRDPS cell over the simulation domain. The top plot shows the variability in the early season (19 November 2018), and the bottom plot shows the variability at the end of the season when the snowpack is at its peak (4 March 2019). Each cell shows a wide spread of SWE values, which indicates that the subgridding of weather parameters is very effective in bringing spatial variability within each cell.

Figure 9 Average SWE aggregated by elevation band and aspect for the season 2018–2019. Atmospheric parameter spatial variability assessment for the season 2018–2019. Red lines correspond to below-treeline elevations (< 1850 m), green is treeline (1850 m < elevation < 2000 m), orange is alpine (2000 m < elevation < 2900 m), and blue is high alpine (> 2900 m).

Figure 10 Boxplot of the SWE modeled by the subgridding framework within each HRDPS cell in the early season and at the end of the season. The labels on the x axis correspond to each HRDPS cell ID. The box spans the interquartile range (IQR), the line represents the median, and the whiskers extend to the minimum and maximum value within 1.5 times the IQR. Outliers have been removed.

5 Discussion

The proposed downscaling framework brings a noteworthy improvement to modeled atmospheric parameters when compared to station values. Bias and MAE for VW and RH are constantly attenuated. For precipitation, in most cases, bias is corrected but MAE increases. Indeed, although the model, on average, underestimates precipitation, major precipitation events are overestimated (Côté et al., 2017). The “one-directional” lapse-rate correction reduces the overall bias by accurately correcting the small and common underestimation errors but accentuates the overall larger and rarer overestimation errors, thus increasing the MAE. We acknowledge that it limits our work, and our lab is currently carrying out research to produce an adaptive bias-correction algorithm. Regarding TA, the mean systematic bias slightly worsens at every station except at Abbott, where bias clearly improves. However, the added error remains under 0.5 °C. MAE reacts similarly, with an exception for Hermit, where the MAE slightly lowers (i.e., improved precision). However, TA presents a positive bias at these three stations (i.e., HRDPS is colder than station measurements), which are all higher in elevation than the nominal elevation of their corresponding HRDPS cell. As a result, a naive inverse distance weighting lapse-rate spatialization scheme would have introduced even more bias and MAE in the system, aggravating the original error. The logarithmic bias correction reduces the bias in TA. However the IDW spatialization depends on the elevation difference between the HRDPS cell and each sub-pixel; IDW surely aggravates a positive TA bias at the parameterization stage if the sub-pixel is higher in elevation. We then argue that the logarithmic lapse-rate parameterization allows the error in the subgridding scheme to be reduced and keeps it within a physically meaningful interval regarding spatial resolution.

Before analyzing the accuracy of the subgridding framework for modeling snow properties, the use of a SNOWPACK run driven by AWS measurements as a validation tool for the subgridding framework needs to be discussed. Madore et al. (2022) performed a detailed parameterization and analysis of the SNOWPACK model at Fidelity station. Results show that SNOWPACK performs very well on modeling HS and SWE and slightly underestimates both parameters (−4.2 % mean error in bulk density). When comparing layer by layer with observed snow profiles, SNOWPACK models the density accurately with a slope of 0.88 and a correlation of 0.97, though the model overestimates low densities and underestimates high densities. Madore et al. (2018) showed that SNOWPACK usually overestimates OGS and introduces more variability than typical on-site observations. However, this lack of variability could be an artifact of the resolution of the sampler used in the field (5 cm). The accuracy of SNOWPACK at Fidelity has been proven, and its biases are clear. Therefore, it seems more than reasonable to consider a SNOWPACK run driven by AWS as a ground truth at Fidelity. The quality of a SNOWPACK run relies almost entirely on its input quality. Abbott and Hermit stations are remote, harder to reach, and hence more difficult to maintain. It was impossible to conduct a similar study at these sites, and we can safely assume that the AWS SNOWPACK runs present substantial bias regarding reality. The results on snow height modeling illustrate this particularly well. When compared against the station SR50 measurements, SGF-SNOWPACK drastically improves HS modeling at Fidelity over HRDPS-SNOWPACK for both seasons (especially in 2019–2020 where SGF-SNOWPACK is almost identical to AWS-SNOWPACK and very close to the SR50 measurements). However, at Abbott, AWS-SNOWPACK is consistently overestimating HS with regards to the SR50 measurements. The SGF is getting the SNOWPACK simulation closer to the AWS-SNOWPACK run but actually further away from reality. Hence, although the SGF improves weather parameters compared to the measurements at the AWSs, the weather inputs and the model parameterization fail to represent reality accurately at this site. Weather measurements at Abbott thus become dubious. Moreover, the SNOWPACK parameterization applied on every cell of the domain has been developed at Fidelity. The site is perfect for studying snow processes as it is sheltered from the wind, but it is not representative of the processes affecting the snow on the whole domain.

Simulated snow parameters are usually improved, depending on the site and period. The most important errors and lowest similarities occur early in the season when the snowpack is starting to build. This is due to differences in snow onset timing between simulations and reality, as well as milder air temperatures early in the season. As a result, a fine margin exists between solid (snow precipitation) and liquid (rain-on-snow), which can strongly alter the microstructure of the snowpack. During this period, the snowpack is also thinner, and discrepancies between layers have a heavier impact on the mean similarity and mean errors. With regards to the microstructure, the subgridding framework allows the OGS overestimation to be decreased by 0.06 mm compared with HRDPS-SNOWPACK simulations at Fidelity and by 0.07 mm at all sites and seasons. Optical diameter usually ranges from 0.1 to 0.4 mm in this region, and the average snow pit OGS is 0.38 mm. Hence, the subgridding framework could improve the modeled OGS by 18 % on average compared to AWS-SNOWPACK-modeled OGS. SWE NSE is improved by 22.5 % at Fidelity over all seasons and by 16 % when averaged over all sites and seasons with the AWS-SNOWPACK as a reference. The SWE bias is improved by 86.5 mm (both season average) at Fidelity and by 56.7 mm on average over all sites and seasons. This represents a major improvement in a remote sensing perspective as SWE is the major driver for SWE inversion algorithms (King et al., 2015; Zhu et al., 2018; Tsang et al., 2022). For instance, King et al. (2015) reported a backscatter increase of 0.82 dB per 1 cm increase in SWE at the Ku-band with Canadian tundra measurements. Similarly, Yueh et al. (2009) found a 0.15 to 0.5 dB increase for every 10 mm increase in SWE at the Ku-band in Colorado. According to these figures, using the subgridding framework in a SWE inversion context could provide an improvement of 1.4 to 7 dB in the Ku-band simulated backscatter.

Finally, the subgridding framework performs well in introducing spatial variability over the simulation domain. SWE distribution across topographic categories respects the elevation gradient, the orientation of dominant winds in the area, and the erosion effect on ridges. Spatial variability is key when considering SAR signal inversion (King et al., 2018), and the subgridding framework should be a highly relevant tool in this context. However, the introduced spatial variability has not been evaluated against distributed snow height measurements, but airborne snow height surveys are expensive and logistically challenging.

6 Conclusions

This study aims at meeting the present need in the snow remote sensing community for both the design and the evaluation of an atmospheric model subgridding framework to perform snow modeling in the context of coupling with a SAR signal inversion routine. To do so, (i) a new NWP downscaling approach was introduced by first parameterizing the 2.5 km HRDPS cells into a virtual weather station array, which was then spatially interpolated using the MeteoIO and Alpine3D models. (ii) Snow simulations were performed using the state-of-the-art model SNOWPACK. Microstructure modeling quality was assessed using the DTW algorithm and an original cost function focusing on density and optical grain size, and SWE modeling improvements were quantified using the Nash–Sutcliffe efficiency coefficient. (iii) The spatial variability in atmospheric parameters and snowpack state variables within one subgridded HRDPS cell was assessed. Finally, we assessed the introduced spatial variability in SWE over the simulation domain, as well as intra-cell variability. The key findings regarding the research questions from the Introduction are as follows:

1. The atmospheric parameter subgridding framework yields an overall good performance, especially for RH and VW. Further research should be carried on to find a better spatialization algorithm for air temperature and an adaptative precipitation rate correction algorithm.

2. The general overestimation of OGS by SNOWPACK when driven by raw HRDPS data decreased by 0.06 mm on average at Fidelity and by 0.07 mm averaged over all sites. This represents an 18 % improvement over raw HRDPS-SNOWPACK-simulated OGS. The Nash–Sutcliffe efficiency coefficient for SWE was improved by 22.5 % at Fidelity and by 16 % on average when compared to HRDPS-SNOWPACK simulations. SWE bias was diminished by 86.5 mm at Fidelity and by 56.7 mm on average. This is a major improvement in a SAR remote sensing context as this could lead to up to a 7 dB improvement in the Ku-band simulated backscatter. In this context, the first 3 months (September to November) of snow simulations should be considered a spin-up phase for the snow model, as discrepancies between reality and simulations are critical before the snowpack is properly established and the similarity stabilizes in December and onward.

3. The subgridding framework introduces a realistic spatial variability in snowpack state variables, respecting altitudinal and orientation gradients as well as ridge effects. The framework brings substantial variability within each HRDPS, reflecting the high spatial variability in the snowpack at the kilometer range in a mountainous environment.

This study shows that downscaling NWP at a 100 m resolution can improve local representation of atmospheric variables and, as a result, improve the modeling of snowpack state variables and spatial variability in the snowpack in complex topography. The modeling of the two key parameters for snow remote sensing, SWE and OGS, was improved. This work is highly relevant in a remote sensing context. The remote sensing community is currently pushing for new SAR satellite missions, and observation system synthetic experiments have proven that satellite SWE measurements would substantially improve SWE products' RMSE (Garnaud et al., 2019; Cho et al., 2023). This study provides finer-spatial-variability forcing data and improved simulated snowpack state variables regarding NWP-driven simulations. As a result, snow simulations performed with such a method can provide a realistic first-guess estimate of the retrieval parameters, bringing a solid basis to overcome the non-unique solution issue in physical retrieval algorithms (Tsang et al., 2022) and steer away from empirical retrieval approaches. In this context, the next logical step is to design a SWE retrieval algorithm, exploiting the vast array of SAR satellites in orbit, such as Sentinel-1 (C-band), TerraSAR-X (X-band), or the SnowSAR mission concept (dual Ku-band) led by Environment and Climate Change Canada and the Canadian Space Agency (Derksen et al., 2021). Future work can thus focus on using these modeled snowpack state variables along with field-inferred distributions as a basis for a SWE SAR retrieval algorithm.

Appendix A: 2019–2020 spatial variability figures

Figure A1 Intra-cell spatial variability in subgridded atmospheric parameters for the season 2019–2020. Each row represents monthly averages for each atmospheric variable and cumulative precipitations for each elevation band. Blue bars correspond to the north aspect, and orange bars correspond to the south aspect.

Figure A2 Snowpack spatial variability assessment for the season 2019–2020. The left column represents snow profiles on the south aspect in the alpine, treeline, and below-treeline elevation bands. The right column represents snow profiles at the same elevation bands but on the north aspect. Colors represent grain type, which are denominated according to the international classification of Fierz et al. (2009).