This paper is of high quality, discussing various implementation strategies for Richards equation in snow cover models. It focusses on an improved treatment of the dry limit of the equation, as well as on a tight coupling with phase changes. This indeed addresses one of the outstanding questions in the field of snow modeling. Some simple test cases are run, to demonstrate the various approaches and setups. The writing is mostly of high quality. I generally think that the paper can be published, after taking my minor considerations into account. 

I have a few issues with somewhat more important feedback. 

1. I would like to see the discussion of existing literature improved. It is too focused on the recent studies, and too focused on the SNOWPACK and CROCUS models. Examples: 

• L31: "it has been proposed in the last decade": I think the earliest implementation of Richards equation in a snow model I am aware of is Jordan 1983 (https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/WR019i004p00979), but there is also work by Colbeck, Illangasekare et al., 1990. I think Daanen and Nieber (2009) were also using Richards equation in their model, before SNOWPACK. I think the paper should provide a bit more historical context, even though there is not a need to do an extensive literature review. It's more that I think it's important to provide a proper historical perspective and context. 

• The sentence in L61 "Rather, snowpack models rely on a ... 1D framework" is too generalized and not acknowledging the work that has been done extending modeling to 2D and 3D. For example the work by researchers Webb, Leroux, Hirashima. I would like to see the work of those researchers, and maybe others I overlooked, included in the discussion. 

2. Section 2.2.1: I think it is critical to discuss that this discrepancy in the dry limit is partly because drying water retention curves are being used (i.e., water retention curves derived from drying snow samples). In a wetting snow sample, the water retention curve may actually nicely approach 0 residual water content. Hysteresis has been used in snow modeling in other studies: Leroux and Pomeroy (2017): https://www.sciencedirect.com/science/article/pii/S0309170817300040. 

3. I think that the paper currently insufficiently discusses that implicit methods can be numerically stable, but still that doesn't mean that they are accurate. I think this is important to convey to the reader. So for example, I would not write in the conclusions: "which can run with timesteps of 15 min and above", without immediately making a remark about the accuracy of the solution. The reason I'm saying this, is because in Fig. 4, it looks like that the most advanced models perform worse at the largest time steps. That is just something to be very careful about. I also wonder if it is not better to show two Figures for Fig. 4: one where each numerical approach is tested for time step sensitivity, by taking 1s simulations for each of the approaches as a benchmark. And the other the current Fig. 4, showing the estimated accuracy under the assumption that Model 1 at 1s timesteps can serve as a benchmark. That would also more clearly deal with the fact that if model 1 run at 1 second time step is considered the reference for determining error-statistics, like RMSE, it is quite obvious that model 1 performs best and creeps closer and closer to the reference simulation when the time step reduces. I wonder if Models 3-5 actually just have another solution, and that for that reason, their performance is relatively constant, compared to the benchmark. 

I also have a few minor comments: 

L68: "effective" may require a definition in this context. 

L79: "neglecting the influence of water vapor": isn't the effect of water vapour included in lambda, as is common for snowcover models? 

Section 2.1: I found this section somewhat overly complex. I don't fully understand why the possibility of phase changes is not directly included in Eq. 1. Now, Eq. 5 is basically incompatible with Eq. 1, because melt is introduced at a later stage. I think this section could be simplified in this regard. 

L162: "as snow below the fusion is considered dry". In fact, an argument that is sometimes used is that snow contains a thin water film, even below zero. And that for that reason, it can be assumed that theta is never truly 0 in snow. I'm not sure what the authors think about this argument. 

L165: If one would set the hydraulic conductivity to 0, it would be possible to suppress any liquid water flow and tiny liquid water amounts. So it is not a given that there is a tendency for percolation, I think. 

L204: "for melting refreezing event in" is a somewhat broken phrasing 

L205-206: This is actually not the case in SNOWPACK. Maybe it was removed at some point, since Bartelt and Lehning 2002 is cited. But I'm quite sure that phase changes translate volumetric contents between the ice and water phase. If models refreeze by increasing volumetric ice content, and melt by reducing element length, repeated melt-freeze cycles basically generate artificial ice layers, because it is an inconsistent approach. The current approach and reasoning in SNOWPACK is exactly what is described in L209-L212. 

L219: "Replacing Eq. 5 ..." is somewhat broken phrasing. 

L228: "we discretize as well" not sure this is proper phrasing 

L330: "when on the of the cell" needs rephrasing 

L346-347: "get stuck in cycles". I would say instead of the algorithm diverges, the solution diverges. And instead of get stuck in cycles, I would say that the solution oscillates. 

L361-362: This sentence is a bit unclear. I assume that this is not shown in the paper, so I would write: "This did not modify the results (not shown)." To signal to the readers that this is not further discussed.