2.1 Modelling geometry

Simulations are performed in semi-distributed geometry. Mountain ranges such as the Alps are discretised into so-called massifs of about 1000 km² to account for regional variability of meteorological conditions. Within each massif, topographically induced variability is taken into account by running the model for a fixed set of topographic classes, e.g. by 300 m elevation bands, for 0, 20 and 40^∘ slopes and eight aspects (see Fig. 1). This set enables us to reproduce the main features of snowpack variability (e.g. Mary et al., 2013).

Figure 13-D schematic view of the semi-distributed geometry, for which the numbers represent the altitudes of the elevation bands (m). From left to right, the three different mountains represent the flat, 20 and 40^∘ slopes.

Download

In this study, we focus on the Grandes Rousses, a single massif in the central French Alps. This area of about 500 km² is represented by N_pts=187 independent topographic classes (see Fig. 1). In the following, specific topographic classes are denoted as elevation_aspect_slope; e.g. 1800_N_40 stands for a 40^∘ slope with a northern aspect at 1800 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$

2.2 CrocO ensemble data assimilation setup

The ensemble data assimilation workflow of CrocO is represented in Fig. 2. In the following, only a short description of the system and its elements is provided. More details on the ensemble modelling setup are available in Cluzet et al. (2020a). Information about its implementation into the Météo-France high-performance computing (HPC) system can be found in Appendix B1.

Figure 2Workflow of CrocO ensemble data assimilation system with four members. $\widehat{x_{\mathrm{0}}}$: initial state at time t₀, F_i: forcing, M_i: ESCROC member, X_b: background state, $x_{\mathrm{b}}^i$: background particles, X_a: analysis, $x_{\mathrm{a}}^i$: analysis particles, y: observation, t₁ and t₂: observation dates.

Download

2.3 The particle filter equations

On a given observation date, we consider a set of observed variables available at several locations, totalling N_y different observations.

• Each member ${\mathit{x}}_{\mathrm{b}}^i$ of the background state X_b is projected into the observation space using the observation operator h. In our case, h is just an orthogonal projection of the N_y observations since HS and reflectance are diagnosed within Crocus (see Sect. 2.2.1). The projection ${\widehat{\mathit{x}}}_{\mathrm{b}}^i=h{\mathit{x}}_{\mathrm{b}}^i=({\widehat{\mathit{x}}}_k^i{)}_{\mathrm{0}<k\le N_y}$ corresponds to the modelled values at each observed variable and/or point.

• These N_y observations are collected in the vector $\mathit{y}=(y_k{)}_{\mathrm{0}<k\le N_y}$. The associated observation error covariance matrix R (Eq. 1) is diagonal (e.g. observation errors are assumed independent):

  $$\begin{array}{}\text{(1)}& \mathbf{R}=\text{diag}({\mathit{\sigma }}_k^{\mathrm{2}},\mathrm{0}<k\le N_y),\end{array}$$

  where ${\mathit{\sigma }}_k^{\mathrm{2}}$ stands for the observation error variance of observation k and depends only on the type of observation y_k (e.g. HS or reflectance).

The PF analysis usually works in two steps.

• 1.

  Compute the particle weights wⁱ as the normalised observation likelihood for each particle (Eq. 2):

  $$\begin{array}{}\text{(2)}& w^i={\displaystyle \frac{e^{-\frac{\mathrm{1}}{\mathrm{2}}(\mathit{y}-{\widehat{\mathit{x}}}_{\mathrm{b}}^i{)}^T{\mathbf{R}}^{-\mathbf{1}}(\mathit{y}-{\widehat{\mathit{x}}}_{\mathrm{b}}^i)}}{{\sum }_{k=\mathrm{1}}^{N_{\mathrm{e}}}{e}^{-\frac{\mathrm{1}}{\mathrm{2}}(\mathit{y}-{\widehat{\mathit{x}}}_{\mathrm{b}}^k{)}^T{\mathbf{R}}^{-\mathbf{1}}(\mathit{y}-{\widehat{\mathit{x}}}_{\mathrm{b}}^k)}}}.\end{array}$$
• 2.

  Resample the particles based on their weights to build the analysis vector X_a. Here, we apply the PF resampling from Kitagawa (1996), which returns $\mathit{s}=(s_i{)}_{\mathrm{0}<i\le N_{\mathrm{e}}}$ ($s_i\in [\mathrm{1}..N_{\mathrm{e}}]$), a sorted vector with duplications representing the particles to replicate.

Figure 3Impact of the inflation ($N_{\text{eff}}^{*}=\mathrm{7}$) vs. no inflation ($N_{\text{eff}}^{*}=\mathrm{1}$) in the 1800_N_40 topographic class (not observed) when assimilating HS of 2015_q80 with the global PF. (a) SWE minimum–maximum envelopes as a function of time, (b) spread and (c) AEM. Dashed lines represent the assimilation dates.

Download

A sample reordering step was added for numerical optimisation with no expected influence on the PF behaviour (see Appendix B2 for more details).

Two simple variants of this algorithm can be identified in a spatialised context:

• for the global approach, perform one analysis over the domain, putting all the available observations in y;

• for the rlocal approach, perform one analysis per model point, assimilating only local observations if any. This corresponds to a localised PF with a block and domain size of 1.

2.3.3 Particle filter and reflectance observations

Assimilating reflectance with the PF requires some adaptations. In Crocus, the TARTES optical scheme (see Sect. 2.2.1) only provides snow reflectance, not all-surface reflectance: no value for the surface reflectance is issued in the absence of snow. Conversely, the weights of the particles are not defined in Eq. (2) if the members are snow-free. These issues were roughly accommodated by setting the reflectances of snow-free members and observations to 0.2 (the value of bare soil broadband albedo in the ISBA model) in the PF for Eq. (2) (Sect. 2.2.3).

Table 1Setup for the height of snow assimilation experiment.

Download Print Version | Download XLSX

3.1 Experiments

3.1.1 Twin experiment setup

In our twin experiment setup, an open-loop ensemble is used as a reference and to generate synthetic observations. Open-loop simulations are carried out with CrocO for four consecutive winters (2013–2017) in the Grandes Rousses (see Sect. 2.1) with 160 members. For each year, the average of SWE over time and space is computed from each member, and members corresponding to the 20th, 40th, 60th and 80th percentiles of the ensemble are extracted to be used as synthetic observations (denoted year_ppercentile, e.g. 2014_p80). This method enables us to evaluate the efficiency of data assimilation experiments under contrasting snow condition scenarios. Before any assimilation experiment, the open-loop member (F_i−M_i couple in Fig. 2) used as the true state is withdrawn and replaced by a new random member.

The spatial coverage of synthetic observations was reduced, mimicking a typical reflectance mask. Synthetic observations were only available above an assumed constant treeline at 1800 m (see Fig. 1) and not available for steep slopes (over 20^∘) and northern aspects (shadows, considering a daily satellite pass around 10:00–11:00 UTC) for the whole snow season. As a result, in this case, only 35 (over 187) topographic classes are observed. Observation dates were chosen corresponding to clear-sky days with a MODIS overpass, resulting in an approximately weekly frequency (e.g. Revuelto et al., 2018; Cluzet et al., 2020a).

Reflectance is sensitive to the surface SSA and LAP (see Sect. 1). A minimal set of two different bands is used, corresponding to MODIS sensor bands 4 (555 nm, sensitive to SSA and LAP) and 5 (1240 nm, usually only sensitive to SSA) (e.g. Fig. 2 of Cluzet et al., 2020a). Observation error variances are set to 1.0 × 10⁻² m² for HS and 5.6 × 10⁻⁴ and 2.0 × 10⁻³ for band 4 and band 5 reflectance, respectively (Wright et al., 2014). These values are only initial values for the inflation in the global and rlocal algorithms. Since the klocal algorithm only uses inflation if k drops to 1 (see Sect. 2.3.2), observation error variances are multiplied by a factor of 5 to enable the klocal algorithm to ingest observations from several points.

Table 2Setup for the first reflectance assimilation experiment.

Download Print Version | Download XLSX

Table 3Setup for the second reflectance assimilation experiment.

Download Print Version | Download XLSX

In order to study the ability of the global, klocal and rlocal algorithms to spread information in space, a first set of experiments is conducted assimilating HS with 40 members (see the setup in Table 1). In order to evaluate the algorithms' ability to assimilate reflectance (band 4 and band 5) a second set of experiments is conducted with other things being equal (Table 2). The ensemble size is increased from 40 to 160 in a third set of experiments assimilating reflectance in order to analyse the influence of a larger ensemble on the algorithm performance (Table 3). Note in Tables 1–3 that $N_{\text{eff}}^{*}$ is adjusted to the ensemble size in order to preserve $N_{\mathrm{e}}/N_{\text{eff}}^{*}\approx$ 5–7 following Larue et al. (2018).

3.2 Evaluation scores

The performance of the assimilation and open-loop run is evaluated against the synthetic truth using several scores. The absolute error of the ensemble mean (AEM) and ensemble spread σ are two common metrics of ensemble modelling. Given an ensemble $E_{m,c,t}$ of N_e members m in topographic class c at time t and the corresponding truth τ_c,t, the ensemble mean is described by Eq. (4):

from which we can compute the AEM (Eq. 5) and the spread (or dispersion) σ (Eq. 6):

where N_t is the number of evaluation time steps.

The continuous ranked probability score (CRPS; Eq. 7; Matheson and Winkler, 1976) evaluates the reliability and resolution of an ensemble based on a verification dataset. An ensemble is reliable when events are forecast with the right probability and has a good resolution when it is able to discriminate distinct observed events. For a reliable system, the resolution is equivalent to the sharpness, which is the spread of the produced forecasts.

If we denote F_c,t the cumulative distribution function (CDF) and 𝒯_c,t the corresponding truth CDF (Heaviside function centred on the truth value), the CRPS is computed at (c,t) following

In this work, the CRPS_c,t value is averaged over time alone or time and space depending on the desired level of aggregation.

The CRPS can be decomposed into two terms following Candille et al. (2015):

where Reli quantifies the reliability of the ensemble. The associated skill scores (CRPSS and ReliS) can be used to compare the performance of an ensemble E to a reference R, here the open-loop run:

A skill score of 1 denotes a perfect score, 0 a neutral performance and −∞ the worst achievable skill score.

4.1 Preliminary results

4.1.2 Correlation patterns

The klocal algorithm relies on background correlation patterns to define localisation domains. To illustrate the potential of using such information in the PF, Fig. 4 shows the correlation patterns of the 40-member open loop in an unobserved topographic class (1800_N_40, red dot) in the midwinter (20 February 2015), several months after the snow season onset. The assimilation variables exhibit strong but contrasting correlation patterns. Band 4 (Fig. 4a) correlations are generally high (0.6–1) and uniform. Many of the observed classes (black dots) are strongly correlated with the considered class. Similar results are obtained for HS (Fig. 4c). Band 5 (Fig. 4b) exhibits substantial correlations, in particular across slopes. However, they are more restricted to the northern aspects, and only a few observed classes in the eastern aspects are substantially correlated with the considered class. Note that negative correlations are evidenced with some lower-altitude south-oriented topographic classes (e.g. 1500_S_40 in Fig. 4b). Finally, these patterns vary with time but remain substantial along the whole season (not shown), and increasing the ensemble size up to 160 leads to identical patterns (not shown).

Figure 420 February 2015 open-loop (40 members) Pearson correlations between the domain points and the 1800_N_40 topographic class (red dot) in band 4 (a), band 5 (b) and HS (c). Left bars show the flat topographic classes in the associated elevation bands, while pie plots show the 20 and 40^∘ slope topographic classes, as depicted in Fig. 1. Black dots denote the observed classes.

Download

4.2 Results of the experiments

4.2.1 Assimilation of the height of snow

In a first step, assimilation of HS from the different synthetic observation scenarios was conducted to serve as a reference for reflectance assimilation. Figure 5 shows the CRPSS (Eq. 9, aggregated over time only) of the HS assimilation with the three PF algorithms considering the synthetic member 2013_q20 as a reference. Results for this specific synthetic member were chosen here as a representative example of the algorithm performance.

Figure 5CRPSS of SWE for the rlocal (a), global (b) and klocal (c) algorithms assimilating the HS of 2013_p20 synthetic observation scenario. The score is computed for the whole snow season for each topographic class. Black dots denote the observed classes.

Download

The rlocal performance compared with the open loop is high (0.7–1) but limited to the observed classes (black dots) since there is no spatial propagation in this algorithm. The global and klocal algorithms have similar overall good performance, managing to strongly reduce modelling uncertainties except at very low altitudes (600–900 m) (skills of −0.2) where snow does not usually last for more than a few weeks.

Figure 6Box plots of SWE CRPS (a, b) and Reli (e, f) for the different algorithms for the 16 different synthetic observation scenarios, separated between observed (a, c, e, g) and not observed (b, d, f, h) classes. Panels (c, d) and (g, h) show the associated skill scores.

Download

This behaviour may vary with the snow conditions, i.e. between the different assimilated synthetic observation scenarios and from one year to another. In order to generalise this result, Fig. 6 shows the CRPS and Reli (aggregated over time and space) of the different algorithms for the 16 synthetic observation scenarios and differentiated between observed and unobserved classes. CRPS and reliability are considerably reduced compared with the open loop (by a factor of 2–3 and 4–5, respectively) for all the algorithms in the observed classes. This suggests that the PF manages to reduce the spread of the ensemble while reducing its errors. In the unobserved classes, the gain is almost as good (CRPSS of 0.6) except for the rlocal algorithm, which is identical to the open loop as expected. No significant difference in skill is obtained between the global and klocal algorithms.

4.2.2 Assimilation of reflectance

Optical reflectance is a promising assimilation variable due to its extended availability in satellite observations, but assimilation of raw reflectance products is not expected to constrain bulk variables like SWE or HS as much as HS assimilation. In order to assess this difference, we conduct assimilation of reflectance only in the same setup as in Sect. 4.2.1 with all other things being equal.

Figure 7Same as Fig. 6 for reflectance with 40 members (filled) and 160 members (hatched).

Download

Figure 7 shows the performance of the reflectance assimilation for the 16 synthetic observation scenarios with 40 members (filled boxes). The different algorithms only lead to moderate improvements in CRPS (median CRPSS of 0–0.2, median ReliS of 0.2–0.4). Moreover, the global and klocal algorithms frequently degrade the performance, suggesting that this configuration is not robust.

Suspecting that 40 members are insufficient to properly represent the multivariate probability density function of reflectance and other model variables, the ensemble size was increased to 160 (hatched boxes), leading to marked improvements in the performance and robustness of the algorithms (median CRPSS of 0.2, median Reli of 0.4–0.6). The reliability of the global algorithm is significantly improved compared to the klocal algorithm.

Figure 8Same as Fig. 5 for the assimilation of the reflectance of the 2016_p60 synthetic observation scenario.

Download

Figure 8 shows the spatial performance of the different algorithms for member 2016_p60. Spatial patterns similar to the HS assimilation are found. The rlocal performance is limited to the observed classes, while global and klocal manage to improve the simulations across aspects and slopes. However, skill scores are lower than for HS (0.2–0.5), and the performance of all algorithms is poor in the classes that are farther away from the observations, i.e. at lower elevations (600–900 m) and in some of the high-altitude steep northern classes (e.g. 2100_N_40 in Fig. 8b and c). Finally, note that slight degradations of performance can sometimes be evidenced, even in the observed classes, for all the algorithms (e.g. in flat conditions at 3300 m in Fig. 8a for the rlocal, not evidenced by this example for the other algorithms).

5.1 Tackling particle filter degeneracy

Because they assimilate several observations at the same time, global and klocal approaches could be prone to PF degeneracy. However, they almost never degrade the performances when assimilating HS in a variety of years and synthetic observation scenario percentiles (Fig. 6). This suggests that inflating the observation errors (as demonstrated by Larue et al., 2018, a result we have generalised in space) and exploiting background correlations to reduce the number of assimilated observations are two efficient approaches to tackle degeneracy.

In several cases though, strong degradation of the score occurs when assimilating reflectance (Fig. 7), which could either be attributed to an algorithmic failure in the PF or an intrinsic lack of informativeness of reflectance in some situations. Based on the good behaviour of the algorithm with HS, and because by construction, the global and klocal algorithms cannot lead to a degenerate PF sample, we consider this to come from the reflectance itself (this point will be further discussed in the following sections).

Beyond tackling degeneracy, the global and klocal algorithms also beat the rlocal approach on Reli and CRPSS (Figs. 7 and 8). This suggests that assimilating multiple observations increases the quality of the PF analysis, even locally. More precisely, most of the improvement is due to the Reli term of the CRPS. This property is crucial for ensemble modelling because it ensures that events are forecasted with the right frequency. However, this is not sufficient; e.g. the climatology has perfect reliability but is not informative at all. Successful assimilation manages to improve general metrics such as the CRPS while improving the reliability. For this aspect, the global and klocal algorithms have a satisfying performance.

5.2 Propagating the observation information

Having sparse observations is one of the most challenging issues for data assimilation systems of snowpack observations (Magnusson et al., 2014; Largeron et al., 2020). In our partially observed synthetic setup, the global and klocal PF variants developed here efficiently propagate the observational information to the unobserved classes with a generally better performance than the open-loop and rlocal approach in the unobserved classes when assimilating HS (Fig. 5).

The algorithms' performance is particularly good across aspects and slopes, with only a few steep northern aspect slopes exhibiting neutral to poor performances (Figs. 5 and 8). This suggests that southern aspect and flat classes are informative for the majority of the simulation domain. Conversely, considering that there are strong background correlations between the western and eastern sides of the domain, we can speculate that observing either side could yield overall good results.

In these figures, propagation of the information is limited towards lower elevation (600–1200 m). At such elevations, the snow cover is usually intermittent and a good discrimination of the precipitation phase is crucial. The PF does this indirectly through HS and reflectance observations because rain causes a decrease in HS through compaction and melting, while band 4 and band 5 reflectances also decrease because of quick isothermal metamorphism (i.e. the surface SSA decreases). However, in our setup, the lowest observed elevation is 1800 m; therefore, indirect observation of the rain–snow line positioning under this level is not possible, potentially explaining the moderate performance of the PF there. In that case, assimilation of snow cover fraction might be the best solution; since the snowpack is intermittent there, the informativeness of this variable is maximal (Aalstad et al., 2018).

The global and klocal algorithms exhibit strong performances when assimilating HS (Fig. 5). HS is closely linked to the SWE (by the bulk density), and the interest of this variable for data assimilation is clear (Margulis et al., 2019). Here, it should be kept in mind that HS assimilation is used as a baseline experiment to evaluate the algorithms and put reflectance assimilation into perspective. The prescribed HS observation errors (σ₀ = 0.1 m) are not necessarily realistic. They should be adapted to the nature of the HS sensor. For example, spaceborne HS observation errors are typically larger (e.g. Eberhard et al., 2020; Deschamps-Berger et al., 2020). The assimilation of such observations would probably yield lower improvements.

Though the performance is lower for reflectance than in our HS experiments, it remains considerable and in line with previous results on point simulations (Charrois et al., 2016), with an average score improvement of 20 %–40 %. This study quite surprisingly suggests that reflectance information can be spread from southern slopes to the northern ones, although in many situations, the snowpack evolves in different ways for these two aspects. For example, in sunny conditions, melt and wet metamorphism will cause a drop in reflectance on southern slopes, while reflectance will not evolve much on northern slopes. Such a phenomenon could explain why low background correlations between southern and northern aspects are exhibited in band 5 (Fig. 4), which is the most sensitive to surface metamorphism through SSA. This example shows that band 5 reflectance observations on southern slopes are not necessarily informative for band 5 reflectance values in the northern aspect per se on every date. On average, however, the positive impact of reflectance observations suggests that they enable the PF to reject the ensemble members with inadequate meteorological forcings (snowfall or cloud cover would lead to wrong reflectance values) or multi-physical parameterisations (influencing e.g. the surface metamorphism), thus correcting the ensemble in the whole domain. These insights are consistent with the study of Winstral et al. (2019), wherein in situ observations are used to correct meteorological forcing parameters across large simulation domains.

Regarding the observations, our study has some methodological limits, however. Observation errors are very roughly prescribed, and the assimilated observations are not corrupted as usually done in synthetic experiments (e.g. Durand and Margulis, 2006). These choices were motivated by the fact that very little is known about the spatial correlation of reflectance observation errors in a semi-distributed setting (e.g. Cluzet et al., 2020a). In a recently submitted paper, the impact of random and systematic errors in reflectance observations on point-scale assimilation experiments is thoroughly investigated (Revuelto et al., 2021). Efforts to better characterise the spatial structure of these observation errors should be conducted in future work.

5.3 Towards the assimilation of real observations of reflectance

Reflectance is an appealing variable for snowpack modelling because of its sensitivity to snowpack surface properties (Dozier et al., 2009) and the abundance of moderate- to high-resolution spaceborne sensors (MODIS, Sentinel-2–3, VIIRS, Landsat) providing us with a handful of observations to assimilate, contrary to HS. The potential for assimilation of SCF, which is retrieved from reflectances, is clear (Margulis et al., 2016; Aalstad et al., 2018; Alonso-Gónzalez et al., 2020). This study demonstrates the potential of the PF to spread information and assimilate raw reflectances with a positive impact (Sect. 5.2). Yet, assimilating real observations of reflectance is another challenge for two reasons.

First, spaceborne reflectance observations are generally noisy and biased (e.g. Cluzet et al., 2020a). Satellite retrievals could be improved in the future (Kokhanovsky et al., 2019; Lamare et al., 2020), and Cluzet et al. (2020a) showed that assimilating ratios of reflectance could be a workaround to tackle this issue. In the near-infrared, the signal-to-noise ratio of reflectance observations might be sufficient to constrain the surface microphysical properties (Durand and Margulis, 2007; Mary et al., 2013), whereas the required accuracy for visible reflectance retrievals to remain informative on snowpack light-absorbing particles content is high (Warren, 2013), and it has yet to be proved whether either approach can achieve this requirement.

Second, in this twin experiment framework, spatial patterns of the synthetic observations are likely compatible with the ensemble since they come from the same modelling system. This may not be the case in reality, therefore making it more difficult to assimilate, and we refer to this issue as model or ensemble realism.

We must assess the strengths and weaknesses of the global and klocal approaches by addressing those two issues. The global algorithm assumes that a global optimum can be found across the whole domain; e.g. the information from the different observations is consistent and can be ingested in one block by the PF. With this strategy, the degeneracy due to the size of the observation vector is efficiently mitigated by the inflation algorithm as discussed in Sect. 5.1. The klocal approach considers only a fraction of the observation information to be relevant to constrain the model state at a given location. This algorithm tries to ingest as much information as possible while rejecting observations coming from snowpack conditions that are too statistically different. As a consequence, because we do not account for the real spatial patterns of observation errors and because we work in a twin experiment setup, a global optimum for the whole domain can exist and can be found by the global algorithm. This might be a reason why it beats the klocal approach (Figs. 6 and 7). In the real world, from the model point of view, there might be contradictory information among the observations that would be difficult to disentangle with a global strategy. The klocal algorithm could be more suited to this situation because it looks for local optima based on the assumption that background correlations are a realistic representation of modelling errors.

These background correlation structures could be overestimated by the ensemble, and tests with real observations are necessary. Strong band 4 correlations (Fig. 4a) might be due to the spatially homogeneous perturbations of LAP fluxes used to force the simulations (see Sect. 2.2.2), a key driver of this variable, and because the same snow model configuration is applied for a given member across the simulation domain. Several studies suggest that LAP fluxes vary with elevation and other topographic parameters (de Magalhães et al., 2019; Sabatier et al., 2020), but to date no reliable model of such processes has been developed for complex terrain. In such a context, assuming uniform LAP forcing seems a reasonable compromise. Strong and almost uniform HS correlations (Fig. 4b) might be caused by the spatial homogeneity of precipitation perturbations and because we do not account for e.g. wind drift, intra-massif heterogeneity of meteorological conditions and gravitational redistribution of snow (Wayand et al., 2018). Despite this semi-distributed framework suffering from obvious limitations, the potential for high-resolution snowpack modelling (Vionnet et al., 2020; Fiddes et al., 2019; Marsh et al., 2020) is hampered by large errors of NWP models in mountainous areas (e.g. Nousu et al., 2019).

In the future, improving the ability of ensemble correlations to represent modelling errors could make the spreading of information an even more challenging task with the klocal algorithm. But significant potential should remain for information propagation, as suggested by results at larger scales (Magnusson et al., 2014; Cantet et al., 2019). The potential decorrelation of topographic classes would also impact the global algorithm. In an unobserved class, constraining the state of the snowpack with information from areas that are not linked to it would likely degrade the forecasting skill, as suggested by the poor performance of the algorithms at low altitudes (Figs. 5 and 8). In contrast, applying CrocO over larger domains (e.g. distributed simulations or a collection of semi-distributed massifs) would probably see the klocal algorithm outperform the global. The increased domain size would make it less plausible to find a global optimum over the domain, whereas spatial flexibility would be an asset of the klocal algorithm. Finally, in the case of modelled coupling between simulation points (e.g. snow drift), which was not the case here, the spatial discontinuities of the klocal analyses (see Sect. 1) might be a drawback compared to the global approach. Spatial discontinuities may also be revealed as impractical for the interpretation of individual simulation outputs by snow forecasters. The klocal approach is likely to reduce these discontinuities compared to the rlocal because similar locations will be treated with similar analyses (i.e. based on similar sets of observations). This issue could be partly mitigated by e.g. state–block–domain approaches (Farchi and Bocquet, 2018).

5.4 Outlook for ensemble modelling and data assimilation

In the snowpack modelling community, ensemble modelling is a powerful tool to represent modelling uncertainties (Vernay et al., 2015; Richter et al., 2020) and for data assimilation (Essery et al., 2013; Lafaysse et al., 2017; Piazzi et al., 2018; Aalstad et al., 2018). This study offers a novel approach to extract valuable information on the snowpack spatial behaviour from spatial correlation patterns of the ensemble. These patterns could be used to diagnose links between locations, transfer information between areas or assess the representativeness of point simulations. More broadly, ensemble background correlations have long been exploited in the NWP and oceanographic communities to refine modelling error representation, which led to significant improvements in the data assimilation systems (Evensen, 2003; Buehner, 2005).

Figure 9Same as Fig. 5 for the assimilation of HS of the 2016_p60 synthetic observation scenario in the 1200–2400 m flat classes.

Download

Ensembles might open a possibility for the assimilation of point-scale observations or sparse remotely sensed observations into spatialised simulations of the snowpack, as suggested by Winstral et al. (2019) and the present work. For instance, there are numerous snow gauges and snow pit observations at ski resorts in the French Alps. These data could be assimilated to correct the ensemble in spatialised simulations (Winstral et al., 2019). The spatial pattern of assimilated observations in the experiments of Sect. 4 does not correspond to the real-life spatial coverage of these kinds of observations. To give insight into their potential, we also applied our methodology to assimilate only five synthetic HS observations with the global PF in the 1200 to 2400 m flat classes. The results are shown in Fig. 9. The assimilation improves the performance in all aspects and slopes. Naturally, this suffers from the same limitation as discussed in Sect. 5.3, not to mention the limited spatial representativeness of in situ observations, but it shows some potential for this idea.

In that way, a more rational use of the available observations could be implemented towards a new ensemble data assimilation system. In the present CrocO system, the SAFRAN reanalysis only assimilates weather station information (precipitation phase, temperature, wind) and makes no use of the numerous snow observations available. Here, snow observations are assimilated by the PF but are not used to correct meteorological forcings (only snow variables; see Fig. 2). In a new ensemble data assimilation system, within CrocO, the SAFRAN meteorological analysis could be bypassed, with the PF directly operating on both the meteorological and snowpack variables through a more comprehensive and coupled strategy.