the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
| 18 Jul 2025
Simulating liquid water distribution at the pore scale in snow: water retention curves and effective transport properties
Abstract. Liquid water flows by gravity and capillarity in snow and modifies drastically its properties. Unlike dry snow, observing wet snow remains a challenge and data from 3D pore-scale imaging are scarce. This limitation hampers our understanding of the water, heat and vapor transport processes in wet snow, as well as their modeling. Here, we explore a simulation-based approach in which 3D images of dry snow were digitally filled and drained with water through imbibition and drainage simulations driven by capillarity. Time series of wet snow images at various water contents were produced. The water retention curves, i.e. the capillary pressure as a function of the water content, are derived for different snow types. New parameters to model the water retention curves of snow are proposed based on the van Genuchten model and compared to existing ones from laboratory experiments. From the 3D images, the hydraulic conductivity, water permeability, effective thermal conductivity, and water vapor diffusivity of wet snow were computed. We analyzed their evolution in relation to water content, density and snow type, compared them with existing regressions and presented new ones when needed. The proposed relationships are a step forward in improving the physical description of wet snow in snowpack scale models.
- Preprint
(2078 KB) - Metadata XML
-
Supplement
(31411 KB) - BibTeX
- EndNote
Status: final response (author comments only)
- RC1: 'Comment on egusphere-2025-2903', Anonymous Referee #1, 12 Aug 2025
-
RC2: 'Comment on egusphere-2025-2903', Michael Lombardo, 13 Aug 2025
The comment was uploaded in the form of a supplement: https://egusphere.copernicus.org/preprints/2025/egusphere-2025-2903/egusphere-2025-2903-RC2-supplement.pdf
-
RC3: 'Comment on egusphere-2025-2903', Anonymous Referee #3, 15 Aug 2025
The study addresses the spatial distribution of liquid water in different snow types and estimates its effect on various transport properties. The authors analyzed in detail 34 scanned dry snow samples and presented the results on water distribution not only in the main text but also in a well-organized and complete Supplementary Information folder (SI). This was well-appreciated and made it very easy (and interesting) to scrutinize the results of the study. Considering that experimental data on liquid water distribution at the pore scale of snow are rare, the authors based their analysis on numerical experiments. This strategy relies on evidence that the computed properties are reasonable and in agreement with experimental findings. More specifically, to trust the outcome of the numerical experiments, I would expect that (i) the capacity of the model is shown for a very few experimental cases (see a proposal below) and that the (ii) general trends are in good agreement with experimental studies. This proof is missing here, and I have some doubts regarding the results of the numerical experiments as I explain below.
Questions on the shape of the wetting curve
Typically, the shapes of the wetting and drainage curves are similar with a shift of the drainage curve to higher (absolute) capillary pressure values. This would be manifested in similar values of the parameter “n” for wetting and drying and in higher “alpha” values for the wetting curves (with “alpha” and “n” as shape parameters of the van Genuchten model to fit the relationship between water content and capillary pressure). While the study shows the expected correlation between the “alpha” values (with alpha_wetting ~1.3 alpha_drainage according to the data in the SI), the values of “n” are very different from the expectations, with more or less constant value “n” for all wetting curves (3 to 6) that differs from the drainage curve with values ranging from 5 to 20 (average 10). Such a clear mismatch of the shapes of the wetting and drainage curves is unusual and also contradicts the experimental outcome presented in Table 1 of Adachi et al. (2020) that reported very similar values of “n” for wetting and drying curves.
In short, I don’t trust the shape of the wetting curve with such smooth shape (different from drainage curve) and converging to the porosity. I would expect a large amount of entrapped air with water saturation values considerably smaller than the porosity.
Application of morphological pore network model (PNM)
I had good experience with the application of morphological pore network models (MPN) for the simulation of drainage processes in porous media, because it is controlled by the size distribution of the pore throats controlling the air invasion. The simulation of the wetting process, however, using spherical structural elements in the MPN, can be more critical: see Figures 5 and 6 in Ahrenholz et al. (2008) - you cited that study - and the sensitivity of the wetting curve with respect to the structural element. See the (i) large fractions of entrapped air in these figures and (ii) the similar shapes of wetting and drainage curves; these results are quite different compared to your findings shown for example in your Figure 4b (and in the data shown in the SI).
I am puzzled by the text on lines 138-140: ”For imbibition, the simulation starts with the 3D image of dry snow. A minimum pore radius is defined according to Eq. (5), which then increases step by step. At each time step, a pore is filled with liquid water if (i) the pore radius is below the minimum pore radius, and (ii) the liquid water phase is connected to the liquid water reservoir.” Such a description is not complete -> the air phase must be continuous as well; it can occur that the invading water encloses air bubbles and then the air must stay entrapped. Was this ignored in this model?
I’m also confused by the insets in Figure 2 that show the fluid distribution during the drainage process. The shape of the interface between water and air does not show the expected curvature ‘with an air sphere invading the water filled pore space’.
Thermal conductivity
The simulated thermal conductivity is a linear function of the water saturation. This linear relationship is unusual compared to findings for example in soils. In analogy to soils, I would assume that there is a fast increase in thermal conductivity with “the addition of a small amount of water” connecting the ice particles at the dry end. How can this strong linear relationship be explained? And what is with the case of theta=0 (after removal of the residual water)? Is there still a linear relationship between thermal conductivity and water content?
RECOMMENDATION
I propose that the authors compare their model with experimental data as follows: Adachi et al. (2020) plotted in Figure 6 measured drainage and wetting curves of three snow samples (“S”, “M”, and “L”). According to Table 1 in Adachi et al. (2020), the dry snow densities were around 500 kg/m^3 and the snow particle diameter was 1.2 and 1.9 mm for the two finer snow samples (“S” and “M”). Considering that the samples of Adachi et al. (2020) are melted forms (but refrozen) and that two of your ‘melted form samples’ (“NH2” and “NH5”) have similar densities and snow particle sizes as in Adachi et al., it would be meaningful to compare samples “S” and “M” with “NH2” and “NH5”, respectively. With such comparison, you may show that the prediction of the drainage curve is close to the experimental values (I doubt a bit that the prediction of the wetting is successful too).
After such experimental proof , the authors could relate the simulated properties to various metrics of the pore space. In the present version, the authors consider the mean curvature but there are other properties of the pore space that should be considered (maximum and width of pore size distribution, and connectivity of the pore space as for example quantified by the Euler-Poincarè-characteristics). After such more detailed quantification of the imaged pore space, new relations between pore space characteristics and transport properties could be deduced with more confidence. You could also consider to fit the 'tortuosity parameter' in the Mualem van Genuchten expression of hydraulic conductivity (now set to ½) to see if it changes with snow structure.
I propose to show the figures on the relationship between water content and capillary pressure using a linear scale of the pressure head, focusing on small absolute capillary pressure values. This allows a better assessment of the shape of the saturation-pressure relationship.
I suggest that the authors also show some cross-sections of the modeled phase distribution (similar to the mean curvature now shown in SI). This would facilitate to assess if the phase distribution is reasonable.
The authors should also provide more information on the morphological pore network model
Specific comments
Line 6: The term “time series” does not fit entirely, because temporal aspects are not considered
Line 39: I would write “shape parameters”, not “hydraulic parameters”
Lines 42: There are other lab methods (evaporation method and dew point methods for dry soils)
Line 53: I think that in Yamaguchi et al. (2012) there are 12 samples that are not melted forms (“rounded grains” in Table 1)
Line 56: Richards equation, not Richard equations
Line 115: I found no specific information on the morphological pore network model in that paper
Line 124: I propose to refer to the Young-Laplace equation because gravity is not included in the model (but the applied pressure)
Lines 128-129: The text “The evolution of the air and liquid water in the ice structure can be seen as a series of erosion and dilation operations that depends on the pore distribution (Ahrenholz et al., 2008)” does not fit here -> erosion and dilation can be used to determine the pore size in complex structure but the evolution of the fluid phases is controlled by the connectivity as well.
Line 131: How is it possible to reach saturation when air will be entrapped (and there must be air entrapment)?
Figure 2: Write in the captions the name of the sample. Is there a melt form sample with such a high porosity?
Line 144: Why the “initial” liquid water? The connectivity of the liquid phase must be tested for each drainage step.
Line 148: Have you tested the REV for the porosity? How is the porosity changing with image size?
Figure 3: this figure is “basic textbook” and is not needed here; please provide always the units of alpha
Line 166: it is the inverse of alpha that is related to the inflection point
Line 167: it is the “width” of the pore size distribution that controls the value of parameter “n”
Lines 168/169: this is not correct considering the presence of entrapped air
Line 169: state that the porosity can be deuced from dry bulk density and the density of ice
Line 229: I am not sure that the smooth shape of the wetting curve is realistic
Figure 6b: I cannot see the drainage regression line
Line 239: The exponential trend is not obvious because you have not enough values on the left side of the x-axis compared to Yamaguchi et al. (2012); in your data, there is a single sample with a x-value around 250000 kg/m^4; and a linear fit is almost as good as the exponential one for x-values ≥250000 kg/m^4.
Line 259: the differences in the findings of Yamaguchi et al. (2012) as function of snow type is important; but in your study you have also a few melted forms; have you checked if you find a trend in your melted-form-subset similar to Yamaguchi?
Lines 262 ff: To characterize the pore size variation, why don’t you use the pore size distribution of the scanned sample? You should have access to the imaged pore size distribution (it is used as basis of the pore network model) and you could easily compute its width. This would be a more direct link to pore sizes compared to the curvature.
Line 272 and Figure 7: the correlation is rather weak, and the nonlinearity is not very strong (similar correlation values for a linear fit); please explain colors in Figure 7
Line 283: this is an interesting result.
Lines 286/287: again-> I don’t think that this is reasonable
Lines 303/304: The correlation with the standard deviation of mean curvature is rather weak and such regression may result in large prediction errors
Line 314: Can you comment on the effect of applying different methods for the quantification of grain size?
Figure 8a and 8b: please use the same intervals on x axis (0.1)
Figure 8: why do the Yamaguchi curves have the same maximum saturation? I thought they used the 90% rule as maximum water content
Line 317: I do not agree with this conclusion considering the shape of the wetting curve, the lack of entrapped air in the model, and not enough samples to proof the non-linearity for small rho/d-values (smaller than 250’000 kg/m^4).
Lines 322/322: the new VG-model reproduces well the simulated values – but it is not clear that the simulations are correct
Figure 10: I’m surprised that there are some samples with conductivity values of ~10cm/day at effective water saturation of 0. How can this be? Do you have water flow along the residual water?
Equation 8 (and figures 10): you use the ‘standard expression’ with tortuosity factor “tau” of 1/2. This is an average value and can change with the porous medium (as discussed in Mualem 1976)
Line 343: the good agreement between the simulations and the ‘standard formulation’ of the Mualem-van Genuchten framework is no strong proof that the simulations are correct. You may try to fit the factor “tau” in Mualem van Genuchten for the different samples and check if the value changes with snow property
Figure 11: What are the predictions when you compute it without residual water? Maybe it is better to plot it as function of water content, not effective saturation
Line 386: how do you know that it is over-estimated?
Line 407/408: That is an interesting finding (similar approaches exist for gas diffusion in soils )
Line 414: you should compare this with other models like Moldrup et al. (2000)
Line 429: you should state this in the captions of figures 10-14
Line 441: the pore size was not presented in this study
Citation: https://doi.org/10.5194/egusphere-2025-2903-RC3
Viewed
| HTML | XML | Total | Supplement | BibTeX | EndNote | |
|---|---|---|---|---|---|---|
| 1,317 | 55 | 24 | 1,396 | 39 | 73 | 71 |
- HTML: 1,317
- PDF: 55
- XML: 24
- Total: 1,396
- Supplement: 39
- BibTeX: 73
- EndNote: 71
Viewed (geographical distribution)
| Country | # | Views | % |
|---|
| Total: | 0 |
| HTML: | 0 |
| PDF: | 0 |
| XML: | 0 |
- 1
The physical properties of wet snow are less well understood than those of dry snow, with fewer measurement data and experimental results available. In this study, the authors reported the relationship between the microstructure of snow and the physical properties of wet snow based on simulations. They also compare the simulation results with existing experimental results and propose new equations for estimating these physical parameters. Their simulation approach is unique, and there is no doubt that their method is useful in situations where experimental results are lacking in wet snow field. For these reasons, this paper makes a significant scientific contribution to the field of wet snow science and is worthy of publication in TC. The paper is well organized and easy for readers to understand.
On the other hand, the discussion in this paper is based on simulation results. Simulation results depend on conditions and parameters, so it is necessary to verify the accuracy of the simulation results based on experimental results or observational results. Until the accuracy of the simulation is guaranteed, it is not possible to evaluate the true scientific contribution of the research. Therefore, I believe that it is necessary to add a discussion on the accuracy of the simulation to the paper before it is accepted for publication.
I will list a few points I have noticed, including these points.
<Major point>
As shown in Equation (5), the pressure head of a liquid depends on the radius of the pore, and this is the basic equation used in the simulations in this study. In order to use Equation (5), it is necessary to calculate the radius of the pores in the snow. However, since the pores in the snow are interconnected, it is necessary to separate them in some way in order to determine their radii. Therefore, the distribution of radii may vary depending on the method used to separate the pores, which in turn may affect the simulation results. For this reason, the authors should add an explanation of the method used to separate the pores and discuss the impact of the pores separation method on the simulation results. Finally, based on comparisons of the measurement results of water retention curve, they need to discuss the accuracy of the simulation results, including the best method for separating each pore.
<Specific points>
L27: I disagree because some studies evaluate results based on observations (e.g., runoff from snow cover) or comparisons with laboratory experiments. If they insist on their claim, they should point out the problems with past studies in more detail.
L117: If gravity is ignored, isn't it impossible to calculate the water retention curve? This is because gravity is included in equation (5).
L126: Since cos12 is 0.97 and cos0 is 1, I don't think this will have a significant impact on the result. However, I think it would normally use 0 degrees, so what is the basis for using 12 degrees? Please provide the reference.
L131: As shown in L117, this model ignores the effects of gravity, so the simulation may not accurately reproduce the actual drainage process. How do they feel about this point?
L249: Because nVG is said to depend on the distribution of pore sizes, the results of this study may be due to the method used to separate the pores. Therefore, this possibility needs to be discussed.
L321: This argument is meaningless. This is because the data used for comparison are simulation results, so it is only natural that the results of this study are the most appropriate results.
L353: As shown in Table 2, the parameter settings used in this study differ from those of Yamaguchi et al. (2012). For example, the αVG value in Yamaguchi et al. (2012) is two orders of magnitude larger than the value used in this study.
On the other hand, as shown in Figure 10, the calculated unsaturated hydraulic conductivity using the parameter settings in this study shows relatively good agreement with the results of Yamaguchi et al. (2012).
This result may suggest that the parameterization of the VG model may not have a significant influence on the unsaturated hydraulic conductivity.
I propose adding further discussion on which parameterization of αVG and nVG is more important in determining the unsaturated hydraulic conductivity. This discussion is considered useful for determining which parameterization requires higher accuracy.