2.1 Surface model

The surface model used in this study, GEOtop, is a physically based model originally developed for hydrological research (Endrizzi et al., 2014). It couples energy and water budgets, represents the energy exchange with the atmosphere and has a multilayer snowpack. Further information is given by Bertoldi et al. (2006), Rigon et al. (2006), Endrizzi (2007) and Dall'Amico et al. (2011). A description of model uncertainty and sensitivity is given by Gubler et al. (2012). Model parameters and soil stratigraphy are set up as defined in Fiddes et al. (2015).

2.2 Downscaling forcing

TopoSCALE is a scheme that generates point-scale model forcing using gridded atmospheric model datasets. It achieves this as follows: (1) data available on pressure levels, including air temperature (T_a), relative humidity (RH), wind speed (U) and wind direction (φ_U) are interpolated to a point of interest in order to provide a dynamic scaling at each time step; (2) incoming longwave radiation (L^↓) is scaled by accounting for downscaled T_a, RH and sky emissivity; (3) we apply a topographic correction to both radiation fields ($S^{\downarrow }/L^{\downarrow }$); and (4) an elevation based lapse rate is applied to precipitation (P). The output is a full set of scaled meteorological fields required to drive a numerical model at hourly time steps.

2.3 Sub-grid scheme

TopoSUB is a scheme that samples land surface heterogeneity at high resolution based on a digital elevation model (DEM) and other surface data (here SRTM 3, 30 m). Input predictors describing dimensions of variability are clustered with a k-means algorithm to reduce computational units in a given simulation domain to a set of clusters. A 1-D surface model is then applied to each cluster using its mean physiographic properties. This approach allows savings in computational effort of multiple orders of magnitude over distributed approaches. For example, a simulation domain represented by an ERA5 grid cell (31 km × 31 km) contains approximately 10⁶ SRTM 3 pixels. This domain can be simulated using 100 TopoSUB clusters, which represents a 10⁴ reduction in the computational load during simulation.

2.4 Data assimilation

We build on previous efforts (e.g. Girotto et al., 2014; Margulis et al., 2015; Aalstad et al., 2018) that focus on the reanalysis of snowpack characteristics (particularly SWE and height of snow, HS) through ensemble-based assimilation of fractional snow covered area (fSCA) retrievals from optical satellite sensors. We choose to use fSCA retrievals because currently only optical satellite sensors can offer the resolution, coverage, accuracy and breadth of information needed to constrain snowpack simulations in complex terrain (see Dozier et al., 2016). We use fSCA retrieved from the MODIS sensors onboard the Aqua and Terra satellites. These retrievals have a sub-kilometric spatial resolution and a near-daily equatorial revisit frequency (in the absence of clouds); thus, the reanalysis we perform could be applied to any mountain range on Earth. By assimilating fSCA observations we exploit the dynamic information content contained in the depletion of the fractional snow cover. The idea is that if one grid cell melts out later than another, there must either have been more snow there to begin with, a slower ablation, or a combination of the two – and vice-versa for an earlier melt out (Martinec and Rango, 1981; Aalstad et al., 2018). This is the essence of traditional snow reconstruction where the snowpack is contrived in reverse from the observed date of disappearance of the snow cover to the day of peak SWE using modelled snowmelt rates (Martinec and Rango, 1981; Dozier et al., 2016). By using ensemble-based DA we can account for uncertainties in the remotely sensed fSCA depletion, the meteorological forcing and the snow model that are ignored in traditional reconstruction (Slater et al., 2013) and arrive at an improved reanalysis (Girotto et al., 2014). Snow reanalysis problems are best approached using batch smoother DA algorithms rather than the more commonly used filters as the snowpack has a long memory (i.e. high temporal autocorrelation) relative to (e.g.) synoptic-scale weather (Margulis et al., 2015; Aalstad et al., 2018). By using a smoother that assimilates all of the fSCA retrievals during the ablation season at once to constrain the ensemble of annual snowpack trajectories, we are able to use the observed ablation to inform the accumulation season which would not be possible with a particle filter.

2.4.1 Generating the prior ensemble

In line with previous studies (e.g. Raleigh et al., 2015), we assume that the main source of uncertainty in modelling the snowpack is in the meteorological forcing and specifically the main variables that control the mass and energy balance, namely air temperature (T_a), precipitation (P), and incoming shortwave (S^↓) and longwave (L^↓) radiation. To generate the prior ensemble we perturb the forcing time series using normally (T_a, S^↓ and L^↓) and lognormally (P) distributed multiplicative perturbation parameters that are fixed throughout the annual integration. Following Navari et al. (2016) we generate a correlated ensemble of perturbation parameters for the different forcing variables. This is to avoid unrealistic perturbations such as a large increase in both precipitation and shortwave radiation. We do this in two steps. First, we generate independent perturbation parameters for each of the forcing variables using normal and lognormal random draws. Secondly, we account for the correlation between the different perturbation parameters by performing a Cholesky decomposition of the covariance matrix. All hyper-parameters used in generating the prior ensemble are given in Table 1 and are based on values from a study in Colorado by De Lannoy et al. (2012) which, in turn, are based on the approach of Reichle et al. (2007), which is a global study.

Table 1Hyper-parameters (means, variances and correlations) defining the joint probability distribution from which the ensemble of multiplicative perturbation parameters are drawn (unitless). These parameters were obtained from De Lannoy et al. (2012) which, in turn, were based on the approach of Reichle et al. (2007).

Download Print Version | Download XLSX

3.1 Meteorological forcing

Driving climate data are obtained from the ERA5 reanalysis from ECMWF. This is the latest reanalysis from ECMWF that updates the ERA-Interim reanalysis. The main improvements are an increase in the spatial resolution to 31 km, hourly temporal resolution and an increase in the vertical model levels to 137. Accumulated values are now from the last time step and not from the last forecast as in ERA-Interim. This means that we can easily obtain the mean rates required to drive our numerical model by simply dividing these accumulations by the hourly time step (Fiddes and Gruber, 2014). Forcing data are detailed in Table 2. For each TopoSUB cluster, defined by the mean physiographic characteristics of a cluster, (Fiddes and Gruber, 2012) the ERA5 meteorological fields are downscaled using TopoSCALE (Fiddes and Gruber, 2014).

Table 2Description of the hourly fields obtained from the ERA5 reanalysis. All of the columns headers are terms defined by ECMWF. “levtype” refers to the level type, “surf” refers to the surface and “pl” denotes the pressure level. The “type” is either “fc”, which refers to forecast, “an”, which refers to analysis, or “inv”, which refers to invariant.

Download Print Version | Download XLSX

3.2 Surface properties

TopoSUB requires topographical parameters as input predictors to the clustering algorithm. We derive the following topographic parameters from the SRTM 3 digital elevation model: elevation, slope, aspect and sky view factor (proportion of visible sky). Surface cover is characterized in a simple three-mode classification in order to approximate subsurface stratigraphies: first a threshold on MODIS normalized difference snow index (NDSI) is used to classify vegetated surfaces, then a simple model further differentiates between steep bedrock and debris slopes. Further details are available in Fiddes et al. (2015).

3.3 Assimilated fSCA observations

We assimilate fSCA retrievals obtained from version 6 of the level 3 daily MODIS snow cover product from the Terra (MOD10A1 product; Hall and Riggs, 2016a) and Aqua (MYD10A1 product; Hall and Riggs, 2016b) satellites. The retrieval algorithm is based on the inversion of a linear regression of the MODIS NDSI on reference fSCA estimated from coincident Landsat imagery and it is given by the “FRA6T” relationship in Salomonson and Appel (2006). The normalized difference snow index exploits the fact that snow is highly reflective in the visible but a good absorber in the shortwave infrared which differentiates it from most other natural surfaces (Painter et al., 2009). If cloud-free retrievals are available from both Terra and Aqua retrievals for a given day then the Terra retrievals are used. Aalstad et al. (2018) compared MODIS fSCA retrievals to reference fSCA estimates obtained from time-lapse photography, imagery from an unmanned aerial vehicle and snow surveys at a site on Svalbard and obtained an RMSE of σ_y=0.13 for the MODIS retrievals. This estimate is in good agreement with those found in the Alps by other studies (e.g. Masson et al., 2018); thus, we use this as the observation error variance (${\mathit{\sigma }}_y^{\mathrm{2}}$) in the assimilation (Sect. 2.3.2).

3.4 Evaluation

3.4.1 Station data

SWE (mm) measurements obtained manually by observers are available at approximate biweekly intervals from snow profiles across Switzerland. Here we use the GCOS dataset which consists of 11 sites (Fig. 2). We call these sites “stations” throughout the paper. The dataset is openly available (Marty, 2017). Automatic HS (cm) measurements performed by sonic ranger (Campbell Scientific SR50) are available from the Intercantonal Measurement and Information System (IMIS) station network at 30 min intervals. This is a high-elevation station network that forms the backbone of the national avalanche service in Switzerland.

Figure 2Experimental set-up: the nine ERA5 grid boxes were selected based on the fact that they contained GCOS SWE monitoring sites (11 stations). All of the IMIS stations in each box are used for evaluation (39 stations). The Weissfluhjoch research station and the flight path of the ADS data are located in the red box, which is also shown at a larger scale in the inset.

4.1 Point DA

Point-scale DA is accomplished by simply mapping the DEM cell corresponding to point of interest (e.g. a validation station) to the corresponding MODIS pixel for that location. The W_p derived from the PBS analysis for that MODIS pixel is then used directly to generate the posterior estimate for that point. Model state results are obtained from the TopoSUB cluster that the DEM cell is a member of. Cumulative distributions of state variables are computed through the ranking of the ensemble of state variables followed by a cumulative summation of the correspondingly sorted weights. These distributions allow for the estimation of quantile values of the posterior model state.

4.2 Spatially distributed DA

The particle batch smoother has typically been applied at the point scale or on regular grids. Here we generalize the method so as to fit the TopoSUB approach. The basic aim is to generate posterior weights for each TopoSUB cluster (which has no location, only physiographic attributes) so that the weights can be used in a highly scalable and powerful manner to generate different products such as a time series of the posterior median aggregated to give basin-level statistics. The key challenge in this aim is how to map the spatial unit of the PBS algorithm, the MODIS pixel, which has a location in space, to a TopoSUB cluster which does not. We achieve this via the following:

where W_p is the N_e×N_p weights matrix and a is a N_p×1 vector containing the fractional abundance (cover) of cluster c represented in each of the MODIS pixels. Here w_c contains the weight of each forcing history for cluster c and is computed for each of the clusters. This yields the weights matrix w_c that contains the weights of the forcing histories for each cluster. As a second step the weights w_c are renormalized to sum to one as that is not guaranteed in Eq. (3).

4.3 Coarse-grid DA

The third DA method addresses bias in forcing at the grid level only, and it is the most efficient and lightweight of the three approaches. It also differs from the previous methods in that the PBS analysis step is computed at the ERA5 grid unit not at the MODIS pixel unit. This makes the analysis step highly efficient and scalable over large areas. It is emphasized that while the two previous methods address both aggregated bias in forcing at the grid level and also correct errors in the sub-grid method (such as physiographic description) and downscaling (such as precipitation distribution), this method only corrects the bias at the grid level. However, it is of interest if we seek a simple and robust way to feedback sub-grid information to large-scale atmospheric grid cells, in this case using ERA5. The method is carried out as follows:

1. MODIS fSCA is computed aggregated to the large-scale atmospheric grid cell (ERA5) while accounting for clouds (max 10 % cloudiness tolerated), and cloud pixels are then filled with the mean fSCA value of the cluster to which the pixel belongs;

2. the predicted observations are computed, i.e. the modelled fSCA, for each cluster and are aggregated to the ERA5 grid cell scale by multiplying by cluster members;

3. PBS is run at the ERA5 grid level to generate a single weight vector for the ensemble.

4.4 Run configurations

All runs are performed using 100 particles, 150 TopoSUB clusters and cover the period from 1 September 2011 to 1 September 2017. The specific temporal period covered by a given result is defined in the text. Throughout the paper a single year refers to the year in which melt occurs, e.g., “2012” refers to the period from 1 September 2011 to 1 September 2012. These “water years” are prefixed with WY (e.g. WY2012) to avoid ambiguity. We measure computational effort through the number of GEOtop model runs (N_r) required per year in an ERA5 grid cell. Recall that the ERA5 grid cell is the fundamental unit on which the downscaling and clustering is performed. In terms of the number of clusters (N_s) and particles (N_e), this effort becomes

In the case of the configuration used in this study (N_e=100, N_s=150), this amounts to 1.5×10⁴ individual model runs. At a 30 m resolution, there are 10⁶ model grid cells within a single 0.25^∘ ERA5 grid cell. Thus, an explicit fully distributed simulation with 100 particles would require N_r=10⁸, which is a 4 order of magnitude increase in computational effort relative to the set-up used in this study.

5.1 Evaluating the forward model

Figure 3 shows performance of the forward model at the Weissfluhjoch (WFJ) research site (see Fig. 2 for position) assessed over the period from WY2012 to WY2017. It illustrates the performance of the downscaling routine with respect to providing an adequate forcing to the model (forcing bias) and the performance of the model in simulating the target variables, SWE and HS, when driven by downscaled ERA5 reanalysis (revealing model and forcing errors) and station observations (revealing model and observation errors). It shows that the TopoSCALE downscaling routine does a reasonable job of providing forcing to the forward model (Fig. 3a–e) with a 0.71 ^∘C RMSE for 2 m air temperatures being particularly low. Conversely, high wind velocities tend to be positively biased, most likely as wind fields representing the free atmosphere on pressure levels have no surface drag that would be present in surface observations. Modelled HS and SWE (Fig. 3f–i) are captured fairly well, reproducing both the onset and melt of the snowpack. However, peak values are generally negatively biased with respect to observations and station-driven model runs. WY2012 is an obvious outlier with large snowfalls not captured by ERA5 precipitation. This can be seen by cumulative precipitation totals computed with and without WY2012 totals (Fig. 3). This is reflected in simulated HS and SWE totals. The performance of the forward model can be analysed by driving the simulations with station measurements to remove most of the uncertainty associated with driving reanalysis data (but with residual observation errors). ERA5-driven simulations are comparable or even outperform station runs in WY2013 and WY2014.

Figure 3Multiyear simulations at the WFJ station (WY2012–2017) in order to show baseline results for the modelling scheme. Panels (a)–(e) assess the downscaling scheme by showing downscaled ERA5 data (ERA5) compared to station measurements (Obs). Panels (f)–(i) assess the simulation of target variables, SWE and HS, in both time series and scatter plots. Here, ERA5 is a simulation driven either by downscaled ERA5 (ERA5) or directly by station measurements (Station). Obs are SWE and HS measurements made at the station. WY2012 is a clear outlier in poor performing ERA5 as shown by cumulative precipitation errors and in the HS and SWE time series. The HS and SWE scatter plots also show this low performance in high values attributed to WY2012. Additionally, ERA5-simulated HS is increasingly biased with depth as errors accumulate over the season to maximum depth values. The same pattern is evident for SWE. It is worth noting that these plots are useful in differentiating sources of error. Obs minus Station approximates model error, whereas Station minus ERA5 approximates the forcing error.

Download

Figure 4DA run at the Truebsee GCOS station, Engelberg. Panel (a) shows the assimilation step with the assimilated fSCA observations represented by black crosses. The shading and solid lines show the 90th percentile range and median of the prior (red) and posterior (blue) estimates. Panel (b) shows the target variable validation (SWE in this case). Posterior/prior are denoted in the same way. Black triangles indicate the measurements used for the validation.

Download

Figure 5Simulated snow depth at the IMIS stations (HS) and snow water equivalent at the GCOS stations (SWE) for both the prior (red) and posterior median (blue) compared to the observation mean. The mean is computed from all values over the entire WY2016. The posterior estimate is markedly improved in both variables. Regression lines compare the fit of the posterior and prior estimates with respect to observations against the 1 : 1 line.

Download

5.2 Point DA

In this experiment we compare the prior and (single-pixel) posterior HS and SWE for WY2016 to the measured values at the respective stations. An example at Truebsee GCOS station (Engelberg) is shown in Fig. 4. This figure demonstrates the effect that the assimilation has not only on the fSCA (which is assimilated), but also on estimates of the other state variables (in this case SWE) which get closer to independent observations. Here one can clearly see how the posterior estimate of SWE (blue shading) is constrained by the assimilation, and the posterior median (blue line) is much closer to validation SWE observations than the prior (red line). We then scaled this up to nine ERA5 grid boxes that span the Swiss Alps and contain 11 GCOS SWE stations. Additionally each box contains multiple IMIS stations measuring HS, which we also looked at in the interest of obtaining more validation data. Significant improvements in the posterior were seen in the estimation of both variables (Fig. 5). We found that the improvement in the SWE was greater than that in the HS. We hypothesize that this is due to representation of snow densities in the snow model. However, the improved representation of the snowpack mass balance as shown by improved SWE estimates is the main variable of interest in our approach. Stations where DA performed worse than the prior can be attributed to poorly characterized melt seasons, lack of MODIS retrievals and/or the MODIS retrievals not being representative at the scale of the observations (cf. Sect. 6).

Figure 6(a–c) Prior/posterior median and observations from ADS sensor flights in the Davos region on 14 April 2014 (see Fig. 2 for location). (d–f) Uncertainty represented by the 90th percentile range of the ensemble and the reduction in the uncertainty in the posterior. The glacier mask is shown in blue. The posterior median is improved with respect to the observations, and the uncertainty is reduced by the DA scheme.

Download

5.3 Spatially distributed DA

We evaluated the performance of the method in improving the spatial patterns and absolute quantities over large areas using data from an airborne digital sensor which has been used to generate high-resolution surfaces of HS in WY2012, WY2013 and WY2014 (Fig. 6, WY2014 only). Both WY2012 and WY2014 show marked improvement in all spatial statistical measures including the mean value, the standard deviation (indicating increased variation) and error statistics such as the RMSE and bias. WY2013 shows little improvement. We would expect a better performance for SWE than for HS due to the previously mentioned issues with the modelled snow density (see Sect. 5.2 and Fig. 5). Figure 6 shows how the 90th percentile range is constrained by the analysis going from the prior to the posterior. Figure 7 shows probability density distributions for observations, prior and posterior in WY2014. The shape and moments, such as the mean, more closely match the observations in the posterior distribution. However, the method fails to capture the very highest accumulations in the distribution (>2.5 m), possibly due to the averaging effects of generalizing weights to TopoSUB clusters.

Figure 7Density distribution plots for HS Obs, prior and posterior within the ADS footprint for 14 April 2014 (see Fig. 6). The observed distribution is better captured by the posterior. Dashed lines give the respective mean values.

Download

5.3.1 Interannual validity of weights

We tested the ability of weights obtained in a given hydrological year to improve results in a different year. We did this by looking at statistics on the Dischma Basin via a cross-validation exercise where each year was forced with results from the 2 other years (WY2012, WY2013 and WY2014). Posteriors forced by the weights of other years improved performance over priors in all cases (Fig. 8). This suggests that the DA method here also works to correct errors that are consistent from year to year. This could be related to spatial patterns of melt, which is a consistent bias in the forcing or errors in the model itself. This is an interesting result that suggests that while this method is primarily a post-processing method, it could be used to improve now/forecasts by using previous year's weights. Additionally an analogue approach could be used to find the years of best fit to the current season in order to select weight sets (Kolberg and Gottschalk, 2010).

Figure 8Interannual validity of weights generated by the DA scheme. Modelled versus observed mean snow depth averaged over the entire ADS zone on ADS acquisition dates WY2012, WY2013 and WY2014 are shown. Posteriors (blue) are generated using the weights of the other 2 years and are compared to the prior (red). For example, posteriors of WY2012 are generated using weights of WY2013 and WY2014.

Download

5.3.2 Large-scale application: seasonal variability

December 2016 was an extremely snow-poor month and start to the winter season. Many ski areas throughout the Alps could not open until late January due to lack of snow. We compare this to December 2011 which was relatively snow rich with above average precipitation and average temperatures for the month (cf. https://www.meteoswiss.admin.ch/home/climate/swiss-climate-in-detail/monthly-and-annual-maps.html, last access: 9 November 2019). Specifically, we investigated how open-loop runs in two contrasting seasons compared to observed spatial patterns of fSCA from MODIS and SLF reports dated 22 December of the respective year (Fig. 9). The model was found to compare well with both spatial patterns of fSCA and SLF snow depth maps, which are operational products created by interpolating station data constrained by the advanced very-high-resolution radiometer (AVHRR)-observed snow extent. Both fSCA and HS show snow-free zones deep into Alpine valleys during December 2016, and the absence of snow over the northern regions and Jura Mountains. The red box indicates the domain of the DA runs for these seasons shown in Fig. 10.

Figure 9Mean December high-resolution (30 m) large area HS simulations (open loop) in two contrasting seasons (c, d) compared to the observed spatial patterns of fSCA from mean December MODIS retrievals (a, b). Modelled HS compares well to the spatial pattern of fSCA. December 2016 was an extremely dry start to the season with many ski resorts unable to open until late January. Both observed fSCA and modelled HS show snow reflecting this fact, with snow-free zones deep into Alpine valleys. The red box indicates the domain of the DA runs for theses seasons shown in Fig. 11. The glacier mask is shown using light blue.

Figure 10The mean HS for December 2016, showing the (a) open loop compared to (b) posterior. (c) The difference plot shows how DA has reduced the low-elevation snow height relatively more than the high-elevation snow height. Variability has been increased. Snow-free valley bottoms in panel (b) show an improved match to (d) MODIS Obs and that the increased spatial variability is realistic.

Download

Next we zoomed in and compared spatial patterns of HS between the deterministic open-loop and posterior run; this demonstrated that DA increased elevation gradients of variability by reducing HS in valley bottoms and increasing it on higher slopes. Thus, DA, as mentioned previously, has the effect of increasing variability in the snow cover distribution. Snow cover extent estimation is also improved by DA with increased snow-free areas in valley bottoms showing an improved fit to MODIS observations.

5.4 Coarse-scale DA

Aggregated series of observed and modelled fSCA are computed at the ERA grid level. This could equally be a hydrological unit such as a basin. The main idea is to correct grid-level biases in the forcing only. If we assume this is the main source of uncertainty, especially with a view to correcting large-scale biases, this is an effective method to apply at the scale of meteorological reanalysis. Data from WFJ are used to illustrate this point. Figure 11 shows two contrasting snow seasons – WY2012 (high) and WY2014 (low) – where mean snow depths and SWE differed by a factor of 2, as recorded at WFJ (Fig. 3). We compare the total ERA5 precipitation (PSUM) over the winter period January–April 2014 (401 mm) and compare it to totals recorded over the same period at the WFJ station (350 mm), ERA5 captures WFJ totals well. It should be added that there is some elevation difference between the ERA5 grid (2024 m a.s.l.) and the WFJ station (2560 m a.s.l.). However, in WY2012 we see quite a different story. The ERA5 grid gives us slightly higher PSUM values of 440 mm, whereas the measured PSUM was almost double this at 826 mm. Figure 11 shows how grid-level biases in the driving forcing from ERA5 have been successfully decreased in WY2012 resulting in increased SWE totals, whereas in WY2014 where ERA5 performance was much better (cf. Fig. 3), DA has had a negligible effect. This simple approach is an extremely cost-effective method of assimilating slope-scale sub-grid information (in this case fSCA) to correct coarse grid-scale forcings (ERA5). It is additionally generic enough that it could be used with various other sub-grid observations, such as soil moisture, to improve grid-level responses.

Figure 11Assimilation of fSCA at the grid level which targets bias in grid-level forcing. Two contrasting seasons, WY2012 (a, b) and WY2014 (c, d), are shown. The shading and solid lines show the 90th percentile range and median of the prior (red) and posterior (blue) estimates respectively. Assimilated observations are indicated in green. Grid-level biases in WY2012 are compensated for by DA. Grid-level forcing was much more accurate in WY2014 and the resulting effect of DA was negligible.

Download

6.1 Sources of error

6.1.4 Scale issues in assimilation

The scale difference between validation data (station or snow profile) and fSCA retrievals from MODIS creates several issues. In this study many of the sites are in valley bottoms so they are accessible on a regular basis. However, this creates a bias caused by how representative the point measurement of the larger MODIS pixel footprint is. In Alpine valley bottoms there tends to be a lot of infrastructure, housing and rivers, which will tend to be snow free earlier (or never snow covered) compared with a station site that is well protected from interference allowing the natural accumulation of snowfall. Therefore, the MODIS footprint will tend to be observed as “snow free” earlier than the validation point. In addition features such as a rivers, road clearance and urban heat islands are not considered in the modelling and will generate bias in the data assimilation. Thus, the most reliable sites for data assimilation, or actually we should say for validating the method, are at high elevations away from effects due to human activity or infrastructure that are not considered in the model. Figure 12 gives an example of a point-scale DA failure that is not due to the DA algorithm (this has worked well), but stems from the representativeness of the fSCA retrievals. The posterior has correctly been pulled in the direction of what has been detected to be the main melt period (red dots), erroneous snowfalls during late spring/early summer are ignored as expected and the end of the winter snowpack (as detected by the fSCA retrievals) has been correctly identified at the beginning of April. However, this site is in the middle of Zermatt town and the MODIS pixel will likely contain a signal from urban effects unaccounted for by the model.

Figure 12An example of poor performance due to the non-representative fSCA retrievals. The posterior has been correctly pulled back to the observed depletion curve. However, it is likely that the depletion curve does not represent the validation station well due to urban effects.

Download

6.2 Applications

Using the methods described in this paper, a range of processing pipelines can be built to address a wide array of both research and operational problems. Specific strengths of the approach are as follows:

• there is slope-scale forcing (climate, reanalysis, forecast) globally;

• the approach explicitly includes the effects of high-resolution topography on surface–atmosphere interactions;

• it is an efficient method to make large ensemble simulations feasible;

• there is data assimilation to correct bias in forcing and quantify uncertainty.

Perhaps most importantly this approach allows applications to be built in remote regions where dense observation networks do not exist, such as High Mountain Asia or parts of North America. These capabilities allow for operational applications such as large area mass movement assessments related to the dynamics of surface and subsurface processes. Driving the system with NWP (forecast) data would allow nowcasting/forecasting applications to be set up with a suitable assimilation framework such as the ensemble Kalman filter (EnKF). While the assimilation of fSCA would be less informative in a sequential method (such as the EnKF), ensemble simulations would still provide a useful quantification of uncertainty.

Transient climate change studies using a combination of reanalysis and climate model data (e.g. CMIP5) would be a valuable research application based on this approach, for example quantifying dynamics of permafrost extent over large areas according to a range of scenarios and models or generating a regional snowpack reanalysis product with projected future changes.

An important operational application and currently a great humanitarian need in many remote regions in Asia (e.g. Afghanistan/Tajikistan) could be an operational avalanche forecast based on a snowpack model (e.g. SNOWPACK, CROCUS), driven by an NWP ensemble to generate a large area probabilistic forecast where few ground stations exist. This would be a relatively cost-effective system to deploy and give first-order hazard assessment where none currently exists.

6.3 Further work

For the moderate (MODIS-like) resolution satellites, we hope that products will emerge from Sentinel-3 and VIIRS to prolong and expand the MODIS record. For high-resolution sensors, there is a strong need for operational products that ideally combine available and emerging sources such as Landsat 8 and Sentinel-2. The French inter-agency initiative THEIA Land Data Centre is starting to produce Sentinel-2-based snow cover products (Gascoin et al., 2019), which – at both high temporal and spatial resolution – will likely prove to be a great product for data assimilation work. Additionally, improved cloud masks are needed as misclassified clouds are potential and significant sources of error in the framework. Both overly strict and overly relaxed cloud masking is problematic, the former leads to valid and potentially important retrievals being discarded, whereas the latter corrupts the signal that we are trying to assimilate (the actual snow cover depletion).

For additional datasets (other than fSCA), land surface temperature can be retrieved from both MODIS and Landsat and provide a means to constrain uncertainty in the surface energy balance. However, the current MODIS products are coarse at 1 km and, therefore, not ideal for mountain regions.

For additional datasets (other than fSCA), land surface temperature (LST) can be retrieved from both MODIS and Landsat and provide a means to constrain uncertainty in the surface energy balance. However, the current MODIS LST products are coarse at 1 km with respect to the expected heterogeneity of LST in mountain regions (Gubler et al., 2011). Snowmelt status (i.e. binary melting/not melting) from synthetic aperture radar (SAR) e.g. Sentinel-1 has the potential to constrain uncertainty in the fSCA during cloudy periods. There is also potential from ICESat-2 which could provide a way to constrain snow depth directly (Treichler and Kääb, 2017).

Assimilation of sparse point data could be an important extension of this work to provide means to assimilate data sources such as ICESat-2 but could also be used to improve TopoSCALE by assimilating point data (stations) to improve the downscaling of reanalysis data. This could be interesting as TopoSCALE performs most poorly is in valleys, where surface effects are poorly represented by the atmospheric model, and this is precisely where stations tend to be most abundant globally. For real-time applications in remote regions, extending the method to assimilate sparse observations is important as fSCA is known to have limited value in sequential (i.e real-time) data assimilation. We have shown that there is some interannual validity of results in our limited test case, suggesting that systematic biases relevant to real-time applications could be addressed through reanalysis. Additionally, creating a long-term library of “best possible” reconstructions/reanalyses and then training an “analogue ensemble” or even a more machine-learning-type approach (like neural nets) could be promising.