Articles | Volume 14, issue 3
Geosci. Model Dev., 14, 1595–1614, 2021
Geosci. Model Dev., 14, 1595–1614, 2021

Development and technical paper 19 Mar 2021

Development and technical paper | 19 Mar 2021

CrocO_v1.0: a particle filter to assimilate snowpack observations in a spatialised framework

CrocO_v1.0: a particle filter to assimilate snowpack observations in a spatialised framework
Bertrand Cluzet1, Matthieu Lafaysse1, Emmanuel Cosme2, Clément Albergel3,a, Louis-François Meunier3, and Marie Dumont1 Bertrand Cluzet et al.
  • 1Univ. Grenoble Alpes, Université de Toulouse, Météo-France, CNRS, CNRM, Centre d'Études de la Neige, Grenoble, France
  • 2Institut des Géosciences de l'Environnement, IGE, UGA-CNRS, Grenoble, France
  • 3CNRM, University of Toulouse, Météo-France, CNRS, 31057 Toulouse, France
  • anow at: European Space Agency Climate Office, ECSAT, Harwell Campus, Oxfordshire, Didcot OX11 0FD, UK

Correspondence: Bertrand Cluzet ( and Matthieu Lafaysse (


Monitoring the evolution of snowpack properties in mountainous areas is crucial for avalanche hazard forecasting and water resources management. In situ and remotely sensed observations provide precious information on the state of the snowpack but usually offer limited spatio-temporal coverage of bulk or surface variables only. In particular, visible–near-infrared (Vis–NIR) reflectance observations can provide information about the snowpack surface properties but are limited by terrain shading and clouds. Snowpack modelling enables the estimation of any physical variable virtually anywhere, but it is affected by large errors and uncertainties. Data assimilation offers a way to combine both sources of information and to propagate information from observed areas to non-observed areas. Here, we present CrocO (Crocus-Observations), an ensemble data assimilation system able to ingest any snowpack observation (applied as a first step to the height of snow (HS) and Vis–NIR reflectances) in a spatialised geometry. CrocO uses an ensemble of snowpack simulations to represent modelling uncertainties and a particle filter (PF) to reduce them. The PF is prone to collapse when assimilating too many observations. Two variants of the PF were specifically implemented to ensure that observational information is propagated in space while tackling this issue. The global algorithm ingests all available observations with an iterative inflation of observation errors, while the klocal algorithm is a localised approach performing a selection of the observations to assimilate based on background correlation patterns. Feasibility testing experiments are carried out in an identical twin experiment setup, with synthetic observations of HS and Vis–NIR reflectances available in only one-sixth of the simulation domain. Results show that compared against runs without assimilation, analyses exhibit an average improvement of the snow water equivalent continuous rank probability score (CRPS) of 60 % when assimilating HS with a 40-member ensemble and an average 20 % CRPS improvement when assimilating reflectance with a 160-member ensemble. Significant improvements are also obtained outside the observation domain. These promising results open a possibility for the assimilation of real observations of reflectance or of any snowpack observations in a spatialised context.

1 Introduction

Seasonal snowpack is an essential element of mountainous areas. Monitoring the evolution of its physical properties is essential to forecasting avalanche hazard (Morin et al.2020) and rain-on-snow-related floods (Pomeroy et al.2016; Würzer et al.2016) as well as monitoring water resources (Mankin et al.2015). Observations alone are too scarce to monitor snowpack conditions. In situ observations provide precise observations of several key variables, but they lack spatial representativeness and have poor spatial coverage. Remote sensing of snowpack variables such as the height of snow (HS; m), snow water equivalent (SWE; kg m−2), visible–near-infrared (Vis–NIR) reflectance and surface temperature provides comprehensive information over large areas but usually has a limited temporal resolution for a small set of variables. Furthermore, these observations are usually available in fractions of simulation domains only, even for spaceborne data (Davaze et al.2018; Veyssière et al.2019; Shaw et al.2019). For instance, snowpack Vis–NIR reflectances from moderate-resolution (250–500 m) satellites such as MODIS and Sentinel-3 can help constrain snowpack surface properties such as microphysical properties (characterised by specific surface area – SSA; m2 kg−1) and light-absorbing particle content (LAP; ggsnow-1) (Durand and Margulis2006; Dozier et al.2009). However, in areas covered by clouds or forests and/or affected by high sub-pixel variability (ridges, roughness, fractional snow cover) and shadows, satellite retrievals are less accurate (Masson et al.2018; Lamare et al.2020), and data should be filtered out (Cluzet et al.2020a). The higher resolution offered by products from Landsat and Sentinel-2 might be an avenue to address this issue (e.g. Masson et al.2018; Aalstad et al.2020), but at these resolutions, reflectance retrievals are quite noisy due to e.g. digital elevation model errors (Cluzet et al.2020a). Finally, note that pixel fractional snow cover (snow cover fraction, SCF) can be accurately retrieved even from noisy reflectances (Sirguey et al.2009; Aalstad et al.2020), but it inherits spatio-temporal limitations. SCF informativeness is also limited in deep snowpack conditions (De Lannoy et al.2012).

Snowpack models of different complexity offer exhaustive spatial and temporal coverage (Krinner et al.2018). They are applied within several spatial configurations, including the collection of points on regular or irregular grids (Morin et al.2020). In this paper, “spatialised” indistinctly refers to any of these configurations. Detailed snowpack models are the only ones able to assess avalanche hazard and monitor water resources (Morin et al.2020), but these applications are limited by their considerable errors and uncertainties (Essery et al.2013; Lafaysse et al.2017). In that context, combining remote sensing observations with models through data assimilation is an appealing solution (Largeron et al.2020). Indeed, data assimilation combines the spatial and temporal coverage of snowpack models with the available information from observations in an optimal way. Assimilation of optical reflectance could reduce modelled SWE errors by up to a factor of 2 (Durand and Margulis2007; Charrois et al.2016), and preliminary studies show its potential for spatialised assimilation (Cluzet et al.2020a). Assimilation of HS is very efficient in reducing modelled SWE errors (Margulis et al.2019). However, the limited spatial coverage of observations stresses the need for data assimilation algorithms able to propagate snowpack observational information into unobserved areas (Winstral et al.2019; Cantet et al.2019; Largeron et al.2020).

A particle filter with sequential importance resampling (PF-SIR; Gordon et al.1993; Van Leeuwen2009) is a Bayesian ensemble data assimilation technique well suited to snowpack modelling (Dechant and Moradkhani2011; Charrois et al.2016; Magnusson et al.2017; Piazzi et al.2018; Larue et al.2018). PF-SIR is a sequential algorithm relying on an ensemble of model runs (particles) which represents the forecast uncertainty. At each observation date, the prior (or background) composed of the particles is evaluated against the observations. The analysis of PF-SIR (later on PF) works in two steps. In a first step, so-called “importance sampling”, the particles are weighted according to their distance to the observations (relative to the observation errors). Then, a resampling of the particles is performed in order to reduce the variance in the weights. The ensemble Kalman filter (EnKF Evensen2003) has also been widely used for snow cover data assimilation (e.g. Slater and Clark2006; De Lannoy et al.2012; Magnusson et al.2014). However, the PF is more adapted to models with a variable number of numerical layers such as detailed snowpack models (Charrois et al.2016).

The PF could be used in a spatialised context to propagate information from observations as suggested by Largeron et al. (2020) and Winstral et al. (2019). Contrary to the EnKF, such applications are rare to date (e.g. Thirel et al.2013; Baba et al.2018; Cantet et al.2019). Indeed, spatialised data assimilation with the PF is not straightforward because of the degeneracy issue, i.e. only a few particles are replicated in the analysis, often resulting in a poor representation of the forecast uncertainties. Degeneracy can be mitigated by increasing the number of particles, but the required population scales exponentially with the number of observations simultaneously assimilated (Snyder et al.2008). Furthermore, an accurate representation of spatial error statistics by the ensemble is essential for the success of the assimilation system. To achieve that, the required ensemble size also scales exponentially with the system dimension, an issue known as the curse of dimensionality (Bengtsson et al.2008). These issues are severe drawbacks when considering applications of the PF to large domains (i.e. implying a large number of observations and/or simulation points) with a reasonable number of particles (Stigter et al.2017).

Several solutions exist to tackle PF degeneracy. A first approach is to inflate the observation errors in the PF. The tolerance of the PF is increased, leading to more particles being replicated. This approach is based on the fact that observation error statistics (including sensor, retrieval and representativeness errors) are usually poorly known and underestimated. It can also be used as a safeguard to prevent the PF from degenerating on specific dates when observations are not compatible with the ensemble. PF inflation was successfully implemented in point-scale simulations of the snowpack (Larue et al.2018). When dealing with a large number of observations, inflation might lead to degeneracy or null analysis (posterior equal to the prior). In this work, we generalise over space the inflation of Larue et al. (2018), trying to ingest all the observations into a single analysis over the domain in a so-called global approach.

PF localisation is a more widespread alternative, tackling degeneracy by reducing the number of observations that are simultaneously assimilated by the PF (Poterjoy2016; Poterjoy and Anderson2016; Penny and Miyoshi2016; Poterjoy et al.2019; italic notations are taken from the review of Farchi and Bocquet2018). In this method, the simulation domain is divided into blocks whereby different PF analyses are performed considering a local subset of observations (domain) based on a localisation radius. This makes it possible to constrain the model in locations that are not directly observed but with nearby observations. Contrary to global approaches, localisation has the disadvantage of producing spatially discontinuous analyses (each point receives a different analysis). This issue can be mitigated in various ways (Poterjoy2016; Farchi and Bocquet2018; Van Leeuwen et al.2019).

The underlying hypothesis of localisation is that model points are independent beyond a certain distance; i.e. constraining one point with the observation from a point that is too distant would be meaningless and likely degrade the analysis performance (Houtekamer and Mitchell1998). However, in the case of small simulation domains or modelled systems driven by large-scale coherent causalities, large-scale correlations (relative to the domain size) may be physically sound, and defining a localisation radius may be a difficult task. In order to address this issue, we developed a new localisation approach called k localisation, whereby localisation domains are based on background correlation patterns.

These developments were implemented into CrocO (Crocus-Observations), an ensemble data assimilation system able to sequentially assimilate snowpack observations with a PF in a spatialised context. CrocO can be implemented in any geometry (e.g. within a distributed (gridded) framework or any irregular spatial discretisation). Here, we apply CrocO in a semi-distributed framework, which is a conceptual spatialised geometry used operationally by Météo-France for avalanche hazard forecasting (Lafaysse et al.2013; Morin et al.2020). This framework is similar to many topographically based discretisations in hydrological models (e.g. Clark et al.2015). This setup enables us to account for the snowpack variability induced by the topography at the scale of a mountain range through meteorological conditions (elevation controls the air temperature and precipitation phase) and the snowpack radiative budget (also dependent on the aspect and slope angle) (Durand et al.1993).

CrocO uses an ensemble of stochastic perturbations from the SAFRAN meteorological analysis (Durand et al.1993; Charrois et al.2016) to force ESCROC (Ensemble System CROCus; Lafaysse et al.2017), the multi-physical version of the Crocus snowpack model (Vionnet et al.2012). The ensemble setup accounts for the major sources of uncertainties in snowpack modelling (Raleigh et al.2015) and was formerly described and evaluated in semi-distributed geometry by Cluzet et al. (2020a).

Inflation and k localisation were implemented into CrocO. Here, we present CrocO and evaluate how it addresses the issues of reflectance observation sparseness and PF degeneracy in the context of snowpack modelling. This problem is divided into two scientific questions. (1) Is CrocO PF able to efficiently spread the information from sparse observations in space without degenerating? (2) Is the spatial information content of reflectance observations valuable for snowpack models? We assess these questions by evaluating the performance of CrocO in modelling the SWE when assimilating synthetic observations of HS and reflectance covering only a portion of the domain.

Section 2 presents the CrocO system, i.e. the ensemble modelling system and the PF algorithms. Section 3 introduces the evaluation methodology. Subsequently, Sect. 4 assesses the performance of CrocO, and Sect. 5 discusses the results. Finally, Sect. 6 provides perspectives and research directions.

2 Material and methods

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 km2 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.


In this study, we focus on the Grandes Rousses, a single massif in the central French Alps. This area of about 500 km2 is represented by Npts=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 ma.s.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. x0^: initial state at time t0, Fi: forcing, Mi: ESCROC member, Xb: background state, xbi: background particles, Xa: analysis, xai: analysis particles, y: observation, t1 and t2: observation dates.


2.2.1 Ensemble of snowpack models

Crocus is a detailed snowpack model coupled with the ground and atmosphere in the ISBA land surface model (Interaction Soil–Biosphere–Atmosphere). It is embedded within the SURFEX_v8.1 modelling platform (SURFace EXternalisée; Masson et al.2013). The TARTES optical scheme (Libois et al.2013, 2015) represents Vis–NIR spectral radiative transfer within the snowpack, driven by snow metamorphism (Carmagnola et al.2014) and light-absorbing particle (LAP; ggsnow-1) deposition fluxes (Tuzet et al.2017). Moreover, TARTES computes the snowpack reflectance with a high spectral resolution, making the model directly comparable to observations. As such, TARTES is both a physical component of Crocus and an observation operator.

ESCROC (Ensemble System CROCus; Lafaysse et al.2017), the multi-physical ensemble version of Crocus, is used to account for snowpack modelling uncertainties. A random draw among 1944 ESCROC multi-physics configurations was performed and used in all the simulations and denoted (Mi)0<iNe, with Ne being the ensemble size (e.g. 40 or 160 members; see Fig. 2). These configurations are considered equiprobable before any data assimilation.

2.2.2 Ensemble of meteorological forcings

Meteorological forcings are taken from the SAFRAN (Durand et al.1993) reanalysis, wherein forecasts from the ARPEGE numerical weather prediction (NWP) model are downscaled and adjusted with surface observations within the massif area. They are combined with MOCAGE LAP fluxes (Josse et al.2004) interpolated at Col du Lautaret (2058 ma.s.l., inside the Grandes Rousses) to constitute the reference forcing dataset. Before the beginning of the simulation, spatially homogeneous stochastic perturbations (e.g. on a given date, the same perturbation parameter is applied across the whole domain) with temporal autocorrelations are applied to this forcing to generate an ensemble of forcings (Fi)0<iNe with the same procedure as described in Cluzet et al. (2020a). More details on the perturbation procedure can be found in Appendix A. At the beginning of the simulation, each forcing Fi is associated with a random Mi ESCROC configuration, and this relation is fixed during the whole simulation.

2.2.3 The particle filter in CrocO

The PF is applied sequentially for each observation date to the background state vectors (soil and snowpack state variables, denoted BG in Fig. 2). Its analysis is an ensemble of initial conditions used to propagate the model forward. The algorithm is implemented into SODA (SURFEX Offline Data Assimilation; Albergel et al.2017), the data assimilation module of SURFEX_v8.1, enabling a continuous execution sequence between ensemble propagation and analysis, as depicted in Fig. 2.

2.3 The particle filter equations

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

  • Each member xbi of the background state Xb is projected into the observation space using the observation operator h. In our case, h is just an orthogonal projection of the Ny observations since HS and reflectance are diagnosed within Crocus (see Sect. 2.2.1). The projection x^bi=hxbi=(x^ki)0<kNy corresponds to the modelled values at each observed variable and/or point.

  • These Ny observations are collected in the vector y=(yk)0<kNy. The associated observation error covariance matrix R (Eq. 1) is diagonal (e.g. observation errors are assumed independent):

    (1) R = diag ( σ k 2 , 0 < k N y ) ,

    where σk2 stands for the observation error variance of observation k and depends only on the type of observation yk (e.g. HS or reflectance).

The PF analysis usually works in two steps.

  • 1.

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

    (2) w i = e - 1 2 ( y - x ^ b i ) T R - 1 ( y - x ^ b i ) k = 1 N e e - 1 2 ( y - x ^ b k ) T R - 1 ( y - x ^ b k ) .
  • 2.

    Resample the particles based on their weights to build the analysis vector Xa. Here, we apply the PF resampling from Kitagawa (1996), which returns s=(si)0<iNe (si[1..Ne]), a sorted vector with duplications representing the particles to replicate.

Figure 3Impact of the inflation (Neff*=7) vs. no inflation (Neff*=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.


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.1 Particle filter degeneracy

Degeneracy occurs when only a small fraction of the particles have non-negligible weights, resulting in a sample s for which only a few different indices are present. It can be diagnosed from the weights using the effective sample size Neff (Liu and Chen1995):

(3) N eff = 1 i = 1 N e ( w i ) 2 .

With a degenerate sample, Neff 1, and with innocuous analysis (all particles are replicated) Neff=Ne.

A first approach to mitigate degeneracy is to use inflation. This heuristic method iteratively inflates R values until the effective sample size is large enough. Here, we develop a variant from the Larue et al. (2018) method, which does not explicitly rely on Neff (Eq. 3). Consider applying an inflation factor 1α to R (0<α1, with α=1 being the value for no inflation) and update Neff (Eqs. 2 and 3): Neff is naturally a decreasing function of α (the more we inflate R the more different particles will be replicated). The idea of our method is to ensure that Neff exceeds a target value, Neff*. If Neff<Neff* (degenerate case), we reduce α (inflate) until Neff=Neff* using Algorithm 1. In the following, inflation is used in the global and rlocal PF (see Sect. 2.2.3).

The core of Algorithm 1 is an hybrid bisection–secant method to find the zero of f:αNeff(α)-Neff* in [0,1]. It is inspired by the rtsafe algorithm (Press et al.1992). The guess function computes a new guess α2 to minimise f. Note that in the unlikely case in which Algorithm 1 does not converge, all the particles are replicated.

2.3.2k localisation

In the k-localisation algorithm, degeneracy is mitigated by reducing the number of observations that are simultaneously assimilated. The PF analysis is applied to each simulation point sequentially. In order to build the analysis at point n, background correlations Bv are computed for each variable v (e.g. HS or reflectance) between n and all the observed points. In a first step, all observations from points exhibiting substantial background correlations (see below for the select_k_biggest function) are used. If the PF degenerates, the number of observations is progressively decreased until degeneracy is mitigated. As earlier, degeneracy is considered mitigated when NeffNeff*. This way, we ensure that a maximal number of observations has been ingested by the PF without degenerating.

In the case of degeneracy, the observation point displaying the lowest correlation is ruled out. The PF weights are computed (Eq. 2), and a new effective sample size is derived (Eq. 3). While the target sample size is not exceeded, this selection proceeds iteratively. The notation k in k localisation refers to the number, k, of retained observations for each variable. This approach is similar to the EnKF localisation algorithm whereby the localisation domain is based on background correlations (Hamill et al.2001). The detailed k-localisation algorithm is described in Algorithm 2, for which the following points apply.

  • For each variable, the select_k_biggest method returns the domain dv of up to k observed points (named p) that are the most correlated (in absolute value) with n and match the following criteria, which were adjusted in preliminary experiments.

    • In xvi, at least 10 % of members are defined in both points. As reflectance is not defined when there is no snow, spuriously high correlation can be obtained when the computation of correlations is based on a very low number of pairs.

    • |Bv(n,p)|>0.3. If the absolute correlation is too low, it is likely that there is poor potential for the distant observation to constrain the ensemble locally. In such a situation, it is better to reject the observation from the local analysis. Negative ensemble correlations can be physically sound, e.g. after a rain-on-snow event between the HS of two points separated by the rain–snow line. In such a situation, an HS observation of either point can hold information on precipitation rates at both locations. At the observed location, the PF will probably select the members with the most appropriate precipitation rates. This sample is likely to perform well at both locations, so it can be used to constrain the unobserved location.

  • The notation d represents the collection of domains dv.

  • The extract_points function extracts d from y,x^i and R.

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 Evaluation strategy

Our strategy is to assess the performance of the analysis by means of twin experiments, i.e. using synthetic observations (e.g. Reichle and Koster2003). The assimilation run is compared to an identical run without assimilation (open-loop run). Synthetic observations are extracted from a model run and assimilated without adding any noise. These observations mimic real observations with perfect knowledge of the true state. Analysis and open-loop experiments can therefore be compared with this true state anywhere for any variable. In a first step, this allows us to get rid of the error and bias issues inherent in real observations (Cluzet et al.2020a), a reason why we did not add any noise to the synthetic observations as commonly done in twin experiments (Lahoz and Menard2010). This way, we can focus on the following two questions (see Sect. 1).

  • Is CrocO PF able to efficiently spread the information from sparse observations into space without degenerating?

  • Is the spatial information content of reflectance a valuable source of information for snowpack models?

In order to disentangle these questions, we run baseline experiments assimilating synthetic observations of HS, which is strongly linked to SWE (Margulis et al.2019). These experiments are used to evaluate the PF algorithm efficiency and as a baseline for synthetic reflectance assimilation experiments evaluating the information content of reflectance.

Three different algorithms are evaluated: the global algorithm (with inflation), the rlocal algorithm (with inflation) and the k-localised algorithm klocal.

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 (FiMi 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−2m2 for HS and 5.6 × 10−4 and 2.0 × 10−3 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 13 that Neff* is adjusted to the ensemble size in order to preserve Ne/Neff* 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 Em,c,t of Ne members m in topographic class c at time t and the corresponding truth τc,t, the ensemble mean is described by Eq. (4):

(4) E c , t = 1 N e m = 1 N e E m , c , t ,

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


where Nt is the number of evaluation time steps.

The continuous ranked probability score (CRPS; Eq. 7; Matheson and Winkler1976) 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 Fc,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

(7) CRPS c , t = R ( F c , t ( x ) - T c , t ( x ) ) 2 d x ( c , t ) [ 1 , N pts ] × [ 1 , N t ] .

In this work, the CRPSc,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):

(8) CRPS = Reli + Resol ,

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:

(9) CRPSS ( E ) = 1 - CRPS ( E ) CRPS ( R ) .

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

4 Results

4.1 Preliminary results

4.1.1 Impact of the inflation

The inflation algorithm was introduced by Larue et al. (2018) in point-scale simulations, but to the best of our knowledge, it has never been applied in a spatialised context. Here we evaluate its impact on the global algorithm by switching it on and off. As an example, Fig. 3 shows the impact of inflation on SWE when assimilating the HS of 2015_p80 (as defined in Sect. 3.1.1) member with the global algorithm in a topographic class which is not observed (1800_N_40, as defined in Sect. 2.1). This choice of member and topographic class is representative of the impact of inflation on the global algorithm.

In this case, both inflation (Neff*=7) and no inflation (Neff*=1) lead to a significant reduction of the ensemble spread compared with the open loop (Fig. 3b). From January 2015 until the peak of SWE in mid-April 2015 (Fig. 3c), the simulation with inflation has significantly lower errors than without inflation and the open loop (10–20 vs. 60–80 and 30–50 kg m−2, respectively), leading to better agreement with the synthetic truth in the melting season (Fig. 3a). During the melting season (mid-April 2015 onwards), the AEM of the assimilation algorithms reaches a peak, coinciding with an absence of observations. In comparison, the open-loop AEM is smaller in the first part of the melting season, but the spread is 3 times larger, making it almost uninformative. For several analyses (21 November 2014 and 30 December 2014, for example) the ensemble spread without inflation drops to 0, while its AEM strongly increases compared to the open loop, suggesting that it is prone to degeneracy.

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.


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.


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.


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).


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.


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 Discussion

In this section, we discuss the performance of CrocO PF algorithms using the assimilation of HS and consider the potential of the assimilation of reflectance in view of assimilating real data.

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= 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 Margulis2006). 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 Margulis2007; Mary et al.2013), whereas the required accuracy for visible reflectance retrievals to remain informative on snowpack light-absorbing particles content is high (Warren2013), 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 Bocquet2018).

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 (Evensen2003; Buehner2005).

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


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.

6 Conclusions

In this study, we introduced CrocO, a new ensemble data assimilation system able to reduce the errors of a spatialised snowpack model in locations that are not observed. The ensemble is built by a combination of meteorological and multi-physical ensembles to represent modelling uncertainties. A particle filter assimilates observations of HS and reflectance. We developed two variants of the PF using inflation or k localisation in order to spread the information from partial observations of the system, without degeneracy of the PF. In the framework of synthetic experiments, we have shown in particular the following:

  1. these variants are able to ingest numerous observations without degeneracy;

  2. an efficient spreading of the observational information towards the unobserved areas is achieved with the global and klocal approaches; and

  3. reflectance assimilation leads to an overall 20 % improvement in CRPS and 60 % in reliability.

We suggest that this approach could be used in any spatialised framework to assimilate sparse observations from e.g. networks of in situ snowpack observations. Beyond the snowpack modelling community, the inflation and k-localisation strategies could help address the problem of partially observed systems. This work is also a first step towards the operational assimilation of reflectance in a semi-distributed context. To reach that goal, biases of reflectance retrievals should be studied and observation error structures duly quantified. Snow cover fraction would be a good companion variable to jointly assimilate with reflectances, requiring the use of an appropriate observation operator. Extending the simulation domain to several massifs would allow the exchange of information between neighbouring massifs with the klocal algorithm.

Appendix A: Stochastic perturbations of the forcings

The stochastic perturbation procedure of the forcings is introduced in Sect. 2.2.2 and is identical to Charrois et al. (2016) for the meteorological parameters and Cluzet et al. (2020a) for the light-absorbing particle (LAP) fluxes. For a given date and forcing variable, perturbation values are the same for all the points in space (no spatial autocorrelation is considered), as SAFRAN semi-distributed massifs have a limited spatial extent (about 1000 km2). Precipitation, incoming radiation, wind speed and air temperature from SAFRAN are perturbed with temporally autocorrelated stochastic parameters. The precipitation, incoming shortwave radiation, and wind speed are perturbed with multiplicative noise. Longwave radiation and air temperature are perturbed with additive noise.

For meteorological variables, the perturbation vector V is built as follows:

(A1) V ( t ) = ϕ V ( t - 1 ) + ε ( t ) ,

where ϕ=e-dt/τ, with dt the forcing time step, τ the decorrelation time (h) and ε a normal law of mean 0 and variance σ2(1−ϕ2). Parameter values for each variable are described in Table A1. The significantly high autocorrelation time of precipitation, 1500 h, was tuned to roughly adjust the ensemble spread to the observed intra-massif variability of yearly accumulated precipitation. Note that the precipitation phase is adjusted with the perturbed air temperature to ensure physical consistency. Further details on the procedure can be found in Charrois et al. (2016).

Regarding LAP fluxes, dry and wet black carbon and mineral dust deposition fluxes from MOCAGE are perturbed with a random factor which is constant throughout the year. Each member has a single multiplicative factor following a log-normal law of mean μ and variance σ (see Table A2). The mean of black carbon random perturbations was adjusted based on comparisons between simulations and field observations at Col du Lautaret, a mountain pass within the considered SAFRAN massif.

Table A1Perturbation parameters for the meteorological variables.

Download Print Version | Download XLSX

Table A2Perturbation parameters for the LAP fluxes.

Download Print Version | Download XLSX

Appendix B: Complements on the implementation

B1 Technical implementation and code performance

CrocO is implemented within the Météo-France HPC (high-performance computing) environment, enabling us to fully parallelise the ensemble (one core per member) and bridge the gap with operational applications (Lafaysse et al.2013; Morin et al.2020). This implementation is strongly parallel. As an example, the execution time of a 1-year assimilation run of 187 model points with 160 members on four nodes of 40 cores each lasts only 2 h. The PF is a lightweight algorithm, and most of the computational burden is due to the propagation of the ensemble and input/output (I/O). Also note that no significant difference in execution time can be noted between the different PF algorithms.

B2 PF sample reordering

As mentioned in Sect. 2.3, a reordering step was implemented after the PF resampling from Kitagawa (1996) for practical reasons.

  • (3) From s, build s̃ such that all elements of the unique values of s lie in the position given by their value. Example with 16 particles:


Indeed, I/O represents a bottleneck in the PF. When building the analysis Xa, the background Xb is already loaded in memory. Since Xa is just a reordering of Xb columns based on s, a reordering of s avoids building a copy of Xb. This way, Xa is built by an online modification of Xb using two pointers. Reordering is a growing consideration in the PF community (Farchi and Bocquet2018).

Code availability

The Crocus snowpack model (including all physical options of the ESCROC system) and the particle filter algorithm are developed in the framework of the open-source SURFEX project. The source files of SURFEX code are provided at (Cluzet et al.2020b) to guarantee the permanent reproducibility of results. However, we recommend that potential future users and developers access the code from its Git repository (, last access: 15 April 2020) to benefit from all tools of code management (history management, bug fixes, documentation, interface for technical support, etc.). This requires a quick registration, and the procedure is described at (last access: 4 January 2021). The version used in this work is tagged as CrocO_v1.0.

Python software called CrocO_toolbox was specifically developed in order to pre-process, post-process and launch CrocO experiments. It is available on GitHub (, release v1.0 of the master branch, last access: 4 May 2020) along with documentation.

The article version of CrocO_toolbox is archived at (Cluzet2020). This software strongly relies on two external Python projects ensuring file management between the different steps of a simulation and the interface with the Météo-France HPC system (including parallelisation and data storage): snowtools and vortex. Their sources are available at (Cluzet et al.2020b) (same archive as SURFEX) to guarantee the permanent reproducibility of results. However, as for the SURFEX project and for the same reasons, it is recommended to access snowtools code from its Git repository (, last access: 4 May 2020). The version used in this work is also tagged as CrocO_v1.0. The vortex project gathers all environment-specific codes of Météo-France modelling systems relative to its HPC system. For this project, only the sources specific to this article's simulations are provided. Common object inheritance is based on vortex version 1.6.1. The version used in this work is also tagged as CrocO_v1.0 in the vortex Git repository.

Because these software programmes could not be applied outside the Météo-France HPC environment, CrocO Python software offers the possibility to run CrocO simulations locally. This functionality was not used here due to the high numerical cost of our simulations, which required the use of the Météo-France HPC environment.

Data availability

Input and output data necessary to reproduce the simulations and figures in this paper are provided at (Cluzet et al.2020c). This archive includes the SAFRAN reanalyses (also available at, Vernay et al.2021), MOCAGE forcings, namelists, configuration files and spin-up files necessary to reproduce the simulations. Raw model outputs can be provided on request, but since they can be up to 500+ GB, only post-processed simulation outputs are provided in this archive, along with scores and scripts to reproduce the figures in the paper.

Author contributions

BC wrote the paper. BC, ML and MD designed the study, and BC developed the code with help from ML, LFM and CA. BC, MD, ML and EC designed the PF variants. All authors contributed to results analysis and discussion.

Competing interests

The authors declare that they have no conflict of interest.


The authors are grateful to François Tuzet, Jesus Revuelto, Rafife Nheili, César Deschamps-Berger, Stéphanie Faroux (SODA) and Matthieu Vernay for their valuable help in collecting input data and code implementation. They would also like to thank Fanny Larue, Joseph Bellier, Pierre De Mey and Guilhem Candille for helpful discussions on the PF inflation and CRPS decomposition.

Review statement

This paper was edited by Richard Mills and reviewed by Kristoffer Aalstad and one anonymous referee.


Aalstad, K., Westermann, S., Schuler, T. V., Boike, J., and Bertino, L.: Ensemble-based assimilation of fractional snow-covered area satellite retrievals to estimate the snow distribution at Arctic sites, The Cryosphere, 12, 247–270,, 2018. a, b, c

Aalstad, K., Westermann, S., and Bertino, L.: Evaluating satellite retrieved fractional snow-covered area at a high-Arctic site using terrestrial photography, Remote Sens. Environ., 239, 111618,, 2020. a, b

Albergel, C., Munier, S., Leroux, D. J., Dewaele, H., Fairbairn, D., Barbu, A. L., Gelati, E., Dorigo, W., Faroux, S., Meurey, C., Le Moigne, P., Decharme, B., Mahfouf, J.-F., and Calvet, J.-C.: Sequential assimilation of satellite-derived vegetation and soil moisture products using SURFEX_v8.0: LDAS-Monde assessment over the Euro-Mediterranean area, Geosci. Model Dev., 10, 3889–3912,, 2017. a

Alonso-González, E., Gutmann, E., Aalstad, K., Fayad, A., and Gascoin, S.: Snowpack dynamics in the Lebanese mountains from quasi-dynamically downscaled ERA5 reanalysis updated by assimilating remotely-sensed fractional snow-covered area, Hydrol. Earth Syst. Sci. Discuss.,, in review, 2020. a

Baba, M., Gascoin, S., and Hanich, L.: Assimilation of Sentinel-2 Data into a Snowpack Model in the High Atlas of Morocco, Remote Sens.-Basel, 10, 1982,, 2018. a

Bengtsson, T., Bickel, P., and Li, B.: Curse-of-dimensionality revisited: Collapse of the particle filter in very large scale systems, in: Probability and statistics: Essays in honor of David A. Freedman, Institute of Mathematical Statistics, 2, 316–334,, 2008. a

Buehner, M.: Ensemble-derived stationary and flow-dependent background-error covariances: Evaluation in a quasi-operational NWP setting, Q. J. Roy. Meteor. Soc., 131, 1013–1043, 2005. a

Candille, G., Brankart, J.-M., and Brasseur, P.: Assessment of an ensemble system that assimilates Jason-1/Envisat altimeter data in a probabilistic model of the North Atlantic ocean circulation, Ocean Sci., 11, 425–438,, 2015. a

Cantet, P., Boucher, M., Lachance-Coutier, S., Turcotte, R., and Fortin, V.: Using a particle filter to estimate the spatial distribution of the snowpack water equivalent, J. Hydrometeorol., 20, 577–594, 2019. a, b, c

Carmagnola, C. M., Morin, S., Lafaysse, M., Domine, F., Lesaffre, B., Lejeune, Y., Picard, G., and Arnaud, L.: Implementation and evaluation of prognostic representations of the optical diameter of snow in the SURFEX/ISBA-Crocus detailed snowpack model, The Cryosphere, 8, 417–437,, 2014. a

Charrois, L., Cosme, E., Dumont, M., Lafaysse, M., Morin, S., Libois, Q., and Picard, G.: On the assimilation of optical reflectances and snow depth observations into a detailed snowpack model, The Cryosphere, 10, 1021–1038,, 2016. a, b, c, d, e, f, g

Clark, M. P., Nijssen, B., Lundquist, J. D., Kavetski, D., Rupp, D. E., Woods, R. A., Freer, J. E., Gutmann, E. D., Wood, A. W., Brekke, L. D., Arnold, J. R., Gochis, D. J., and Rasmussen, R. M.: A unified approach for process-based hydrologic modeling: 1. Modeling concept, Water Resour. Res., 51, 2498–2514, 2015. a

Cluzet, B.: bertrandcz/CrocO_toolbox: CrocO_toolbox: GMD mauscript version (Version 1.0), Zenodo,, 2020. a

Cluzet, B., Revuelto, J., Lafaysse, M., Tuzet, F., Cosme, E., Picard, G., Arnaud, L., and Dumont, M.: Towards the assimilation of satellite reflectance into semi-distributed ensemble snowpack simulations, Cold Reg. Sci. Technol., 170, 102918,, 2020a. a, b, c, d, e, f, g, h, i, j, k, l, m

Cluzet, B., Lafaysse, M., Cosme, E., Albergel, C., Meunier, L.-F., and Dumont, M.: CrocO: model source code and external librairies (Version v1.0), Zenodo,, 2020b. a, b

Cluzet, B., Lafaysse, M., Cosme, E., Albergel, C., Meunier, L.-F., and Dumont, M.: CrocO: manuscript Dataset (Version v1.0) [Data set], Zenodo,, 2020c. a

Davaze, L., Rabatel, A., Arnaud, Y., Sirguey, P., Six, D., Letreguilly, A., and Dumont, M.: Monitoring glacier albedo as a proxy to derive summer and annual surface mass balances from optical remote-sensing data, The Cryosphere, 12, 271–286,, 2018. a

De Lannoy, G. J., Reichle, R. H., Arsenault, K. R., Houser, P. R., Kumar, S., Verhoest, N. E., and Pauwels, V. R.: Multiscale assimilation of Advanced Microwave Scanning Radiometer–EOS snow water equivalent and Moderate Resolution Imaging Spectroradiometer snow cover fraction observations in northern Colorado, Water Resour. Res., 48, W01522,, 2012. a, b

de Magalhães, N., Evangelista, H., Condom, T., Rabatel, A., and Ginot, P.: Amazonian Biomass Burning Enhances Tropical Andean Glaciers Melting, Sci. Rep.-UK, 9, 1–12, 2019. a

Dechant, C. and Moradkhani, H.: Radiance data assimilation for operational snow and streamflow forecasting, Adv. Water Resour., 34, 351–364,, 2011. a

Deschamps-Berger, C., Gascoin, S., Berthier, E., Deems, J., Gutmann, E., Dehecq, A., Shean, D., and Dumont, M.: Snow depth mapping from stereo satellite imagery in mountainous terrain: evaluation using airborne laser-scanning data, The Cryosphere, 14, 2925–2940,, 2020. a

Dozier, J., Green, R. O., Nolin, A. W., and Painter, T. H.: Interpretation of snow properties from imaging spectrometry, Remote Sens. Environ., 113, S25–S37, 2009. a, b

Durand, M. and Margulis, S. A.: Feasibility test of multifrequency radiometric data assimilation to estimate snow water equivalent, J. Hydrometeorol., 7, 443–457, 2006. a, b

Durand, M. and Margulis, S. A.: Correcting first-order errors in snow water equivalent estimates using a multifrequency, multiscale radiometric data assimilation scheme, J. Geophys. Res.-Atmos., 112, D13121,, 2007. a, b

Durand, Y., Brun, E., Mérindol, L., Guyomarc'h, G., Lesaffre, B., and Martin, E.: A meteorological estimation of relevant parameters for snow models, Ann. Glaciol., 18, 65–71, 1993. a, b, c

Eberhard, L. A., Sirguey, P., Miller, A., Marty, M., Schindler, K., Stoffel, A., and Bühler, Y.: Intercomparison of photogrammetric platforms for spatially continuous snow depth mapping, The Cryosphere Discuss.,, in review, 2020. a

Essery, R., Morin, S., Lejeune, Y., and Bauduin-Ménard, C.: A comparison of 1701 snow models using observations from an alpine site, Adv. Water Res., 55, 131–148,, 2013. a, b

Evensen, G.: The ensemble Kalman filter: Theoretical formulation and practical implementation, Ocean Dynam., 53, 343–367, 2003. a, b

Farchi, A. and Bocquet, M.: Review article: Comparison of local particle filters and new implementations, Nonlin. Processes Geophys., 25, 765–807,, 2018. a, b, c, d

Fiddes, J., Aalstad, K., and Westermann, S.: Hyper-resolution ensemble-based snow reanalysis in mountain regions using clustering, Hydrol. Earth Syst. Sci., 23, 4717–4736,, 2019. a

Gordon, N. J., Salmond, D. J., and Smith, A. F.: Novel approach to nonlinear/non-Gaussian Bayesian state estimation, IEE Proc.-F, 140, 107–113, 1993. a

Hamill, T. M., Whitaker, J. S., and Snyder, C.: Distance-dependent filtering of background error covariance estimates in an ensemble Kalman filter, Mon. Weather Rev., 129, 2776–2790, 2001. a

Houtekamer, P. L. and Mitchell, H. L.: Data assimilation using an ensemble Kalman filter technique, Mon. Weather Rev., 126, 796–811, 1998. a

Josse, B., Simon, P., and Peuch, V.-H.: Radon global simulations with the multiscale chemistry and transport model MOCAGE, Tellus B, 56, 339–356, 2004. a

Kitagawa, G.: Monte Carlo filter and smoother for non-Gaussian nonlinear state space models, J. Comput. Graph. Stat., 5, 1–25, 1996. a, b

Kokhanovsky, A., Lamare, M., Danne, O., Brockmann, C., Dumont, M., Picard, G., Arnaud, L., Favier, V., Jourdain, B., Le Meur, E., Di Mauro, B., Aoki, T., Niwano, M., Rozanov, V., Korkin, S., Kipfstuhl, S., Freitag, J., Hoerhold, M., Zuhr, A., Vladimirova, D., Faber, A.-K., Steen-Larsen, H. C., Wahl, S., Andersen, J. K., Vandecrux, B., van As, D., Mankoff, K. D., Kern, M., Zege, E., and Box, J. E.: Retrieval of snow properties from the Sentinel-3 Ocean and Land Colour Instrument, Remote Sens.-Basel, 11, 2280,, 2019. a

Krinner, G., Derksen, C., Essery, R., Flanner, M., Hagemann, S., Clark, M., Hall, A., Rott, H., Brutel-Vuilmet, C., Kim, H., Ménard, C. B., Mudryk, L., Thackeray, C., Wang, L., Arduini, G., Balsamo, G., Bartlett, P., Boike, J., Boone, A., Chéruy, F., Colin, J., Cuntz, M., Dai, Y., Decharme, B., Derry, J., Ducharne, A., Dutra, E., Fang, X., Fierz, C., Ghattas, J., Gusev, Y., Haverd, V., Kontu, A., Lafaysse, M., Law, R., Lawrence, D., Li, W., Marke, T., Marks, D., Ménégoz, M., Nasonova, O., Nitta, T., Niwano, M., Pomeroy, J., Raleigh, M. S., Schaedler, G., Semenov, V., Smirnova, T. G., Stacke, T., Strasser, U., Svenson, S., Turkov, D., Wang, T., Wever, N., Yuan, H., Zhou, W., and Zhu, D.: ESM-SnowMIP: assessing snow models and quantifying snow-related climate feedbacks, Geosci. Model Dev., 11, 5027–5049,, 2018. a

Lafaysse, M., Morin, S., Coléou, C., Vernay, M., Serça, D., Besson, F., Willemet, J.-M., Giraud, G., and Durand, Y.: Toward a new chain of models for avalanche hazard forecasting in French mountain ranges, including low altitude mountains, in: Proceedings of the International Snow Science Workshop – Grenoble and Chamonix, 162–166, 2013. a, b

Lafaysse, M., Cluzet, B., Dumont, M., Lejeune, Y., Vionnet, V., and Morin, S.: A multiphysical ensemble system of numerical snow modelling, The Cryosphere, 11, 1173–1198,, 2017. a, b, c, d

Lahoz, B. K. W. and Menard, R.: Data assimilation, Springer, Berlin, Heidelberg,, 2010. a

Lamare, M., Dumont, M., Picard, G., Larue, F., Tuzet, F., Delcourt, C., and Arnaud, L.: Simulating optical top-of-atmosphere radiance satellite images over snow-covered rugged terrain, The Cryosphere, 14, 3995–4020,, 2020. a, b

Largeron, C., Dumont, M., Morin, S., Boone, A., Lafaysse, M., Metref, S., Cosme, E., Jonas, T., Winstral, A., and Margulis, S. A.: Towards snow cover estimation in mountainous areas using modern data assimilation methods: A review, Front. Earth Sci., 8, 325,, 2020. a, b, c, d

Larue, F., Royer, A., De Sève, D., Roy, A., and Cosme, E.: Assimilation of passive microwave AMSR-2 satellite observations in a snowpack evolution model over northeastern Canada, Hydrol. Earth Syst. Sci., 22, 5711–5734,, 2018. a, b, c, d, e, f, g

Libois, Q., Picard, G., France, J. L., Arnaud, L., Dumont, M., Carmagnola, C. M., and King, M. D.: Influence of grain shape on light penetration in snow, The Cryosphere, 7, 1803–1818,, 2013. a

Libois, Q., Picard, G., Arnaud, L., Dumont, M., Lafaysse, M., Morin, S., and Lefebvre, E.: Summertime evolution of snow specific surface area close to the surface on the Antarctic Plateau, The Cryosphere, 9, 2383–2398,, 2015. a

Liu, J. S. and Chen, R.: Blind deconvolution via sequential imputations, J. Am. Stat. Assoc., 90, 567–576, 1995. a

Magnusson, J., Gustafsson, D., Hüsler, F., and Jonas, T.: Assimilation of point SWE data into a distributed snow cover model comparing two contrasting methods, Water Resour. Res., 50, 7816–7835, 2014. a, b, c

Magnusson, J., Winstral, A., Stordal, A. S., Essery, R., and Jonas, T.: Improving physically based snow simulations by assimilating snow depths using the particle filter, Water Resour. Res., 53, 1125–1143, 2017. a

Mankin, J. S., Viviroli, D., Singh, D., Hoekstra, A. Y., and Diffenbaugh, N. S.: The potential for snow to supply human water demand in the present and future, Environ. Res. Lett., 10, 114016,, 2015. a

Margulis, A., Cortés, G., Girotto, M., and Durand, M.: A Landsat-Era Sierra Nevada Snow Reanalysis (1985–2015), J. Hydrometeorol., 17, 1203–1221,, 2016. a

Margulis, S. A., Fang, Y., Li, D., Lettenmaier, D. P., and Andreadis, K.: The Utility of Infrequent Snow Depth Images for Deriving Continuous Space-Time Estimates of Seasonal Snow Water Equivalent, Geophys. Res. Lett., 46, 5331–5340, 2019. a, b, c

Marsh, C. B., Pomeroy, J. W., and Wheater, H. S.: The Canadian Hydrological Model (CHM) v1.0: a multi-scale, multi-extent, variable-complexity hydrological model – design and overview, Geosci. Model Dev., 13, 225–247,, 2020. a

Mary, A., Dumont, M., Dedieu, J.-P., Durand, Y., Sirguey, P., Milhem, H., Mestre, O., Negi, H. S., Kokhanovsky, A. A., Lafaysse, M., and Morin, S.: Intercomparison of retrieval algorithms for the specific surface area of snow from near-infrared satellite data in mountainous terrain, and comparison with the output of a semi-distributed snowpack model, The Cryosphere, 7, 741–761,, 2013. a, b

Masson, T., Dumont, M., Mura, M. D., Sirguey, P., Gascoin, S., Dedieu, J.-P., and Chanussot, J.: An Assessment of Existing Methodologies to Retrieve Snow Cover Fraction from MODIS Data, Remote Sens.-Basel, 10, 619,, 2018. a, b

Masson, V., Le Moigne, P., Martin, E., Faroux, S., Alias, A., Alkama, R., Belamari, S., Barbu, A., Boone, A., Bouyssel, F., Brousseau, P., Brun, E., Calvet, J.-C., Carrer, D., Decharme, B., Delire, C., Donier, S., Essaouini, K., Gibelin, A.-L., Giordani, H., Habets, F., Jidane, M., Kerdraon, G., Kourzeneva, E., Lafaysse, M., Lafont, S., Lebeaupin Brossier, C., Lemonsu, A., Mahfouf, J.-F., Marguinaud, P., Mokhtari, M., Morin, S., Pigeon, G., Salgado, R., Seity, Y., Taillefer, F., Tanguy, G., Tulet, P., Vincendon, B., Vionnet, V., and Voldoire, A.: The SURFEXv7.2 land and ocean surface platform for coupled or offline simulation of earth surface variables and fluxes, Geosci. Model Dev., 6, 929–960,, 2013. a

Matheson, J. E. and Winkler, R. L.: Scoring rules for continuous probability distributions, Manage. Sci., 22, 1087–1096, 1976. a

Morin, S., Horton, S., Techel, F., Bavay, M., Coléou, C., Fierz, C., Gobiet, A., Hagenmuller, P., Lafaysse, M., Ližar, M., Mitterer, C., Monti, F., Müller, K., Olefs, M., Snook, J. S., van Herwijnen, A., and Vionnet, V.: Application of physical snowpack models in support of operational avalanche hazard forecasting: A status report on current implementations and prospects for the future, Cold Reg. Sci. Technol., 170, 102910,, 2020. a, b, c, d, e

Nousu, J.-P., Lafaysse, M., Vernay, M., Bellier, J., Evin, G., and Joly, B.: Statistical post-processing of ensemble forecasts of the height of new snow, Nonlin. Processes Geophys., 26, 339–357,, 2019. a

Penny, S. G. and Miyoshi, T.: A local particle filter for high-dimensional geophysical systems, Nonlin. Processes Geophys., 23, 391–405,, 2016. a

Piazzi, G., Thirel, G., Campo, L., and Gabellani, S.: A particle filter scheme for multivariate data assimilation into a point-scale snowpack model in an Alpine environment, The Cryosphere, 12, 2287–2306,, 2018. a, b

Pomeroy, J. W., Fang, X., and Marks, D. G.: The cold rain-on-snow event of June 2013 in the Canadian Rockies–Characteristics and diagnosis, Hydrol. Process., 30, 2899–2914, 2016. a

Poterjoy, J.: A localized particle filter for high-dimensional nonlinear systems, Mon. Weather Rev., 144, 59–76, 2016. a, b

Poterjoy, J. and Anderson, J. L.: Efficient assimilation of simulated observations in a high-dimensional geophysical system using a localized particle filter, Mon. Weather Rev., 144, 2007–2020, 2016. a

Poterjoy, J., Wicker, L., and Buehner, M.: Progress toward the application of a localized particle filter for numerical weather prediction, Mon. Weather Rev., 147, 1107–1126, 2019. a

Press, W. H., Flannery, B. P., Teukolsky, S. A., and Vetterling, W. T.: Numerical recipes in Fortran the art of scientific computing, 2nd Edn., Cambridge University Press, Cambridge, 1992. a

Raleigh, M. S., Lundquist, J. D., and Clark, M. P.: Exploring the impact of forcing error characteristics on physically based snow simulations within a global sensitivity analysis framework, Hydrol. Earth Syst. Sci., 19, 3153–3179,, 2015. a

Reichle, R. H. and Koster, R. D.: Assessing the impact of horizontal error correlations in background fields on soil moisture estimation, J. Hydrometeorol., 4, 1229–1242, 2003. a

Revuelto, J., Lecourt, G., Lafaysse, M., Zin, I., Charrois, L., Vionnet, V., Dumont, M., Rabatel, A., Six, D., Condom, T., Morin, S., Viani, A., and Sirguey, P.: Multi-Criteria Evaluation of Snowpack Simulations in Complex Alpine Terrain Using Satellite and In Situ Observations, Remote Sens.-Basel, 10, 1171,, 2018. a

Revuelto, J., Cluzet, B., Duran, N., Fructus, M., Lafaysse, M., Cosme, E., and Dumont, M.: Assimilation of surface reflectance in snow simulations; impact on bulk snow variables, in preparation, 2021. a

Richter, B., van Herwijnen, A., Rotach, M. W., and Schweizer, J.: Sensitivity of modeled snow stability data to meteorological input uncertainty, Nat. Hazards Earth Syst. Sci., 20, 2873–2888,, 2020. a

Sabatier, T., Largeron, Y., Paci, A., Lac, C., Rodier, Q., Canut, G., and Masson, V.: Semi-idealized simulations of wintertime flows and pollutant transport in an alpine valley: Passive tracer tracking (Part II), Q. J. Roy. Meteor. Soc., 146, 827–845,, 2020. a

Shaw, T. E., Gascoin, S., Mendoza, P. A., Pellicciotti, F., and McPhee, J.: Snow depth patterns in a high mountain Andean catchment from satellite optical tri-stereoscopic remote sensing, Water Resour. Res., 56, e2019WR024880,, 2019. a

Sirguey, P., Mathieu, R., and Arnaud, Y.: Subpixel monitoring of the seasonal snow cover with MODIS at 250 m spatial resolution in the Southern Alps of New Zealand: methodology and accuracy assessment, Remote Sens. Environ., 113, 160–181,, 2009. a

Slater, A. G. and Clark, M. P.: Snow data assimilation via an ensemble Kalman filter, J. Hydrometeorol., 7, 478–493, 2006. a

Snyder, C., Bengtsson, T., Bickel, P., and Anderson, J.: Obstacles to high-dimensional particle filtering, Mon. Weather Rev., 136, 4629–4640, 2008. a

Stigter, E. E., Wanders, N., Saloranta, T. M., Shea, J. M., Bierkens, M. F. P., and Immerzeel, W. W.: Assimilation of snow cover and snow depth into a snow model to estimate snow water equivalent and snowmelt runoff in a Himalayan catchment, The Cryosphere, 11, 1647–1664,, 2017. a

Thirel, G., Salamon, P., Burek, P., and Kalas, M.: Assimilation of MODIS snow cover area data in a distributed hydrological model using the particle filter, Remote Sens.-Basel, 5, 5825–5850, 2013. a

Tuzet, F., Dumont, M., Lafaysse, M., Picard, G., Arnaud, L., Voisin, D., Lejeune, Y., Charrois, L., Nabat, P., and Morin, S.: A multilayer physically based snowpack model simulating direct and indirect radiative impacts of light-absorbing impurities in snow, The Cryosphere, 11, 2633–2653,, 2017. a

Van Leeuwen, P. J.: Particle filtering in geophysical systems, Mon. Weather Rev., 137, 4089–4114, 2009. a

Van Leeuwen, P. J., Künsch, H. R., Nerger, L., Potthast, R., and Reich, S.: Particle filters for high-dimensional geoscience applications: A review, Q. J. Roy. Meteor. Soc., 145, 2335–2365, 2019. a

Vernay, M., Lafaysse, M., Merindol, L., Giraud, G., and Morin, S.: Ensemble Forecasting of snowpack conditions and avalanche hazard, Cold Reg. Sci. Technol., 120, 251–262,, 2015. a

Vernay, M., Lafaysse, M., Mérindol, L., Hagenmuller, P., Verfaillie, D., Morin., S., and Monteiro, D.: A 61-years meteorological and snow conditions re-analysis over the French mountainous areas (1958–2019) [Data set], AERIS,, 2021. a

Veyssière, G., Karbou, F., Morin, S., Lafaysse, M., and Vionnet, V.: Evaluation of Sub-Kilometric Numerical Simulations of C-Band Radar Backscatter over the French Alps against Sentinel-1 Observations, Remote Sens.-Basel, 11, 8,, 2019. a

Vionnet, V., Brun, E., Morin, S., Boone, A., Faroux, S., Le Moigne, P., Martin, E., and Willemet, J.-M.: The detailed snowpack scheme Crocus and its implementation in SURFEX v7.2, Geosci. Model Dev., 5, 773–791,, 2012. a

Vionnet, V., Marsh, C. B., Menounos, B., Gascoin, S., Wayand, N. E., Shea, J., Mukherjee, K., and Pomeroy, J. W.: Multi-scale snowdrift-permitting modelling of mountain snowpack, The Cryosphere Discuss.,, in review, 2020. a

Warren, S. G.: Can black carbon in snow be detected by remote sensing?, J. Geophys. Res., 118, 779–786,, 2013. a

Wayand, N. E., Marsh, C. B., Shea, J. M., and Pomeroy, J. W.: Globally scalable alpine snow metrics, Remote Sens. Environ., 213, 61–72, 2018. a

Winstral, A., Magnusson, J., Schirmer, M., and Jonas, T.: The Bias-Detecting Ensemble: A New and Efficient Technique for Dynamically Incorporating Observations Into Physics-Based, Multilayer Snow Models, Water Resour. Res., 55, 613–631, 2019.  a, b, c, d, e

Wright, P., Bergin, M., Dibb, J., Lefer, B., Domine, F., Carman, T., Carmagnola, C., Dumont, M., Courville, Z., Schaaf, C., and Wang, Z.: Comparing MODIS daily snow albedo to spectral albedo field measurements in Central Greenland, Remote Sens. Environ., 140, 118–129, 2014. a

Würzer, S., Jonas, T., Wever, N., and Lehning, M.: Influence of initial snowpack properties on runoff formation during rain-on-snow events, J. Hydrometeorol., 17, 1801–1815, 2016. a

Short summary
In the mountains, the combination of large model error and observation sparseness is a challenge for data assimilation. Here, we develop two variants of the particle filter (PF) in order to propagate the information content of observations into unobserved areas. By adjusting observation errors or exploiting background correlation patterns, we demonstrate the potential for partial observations of snow depth and surface reflectance to improve model accuracy with the PF in an idealised setting.