The manuscript deals with a classical 1D analysis of the unsaturated flow stability during a wetting event. The numerical treatment of the problem at stake, although quite standard and therefore with no new insight, appears to be correct (I went through it, and I didn’t find any error). Instead, the English usage requires a very solid and sounding proofreading. Besides this (marginal) aspect, the main issue which I see with the manus is of methodological nature. In particular, my skepticism is two-fold. First, accounting for the gravity solely (thus neglecting the impact of retention) may be unrealistic especially if one is interested (as it is usually happens in the applications) on the “onset” of the stability vsinstability”. Second, the authors have carried out a long and intensive analysis of the flow rates which make (or not) stable the flow, but what and where is the stability analysis? The very new and innovative insight could have been a Touring analysis of the stable/unstable flow patterns in order to highlight which ones are those parameters (and perhaps the infiltration rate is the most important one) that regulate such a stability. Instead, the manuscript, as it is, is nothing more a numerical analyses (followed by an experimental benchmark), quite similar to many others, already existing in the literature.

The paper introduces a hysteretic model that can generate unstable infiltration fronts and the resulting fingering patterns. The model is tested by reproducing highly idealized experiments in two-dimensional sand tanks, with very limited additional analysis.

The model and its equations are not terribly well explained, which creates some confusion. The model test on experimental data is convincing, but there is no follow up. There is no indication if the model is of any use to apply to natural conditions or to address real world problems. This gives the impression of only half a paper: the presentation of a model and its validation, but no meaningful application.

Major comments (also at the top of the attached annotated manuscript):

The model is not new, and the simulations cover a very limited scope, merely attempting to reproduce laboratory experiments under conditions that have no relevance for the field. This strongly constricts the overall contribution of the paper.

Although the authors discuss the literature at length throughout the paper, they seem to have missed many relevant papers. In the detailed comments I provide numerous references that may be useful to turn this paper into material that deserves publication.

The English needs to be improved a lot. I had a hard time understanding important sections of the paper. In combination with the inadequate explanation of the equations and their variables, it is not really possible to grasp the theoretical elements of the paper.

The model as such is interesting, and could be part of an interesting paper if it is better explained. To make the paper complete, the model would have to be used to carry out simulations that have some relevance for the real world, not the world of 2D tanks filled with air-dry artificial porous media.

The Results are tedious to read.

Too much literature is discussed in the Results section. It slows down the pace and increases the wordiness of the section. The literature review belongs in the Introduction. 

  

You only appear to pursue to reproduce experiments in 2D sand tanks under highly idealized conditions that have no relevance for field scale problems related to preferential flow.

What is the point of developing a model that cannot be applied to address problems in natural and agricultural soils, but is entirely focused on experiments with artificial porous media that were carried out to reveal the physical mechanism driving wetting front instability, so that the science could move on and make this knowledge operational for field problems?

You do not convince the reader why your model matters. Confirming experimental results only serves to show that you model works, so you can apply it to problems beyond experimental reach. Instead, you made it the main point of the paper.

Figure out what the real contribution of the paper is, pick the most interesting results, and discuss these intelligently, with an eye to theoretical advances as well as relevance for the real world. Rigourously cut away the obvious stuff. This will require a whole suite of new model runs, new figures, and a more lively discussion. 

You present a state-of-the-art hysteretic model, but use it with outdated soil hydraulic functions. Preferential flow under natural conditions will manifest itself during a season with many cycles of rain and evapotranspiration. This is a wonderful playground for a model that can handle hysteresis, but you have not considered even one cylce. This is one of many missed opportunities. See my comments on the Results Section for more thoughts on this.

I read with great interest this manuscript which reports about the attempt at describing the fingering process in a two dimensional porous medium by means of a semicontinuum approach. The article is rich in the literature and in the analyses. As a general recommendation I suggest:

1. To describe with some more details the background assumptions and the model;

2. To perform a more detailed comparison (even quantitative, if possible) between previous experiments and the numerical findings

In the followings I detail a little more the questions arised from the reading of the paper, that I recommend to address in the review.

The Authors represent a soil section with three main assumptions:

1. The SWRC is represented by a mixed form of a classical van Genuchten's SWRC (in an hysteretic form) with a minimal Prandtl's hysteresis operator. The merge of the two curves is ruled by the dimension of the simulation cell. This idea is close to that the experimental SWRCs are sensitive to the dimension of the laboratory sample, being more flattened as soon as the sample dimension is reduced. In this sense the Prandtl operator is seen as the minimal behaviour of a single capillary. Then they fix the cell dimension and accordingly a SWRC for the soil is found.

I recommend to better evidence the reason for the choice of the cell dimension, and to express with more detail whether the simulations are expected to be sensitive to the cell dimension (as the SWRC seems conditioned to it).

2. The intrinsic permeability is represented by means of a stocastic simulation. It is really needed? Probably it seems not, as the Authors report in the discussion, but it anyway could mix the effects of hysteresis with the effect of the permeability field on affecting the formation of preferential pathways. Have the Authors considered the possibility of performing separate simulations to evidence the relative importance of hysteresis vs permeability field?

3. The relative permeability is represented by means of a classical power law / or van Genuchten Mualem shape.

With these hypotheses, it is assumed that the conservation of mass with the Darcy--Buckingham law, admits the formation of a saturation overshoot in the 1d form, and the formation of front instability in the 2d form (as it does in the simulations). The obtained equation is in my opinion a form of the Richards equation with hyteretic SWRC, because, as far as the soil is discontinuous and we aim at representing its properties by means of descriptive index properties (as the saturation or the porosity are), we implicitly admit that the Richards equation is locally defined by means of average values on a certain small domain, which is commonly referred to as the Representative Elementary Volume -- even if not defining its dimension. What moves the present model from the Richards equation is the choice of the dimension of the cell, which rules the behaviour of the SWRC.

Regarding the formation of an overshoot the Authors refer together to two kinds of overshoots which are very different in their meaning. In fact an overshoot related to a pulsation in the Dirichlet boundary condition is in any case in agreement with the parabolic behaviour of the Richards equation, as it can be in any case framed within the maximum principle which guarantees the solution of the parabolic operators. Consider, for example, the Stokes problem of the velocity profile in a seminifinte steady fluid with pulsating wall.

Another is the case for which the maximum principle is not proven, as it can happen in some cases in which the soil is not homogeneous (see Barontini et al, 2007, WRR). This is the intriguing case which, according to the Authors, may lead to the instability. Is it possible to check the applicability of the maximum principle for the investigated case? I mean: the maximum principle is unprovable in the equation? or it is provable in the equation but it can be put into discussion as a consequence of the unhomogeneities of the soil permeability?

The results are interesting, particularly as the model describes a concentration of flow in the most saturated soil, even in the case in which the fingers are not developed. This is in agreement with the fact that the relative conductivity drops down as soon as the soil is not completely saturated. This behaviour is most evident in organic soils, where van Gencuchten's n is small (smaller than 2), but it is evident in any porous medium, also in this case for n between 6 and 8.

Yet one may argue whether the results account for a physical behavious or for a mathematical description intrinsic to the model. This is why I recommend to the Authors (1) to better focus on which is in their opinion the physical source fo the instaibility, whether it is the hysteresis or the variability of the permeability. (in this case it cannot be the presence of macropores, as -- if I properly understand -- macropores are not describerd in the model) and (2) to provide closer comparison with literature experimental results.

Before closing, I add some minimal notes:

-- The Authors express the conductivity of the cell as the geometric average of two conductivities at different water content (the minimum and the maximum, it seems, but it shoud be probably better enlightened). Where does this scheme come from?

-- Mualem's parameter \ell is usually set at 0.5 (despite it can be changed, if needed): why did the Authors set at 0.8? Does it come form a fit of an experimetnal conductivity curve?

-- Figure 1 is not very clear, I recomment to add the direction of the potential axis

-- eq. 3: why 0.5?

-- l.226: please better detail the meaning of residual and initial saturation as residual is greater than initial

-- How the bottom boundary condition was chosen? If the column is suspended a seepage face condition would have been more realistic, yet it requires soil saturation before the leakage starts;

-- l.443: it is referred to stabilized flow as a percolation flow (i.e. with null spatial rate of change of the tensiometer--pressure potential), but due to the boundary conditions it might not be a pure percolation flow.

Thank you for the attention.