Dear editor,

thanks for having given me the opportunity to review this interesting manuscript. I believe this is a very well written paper that presents a significant advance of IMEX_SfloW2D.

General comments

This new version now allows simulating dilute pyroclastic density currents (DPDCs) with a good balance between the degree of simplification, necessary if one wants to carry out probabilistic volcanic hazard assessment with the currently available computational resources, and accuracy of the physical treatment of the fundamental processes characterizing DPDCs.

The model is well presented, with a robust treatment of the equations and the numerical methodology and several very well-documented test cases.

I therefore recommend this paper to be accepted for publication in GMD, provided the authors conduct minor/moderate revisions addressing the following points. I also attach a commented version of the pdf document, where the authors can find my specific comments at the position they refer to and some other minor comments.

Specific comments

• The authors should better explain how a shallow-water model, which by definition cannot solve for vertical gradients of the flow properties, can be employed to simulate DPDCs, which are strongly vertically-stratified (particle concentration, flow density, flow velocity). The authors clearly states in the Conclusions section that the vertical stratification is not taken into account. However, a couple of years ago Biagioli et al. (2021, https://doi.org/10.1016/j.apm.2020.12.036) published a paper on a method to take vertical profiles into account in shallow-water models. Hence, unless I misunderstood this paper completely, I expected this method to be implemented into IMEX v2 in view of its application to DPDCs.
• In the Introduction section a significant space is dedicated to the effect of entrainment of ambient air into the DPDC. Indeed, entrainment is a key factor controlling DPDCs, however the role of sedimentation should be equally emphasized in this introduction since, like entrainment, it controls the rate of change of the flow bulk density and, hence, its existence. Indeed, some words should be spent in emphasizing the concept that a DPDC, as any other density current (e.g., turbidity currents), exists until there is a density contrast with the ambient in which it flows.
• The Github repository lacks of a User manual. This is a very important document for users.
• Sometimes the simplifications and assumptions of the models should be better explained, in particular their consequences. For example, in paragraph 154 it is stated that the conservation of total energy is not taken into account. What are the consequences of this simplification?
• Eq. 14. Here the assumption is made that, everything that is sedimenting is also creating a deposit. This is not always necessarily the case, since sedimenting particles (in fact, I would use a sedimentation term in the conservation equations, instead of a deposition term) can be moved as a traction carpet or even re-entrained. This assumption is not explained in the manuscript.
• IMEX should not be restricted to radially evolving flows, however I think this should be stated more clearly in the manuscript.
• Entrainment model. Eq. (17), quantifying the coefficient of rate of entrainment of ambient air into the flow, is from Parker et al. (1987). This paper is not cited here in the manuscript (https://doi.org/10.1080/00221688709499292). Furthermore, there are other relationships available in the literature, e.g.  https://doi.org/10.1029/2003JF000052,  https://doi.org/10.1029/2019GL084776. While these should not necessarily be implemented into the code, they can be at least cited perhaps in the introduction when you first introduce the entainment.
• Paragraph 225. Here the authors make a first assumption that 10% of the particles' thermal energy is available to produce steam. A reference or more grounds on this assumption would be helpful. Furthermore, the larger fraction for smaller particles is reasonable and can perhaps be explained by the larger specific surface (hence, also the type of particles (e.g., pumice vs. lithic fragments) via the shape).
• Sedimentation rate vs. deposition rate. The model equations contains a sedimentation term which is represented by a "D", which in turn recalls deposition. The authors refer to this term as "volumetric rate of deposition". In my opinion, this is actually the volumetric rate of sedimentation, i.e. the rate at which particles leave the flow. The latter, though, does not necessarily coincide with the deposition rate, which quantifies the rate at which the deposit “grows”. I suggest the authors to review the use of this terminology. The model like the one used in IMEX (eq. 19) and the sink term in the conservation equations are sedimentation rate terms, deposition rates may differ because there may be processes at the deposition interface that mobilize (e.g., traction carpet) or even re-entrain the particles settling onto the ground.
  
  Furthermore, the use of a model like eq. 19 has the limitation that it does not take into account the effect of the flow itself on the sedimentation of particles, specifically the turbulence. In fact, models like eq. 19 assume a constant sedimentation rate that does not take into account the turbulence, which can keep particle in suspension and this is regulated by the Rouse number as well established in classic sedimentology. See for example the recent works of Dellino et al. (https://doi.org/10.1111/sed.12485,  https://doi.org/10.1111/sed.12693) and cited literature therein. I am not asking the authors to implement these models (which may add a further complexity) nor to cite them, but suggesting to make these assumptions and their consequences (i.e., a potential overestimation of the sedimentation rate, hence underestimation of the flow runout) clear in the text.
• Settling velocity of particles. There are more recent and accurate treatments of the particle settling velocity calculations. For simple spheres (which, in turn, is a strong and limiting assumption, see my comments in the following), there are drag laws that avoid the switch condition in the Lun and Gidaspow drag law (e.g. https://onlinelibrary.wiley.com/doi/10.1002/cjce.5450490403; https://doi.org/10.1016/0032-5910(89)80008-7), which in some situations may cause problems. 
  
  Then there are the problems arising from using a drag law for spheres in the case of natural sediments/volcanic particles, which are (often very) far from spheres, with strong implications on the drag coefficient and therefore the terminal (settling) velocity. Even if the authors do not want to implement such models (e.g., https://doi.org/10.1016/j.powtec.2018.12.040, https://doi.org/10.1002/2017JB014926, https://doi.org/10.1016/0032-5910(93)80051-B), they should clearly state that their approach is a simplification and that can lead to overestimations of the terminal velocity, hence of the sedimentation rates in the flow. 
• "Realistic" test case. The authors apply IMEX to the Krakatau 2018 eruption test case, but they state that "the simulations are not aimed at replicating the eruption..." (I assume, they refer to the DPDCs of the eruption). Whilst I acknowledge the usefulness of the exercise they present, I think that a complete comparison with a real case, for which data are available, would add more value to the manuscript and the software. The potential user would be more keen to use IMEX for realistic applications if they saw a successful benchmarking against a real case. It would be also interesting to see benchmarking applications against large-scale experiments, like the experiments carried out recently in New Zealand (PELEE) or the older ones by the University of Bari group. Maybe the former are more complicated for this code, since they produce channelised DPDCs; the latter, though, should be straightforward to reproduce.
• Finally, since this code should guarantee run times compatible with probabilistic volcanic hazard assessments, it would be useful to have some data on the run times vs. used computational resources of a realistic test case (e.g. the presented Krakatau test case).

Technical corrections and minor points

The author can find other few comments on minor points and technical corrections in the attached commented manuscript.

I hope my comments will represent a well-received contribution to improving the manuscript.

With kind regards.

Fabio Dioguardi, PhD

Review of the manuscript: IMEX_SfloW2D v2: a depth-averaged numerical flow model for volcanic gas-particle flows over complex topographies and water, by Mattia de’ Michieli Vitturi, Tomaso Esposti Ongaro, and Samantha Engwell.

The manuscript presents new developments to the numerical code IMEX_SfloW2D. Improvements include sedimentation and entrainment, transport equations for solid particles of different sizes, transport equations for different components of the carrier phase, an equation for temperature/energy, vaporization and entrainment of water and the simulation of sub to supercritical regimes.

The paper is clear and well presented. The systems of equations seem correct. However, as I am neither a mathematician nor a numericist, I cannot say that there are no errors and I recommend the advice of a specialist in the mathematical and numerical methods presented.

This is a very mathematical and methodological article describing the capabilities of code. The purpose of the code is not to reproduce or explain volcanic eruptions, even if the Krakatau application shows impressive potentialities. The manuscript is probably difficult to access for a more geological or volcanological audience. However, the scope of the paper is totally suitable to the journal, the article lays the foundations for a promising code for future studies in volcanology and I think it deserves to be published.

I add that I broadly agree with the comments posted by the other reviewer on 16 June.

Please find below minor comments and questions:

Equation 1-4: the notation must be similar. s,is in the equation, sis in the text.

Line 118: if possible try to avoid the exponent that can be mistaken for “to the b power”

Line 143: it's difficult to know whether this simplification has little or a major impact on the results. Under what conditions could induced errors alter the model's results in relation to the natural phenomenon?

“This can lead to numerical errors associated with the mixture temperature obtained from the total mixture specific energy and the kinetic energy computed from mass and momentum equations. For this reason, in some cases, instead of the full energy equation as presented above, it is preferable to solve a simpler transport equation for the specific thermal energy CvT”

175: Froude number

178: “The code we present is mostly aimed at simulating 2D spreading flows, for which we do not simulate the initial phase”. Once published and distributed, this code will most likely be used to simulate complete eruptions, including initial phases. Can it be used? Is it possible to quantify the errors if it is used to reproduce the initial phase? Errors of a few percent can be a problem from a mathematical point of view, but are perfectly acceptable in volcanology, where knowledge of source conditions is often limited.

Equation 16 / line 125: Something is not clear to me. Is the entrainment calculated only on the edges or in the whole flow? From equation 16, it seems that entrainment also takes place within the flow itself. In that case, the hypothesis of line 125 does not seem right to me because the flow carries momentum: “There are no terms associated with air and water vapour entrainment, because they do not carry any horizontal momentum into the flow.”

They way chosen to define the Ridchardson number is close to a Froude number ( Ri = 1/Fr² ). I am a bit lost to see the link with the ratio of the stabilizing stratification of the current to destabilizing velocity shear and, consequently, with the entrainment.

289: Could the authors explain what is an “opportune slope limiter”? On what criteria are they chosen? This choice seems to have strong consequences on edge values and therefore on the fluxes.

Figure 1: The authors use two flux values, Q(w, j-1, k) and Q(E, j, k) for example. What is the physical reality of these two flows at the same location? I didn't understand whether there was a procedure to make them compatible.

Line 311: same question here. I wonder if it's physically correct to consider 2 speeds, one to the left, the other to the right, in the same place? Is this choice linked to the physics of the phenomenon or to numerical stability?

Lines 325-327: I'm not able to understand what these choices imply for the model's results.

Figure 2: captions d and c are inversed.

Figure 2, 3, 4… : if I'm not mistaken, the proposed solutions are only those of the model. Is there no analytical solution to prove that the model reproduces them correctly or to help the reader visualize the model's precision? How do the reader know if the model is working properly?

Figure 10 (13 and the others): if possible, change the curve lines to adapt them to black and white printing

Figure 13: why at the source, the entrained air represents 100% of the gases? Can't you distinguish between air initially present and entrained air? The scale does not seem right for the low particle concentration. Perhaps a more specific scale on the right-hand axis would be more appropriate.

Line 500 and 365: kg m-3 (with a space)