Referee 1 :

General Comments

The article identifies relative changes in subsoil stress caused by the snow cover in its fresh and dry state (when it is melting). The authors use ambient seismic noise to calculate these changes using coda DV/V wave interferometry. The hypothesis is that the melting snow can percolate through the soil surface and increase the pore pressure and density, leading to possible mass slips.

The article is well structured and adequately written. A significant contribution is that experimental results can be correlated with numerical simulations, which show that relative stress changes can be reproduced for the two physical states of the snow cover.

The article can be accepted with minor corrections.

I suggest a discussion of variations in dV/V estimation if atmospheric effects are taken into account.

• Such atmospheric effect might probably influence the measured dV/V where atmospheric change occurs, but we expect its amplitude negligible for most cases, compared to other environmental influences (Le Breton et al., 2021; Hotovec‐Ellis et al., 2014). The additional loading from the snow cover is much more important than atmospheric pressure variations. Then, following previous literature, we can argue that a dV/V variation less than 0,1 % for atmospheric changes (few kPa) is expected. We can add this discussion into the text.

Also, the authors should include a figure showing the correlations obtained and indicating the part of the waveform in which the dV/V estimation is made.

• Yes, we agree with this suggestion. Find the new figure in supplement materials (Figure 1) showing the entire correlogram from seismic data, including the part of the waveform from which dV/V is estimated.

The authors assume that the coda is mainly composed of surface waves (not exist wave scattering). If this is not entirely true, then it should be discussed, as is the confrontation conducted in the modeling (section 4.1), where results are generated for Rayleigh waves.

• Although the coda is composed of both diffused surface and diffused body waves (see (Obermann et al., 2013)), we assumed that coda that we used for computing dV/V is mainly composed of diffused surface waves (early times of the coda, see Figure 1), that are most sensitive to shallow changes. This is an assumption, which we will explicitly mention in the revised version, that we used for modelling : actually we modelled Rayleigh wave velocity (diffusion is here not an issue). In fact, the energy partitioning dynamics favors Rayleigh waves in the early part of coda, when considering sensors and most of seismic noise sources being at (or almost) the surface (Obermann et al., 2013).

Specific comments 

Lines 80-90. For the reader to visualize how the influence of the snowpack will be modeled, the authors must make a scheme that illustrates this procedure.

• Yes, we agree with this statement. We propose to add a scheme (shown Figure 2 in supplement) in the new version of manuscript (part Modelling) in order to explain how both snowpack and ground are modeled, for one snowfall event (Snowfall 2), as well as the schematic location of the seismometers.

Line 150.  Here's an error, the depth of research in the refraction study can't be the same as the length of the line.

• The depth of refraction survey is of the order of 1/3 of the array aperture, say about 20m. The text was not clear and is now clarified -> “… a three layers model down to a depth of about 20m”

Line 200. In this part, the authors consider dV/V calculations on surface waves, however. The previous section does not indicate if the analysis of dV / V involves direct surface waves or assumes that the correlograms' coda is composed purely of surface waves.

• Direct surface waves are not used in the time window considered here : we process only the early part of the coda. In fact we assume that the part of the correlogram that we used for computing dV/V is mostly (not purely) dominated by diffused (not direct) surface waves. Since we use only vertical component geophones, we assimilate measured seismic wave velocity changes as Rayleigh wave velocities, rather than Love wave ones. We will indicate it in the new version.

References

 Hotovec‐Ellis, A. J., Gomberg, J., Vidale, J. E., and Creager, K. C.: A continuous record of intereruption velocity change at Mount St. Helens from coda wave interferometry, 119, 2199–2214, https://doi.org/10.1002/2013JB010742, 2014.

Le Breton, M., Bontemps, N., Guillemot, A., Baillet, L., and Larose, É.: Landslide monitoring using seismic ambient noise correlation: challenges and applications, Earth-Science Reviews, 103518, https://doi.org/10.1016/j.earscirev.2021.103518, 2021.

Obermann, A., Planès, T., Larose, E., Sens-Schönfelder, C., and Campillo, M.: Depth sensitivity of seismic coda waves to velocity perturbations in an elastic heterogeneous medium, Geophys J Int, 194, 372–382, https://doi.org/10.1093/gji/ggt043, 2013.

Referee 2 :

This study presents an interesting topic on the potential of using ambient noise seismology to monitor changes in the snow layer through precipitation/melting through quantification of velocity changes. It presents findings from a well-designed study of noise recordings over a period of a few months along with meteorological measurements and cameras for visual snow thickness estimation. The findings show that the noise recordings can appropriately pick up velocity changes linked to the snow fall/melt. A simulation of the snowpack is used to assess ; a simulation of the snowpack is used to assess the snow elevation. Finally, a simplified profile is used to forward model the velocity change as a function of frequency of Rayleigh waves, this profile assumes the soil layer is unchanged and assumes a temperature relation for density and P-wave velocity with unchanged poisson ratio, which assumes dry snow (no liquid water). The authors show their model to explain the decrease in velocity during the snowfall period, but cannot explain the snowmelt events.

In my opinion the paper deserves publication after major review.

I suggest several edits:

- introduction: you mention previous delta V/V measurements in snow that are ambiguous or contradicting. It would be worthwhile to comment yourselves on the previous studies, as this is exactly what you are trying to solve.

• We will add this following part in the new introduction (after l. 42): “Some observations show a positive correlation between snow depth and dV/V measurements at seasonal scale (Hotovec‐Ellis et al., 2014; Wang et al., 2017), whereas Wang et al. (2017) mentioned a negative correlation during intense snowfalls. In ice sheets, Mordret et al. (2016) modelled the effect of snow accumulation in by using poroelasticity and viscoelasticity at seasonal scale. But to the best of our knowledge, the  effect of snow on dV/V in snowy temperate regions has not been properly studied with high resolution (Larose et al., 2015, Fig. 10).”

-introduction: again the readers would benefit from a short intro on what the snow cover consists of (definitions of dry vs. fresh snow, phenomena of compaction, where does melting take place), how it interacts with the subsurface, and consequently how any snow changes affect velocity (is it all temperature related, phase changes, and if at all geomechanical changes due to loading/unloading?

• We will add the following part to the introduction: “Snow is a highly porous material with low density and low elastic modulus (Gerling et al., 2017) Typical densities for a seasonal snow cover range from 50 to 500 kg/m3 (Schweizer and Jamieson, 2003). Fresh snow generally has a density between 50 and 150 kg/m3, yet due to snow settlement (compaction), density rapidly increases. Snow is a material that exists very close to its melting point, causing rapid micro-structural changes (e.g. Herwijnen and Miller, 2013). During the winter season, when air temperature mostly remains below freezing, there is no liquid water in the snowpack and snow temperatures are below zero. This is called a dry snowpack. In spring, warm temperatures and solar radiation cause daily surface melting. As a result, snowpack temperatures gradually increase to zero degrees, and the liquid water content increases. This is called a wet snowpack.  Elastic wave velocities in snow, like most of its mechanical properties, including the elastic modulus, are highly dependent on snow density, temperature and liquid water content. While the effect of snow density and temperature are well documented (e.g. Schweizer and Camponovo, 2002; Sayers, 2021), the influence of liquid water content is still poorly understood.”

-field site and instrumentation: Somewhere in this section or the next, you could include a schematic cross-section that defines the main layers at your field site would be beneficial. Also perhaps a time-lapse cross-section to show which layer from the snow changes (do both bottom and top snow layer change velocity, or do we define the bottom layer of snow as the one without changes in velocity, merely changes in thickness).

• Yes, we agree with this statement. We propose to add a scheme (see Figure 1 in supplement) in the new version of manuscript (part Modelling) in order to explain how both snowpack and ground layers are changed, during one snowfall event (Snowfall 2), as well as the schematic location of the seismometers on the field.

-field site: are sensors installed 30 cm to 50cm into the soil? Is there no snow at the time of burial? 

• Yes, the sensors are set up into the soil, without snow at the time of burial. Here it is worth noticing that our study doesn’t depend on the depth of sensors that we used, since we addressed here only surface wave velocities, that are not depth-dependent (contrarily to the bulk or surface wave amplitude). Actually, phase and velocity of surface waves do not depend on depth.

-field site: sensors installed at a depth with snow falling on them. Does this imply that when you model Rayleigh waves, you should extract Rayleigh wave solution not at the free surface but at a certain depth. Is this effect negligible?

• Again, we will indicate that our study doesn’t depend on the depth of sensors that we used, since we studied only surface wave phase velocities that are not depth-dependent (contrarily to the wave amplitude).

-CC and delta/V estimation in Figures 3-7. Why do you use a different CC cut-off? I understand that a smaller CC would decreases your confidence in estimated delta V/V. Is there a remedy to increase the CC by choosing a moving average for the reference? And then accumulate the delta v/v from this moving reference to a zero reference of your choice?

• We tried to increase the CC by using a moving-averaged reference period, but the results weren’t much different actually, as well as for an averaging over the whole data period. For all cases, the response of dV/V just after snowfall events are visible. Thus, since the type of reference period doesn’t seem to play a role in our study, we chose to only show the results of dV/V and CC over the entire winter season (see Figure 2 in our article) for a fixed reference period (from January to February). But we used such moving average for the modelling step, in order to compare the elastic state of the reference (few hours before the snowfall event) to the one during the snowfall peak (few hours just after the snowfall).
• Also, we used a different CC cut-off for preventing dV/V biaises (outliers due to cycle skipping or measurement error) that are not interpretable to our knowledge (one peak of +10% increased dV/V in 6 hours is impossible). CC cut-offs are different over the snowfall events, but it can be physically understandable, as the moving average for the CC computing is time-dependant, by definition.

-snowpack simulation: please give more details in the physics of snowpack (ie FEM, what problem it solves mechanical, thermal, fluid flow?) State here your outputs  snow elevation, density, temperature as a function of depth etc.

• We will add the following to provide more details about SNOWPACK: “SNOWPACK simulates snow microstructure and the layering of the snowpack based on weather data. It is based on a Lagrangian finite element implementation and solves the non-stationary heat transfer and settlement equations. It encompasses phase transitions and the transport of liquid water. SNOWPACK provides detailed information on the mechanical and physical properties of each snow layer, including temperature, density, liquid water content and snow microstructural descriptors.”

-I do not understand since you are getting a temperature profile from the snowpack simulation why you then choose to have a two layer model with average temperatures to use in the numerical model for Young’s modulus E? The Young’s modulus model appears to have an exponential relation with temperature, averaging temperature prior to estimating the Young’s modulus will not give you the average Young’s modulus. I see no merit in predicting a two-layer model for the snow’s elastic constants in the forward model, as you never quite interpret (invert) your dispersion curve for two layers. It would make some sense to at least plot the continuous depth temperature and consequently continuous modelled velocity profiles to understand where the velocity changes are occurring. Then you may have an argument on modelling two-layer snow model in your Rayleigh velocity modeller, even though Geopsy I am sure accepts a velocity profile (if this is not possible could input many small layers). I believe that simplifying without reason your model may be one of the reasons for some discrepancies noted in your figures 9-12.

• We agree with this statement. We primarily chose to model the snowpack with only two layers for the sake of simplicity. Depth-averaging temperature and density profiles along only two sub-layers (one composed of fresh snow, the other composed of older and more compacted snow, as depicted in part 3.1) is a common procedure to address the complexity of a snowpack in snow physics. In our case, we can improve our depth resolution by considering much more sub-layers than only two. Thus, we agree with showing the results of modelled dV/V from the forward problem solved by Geopsy, and discuss if the depth resolution would actually be improved and if our simplifying is consistent. We test this procedure for the most significant snowfall event (SF2), by modelling snowpack properties with 10 cm thick sub-layers. The corresponding results are shown in Figure 2 (see Figure 2 in supplement) : whereas the modelled dV/V (blue curve) over-estimates the observations (red squares), the order of magnitude of the decreasing is still right over the frequencies. This findings shows that we do not precisely control the qualitative parameters of the study, but also that our physical model and its interpretation are valid. In this sense, our depth-averaging is then consistent, if considering this work as a first step towards a comprehensive study of the relation between dV/V and snow.

-It seems really low-hanging fruit to properly invert for your velocity curve from your experimental data. This would then properly give an understanding of where deltav/v varies more/less and why. We do not have this information readily available in the delta V vs. Freq. domain. This would also give you an opportunity to check constraints – having constant soil layers, or allowing the velocity in the first soil layer to change etc. Based on your Vp/Vs/density/temperature relationships, you could have even inverted for the temperature and seen if it agrees with the snowpack model.

• In the first version of our study, we chose not to invert dispersion curves from our dV/V measurements in order to solve the Vs and Vs changes along depth, because of the lack of precision for dV/V measurements in broad frequency ranges. Furthermore, the solution of this inverse problem would be non-unique, without important assumptions (dealing with at which depth no change occurs).
• To our mind, the inversion of temperature and density profile from the results of the inversion of our experimental dispersion curve is very challenging. Indeed, the relations between temperature, density and elastic properties of the snow are not well constrained, making the inversion of this model quite impractical. But it is true that seismic dV/V retrieved from ambient noise will probably bring additional constraints on the mechanical structure of the snowpack ; this present work is the first step in this direction.
• Although these cautions, we can try to invert for the dispersion curve (interpolated from our measured dV/V along frequency) for one example of snowfall event, assuming no changes in deep layers (bedrock and soil) for testing. And we will add a discussion dealing with this inversion in the new version.

-State early on when you present your model in 4.2 that it is for dry snow. It should be made clear that it cannot address the melting of the snow.

• Yes, we will indicate earlier and clearer this statement in the new version (part 4.2).

-While you state that a 3-phase modelling is outside of your scope, at least what is the qualitative understanding on the effect of having water on your delta v/v. How does this qualitative understanding translate to figure 12.

• The influence of water on dV/V in the soil (shallow subsurface) is actually well documented by observations and experiments (monitoring of groundwater level by coda wave interferometry, by example (Voisin et al., 2017; Le Breton et al., 2021 for a review). All studies show an negative correlation between the water content and dV/V : when water is increasing, dV/V is decreasing. But this effect is rarely modelled. In our case, one qualitative understanding may suppose that melting water percolation in late winter can saturate the soil in shallow porous layers, leading to a sudden drop of shear modulus, and thus a decrease in S-wave velocity. While Rayleigh wave velocity is mainly controlled by S-wave velocity (Grêt et al., 2006), we can then assume a negative drop of dV/V during melting periods (as observed actually). The quantitative understanding of the influence of water content in porous media is often studied by using Biot-Gassmann theory (for example, see (Voisin et al., 2016; Sidler, 2015), providing a drop of several percent for similar study sites. But in our study, the influence of snowpack and the melting water content are both unknown, making a 3-phase modelling very challenging. There is currently no literature on seismic waves in wet snow. We can expect that water will increase the density, and melting will decrease the rigidity (contacts between grains), all together decreasing the shear wave velocity, and thus decreasing dV/V.

References

Gerling, B., Löwe, H., and van Herwijnen, A.: Measuring the Elastic Modulus of Snow, 44, 11,088-11,096, https://doi.org/10.1002/2017GL075110, 2017.

Grêt, A., Snieder, R., and Scales, J.: Time-lapse monitoring of rock properties with coda wave interferometry, 111, https://doi.org/10.1029/2004JB003354, 2006.

Herwijnen, A. V. and Miller, D. A.: Experimental and numerical investigation of the sintering rate of snow, 59, 269–274, https://doi.org/10.3189/2013JoG12J094, 2013.

Hotovec‐Ellis, A. J., Gomberg, J., Vidale, J. E., and Creager, K. C.: A continuous record of intereruption velocity change at Mount St. Helens from coda wave interferometry, 119, 2199–2214, https://doi.org/10.1002/2013JB010742, 2014.

Larose, E., Carrière, S., Voisin, C., Bottelin, P., Baillet, L., Guéguen, P., Walter, F., Jongmans, D., Guillier, B., Garambois, S., Gimbert, F., and Massey, C.: Environmental seismology: What can we learn on earth surface processes with ambient noise?, 116, 62–74, https://doi.org/10.1016/j.jappgeo.2015.02.001, 2015.

Le Breton, M., Bontemps, N., Guillemot, A., Baillet, L., and Larose, É.: Landslide monitoring using seismic ambient noise correlation: challenges and applications, Earth-Science Reviews, 103518, https://doi.org/10.1016/j.earscirev.2021.103518, 2021.

Mordret, A., Mikesell, T. D., Harig, C., Lipovsky, B. P., and Prieto, G. A.: Monitoring southwest Greenland’s ice sheet melt with ambient seismic noise, 2, e1501538, https://doi.org/10.1126/sciadv.1501538, 2016.

Sayers, C. M.: Porosity dependence of elastic moduli of snow and firn, 1–9, https://doi.org/10.1017/jog.2021.25, undefined/ed.

Schweizer, J. and Camponovo, C.: The temperature dependence of the effective elastic shear modulus of snow, 35, 55–64, https://doi.org/10.1016/S0165-232x(02)00030-7, 2002.

Schweizer, J. and Jamieson, J. B.: Snowpack properties for snow profile analysis, Cold Regions Science and Technology, 37, 233–241, https://doi.org/10.1016/S0165-232X(03)00067-3, 2003.

Sidler, R.: A porosity-based Biot model for acoustic waves in snow, 61, 789–798, https://doi.org/10.3189/2015JoG15J040, 2015.

Voisin, C., Garambois, S., Massey, C., and Brossier, R.: Seismic noise monitoring of the water table in a deep-seated, slow-moving landslide, Interpretation, 4, SJ67–SJ76, https://doi.org/10.1190/INT-2016-0010.1, 2016.

Voisin, C., Guzmán, M. A. R., Réfloch, A., Taruselli, M., and Garambois, S.: Groundwater Monitoring with Passive Seismic Interferometry, 09, 1414, https://doi.org/10.4236/jwarp.2017.912091, 2017.

Wang, Q.-Y., Brenguier, F., Campillo, M., Lecointre, A., Takeda, T., and Aoki, Y.: Seasonal Crustal Seismic Velocity Changes Throughout Japan, Journal of Geophysical Research: Solid Earth, 122, 7987–8002, https://doi.org/10.1002/2017JB014307, 2017.