General Comments:

Overall, this is an excellent paper worthy of publication in GMD. The topic of snow data assimilation is of high scientific importance and providing a unifying framework for implementing such methods is to be commended and should be of high value to the community. The structure of the manuscript and presentation are generally clear. There are many specific comments listed below which aim to improve the communication of the proposed framework and aid users in its implementation and use. In particular, more details on the sample problems would aid in the reproducibility and extension of the work to other problems. In testing the code, it appears that the github repository code works, but the sample provided on zenodo has a bug. More details are provided below.

Specific Comments:

The comments provide herein represent a list of relatively minor additions and/or corrections that would improve the paper.

1. Title: The usage of “Multiscale” in the title does not seem particularly warranted. My expectation based on the title was that the implementation would be flexible enough to model snow at multiple scales (resolutions) and/or assimilate data at multiple scales (resolutions). It is not clear from the presentation whether either is the case. Or is the meaning meant to convey multiple temporal scales? The authors should consider whether the title should be changed for clarity. If there is an aspect of what you are proposing that is indeed “multiscale” you should emphasize that more for the reader’s benefit.

2. Line 85: In mentioning the “posterior mean snow simulation from FSM2” it would be useful to know what variables that contains. Maybe the FSM2 variables could be shown in a Table?

3. Line 132: In mentioning the additive/multiplicative perturbations there is no description as to whether they are perfectly independent or perfectly correlated or something in between. In other words, do all pixels get the same perturbation (i.e. from the same random number) or do they get fully independent perturbations (i.e. each sampled independently). Mentioning here or elsewhere that this neglects spatially-correlated errors/uncertainties would be appropriate. It is mentioned earlier that the model structure is fully independent, but saying whether the perturbations are as well would clarify the setup.

4. Line 413: In describing the time-invariant perturbations it may be worth mentioning what the implications of that are vs. other options (independent in time or correlated in time).

5. Line 654: I think a couple of sentences describing the mechanics of the iterative nature of the method and why it outperforms other methods would be warranted.

6. Lines 696-703: More explanation of what it meant by “inflated observation errors” and how it fits into the method would be helpful to the reader. An elaboration on the note about the “multiple data assimilation approach does not actually violate” would also be helpful. Since this particular method is less standard than others, I assume most readers would benefit from more detail here.

7. Line 741: Can you provide some justification of the choice of four (4) for the number of assimilation cycles?

8. Section 4: I would urge that more consistent (and maybe simpler) language be used throughout to refer to the three experiments being done so that readers can follow more easily. It doesn’t seem that “drone data” or “satellite data” are as relevant to the first two experiments compared to the first being a spatially-distributed (snow depth) data assimilation experiment and the second being a point-scale joint (FSCA+LST) data assimilation experiment. For the benchmark case, it is not clearly defined what experiment is actually being done. Is it snow depth or LST+FSCA assimilation? There is a mention of “single cell” which may imply it is the same setup as the second experiment, but this is not clear as currently presented.

9. Section 4: In an effort to make the sample experiments more reproducible for the readers, I would suggest tabulating any key parameter differences (beyond default values) in the config.py and/or constants.py input files that are specific to each experiment being done. It would also be helpful to connect the individual experiments to the theory provided earlier in the paper, i.e. description of the states, measurement, etc. In particular, if transforms are used with respect to the measurements (as referred to on Lines 359-370), it would be useful to see the form of the those transforms in the experimental setup in Section 4.

10. Section 4: Perhaps in each case you can explain what the measurement model is for that experiment, i.e., is it just an internal model state (snowdepth, LST?) or a prescribed diagnostic relationship (FSCA?). In cases where it is a prescribed diagnostic relationship, how is that handled within the framework? I imagine that the current FSCA is built-in to FSM2, but what if an alternative representation was desired. Would that be handled via modification of the FSM2 snow model, or via another method.

11. Section 5: I found the organization of Section 4 flowing into Section 5.1 hard to follow. As mentioned above, I would suggest using the same language to refer to the three experiments throughout to help in this regard. It wasn’t clear to my why the benchmark (single cell) and distributed snow depth results were presented together in Section 5.1. They are described as two different experiments in Section 4 and so I think it would be easier to follow if they were treated as such in Section 5. Lumping them together in 5.1 seems a bit disjointed. Or maybe “single cell” here does not refer to the benchmark case (although “single cell” is used in that context too)?

12. With respect to Table 1, it is not clear what the reference data being used to compute RMSE is. It implies snow depth, but the description of what data was assimilated in the benchmark experiments is unclear (see comment above). The notation used for each scheme is also not defined. Perhaps define PF-c, PF-r, in caption?

13. In the context of Figure 3, it would be helpful to explain the meaning of “MDA” when only snow depth is being assimilated. I believe this method actually differs in this case due to its iterative nature rather than multi-data? This comes into play later where different notation is used to refer to iterative versions of method. Perhaps you can harmonize how you refer to iterative methods across the manuscript.

14. Figure 3. Refer to which experiment this corresponds to. And is this a particular cell? Is it the one shown in Figure 2?

15. Line 797: The reader would benefit from more description of how the prior forcing perturbations are generated in this context and how the posterior emerges from that. Can you clarify whether prior was identical across space and why patterns in the posterior emerge. Is there anything to be learned from the posterior uncertainty of these, i.e. is one more certain than the other (i.e. precip. vs. temperature? And why are the posterior patterns between the two fields seemingly so highly correlated. More discussion either here or in Section 6 would benefit the reader.

16. Discussion associated with Figure 6. Indicate that the fields in Figure 6 are the posterior mean. Units should be associated with temperature. Could more discussion be provided to hypothesize why the patterns are what they show.

17. Figure 7: There are inconsistencies (and typos.) between the use of what should be “LST” in the caption and “SST” in the figure. Is SST meant to be “snow surface temperature”. If that is preferred, SST should be used throughout instead of LST. The acronym “IKS” should be defined in the caption.

18. Figure 8: Acronyms need to be defined in the caption and reconciled with earlier ones. How does the Ensemble Smoother – MDA compare to any of these? Is it the same as I-ES?

19. Line 869: It is not clear what is meant by: “The assimilation of the FSCA provides information when FSCA saturates at 1, …”. Should this read “… does not provide information”?

20. Code and data availability: It seems that the MuSA code from the original github repository vs. the version provided on zenodo are different. In particular, when run on a mac, the github version worked, while the zenodo version did not. It appears to center on differences in the code, where the latter crashed out when finding the OS to be ‘darwin’ (macOS) instead of ‘linux’. I suggest making sure to reconcile the two so that the one posted on zenodo works. It would also be helpful for reproducing the results to 1) tabulate key parameters specific to each experiment (as suggested above) and 2) providing the actual input files for each experiment with the code distribution. This would make it much easier to reproduce the results from the paper and extend the framework to other cases rather than having to interpret which parameters to change.

Technical Corrections:

This is not an exhaustive list of typos., but ones that jumped out:

1. In Figure 1 there is a typo., where “weigths” should instead be “weights”.
2. Line 205: Typo. in the phrase “the are usually”.
3. Line 383: Typo. in the phrase “we will let denote anamorphosed”.
4. Line 391: “SM2” should be “FSM2”.
5. Line 491: “converege” should be “converge”.
6. Line 790: “smothers” should be “smoothers”.

In this paper, the authors developed MuSA (v1.0), a standalone snow data assimilation system. MuSA provides a comprehensive framework encapsulating FSM2 (a widespread opensource intermediate complexity snowpack model) with six variants of ensemble bayesian data assimilation algorithms revolving around the Ensemble Kalman Filter (EnKF), the Particle Filter (PF) and their smoothing counterparts. As such, this framework is a great contribution to the snow data assimilation community, as it allows to seamlessly compare different assimilation strategies. It also seems accessible as an educational tool and outreach. Furthermore, the literature only offers data assimilation implementations that are strongly bound to operational constraints (such as country-specific numerical weather prediction models or high-performance computing infrastructures), and are therefore not transferable in space, while MuSA could conceptually be used anywhere, anytime, something which is definitely missing at the moment. Indeed, it is sometimes virtually impossible to compare the performance of algorithms produced by different teams worldwide. The authors present the potential of MuSA in two different assimilation experiments: assimilating drone-base high-resolution maps over a small catchment in the Pyrenees (with a focus on a single “pixel”), and a combination of snow cover fraction and Land Surface temperature from MODIS, which allow to exhibit the different behaviour of the respective algorithms over two snow seasons and make a good entry point for new users.

The data assimilation algorithms implemented in MuSA are all already available in the literature. Only a tiny level of novelty resides in the use of the heuristic “redrawing from a normal approximation of the posterior” which is introduced in l. 99-107 to fight PF degeneracy but would deserve more stance and details. The idea is similar to approaches where the perturbation parameters are not resampled, only resampling the model states (Cluzet et al., 2021), or where perturbation are added to the posterior resampled parameters (Piazzi et al., 2018).

Then comes a very long Section 3 (17 pages) of theoretical developments and digressions which is rather ambitious as it tries to bring all the algorithms together into the framework of Bayesian inference almost in a review-like style. This comes at the cost of lengthy digressions and theoretical developments that have little interest in explaining the interest of MuSA and the presented results but have the merit of putting everything together in one single piece of literature, which will be really appreciated by the community outside of data assimilation experts and is in the spirit of GMD. I strongly recommend the authors to try and make this section more concise and accessible to such a public and am pleased to provide some suggestions in the attached annotated manuscript.

The results represent a nice illustration of the capabilities of MuSA. They have a limited level of scientific novelty and lack methodological depth for any significant conclusion to be drawn. At the light of this consideration, I recommend the authors to temper or even remove most statements regarding the conclusions drawn from their results in the discussions (see the notes below and the attached PDF). I don’t think, for instance, that this study is substantiated enough to discuss the relative advantages of EnKF vs. PF approaches, but it definitely lays a nice context to do so in further work. In this spirit, I recommend reorganizing their discussion with the aim of highlighting the novelty and versatility of the method, detailing potential use cases including educational purposes. More references would also be needed in that section. I think that this would ultimately sound more GMD-like and would increase the impact of the paper.

The code is well documented and seems written in a neat way that will guarantee its accessibility to the community and beyond.

Overall, I think that this paper is relevant for publication in GMD, as it provides a significantly new modelling tool to the snow data assimilation community and beyond. I am confident that the authors can address the minor recommendations detailed above. Below are listed some technical points. Please also make sure to address all the comments and suggestions in the attached commented manuscript.

Technical notes

• Several code performance/implementation features are disseminated throughout the manuscript. This is definitively an important added value from MuSA, that would deserve more stance. Consider grouping them?
• - l. 84: a thorough description of the direct insertion is required, as this is not as trivial as it seems: when assimilating HS, are you adding/removing mass, or do you just squeeze the snowpack layers? How do you handle the relayering?
• 125-135: please improve the style, this could be written in a more concise/impacting way.
• 157-163: same comment.
• 168-171: While assuming independent observation errors seems a reasonable assumption, attributing the same value for all observations of a given variable does not seem state of the art (take the example of the snow depth model error of e.g. Magnusson et al., 2014). Consider commenting on that since this is not much of a technical hurdle for further developments.
• - l. 172-174: while there is no conceptual shortcoming in the case of smoothers, please give more implementation details as to how time-irregular observations are assimilated in the sequential filter algorithms: are they assimilated ‘on the fly’, or grouped by daily/hourly batches? For instance, I can imagine some IO bottlenecks or PF degeneracy problems in the situation where an observation would come every 15mins and would be assimilated on the fly.

• 494-512: could be removed
• 577-579: Numerical instabilities: please expand as the reference is not open source. And what if several/all particles obtain likelihoods below the machine precision? How do you ensure that the likelihoods are still “sorted” in a proper way?
• 4: consider enclosing each individual algorithm into a numbered algorithm block
• 3.5: consider condensing the bibliography on the EnKF
• 6:

• In the first section (l. 820-840), the authors consider that the literature often lacks sufficient justification for their DA algorithmic choice, arguing that EnKF-based algorithms (both in the filtering and smoothing fashion), have been overlooked by the community for bad reasons (l. 827-828) (“due to the Gaussian linear assumption”). This statement is hazardous, and in the case of Largeron et al., (2020), just false. In Sec. 4.3 of Largeron et al., (2020), they list several additional advantages of PF vs. EnKF in the case of detailed snowpack models. Piazzi et al., (2018) as well as Cantet et al., (2019) also discuss the advantages of PF vs. EnKF in similar terms. These reasons may be bad or weak, but if the authors want to contradict them, they should be cited and contradicted extensively, otherwise this discussion is just not fair. Nevertheless, while I consider that the authors might be right, even with more literature substance I am not sure that this fits in the scope of this paper, and I would rather see this in a following, dedicated publication.
• I also find that that section is very badly structured, jumping back and forth between filters and smoothers. Snyder (2008) only treats the question of Particle filtering. While I agree that their statements can be transposed to particle smoothers, it is necessary to acknowledge that jump. More generally, discussing the advantages of smoothers vs. filters does not seem very interesting, as they are different solutions for different problems: reanalysis and operational forecast, respectively.
• 848-849: “The assimilation of snow depth products has been shown to be a very robust approach to SWE estimation.” This statement is not substantiated by any results from the paper.
• L849-851: “The posterior maps of the perturbation parameters showed intricately detailed and physically-consistent spatial patterns, especially considering the fact that there is not any cell 850 intercommunication in MuSA”. Please use another term than “physically consistent”: although that’s the best we can do, it is not “physically consistent” to compensate an unrepresented process (wind drift) by other processes (precipitation and temperature biases).
• 854: “Thus, MuSA can be used to study the importance of missing snow processes in the FSM2 model”. When looking closer at Fig. 6, the acute reader can realize that the assimilation is not only compensating for unrepresented wind redistribution, but also from strong precipitation and temperature biases coming from the downscaled ERA input: is it really possible to disentangle “missing snow processes” from input errors in such a setting? This is highly debatable. And a reference to Günther et al., (2019) would be required here.

Figures

• Fig3: In l 397-400, the authors assume that the snow can be reset to zero during the summer (seasonal snow). Their theoretical developments in the corresponding section seem to rely on this assumption. However, in the top-right and top-bottom panels of Fig3., this is clearly not the case, with persistent snowpack through the 2019 summer. Can the authors comment on that?
• 4: the “spaghetti” representation of ensemble members with imbalanced weights must account for the weight of the ensemble in the transparency parameter. Otherwise, this visually gives “weight” to particle with negligible weight and can lead to very misleading conclusions. For example, here the PBS could be completely degenerate, with only one scenario carrying all the weights (which is actually the case, we see in the 2019 summer that the mean drops to strictly zero (without any “tail” as with e.g. the ES-MDA) while some members have still non-negligible snow amounts. Alternatively, one can only display ensemble statistics (mean, quantiles) in such a situation.
• Fig4: on the bottom-right panel, it is difficult to distinguish the timeseries at 100% zoom factor.
• 1 + Fig.5: (spatial validation with independent observation) We are missing crucial pieces of information here: what are the precise dates of the assimilated maps and when is the independent validation date? It obviously makes a big difference in the expected performance, especially with a smoother approach, whether the validation date is “closely surrounded” by assimilated observations or not.
• 5: the scale is not consistent with the scale in Fig. 2b
• 6: Please use fair color scales for temperature and precipitation. In saturation/greyscale, the +1.5 K values should be indistinguishable from the -1.5K values which is not the case here. If you’d use such a fair color scale, you’d make the strong positive bias in the prior much more obvious. A significant prior underestimation of the precipitation would also be visible. Please also comment on both in the results section.
• Fig 7: please introduce month/year in the time axis labels. For the sake of interpretation, I would put the SWE timeseries panel below the SST and FSCA panels, and with the same horizontal extent.

Code and Data policy

• No statement on the drone observations, I think.
• Adding MODIS data to the zenodo archive is necessary to facilitate full reproducibility of assimilation experiments
• A statement should be made regarding the availability of the Micromet code which was used to downscale ERA forcings.
• Portability: state somewhere (as in the github repo) MuSA does not support Windows yet (https://github.com/ealonsogzl/MuSA/blob/94a4a0377e84dbe5a8e1441ed28406bd2769e846/main.py#L140-L145), but that it seems definitely possible.

References:

Cantet et al., JoH, (2019), doi:10.1175/JHM-D-18-0140.1

Cluzet et al., TC, (2022), doi:10.5194/tc-16-1281-2022

Odry et al.,TC,  (2022), doi:10.5194/tc-16-3489-2022

Winstral et al., (2019), doi: 10.1029/2018wr024521