The authors present, describe, and evaluate a simulation module for dense flow avalanches (com1DFA) implemented within the newly developed, flexible Avaframe mass flow simulation framework. After a detailed presentation of the physical basis and numerical implementation of the tool, it is verified through two benchmark tests against semi-analytical solutions on a simple topography and a dam break scenario. The work is highly relevant from both a scientific and a practical perspective, and is certainly suitable for the journal. The manuscript is well-organized and well-written, with appropriate and informative illustrations. The methodical part targets very much at an audience with a strong mathematical background (which is fine for this type of model description paper) but will be hard to follow for anyone else. I would certainly like to see the work published. Before, I recommend some minor revisions. Please find below my comments and suggestions.

L23f: “… combining over a decade of operational application …”: This formulation is not very clear, maybe better mention that this operational application refers to SamosAT?

L42: “… proposed in Mangeney-Castelnau et al. (2003), …”

L64: It is not clear to me how the word “AvaFrame” fits in this sentence. Further, the “F”in “AvaFrame” is written in upper case here and in other places in the manuscript, but in lower case in the title of the paper. Please ensure consistency, or explain why there is this difference.

L71: Maybe I overlooked it – did you introduce the abbreviation “DFA” before? It is very clear what it means, but still it has to be defined.

L144f: Please check brackets in the reference to the Voellmy model.

L255: I am not a native speaker, but “… using such an SPH method …” might be better.

L256: “… previously described grid method …”

L311: “… used in the com1DFA code …”

L321f: Shouldn’t the DEM also be described by the number of cells in x and y direction?

L400: Maybe better: “… a (semi-)analytical solution …”. But I am not sure, it is a tricky issue.

L411: The similarity solution and dam break tests come very sudden here (not mentioned before). Maybe better reformulate this sentence a little bit.

L437: Is this really a dam break? I would rather describe it as the collapse of a cliff. Shouldn’t a dam break imply that there is some fluid behind, which is released?

L474: “… between analytical and numerical solution …”

L480: “… with an α exponent …”

L481f: “… α = 0 for the dam break test …”

Figure 4, header: “Similarity solution test …”

L493: “… only Coulomb friction …”

L541: “… and questions like …”

L562: “The results of these tests …”

L569: “at the date of publication of this manuscript”: this formulation sounds somehow strange to me, better reformulate.

L590: Better: “… throughout the development.”

L591: “Note that the computational efficiency …”

L596f: “feed back” --> “feedback”

Content, suitability for journal, interest to audience

The manuscript describes the set-up and numerical solution techniques of an open-source model for simulating the flow of dense-flow avalanches geared at use in practical problems (hazard and risk mapping, design of protection structures), and it also presents validation cases. Comprehensive documentation of such a code is a pre-requisite for applying it in contexts where the simulation results can have legal consequences, such as issued or declined building permits, decisions on large investments in protection structures, dams or extra reinforcements in buildings. Having such documentation undergo peer-review is valuable and in the best interest of potential users of the code. Given the focus on code implementation and performance, this journal appears to be a good match for this manuscript. Conversely, the manuscript is suitable for Geoscientific Model Development and will attract the interest of a part of GDM's audience, even though there is no scientific break-through and there have been many papers on similar models in other journals over the past decade.

Presentation, language, referencing

The paper is clearly structured, and with a few exceptions indicated below (Minor remarks) or marked in the manuscript, all points are suitably explained. At a few places, the mathematical derivations are on the explicit side and could be shortened. The manuscript contains a fair number of formulas, which employ suitable and intuitive notation and are typeset carefully. Generally, it is apparent that the authors have used a considerable amount of care in preparing the manuscript.

In a large part of the manuscript, the English is good. There are, however, passages with grammar errors and at times rather clumsy and even ambiguous formulations. The annotations contain some suggestions on how to improve this. Nevertheless, the authors need to critically re-read the text and bring these passages up to the level of the rest.

The referencing appears to be adequate, and there is no excessive self-citing. I also appreciate that the authors are aware of the pioneering work of the group at Moscow State University as early as the 1960s. It could be useful for the readers if the authors cited the most relevant open-source models for geophysical mass flows, like Titan2D, Clawpack and its derivatives or MassMov2D, and perhaps summarized the salient differences from their approach (from a user perspective).

Major remarks

1. To the best of my understanding, the geometric effects are rendered correctly in the formulas. However, the explanations given in the text appear confusing if not misleading, particularly in Appendix B2.
   
   It is crucial to distinguish between the terrain surface, which is a 2D manifold S embedded in 3D Euclidian space E₃, and the 2D Euclidian tangent spaces T_x, one associated with each point x on the manifold at a given time. In some contexts, these 2D tangent spaces are extended to 3D (Euclidian) spaces T_x3 in the manuscript; this can be done naturally, but these spaces are not the same conceptually as the embedding space E₃. The fields h, u (and also the stress tensor σ) are not functions from the tangent plane, as stated in Appendix B2, to Euclidian space or R but from the surface S to the non-negative real numbers (h), the T_x or T_x3 (u) and a 9D space for the stress tensor. The appropriate modern mathematical concept to use would be the fiber bundle. I think that the fiber bundles in the present situation are essentially trivial, provided the fibers for the field h are limited to the interval [0,1/κ_max[, where κ_max is the maximum of curvature (in any direction) of the base manifold S (the topography). The same results can, however, also be obtained with classical differential geometry. A fundamental result is that the whole theory can be developed without reference to the Euclidian space E in which S is embedded, but the price to pay is using the metric tensor and its derivatives, from which the connection coefficients (Christoffel symbols) can be derived. Most models known to me that tackle at least some of the intricacies of flows over curved surfaces use this approach, defining a coordinate system on S and letting the tangent space at a given point be spanned by the tangent vectors to the coordinate lines through that point. These coordinate systems cannot be orthonormal unless S is flat.
   
   In contrast to this, the authors explicitly use the embedding Euclidian space E₃, which allows to give absolute orientation to all tangent spaces (and their 3D extensions) associated with the points on S. The "natural coordinate system" (NCS) used by the authors is not a coordinate system on S or in _E3 but a separate coordinate system in each of the T_x3; moreover, the NCSs are time-dependent, being advected along with the "particles". The choice of spatially and temporally varying coordinate systems in each of the extended tangent planes would be a computational nightmare. However, the authors cleverly exploit the fact that, in essence, only the basis vector v₁(x,t) is needed to compute the gradient of the kernel W.
   
   To be clear: I do not criticize the authors' choice of dealing with the geometry—it does indeed have several attractive properties. However, for the readers' sake I wish they clearly distinguished between base space and target space of the functions and briefly explained the difference with traditional approaches.
   
   There are also some points that are important for understanding the method but are not explained as far as I can see. In the mass-balance equation (1), dV(x,t)/dt = q(x,t)/ρ dA, and Eq. (5) implies dV = h(x,t) dA. This shows that dV/dt = dh/dt dA + h dA/dt, i.e., in order to find the flow depth one must know how the basal area of the SPH particle changes over time. To this end, one will need the 2D divergence of the velocity field. I think I can see how this can be obtained by expressing the basis vectors at all points in terms of the global Cartesian coordinate system, but I cannot find this explained in the manuscript.
   
   In this context, I am also surprised that the terrain curvature (more precisely, the Gaussian curvature) only briefly appears in a formula after Eq. (4), even though it plays a significant role in the equation of motion later on, where it is masked as dv₃/dt. I think the mathematical development would gain clarity if the connection between dv₃/dt, ∇v₃ and curvature κ were stated explicitly.
2. Given there are many good Eulerian mesh-based solvers for hyperbolic equations and that com1DFA seems to be fairly complex and not in the top league when it comes to speed, the authors ought to explain why they chose this approach. Moreover, it is known that boundary conditions are not straightforward to implement in SPH codes, thus at least a brief discussion of how this is done here (or avoided if the particles' kernel areas never approach the domain boundary) should be included.
3. There is one test of the model I would like to see, namely impact on a wall or, if the code does not implement impenetrable boundaries, run-up on a steep counter-slope. Such a situation occurs in many practical applications. This serves to test the shock-capturing capabilities of the code, as it has to create a backward traveling shock. The simplest form of this test is for an inviscid jet traveling horizontally with a prescribed speed and flow depth before hitting a vertical wall. In this case, the shock speed and height can be obtained analytically. In a test with friction on a slope, the shape of the final deposit would also be of interest because many practically used models of Voellmy-type (of which the Coulomb model presented here is a special case) tend to creep (sometimes intermittently) for very long times until the surface slope of the deposit nowhere exceeds the bed friction angle.
4. Section 5.2 (except Sec. 5.2.3) takes a rather large fraction of the entire manuscript, even though the theory behind these tests has in its essence been known for 40 years or perhaps even for a century. Given the results of the tests presented in Sec. 5.1, I would have been surprised if the energy-line tests had shown poor results. This part could easily be shortened.
5. Section 6 is titled "Conclusions", but except for lines 591–595 it is simply a brief summary of the preceding five sections. Readers who have read the paper, will be bored, and those who just read the abstract and skipped the rest will not be enlightened either. For this section to be interesting, the authors could try to answer some of the following questions: 1) For which type of application is com1DFA particularly suitable, for which is it not suitable (at least presently)? 2) How does com1DFA compare to other codes "on the market"? 3) What are the most interesting, useful and feasible extensions of com1DFA?
   
   Lines 591–595 contain relevant information with regard to practical applications. The authors state that "simulations compute ... within minutes on an 8-core laptop." This suggests that the code uses some degree of parallel processing (for which SPH is well suited). From my own and my colleagues' experience, I consider this to be a bottleneck because consultants often run perhaps ten simulations per path to test sensitivity to the choice of parameters and hate to wait for an hour to get the results. Other codes like NGI's MoT-Voellmy achieve such simulations on a single core typically within seconds to tens of seconds (admittedly at less precision, but this is of little concern if the initial conditions are as poorly known as usually is the case).
6. In Appendix A, the authors state the mathematical formulation of the entraiment models that have been implemented in com1DFA, even though they are not used in this manuscript. This is valuable as a documentation of the present state of the code. For this reason, a few critical remarks concerning these formulas may be in order:

• In Eq. (A2), neither ρ_ent nor e_b are defined.
• The expression (m_k/(ρ₀h_k))^1/2 seems to be used as the width of the SPH particle. This would seem to apply only to a square-shaped particle.
• The formula for the entrainment rate due to plowing uses the flow velocity. It has been pointed out already in the 1960s (Eglit, Grigorian, Yakimov) that the entrainment front propagates faster than the flow. It would seem more physically justifiable to create new particles from the plowed snow cover, with the rate of particle creation accounting for the shock propagation velocity whereas the new partcles themselves move at the lower speed of the flow front. Another weakness of the proposed formula is that the erosion depth h_ent nust be specified beforehand; on physical grounds, it should be determined dynamically. The formula does not account for the shear/compression resistance of the snow cover either (see (Eglit et al., 2020) for more details).
• The formula for scouring erosion has similar shortcomings in that it does not account for the shear strength of the snow—unless this is somehow included in e_b (if so, this needs to be explained!). I do not understand the physical reasoning behind the factor ||u||—based on the arguments in (Issler and Pastor, 2011) or (Issler, 2014) I would expect 1/u.

Minor remarks

In the attached commented manuscript, there are about two-hundred remarks that the authors should look at when revising the manuscript. I would, however, ask them not to reply individually to each of them, irrespective of what the journal's policies are. Instead, please comment only on the most important remarks you do not agree with! The vast majority of the remarks concerns language issues. Among the remainder, the most relevant are the following:

• The authors state that they use prismatic control volumes (or particles). This would lead to "crevasses" in the flow body wherever the terrain curvture is negative, and to "caves" (or, alternatively, interpenetration of parts of the control volumes) in areas of positive cutrvature. Perhaps they mean truncated pyramidal volumes instead?
• The description of the SPH method is rather confusing in my eyes, as it is being said that the value of some quantity at the location of particle j is the weighted sum of that quantity over all other particles. But the same quantity at the location of another particle depends, among others, on the quantity of the first particle. In the extreme case with only 2 particles, we get
  
          f₁ = w₁₂f₂   and  f₂ = w₂₁f₁ = w₁₂f₁ = w₁₂²f₁   ⇒   w₁₂ = ±1   ⇒   f₂ = ±f₁,
  
  which is hardly what is intended. The way I understand SPH is that it provides a way of interpolating a field value at any location from the corresponding quantity associated with a finite number of locations (i.e., particles). This seems to be supported, e.g., by the documentation of Abaqus, which states:
  
  "At its core, the method is not based on discrete particles (spheres) colliding with each other in compression or exhibiting cohesive-like behavior in tension as the word particle might suggest. Rather, it is simply a clever discretization method of continuum partial differential equations. In that respect, smoothed particle hydrodynamics is quite similar to the finite element method. SPH uses an evolving interpolation scheme to approximate a field variable at any point in a domain. The value of a variable at a particle of interest can be approximated by summing the contributions from a set of neighboring particles, denoted by subscript j, for which the “kernel” function, W, is not zero:
  
          ⟨f(x)⟩≃∑_j m_j ρ_j f_j W(∣∣x−x_j∣∣,h)."
  
  I am not very familiar with the SPH method, but I think other readers might stumble over the same issue; therefore, I ask the authors to consider this naïve objection and to make the description as intuitive as possible for a non-specialist audience. In particular, it would be important to explain at which stage this averaging over neighboring particles is carried out—I could not find any mentioning of it in the text.
• The authors should explain why they think that the convergence criteria developed by Moussa & Vila (2000) also apply in the case at hand.
• I think that the dam-break solution (Fig. 6) can be applied at the upstream end of the released mass as well (the solution consists of a dam break towards the left on a horizontal plane with gravity reduced by a factor cos θ, superposed with a uniformly accelerated downstream motion). Doing so would provide yet another test of the code.
• In the reference list, please include DOIs where available.

Assessment and recommendation to the editor

Even though the manuscript does not score particularly high in terms of scientific innovation, it serves an important purpose for the community of avalanche researchers and practitioners as it clearly explains the foundations, and demonstrates the performance, of an open-source model meant for practical use. Therefore, I recommend acceptance for publication.

While I have not found any serious errors, there are a few conceptual points that need to be straightened out, and the manuscript needs small adjustments in a multitude of places. This revision needs to be done comprehensively and carefully. If the editor checks the revised manuscript in this regard, a minor revision without another round of reviewing will suffice.

April 11, 2023

Dieter Issler