This well-written paper combines theory, computational methods, analogue laboratory experiments and field observations of cohesive debris flows.

The paper presents a depth-integrated model for cohesive-frictional debris flows. The model resembles lubrication theory, in that inertia of the fluid is neglected, and the depth-integrated flux is therefore an instantaneous function of the local surface gradient and thickness.

These model equations are solved using a control volume finite-element technique, which is validated by comparison to exact solutions. The model predicts very well the deposits of laboratory experiments of sand/kaolinite/water flows, and is compared with with field observations from Coussot et al. (1996).

The agreement between the model predictions and laboratory experiments is striking, and I have no doubt that the modelling approximations made (shallow inertia-free flow with homogeneous cohesive-frictional rheology) are very well suited to the flows in the lab. As a description of laboratory experiments it is therefore a strong piece of work.

However, it is much less clear to me that the physics studied here is relevant to many natural debris flows. There is some acknowledgement of the differences between the modelling here and natural debris flows (lines 52-58 and 441-446). But in my view there should be a much more comprehensive discussion and justification of which predictions from this paper (obtained at laboratory scale) would be expected to hold true at field scale, and which would not. For example:

1. It is not clear that yield stress is responsible for the blunt snouts of natural debris flows. In the model used by this paper, the yield stress required to produce blunt snouts at field scale is very large, evidenced by a fit of τ_Y/(ρ g) ≈ 0.5m in the field observations, compared to ~0.001m in the experiments. How can this difference of a factor of ~500 in yield stress be explained?
   
   It seems likely to me that the formation of blunt snouts at field scale could due to a different process, for example the loss of excess pore pressure at the flow front, resulting in a substantial increase in the frictional part of the stress here. That is, the blunt snout could be due to a rheology that is inhomogeneous but not cohesive. Section 8.1 would benefit from some discussion of these points.
2. More generally, the difference in physics between small-scale analogue experiments and natural-scale flows has been raised by several authors, including Iverson (e.g. https://doi.org/10.1016/j.geomorph.2015.02.033). The paper would benefit from a discussion of the parameter regime realised in experiments (Froude number, Reynolds number, Savage/inertial number etc.) and a comparison of this with natural examples.
3. An important feature of the model is that the inertia of the flow is neglected (line 76). Is this really valid for natural debris flows? Many debris flows are supercritical and exhibit features such as shocks, roll waves and superelevation in curved channels, which require inertia. Is there evidence that inertia can be neglected at field scale?
4. It is tempting to attribute the similar conical shapes of the natural debris flow fan (figure 1) and experimental deposits (figure 9) to a similar formation mechanism. Though I do not know the 2009 Xinfa debris flow shown in figure 1, the inundation of houses in this figure suggests a flow of perhaps ~2m deep occurring on a much larger (perhaps ~30m tall?) pre-existing debris-flow fan. If so, this is clearly a completely different mechanism from the en masse deposition of the entire fan in the experiments. There is some acknowledgement of this around line 55, but in my view a much clearer statement is needed as to the differences between the modelling/experiments in this paper and the natural deposit in figure 1.

Minor points:

•  Equation (1): There is no source term in this equation corresponding to the inflow. Is the inflow flux Q_in modelled as a source term of limited spatial extent on the right hand side of (1)?
• Equation (3) / line 83: I initially misunderstood the statement on line 83 and believed that S_c was a constant for a given material. It may be useful to make it clearer that S_c is dependent on local instantaneous flow depth and free-surface slope, and that this dependence is derived in section 4.
• Equation (4): how is this equation derived? From (3), it is clear that the free surface slope is no greater than S_c, but not clear to me why it is exactly equal to S_c. (Derivation of (4) must require some constraints on the initial conditions or inflow functions Q_in, as for certain choices of these are counterexamples to (4). For example, if z_b(x,y) = 0 and Q_in(x,y) = k and the initial conditions are \tilde z = 0 at t = 0, then the exact solution is \tilde z = k×t, which does not satisfy equation 4.) In a recent paper (https://doi.org/10.1017/jfm.2021.1074, section 7) we referred to regions of a dry granular flow deposit that do satisfy (4) as "maximal", but this was not true of the entire deposit.
• Figure 5 / line 209: What is the source function Q_in for these solutions? Presumably the source is at a different value of x for each flux?
• Line 235: I don't fully understand "an approximate analytical solution obtained by setting tan β = 0": why does tan β=0 allow an analytical solution, and what exactly is being compared (is a numerical solution with β ≠ 0 compared to an exact solution with β = 0?)
• Figure 6: Should "(b,f,j) transverse deposit profiles" read "(d,h,l) transverse deposit profiles"?
• Table 1: What is the order of convergence of the numerical scheme in space and time? From the time discretisation in equation (8), it appears to be first order in time. Is it also first order in space?
• Figure 10 and 11: these contour plots are noticeably slow to plot in my PDF viewer: are they particularly large figures that could be reduced in resolution?

Chris Johnson, University of Manchester

In this paper, a theoretical and numerical model for the morphology of the deposits of debris flows is presented. As a main simplification compared to existing models, effects of inertia are neglected. While existing models are based on shallow-water type (Savage-Hutter) equations, this approach arrives at a nonlinear diffusion equation with a threshold slope. For validation, analytical solutions, topographies of real debris flow deposits, and laboratory experiments (being a part of the study) are used.



First, I would like to emphasize that both the theory and the numerical implementation are described very well and in great detail. Since the  diffusion equation is numerically not very challenging, one might even ask whether such a detailed and basic level is necessary. However, I do not complain about this.



My main criticism concerns the simplification by neglecting effects of inertia. As stated by the authors, this limits the applicability to low  velocities. The question whether this is a serious limitation for the application to real debris flow is not addressed sufficiently. All results used for validation are solely based on the final final topography and thus on the very end of the movement when the velocities should indeed be small. On the other hand, the introduction starts from the hazard of debris flow, where the runout length is more important than the morphology of the deposits. So the authors should point out more clearly that the referenced existing models also attempt to predict the runout also at high velocities, while this is not tested for the new model. It even looks as if the new model mainly constructs a final deposit topography that obeys a predefined relation between slope and thickness.

As a second point concerning neglecting effects of friction, I am not fully convinced that it makes things simpler or more efficient. It is stated that  the existing models require a large amount of input data. However, can go back to the original Savage-Hutter equations with a simple static friction term and nothing else. Then the coefficient of friction would be the only model parameter. We could also go a step further and use the Mohr-Coulomb criterion as proposed in the recent manuscript. The number of parameters and their meaning would be almost the same in both models then. This scenario would allow for an assessment of how much we lose by neglecting effects of inertia and how much we save. Theoretically, we save much because the equations become simpler. However, the results about the computational performance given in Table 1 are disappointing. It seems that the diffusion model model with the explicit time step requires very small time increments. Without having data for comparison available, it looks to me as if the new model was quite inefficient compared to existing models.



To summarize these points, it would be essential for me to see a thorough analysis of what we lose with regard to real debris flow with the new model and whether there is any increase in numerical efficiency.



Provided that this can be done, I would also suggest to consider the following aspects:

Section 2: If I got it correctly, the flux is only dependent on the slope, but not on the thickness (above a minimum thickness). This means that the flow velocity increases with decreasing thickness. I would have rather expected a flow velocity that depends on the slope only. I guess that the rather high fluxes at low thickness arising from the approach used here are not very good for the numerical performance. Is there a specific reason for this approach?



Section 3: Rather for curiosity (since I am not an expert on this): Why did you not use a standard Delaunay triangulation in combination with Voronoi polygons as control volumes?



Equation 10: How did Qin come in here compared to Eq. 8, and what is it used for? I thought you start the simulation with a given thickness distribution. Or is it just the source term for reproducing the laboratory experiments?



Figure 5: If the deposit thickness H is measured at the apex, I have some difficulties in relating the values to the legend.



Section 7: I am not convinced that the comparison with analytical solutions should be considered so extensively. These comparisons only illustrate that the numerical implementation of the model works and have nothing to do with the applicability of the model. So the excellent agreement should not be stressed too much.



Anyway, I enjoyed reading the manuscript and like the approach in principle, despite my criticism.



Best regards,

Stefan Hergarten