Ensemble-based snow depth data assimilation for a multi-layer snow scheme over the European Arctic

Bakketun, Åsmund; Blyverket, Jostein; Müller, Malte

doi:https://doi.org/10.5194/egusphere-2025-1693

Preprints

https://doi.org/10.5194/egusphere-2025-1693

Preprints

09 May 2025

| 09 May 2025

Ensemble-based snow depth data assimilation for a multi-layer snow scheme over the European Arctic

Åsmund Bakketun, Jostein Blyverket, and Malte Müller

Abstract. Reliable estimates of Earth system conditions are important for weather forecasting, hydrological modelling and their downstream applications. Both real-time prediction systems and historical reanalyses use a combination of observations and physical laws embedded in numerical models to generate gapless and accurate estimates of weather, climate and hydrological conditions. Data assimilation systems merge information from the two sources in an objective way, accounting for their respective uncertainties. In this work we present a regional reanalysis system, focusing on the land surface component. The system uses a multi-layer snow model together with the ensemble-based Local Ensemble Transform Kalman Filter (LETKF) data assimilation scheme. The system is run for a 4 year period over the European Arctic, assimilating in situ snow depth observations. Evaluation of the snow depth estimates showed reduced errors compared to existing products and positive impact of the data assimilation over the domain. Furthermore, a significant difference in total accumulated snow water was seen over the domain, implying a potential impact on downstream hydrological applications. The ensemble correlations between the total snow depth and the relatively large control vector indicated that the ensemble was able to represent snow compaction processes. The LETKF is thus able to account for these processes, which are often neglected in snow depth data assimilation. The system presented in this study allows for future extensions, including other types of observations and analyses beyond snow variables.

Received: 09 Apr 2025 – Discussion started: 09 May 2025

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors. Views expressed in the text are those of the authors and do not necessarily reflect the views of the publisher.

Download & links

Åsmund Bakketun, Jostein Blyverket, and Malte Müller

Status: final response (author comments only)

Subscribe to comment alert

CC1: 'Comment on egusphere-2025-1693', Nima Zafarmomen, 15 May 2025

The paper addresses an important gap: bringing a flow-dependent ensemble Kalman framework to a multi-layer snow scheme for a high‐latitude regional reanalysis. The topic is timely and the modelling chain (SURFEX–ISBA explicit snow + LETKF, driven by CARRA) is potentially valuable for cryospheric and hydrologic communities.

Forcing and domain
The manuscript states that CARRA forcing is “interpolated to the model grid” but omits grid spacing for both driver and land model. Domain limits are described in text but not in coordinates; include bounding box and spatial resolution.
Ensemble generation
Only perturbing atmospheric forcing inevitably under-represents uncertainty in snow compaction, albedo metamorphism, and interception. You should quantify how ensemble spread compares to innovation statistics (e.g., spread-skill ratio) to demonstrate sufficiency of the perturbation strategy. If spread is systematically low, adding multiplicative inflation alone is insufficient; process perturbations or parameter perturbations may be needed.
The “remapping” approach for precipitation displacement is innovative, yet Appendix A lacks diagnostic evidence that the scheme produces realistic error structures. Provide at minimum a variogram or visual comparison between perturbed and reference precipitation fields.
Ensemble size
Ten members is very small for a 36-variable profile. You report that 20 members offered “no considerable degradation”, but give no metrics. Include a sensitivity figure (e.g., CRPS vs. ensemble size) to justify the final choice.
Increment analysis (Sect. 3.1)
Figure 4 shows domain-mean increments of several millimetres water equivalent per day—this is large. Provide histograms or spatial standard deviations to make clear whether these increments are isolated to specific subregions or pervasive. Without that context, the reader cannot judge if the LETKF is “adding missing precipitation” or merely compensating biased forcing.
Skill metrics
– CRPS and MAE are reported, but no sampling uncertainty is provided. Bootstrap confidence intervals would show whether the apparent improvements are statistically robust.

– Station splits (OBS-ONLY-Pv1, etc.) prove useful, yet the sample sizes differ dramatically. Present RMSE normalised by climatological variance to avoid overweighting dense station clusters.
SWE validation
Only six pillow sites are available, but you can still compute Kling-Gupta or Nash–Sutcliffe across time to give hydrologists a sense of hydro-logical skill. Also, the negative bias at one degraded site coincides with orographic precipitation maxima; examine whether forcing under-catch is the root cause.

I strongly recommend that the authors expand their discussion, as data assimilation is not only applicable to snow schemes but is also widely used in other areas such as streamflow prediction. I also recommend citing the following papers: Optimising ensemble streamflow predictions with bias correction and data assimilation techniques; Assimilation of Sentinel-Based Leaf Area Index for Modeling Surface-Ground Water Interactions in Irrigation Districts

Citation: https://doi.org/10.5194/egusphere-2025-1693-CC1
RC1: 'Comment on egusphere-2025-1693', Matthieu Lafaysse, 24 Jul 2025

General comments
Bakketun et al. present the results of a new land surface reanalysis refining the snow scheme compared to the forcing CARRA reanalysis and focusing on the assimilation of in-situ snow depths with a new ensemble assimilation framework derivated from the well-known Ensemble Kalman Filter. The introduction presents clearly the context and objectives of this work and the production of an Arctic land surface reanalysis with improved snow modelling and snow data assimilation has certainly a great potential for various related scientific fields.
However, in general, the methods are not described with a sufficient accuracy level to allow their reproductibility (parameters used for perturbation, interpolation, observation errors, inflation, etc.). I believe the manuscript shoud be revised in that spirit. As a result, it is sometimes difficult to evaluate the relevance of the proposed methodology, see in particular my detailed comments relative to L149-151, L153-163 ; L214-217 ; L250. My understanding after reading the Methods section was that the main advantage of the EnKF was lost by the local assumption for application while the estimation of the snow model states are much more challenging to obtain than in the context of particle filters or particle batch smoothers. More explanations are necessary to understand this key methodological choices, because although I may have misunderstood something, at this stage, my feeling was that the chosen methodology might not be the most appropriate for this application. When reaching the results Section, I understand that LETKF is actually more complex than described, probably including neighbour pixels in the assimilation process but the readers should have fully understood that after reading the Methods even for readers not aware of all the variants of ensemble data assimilation algorithms.
Another critical aspect of this manuscript is the fact that the independence of observations is not considered as a mandatory property for evaluating the added value of data assimilation. Although some datasets are more or less independent, the results do not focus on this part. More critically, when considering only this dataset, I am afraid that the results are not really convincing to demonstrate the added value of their land surface reanalysis compared to the CARRA reference. As a result, I stay with a mitigated opinion where I see a really nice attempt to produce snow reanalyses with state-of-the-art methods, but leading to not really convincing results, questioning whether the chosen assimilation methodology was really appropriate, and with an unclear understanding of the behaviour of the assimilation algorithm due to unsufficient methodology description and maybe lack of examples on specific cases.
Therefore, I would recommend a major review for this manuscript to give the opportunity to authors to improve their description of methods, refine their evaluation process focusing on independent dataset, and maybe find new results emphasing better the added value of their methodology compared to their benchmark.
This being said, I am not fully convinced that assimilating only in situ snow depth observations has really a significant potential in mountainous areas with sparse observations with low spatial representativity, and I am really wondering why the authors do not prefer to start by assimilating snow cover fractions already assimilated in CARRA.
Detailed comments
L52 Note that this reference presents twin experiments (with synthetic observations). To the best of my knowledge, real observations of surface temperature have not yet been assimilated in a snow cover modelling system. I suggest to be more accurate on that point.
L54-56 The authors classify the DA algorithms according to their management of uncertainties but the algorithms also differ in their ability to preserve the multivariate consistency of model variables. It would also be useful to emphasize this point which is also a strength for instance of particle filters.
L59-60 Variable transformation is also a common practice to fit with the Gaussian assumption.
L61 Unclear at this step of the introduction which model refers to « the multi-layer snow model ».
L85-86 How was estimated this lookup table and is it variable in space ?
L87 and elsewhere : The display of units should be improved (space between kg and m^-2)
L88 I guess there is a typo in the unit (kg m^-2 and not kg m^-3). Can you explain the philosophy behind excluding observations « if the model exceeds 25 kg m^-2 » ?
L90-93 Does this mean that pseudo-observations and in-situ observations are considered to have the same observation error ?
L93-94 Are there other variables to describe the snow cover state in the model ? Are they kept at the same value during this process ?
L94-95 Can you explain the reason behind that ? I guess this is because surface temperature can not be updated ?
L96 « For modelling the land surface » → I understand that the idea is to produce a complementary offline simulation with an improved surface model compared to the CARRA coupled system, but this could be more explicitely explained in this transition.
L101 « the reduced heat flux between soil and snow » As a coauthor of Monteiro et al, 2024, I know what this is about, but I can imagine than most other readers would need more explanations.
L112 To encourage the reproductibility of the results, could you be more accurate on the spatial and temporal interpolation methods used for the different forcing variables ?
L129-130 Snow albedo and snow age do also need to be updated
L128-132 Direct insertion of snow depth has been applied in some studies with multilayer snow models by preserving the ratio between layers and other snow properties. The assumptions are strong in that case, but it could be mentioned to emphasize the added value of more advanced ensemble approaches.
L149-151 If EnKF is applied independently at each grid point, how can the increments be propagated in space ? This is rather counter-intuitive in the context of assimilation of in-situ snow depths in a spatialized system. Please provide explanations.
L153-163 Could you mention if vector x includes the 36 prognostic variables described line 129 ? If yes, how do you manage snow age and snow albedo ? What is the implication of computing error covariances with such a large vector ? What guarantees that the resulting analysed vector fills the discretization rules of the snow scheme ? Considering that for a given snow depth, the vertical discretization of ISBA-ES is fixed, I am not sure what is the interest to consider the 12 layer snow depths as independent variables in the state vector. I imagine it helps constraining unobserved density and heat contents. But could you justify this choice and provide more details to help the understanding beyond the mathematical formalism ?
Then, does vector x only include the prognostic variables of one single point ? If yes, how do you propagate information in space ? I see further from the results, that increments are spatialized but I do not see from this formalism how this happens.
L165 This assumption would mean that observations are as representative as the open patch than as the forest patch. However, in a large majority, observations are operated on open areas. Could you comment on that ?
L182-183 This is true but monovariate perturbations may produce unrealistic meteorological conditions with physical inconsistencies between meteorological variables, potentially leading to unrealistic model states.
L185 As mentioned before, the length of the control vector might be reduced thanks to the bijectivity between total snow depth and layer snow depths.
L187-189 The description of perturbations refers to several references which is not accurate enough for reproductibility. Please provide a table with perturbation methods and parameters (variance, auto-correlation time constants, etc.) for each variable.
L190-203 This is an interesting idea but even in Appendix A, the perturbation parameters to prescribe are unclear and not provided (I guess at least some horizontal distance defining the statistical properties of the random advection ?). Could you think again the description of the method as a basis allowing reproductibility ?
L212 « we initialized zero snow members with the ensemble-mean of members with snow. » It is unclear how is computed the mean. The different members have different layer thicknesses. Do you mean that the properties (density, heat) are averaged for each numerical layer regardless its depth in the snowpack ? Is it physically realistic ? Would not it be sufficient to extract the values of one single member with low snow depth rather than potentially mixing the properties of thin snowpacks with properties of thick snowpacks ?
L214-217 If the algorithm provides unrealistic states and crashes happen, I am wondering if applying thresholds is really the appropriate option compared to improvements in the assimilation formalism. To my mind, the main added value of the EnKF compared to particle filters is lost if it is applied independently at each pixel point. At the local scale and to assimilate only snow depths, particle filters would be best suited algorithms to retrieve all variable states without any possibility to obtain unrealistic model states, and they have already been largely applied within SURFEX even with more complex snow schemes. The added value of the proposed approach is unclear for me.
Table 1 : I do not think a threshold on temperature is necessary as snow heat, snow density and snow mass are sufficient to diagnose snow temperature. Or do you mean the threshold on temperature is applied to constrain the heat content ?
Section 2.6 Could you provide the elevation distribution of the different observation datasets ? Is a threshold applied on elevation as commonly done in NWP ?
I also understand that optical snow cover maps were not used here while they were assimilated in CARRA. This should be more explicit and justified as it could be questionable that an offline land surface reanalysis assimilates less observations than the original coupled reanalysis.
L231 How are initialized the simulations ? Do you apply any spinup to initialize the soil temperature profiles ? What is the assumption for snow depth ? I guess some pixels are affected by permanent snow ?
L242 Do you use independent observation datasets (not used in the assimilation process) to evaluate the system ?
L250 The fact that Figure 3 shows spatial patterns of increments means that something is definitely missing in Section 2.3 to produce the spatial propragation of increments when applying EnKF at the pixel scale. I guess this is associated with « the Gaussian localization function » mentioned L256 but this should be explained in Section 2.3 (providing methodological details as well as control parameters) as after my first reading I thought only one pixel was considered in vector x. The paper should be self-sufficient to understand everything without reading other papers describing LETKF. This raises numerous questions by the way about the localization radius. Is it constant or spatially / temporally variable ? Does it depend on the density of the observation network ? Etc.
L256 « The patterns of these increments are similar to the Gaussian localization function ». Without any knowledge about this localization function, I do not understand this comment.
L259-260 « The stations with mean negative innovations are situated in mountainous areas and have small impact on the analysis. » Why? Do you reduce some localization radius as a function of elevation ? Methodological details are missing to understand.
L261 « domain-averaged ». There are two different domains between Figure 2 and Figure 3. To which domain is applied this average ?
Figure 5 presents ensemble correlations between between layer variables and total snow depth. The results are useful to understand the assimilation process. However, it is not explained how this information is represented at a given depth in the snowpack as the different ensemble members have different snow depths and layer thicknesses. Do you choose a member ? Or compute median thicknesses ? Please provide the details in the text.
L267-270 The correlations between total snow depth and mass of each layer could be commented at the light of the vertical discretization rules of the snow scheme, explaining easily the weak correlation between the total snow depth and the mass of the surface and bottom layers.
L284-291 I would have expected that the effect of assimilation on density would have been described here (comparing densities obtained with CARRA-Land-Pv1 and CTRL), as done before with snow depth.
L292-296 The observations used to compute errors are not described. Are they independent from the dataset used for assimilation ? If not, I think it is unsufficient to evaluate a data assimilation algorithm with the assimilated data themselves. Then, the fact that assimilation deteriorates snow water equivalent is really problematic as the main advantage the authors emphasize in their method is its ability to attribute snow depth errors in both mass and density. Unfortunately, the results suggest that this idea fails. A much more simple algorithm assuming that the simulated snow density is correct may have been a better assumption ?
L300-301 If I understand well, OBS-ONLY-CARRA is the only relatively independent dataset to evaluate the added value of assimilation. So this is of course more challenging, but the description of results in the whole section should be mainly focused on this dataset, as imrproving snow depth at the assimilated observations is definitely not a proof that the system is valid at large scale. The word ‘validation’ is inappropriate in the title of the subsection (prefer evaluation). Unfortunately, the fact that the whole method deteriorates scores compared to CARRA (Fig 8c, value on the right of each column in Tab 3) is rather problematic while the goal of the study was to improve the simulation of snow cover with a dedicated offline reanalysis.
L307 It should be explained that CRPS is identical to MAE in the case of a deterministic forecast.
L333-334 I would like to be optimistic but compared to independent observations (OBS-ONLY-CARRA) not assimilated in the offline reanalysis, the results do not exhibit improvements of snow depths (Fig 8c, Table 3) and exhibit significant a deteriotation of density (Fig 7). I think a more qualified discussion of results is needed allowing a better questioning of the choices for data assimilation.
L343-357 I think the discussion about the link between density and mass could be improved.
« We emphasize that these stations are not collocated with assimilated snow depth observations ». In a spatialized assimilation system, I think the goal is to improve snow cover at large scale, not only at the assimilated observations.
L358 « By exploring the different observation sets used in CARRA and CARRA-Land-Pv1, we evaluate the analysis performance where observations are not available for assimilation ». Actually, only one dataset is independent and the results do not really focus on this one. I think the evaluation dataset used in this paper should be questioned and probably extended. A leave-one-out approach as in Cluzet et al., 2022 could be a valuable approach to consider all the snow depth datasets. Independent observations of satellite snow cover fractions could also be considered for a spatialized evaluation of the added value of data assimilation, especially in a context where such observations were assimilated by CARRA and are no longer assimilated in CARRA-Land-PV1.
L361 It is difficult to have a critical reading of this discussion due to the lacks in the methodology description. For instance, authors say that « the distance based localization used in CARRA-Land-Pv1 cause several stations in OBS-ONLY-CARRA to be unreachable in the analysis. C ». But the localization distance has not been defined neither provided, and the typical observation density is not known. I think a discussion should more rely on quantitative considerations provided to the reader.
L384-391 Authors explain again the theoretical added value of their framework compared to the simple assimilation process of CARRA. However, my feeling is that they have not found the results illustrating how considering variable background uncertainties is able to improve assimilation compared to their reference benchmark. If authors have illustrations on specific cases, as suggested, I would encourage to present the material in the paper to better emphasize the potential added value, if it can not be demonstrated with finalized scores.
L393 implemented
L396-397 Again, this conclusion seems rather optimistic considering only independent evaluation datasets.
L399-401 Again, my feeling is that this theorical advantage is not completely supported by results.

Citation: https://doi.org/10.5194/egusphere-2025-1693-RC1
RC2: 'Comment on egusphere-2025-1693', Anonymous Referee #2, 12 Sep 2025

My comments are in the attached pdf file

Citation: https://doi.org/10.5194/egusphere-2025-1693-RC2

Åsmund Bakketun, Jostein Blyverket, and Malte Müller

Viewed

Total article views: 427 (including HTML, PDF, and XML)

HTML	PDF	XML	Total	BibTeX	EndNote
365	48	14	427	16	35

HTML: 365
PDF: 48
XML: 14
Total: 427
BibTeX: 16
EndNote: 35

Views and downloads (calculated since 09 May 2025)

Month	HTML	PDF	XML	Total
May 2025	79	23	3	105
Jun 2025	35	5	7	47
Jul 2025	28	10	4	42
Aug 2025	105	7	0	112
Sep 2025	118	3	0	121

Cumulative views and downloads (calculated since 09 May 2025)

Month	HTML	PDF	XML	Total
May 2025	79	23	3	105
Jun 2025	35	5	7	47
Jul 2025	28	10	4	42
Aug 2025	105	7	0	112
Sep 2025	118	3	0	121

Viewed (geographical distribution)

Total article views: 414 (including HTML, PDF, and XML) Thereof 414 with geography defined and 0 with unknown origin.

Country	#	Views	%

Latest update: 12 Sep 2025

Short summary

Obtaining accurate estimates of seasonal snow conditions requires a combination of observations and numerical models. We use a model accounting for the vertical structure of the snow, and a data assimilation method representing varying uncertainty of the model in time and space. Compared to existing products, neglecting these considerations, our system produced improved estimates of seasonal snow conditions. Snow mass estimates suggest a potential impact on derived hydrological applications.


Total:	0
HTML:	0
PDF:	0
XML:	0