General comments:

The study ‘Impact of rain-on-snow events on snowpack structure and runoff under a boreal canopy’ by Bouchard et al. explores how rain-on-snow (ROS) events impact sub-canopy snowpack structure based on experimental findings, and further assesses whether the physics-based, multi-layer model SNOWPACK succeeds in capturing these effects. Improvements to the representation of canopy snow in SNOWPACK are also suggested. I really enjoyed reading this manuscript. It addresses multiple topics that are known to be important and understudied in the forest snow (modelling) community, and I commend the authors for the laborious and extensive field campaigns – such datasets are unique and highly valuable for the community. The manuscript is well written and generally well structured. The methods are very clear and easy to read (although some additional details are needed, see below). The results include many interesting findings but would benefit from some streamlining (potentially just an introductory paragraph). In terms of content and logic, I identified some issues that should be addressed before the study can be considered for publication, see specific comments below. The suggested changes mainly aim at making the argumentation easier to follow, do not question the novelty and value of the study, and I don’t see any of them being problematic to tackle. I am confident that this article will make for a great contribution to the literature following revisions. I encourage the authors to reach out if they have questions to my comments and wish to discuss any of them further (giulia.mazzotti@meteo.fr).

Specific comments:

1. Study goals (L 94 ff): I feel that the current way of introducing the study objectives does not do justice to the interesting research presented thereafter. The statement ‘this study aims to help fill the research gap on ROS events in the boreal forest’ is very vague. Then, it is suggested that modelling is used as a tool to help field data interpretation. Then again, the third objective addresses a model improvement with a non-obvious link to ROS (and remains somewhat disconnected from the two former objectives). And finally, ‘improved understanding’ is stated as overall goal, which disregards the model development work.

You could make a much stronger case by better connecting the three aspects of the study and by highlighting the model testing/development as part of the overall objective. I suggest arguing somewhat along these lines: ‘This study addresses observation and modelling of ROS events in the boreal forest. Specifically, we first use field observations […]. These field observations provide the unique opportunity to evaluate the multi-layer, microstructure resolving model SNOWPACK […]. Motivated by model shortcomings identified in this evaluation (or: motivated by the hypothesis that unrealistic properties of unloaded snow hamper the performance of SNOWPACK in simulating how sub-canopy snowpacks react to ROS events), we finally suggest an alternative representation of canopy snow properties and assess its impact on SNOWPACK simulations […]. Altogether, this work improves both our understanding of ROS events in forests and our ability to capture their impacts on snowpack structure and runoff by physically-based models (which is important as such events are expected to become more frequent).’ That’s just a suggestion, I hope you get the idea…

2. Logic of linking ROS and canopy snow: While the title suggests a focus on rain on snow events, the modelling part focuses on the representation of canopy snow. The authors seem to implicitly assume that the two processes are interlinked., i.e., that canopy snow and its unloading are what makes the impact of ROS in forests different to their impact outside of forests. While this might be correct to some extent, it is not intuitive/self-explanatory, and not yet sufficiently justified in the study. This aspect needs to be worked out better, and I think there are multiple components to it.

• Experimental evidence: it needs to be stated clearly that the snowpack features of the sub-canopy snowpack that favor runoff from ROS events are not present in canopy gaps (or open areas) and result from interception and unloading. This is implicitly alluded to in the paragraph starting L 54 but including a short summary of the main differences between sub-canopy snowpacks and the snowpack in canopy gaps found by Bouchard et al. (2022) would make this link more explicit.
• Causality: The title suggests that the study addresses how ROS events affect the snowpack structure, but a big part of the study is about how a snowpack structure shaped by interception and unloading affects the snowpack’s reaction to ROS. By not acknowledging this distinction, the two processes become somewhat entangled. It might be worthwhile to reconsider the title (although I don’t have a concrete suggestion).
• Model choices and structure in SNOWPACK: one of the reasons why the above-mentioned causality is confusing is that in SNOWPACK, rain on snow and unloading events may indeed be entangled, but this may be a consequence of some model choices that are quite unusual; in particular, the fact that unloading occurs only when interception capacity is exceeded (as far as I know, most common unloading parametrizations include some sort of time decay function to remove canopy snow even if interception is below capacity). In Gouttevin et al., it is stated that this choice allows for a gradual unloading. However, a maybe less evident consequence of this is that ROS events induce big unloading events because the onset of rain entails a large drop in interception capacity (I hope I am interpreting this correctly, but it seems to follow from what you write in L 657 ‘when canopy interception is in the liquid phase’, and it agrees with Fig 10). This could be why ROS and unloading become linked. I wonder how much your results would change if unloading was triggered differently in SNOWPACK (e.g. by wind, when capacity is not yet reached, or using other approaches as summarized e.g. by Lumbrazo et al. 2022). I am not saying that you need to test alternative parametrizations, but I think that this particularity of SNOWPACK should be addressed in the discussion to facilitate interpretation of your results, and you should acknowledge that your results are somewhat specific to SNOWPACK (probably in Section 5.3). In the same section, I suggest you also comment on the fact that the unloading implementation in SNOWPACK still lacks some processes that would likely impact the subcanopy snowpack (e.g. drip unloading?).
• Comparison to non-forested sites: Since Bouchard et al. 2022 and 2023 have observations in forest gaps, it would have been really cool to see what the model does there (i.e. to collect further evidence that it is indeed the canopy snow that makes the difference). The study in its present form is already rich enough and I am aware that conducting such simulations would be a lot of additional work, so I am not asking to add this to the revision, but if you plan to expand on this work in the future, I think this would be the perfect next step (consider mentioning in outlook?). I know that SNOWPACK is not yet suited to simulate snow in canopy gaps, but I think that the approach used here: https://doi.org/10.5194/egusphere-2023-2781 could easily be applied to SNOWPACK as well.

3. ISD implementation – motivation, interpretation of results, discussion: The implementation of canopy snow densificaty (and microstructure) evolution is doubtlessly a nice contribution addressing a need that has been highlighted multiple times in the literature and could be valorized more. At the beginning of sectionn 3.2.1, instead of stating ‘as part of our third research objective’, I would recall the reasoning behind this development (i.e, currently unrealistic representation of unloading snow properties). In particular, I would note that we expect a microstructure resolving, multi-layer model like SNOWPACK to be more sensitive to the simplified representation of canopy snow than ‘simpler’ snow models these parametrizations were originally intended for.

The results section left me with some open questions that should be addressed (or maybe you just need to guide the reader a bit more…):

• Figures 10 and 12: I am a bit surprised that the density of the unloaded snow with the IM version is always so low. After all, density does depend on temperature, and the Appendix states that it can reach up to 250kg/m3. Since most unloading (and especially the example shown in Figure 12) occurs during ROS, I assume this to coincide with a rather high air temperature so I would have expected a higher density, too. Do you have an explanation for this? Maybe it would help to describe more in detail what happens in the model when snow unloads at the onset of a ROS event.
• In Section 4.3.2: Is the formation of MF layers that are detected by ISD but not by IM always preceded by an unloading event? I would specify this. Related: did you look at the impact of unloading snow at the snow surface also when there is no ROS?
• In Section 4.3.3: You suggest that more preferential flow and snowpack runoff occurs with IM than with ISD. Wouldn’t this match your observation of preferential flow channels better? The decreased runoff with ISD seems contradictory to Bouchard 2022 who states that the subcanopy snowpack structure should favor rapid runoff?
• Section 4.4: The large sensitivity to dg left me wondering whether the real impact of your model adaptation really comes from the densification, or whether it is from the treatment of snow microstructure descriptors. My understanding is that while densification is a continuous function (actually, also in IM), the treatment of canopy snow grain type/size is binary (in ISD, while in IM they always get the same properties, see Fig 2?). I am not familiar with the relationship between water content, saturation threshold for preferential vs. matrix flow, and snow properties (density, grain size?). From a quick look at Wever et al. I get that grain radius is the key parameter. Now my question: does the modification of density relative to IM really make the difference, or was the real problem in IM that the unloading snow was attributed fresh snow properties? I think it’s crucial to explore and clarify this (either way it’s an interesting finding – you identified an important shortcoming of the original implementation!).

Technical comment:

Most of these questions/suggestions are aimed at improving readability / eliminate unclear phrasing and typos.

L 20: lead to the formation (without s)

L 28: Use of ‘indeed’ is awkward. ‘Specifically’? ‘In fact’?

L 30: ‘Our results show …’ – this statement will need revision. To really show this, you would need to present data of a snowpack outside of a forest as well.

L 66: 2015 is no longer ‘recently’ (unfortunately! ;-) )

L 76: ‘This is partly due to limited field data’ – a bit too vague. What data? And do you mean ‘data to inform model improvements’?

L 84: ‘can undergo metamorphism before being unloaded’ – this would be a good place to note that these interception and unloading parametrizations were originally intended for models that are much less complex than SNOWPACK and do not resolve snow microstructure, see my comment above.

L 119: ‘when estimated’ – maybe worth mentioning the method (just say if that was laser altimetry, visual inspection, photogrammetry?’). Same for BRV site.

Table 1: I would move this to the Appendix, it’s very technical. But that’s a matter of taste.

L 138: 200km is a lot!

L 143: The big difference in LAI is surprising. The hemispherical images look quite similar, and the trees are taller in the case of BRV. Could it be due to the method used to derive these values? What does this imply for the use of detailed snowpack models?

L 186: It is unclear to me whether this was done based on the data measured at the sites (i.e., the automated snow under the canopy measurements), the time-lapse cameras, the data at the stations, or a combination of all. Please be more specific.

L 187ff: Here, I think you would do your storyline a favor by justifying why you chose to use SNOWPACK (and highlighting that your data offers an excellent opportunity to evaluate and further improve it for forest applications).

L 202ff.: It would be helpful to add a few words on what snow properties are used to determine this saturation threshold and therefore the transition from preferential to matrix flow, see my major comment #3. I think this will facilitate interpretation of Fig. 12 and make the justification of why you use a higher threshold than Wever et al. more convincing.

L 209: ‘not shown’: consider putting in supplementary material if you feel this is useful for readers (not a requirement)

L 210: ‘geometric mean’: it is unclear of what, and I am unclear about the purpose of providing this information. Of the hydraulic conductivity of the two layers? To get hydraulic conductivity at the nodes? Why would you need that? Please rephrase to make this clearer.

L228: This sounds as if Koch et al. derived this formula for this exact purpose, but that’s not the case. I think you should give some context to avoid misunderstandings (something like: ‘we use the formula from Koch et al., originally developed for … based on data from …’)

Table 2: Heading ‘For’ seems a typo, please revisit. It would further be useful to add symbols to the variable names (e.g., I think that the ‘direct throughfall fraction’ is the gap fraction cf (–) in Eq. A2, but it would be good to make this unequivocal).

L294: Does this mean that the phase determined from the images replaces the phase determined with the dual threshold? If yes, why do you need the dual threshold? Please clarify.

Figure 2:  For IM, is rho_new the same as rho_s,int in L 661 / A4? I guess so, but I would advise you to be consistent with the symbol because it can get confusing really quickly…

Table 3 is not referenced anywhere in the text. I would mention here that you do a sensitivity analysis, so that Section 4.4 comes as less of a surprise.

Section 4: consider giving a short overview of the result sections to follow, I got a bit lost in the many subsections…

Figure 3 is nice, but it would help a lot to add the acquisition times of snow pits. This would also help interpret the information in Table 4, because I guess you are much more certain about the formation date of layers where you have more frequent observations? And the number of observations is also directly related to number of field visits.

Table 4: ‘The formation date is assumed to be the first rain-on-snow event since the previous layer was observed for the first time.’ This took me multiple reads to digest. Suggest changing to: ‘The formation date is assumed to be the date of the first rain-on-snow event that occurred after the date of the snow pit acquisition during which the layer beneath the melt-freeze layer was observed for the first time.’ (I know it’s bulky…)

Figure 5: Really cool pictures. For d): what made you think that these water pockets formed in the snowpack, rather than being chunks of unloading snow? Same for e)/f)?

Figure 6: ‘The observed temperature was forced to a maximum of 0°C’: Do you mean that positive values were set to 0? Are these noise or a constant bias?

Figure 7 and potentially 8: Why did you decide not to include the ISD simulation in these Figures? It should be straightforward and would help contextualize the later results (without making the figures too dense)

Figure 10a: It seems that between November and mid-March, unloading in the model happens mainly in three big events, the two first of which match time periods in which the time lapse cameras detect the canopy becoming snow free. If I understand the way unloading is treated in the model, these large events can only be triggered by sudden large change in interception capacity and would coincide with ROS events. If that’s correct, please mention it explicitly, it would really help interpretation of the results.

Figure 10c: Looking at the other density profiles in the Supplementary Material, this seems to be the only case where the impact of ISD on surface density was observed. So it might be useful to explicitly point at the ‘layers resulting from unloading’ mentioned in L 442.

L 444 ff (to end of section): Really nice result!

L 505: Strictly speaking, this is just the model sensitivity, you don’t have experimental evidence. I would remove this statement / integrate with the previous sentence.

L 517: Awkward sentence. Do you mean ‘simulated runoff is more sensitive to canopy snow parameters than the number of melt-freeze layers’? Please revisit.

L 539: ‘Air temperature and amount of precipitation during a ROS event do not appear to be good predictors of melt–freeze layer thickness’: Please link to the result section / figure that backs up this statement, I wasn’t sure what you are alluding to.

L 558: The ‘rapid runoff response’ would better match the results you obtained with IM – so this statement is a bit confusing.

L 561: ‘General overestimation of snow density’ – where do you see that? The information from S8 to S15 is somewhat hard to digest...

L 564: I am not sure that this is because storage capacity is reached, it might just be too small interception (maybe I am missing something?).

L 570ff: Confusing; if SNOWPACK simulates too little radiation reaching the ground, it would rather UNDERestimate surface temperature? And why would the underestimation of snow surface layer density only impact the nighttime? To me this looks more like a longwave radiation effect, but I might be misunderstanding your statement. Please doublecheck.

L 619: Here it would be fair to acknowledge that some multi-layer snow models ARE already able to resolve tree-scale processes (SNOWPALM, ; FSM2, https://doi.org/10.1029/2020WR027572) – these models do not include microstructure, but concepts used therein could in principle be applied to SNOWPACK as well in the future. In this context, this preprint might be of interest as well: https://doi.org/10.5194/egusphere-2023-2781.

L 670: ‘Solid or liquid unload’: shouldn’t this be ‘unloading’?

Supplementary Material: I noticed that the ISD simulation in figure S20 onwards is labelled ‘Snow Tracking’ – maybe from an older version? It’s a detail and only the Supplementary Material, but consider correcting this for consistency.  

This study presents a detailed assessment of the reproduction of rain-on-snow events by an improved version of the SNOWPACK model, using multi-source information collected at two sites in Québec, CANADA. The paper is very well-written and documented, provides clear and illustrative figures, and interesting results and conclusions. I recommend publication and I will just provide some comments that could help to refine the presentation of this study.

General comments

In the literature, the ERA5-Land reanalysis has been shown to be globally reliable concerning mean precipitation characteristics, e.g. annual/seasonal accumulation but is quite limited for the reproduction of local/intense events (Reder et al., 2022, https://doi.org/10.1016/j.wace.2022.100407). It is also known to produce too many wet days (Bandhauer et al., 2022, https://doi.org/10.1002/joc.7269). I was a bit surprised by the motivation for the use of  ERA5-Land instead of the Hydro-Quebec gauged data (because of the mm resolution). It seems that ERA5-Land caused some issues (l. 374-375, l.385, l.564) and that there are important discrepancies in Fig. S3. Have you tried to use the gauged data instead of ERA5-Land data as forcing data even if the measurements have a limited precision? In any case, I find that the arguments against the gauged data should be strengthened. 

Minor comments

1. 27-28: I did not understand where the numbers 27 and 18 were coming from. In any case, I find that this is a very specific result that should not appear in the abstract since it is difficult to appreciate this improvement without knowing what they represent.
2. 29: “the unloading of dense unloaded …” could be rephrased to avoid the repetition.
3. 31: “modulates the sub-canopy snowpack structure and runoff from rain-on-snow events”: I think it is a very important result and a sentence could be added to describe in what way the snowpack structure and runoff are modified, as indicated in the conclusions (l. 635).
4. 61-62: At this stage of the manuscript, I could not understand the meaning of this sentence. After reading the whole manuscript, I guess it refers to l. 611-612. This comment belongs to the discussion/conclusion parts in my opinion.
5. 166-167: it could be helpful to add just one sentence to explain how the precipitation phase is obtained (without looking at the paper by Floyd and Weiler, 2008).
6. 194: SNOWAPCK -> SNOWPACK.
7. 231: “it could to be” -> it could be.
8. 246: are is in italic.

Table 2: missing point at the end of the caption.

1. 393: The point after “scale” should be removed.

Figure 8: The point after “Observations” in the legend should be removed.

1. 478: for both ROS events?