The footprint of column collapse regimes on pyroclastic flow temperatures and plume heights

The gravitational collapse of eruption columns generates ground-hugging pyroclastic density currents (PDCs) with highly variable temperatures, high enough to be a threat for communities surrounding volcanoes. The reasons for such great temperature variability are debated in terms of eruptive versus transport and emplacement processes. Here, using a three-dimensional multiphase model, we show that the initial temperature of PDCs linearly correlates to the percentage of collapsing mass, with a maximum temperature decrease of 45% in the case of low percentages of collapse (10%), owing to an efficient entrainment of air into the jet structure. Analyses also demonstrate that column collapse limits the dispersal capabilities of volcanic plumes, reducing their maximum height by up to 45%. Our findings provide quantitative insights into the mechanism of turbulent mixing, and suggest that temperatures of PDC deposits may serve as a marker for determining column collapse conditions, which are of primarily importance in hazard studies.


Review report on "Thermal regime of eruption column collapses and generation of hot pyroclastic flows" by Trolese et al.
This paper provides interesting and thoughtprovoking numerical results focusing on collapsing eruption columns, and discusses their applications to field data concering temperatures of pyroclastic density currents (PDCs). The subject is important and some of the present results are novel. The present numerical results support the column collapse condition which has been recently proposed by Koyaguchi and Suzuki (2018) and Koyaguchi et al., (2018) (see Appendix of this review report). In addition, the results have quantitatively shown how the partial collapse develops in the transitional state of column collapse, which is certainly novel. I agree with the authors' conclusion that the present study provides a new insight into the mechanism of partial collapse of eruption column. On the other hand, there seem to be several points to be improved in this paper; some of the analyses of the simulation results look still preliminary and do not fully explain the observations in nature (e.g., those in Figure 3). Accordingly, I recommend the publication of this paper after some revisions. The points that can be improved are listed below. Figure 2c shows that the proportion of the collapsed flow, Q c , is correlated with the temperature of collapsing downflow, (T c −T a )/(T 0 −T a ). This is an interesting observation, and it is the key idea of this paper to explain the variation in temperature of PDCs in nature (i.e., Figure 3). However, the physics behind the relationship in Figure 2c is not clearly presented in the paper. As a consequence, a part of the geological implication using Figure 3 is not convincing as it stands.

Major comment
In fact, there is one important unresolved problem in Figure 2c. That is, how can the mixture with low (T c − T a )/(T 0 − T a ) be denser than air to form a PDC? The variations in Q c and (T c − T a )/(T 0 − T a ) observed in Figure 2c reflect the variation in the mass fraction of ejected magma (i.e., mixing ratio of ejected magma and air) in the collapsed flow (let us define this mass fraction as ξ); Q c and (T c − T a )/(T 0 − T a ) decrease as ξ decreases. The problem is that the mixture with low ξ is inevitably buoyant and it cannot form a PDC.
Let us consider a homogeneous mixture with ξ as a reference state. For simplicity, it is assumed that the effect of particle separation is negligible. In this reference case, the temperature of the eruption cloud, T , is calculated from where C Pair and C Pm are the specific heats of the air and the magma at constant pressures, respectively, and T a and T 0 are the temperatures of the ambient air and the magma, respectively. The specific heat of the magma is calculated from those of pyroclasts (C s ) and volcanic gas (C Pg ) as C Pm = y w C Pg + (1 − y w )C s . This equation means that the vertical axis of Figure 2c is closely correlated with ξ (see the black curve in Figure 1 in this review report). On the other hand, substituting T in Eq. (1) into the equation of state, we can calculate the density of the mixture ρ as  (1) and (2). T 0 = 1123 K and y w = 0.02 are assumed. The horizontal thin red line indicates the density of air at p = 1.013 × 10 5 Pa (ρ air ).
where P is the pressure, R air is the gas constant of the air. Here, we used the approximation that the volume fraction of pyroclasts is negligibly small in the eruption cloud. From Eq.(1) and Eq.
(2) we obtain the relationship between the mass fraction of the ejected material ξ and the density ρ as the red curve in Figure 1. My preliminary calculation in Figure 1 suggests that the mixture with (T c − T a )/(T 0 − T a ) < 0.7 is expected to be buoyant ( Figure 1). This simple consideration based on the reference state leads to a conclusion to explain the occurrence of low temperature PDC in Figure 3. The values of T and ξ for ρ = ρ air (the star in Figure 1) primarily depend on T 0 ; they decrease as T 0 becomes low, as is the case of phreatomagmatic eruptions (see Koyaguchi and Suzuki (2018) for more comprehensive discussions). Thus, it is inferred that very low temperature PDCs in Figure 3 resulted from phreatomagmatic eruptions with low T 0 . Figure 2c in this paper, on the other hand, indicates that the variation in PDCs' temperature can result from variation in U 0 (i.e., initial momentum) for given mass eruption rate Q 0 . As U 0 increases for given Q 0 , Q c decreases (see Appendix below), which generates a PDC with a low (T c − T a )/(T 0 − T a ). To apply this result to the observations in nature (e.g., the observations in Figure 3), we must specify the minimum value of (T c − T a )/(T 0 − T a ) for the mixture to be denser than air to form a PDC (i.e., whether it is as low as 0.5 like the case of Figure 2c or whether it can be even lower). For this purpose, we must understand the mechanism how the mixtures with (T c − T a )/(T 0 − T a ) < 0.7 can be denser than air. There may be some possible mechanism to form a dense mixture with low (T c − T a )/(T 0 − T a ) as follows.
• The mixture of the erupted material and air is not necessarily a consequence of simple binary mixing, which is assumed in the calculations of Eqs. (1) and (2), but may be heterogeneous. If the thermal equilibrium in the mixture is incomplete, the density of such heterogeneous mixture may deviate from the red curve in Figure 1 to considerable extent (see the light blue region in Figure 1 in this report).
• The above calculations of Eqs. (1) and (2) is based on the assumption that the effect of the particle separation is negligible, whereas Supplementary Figure 3 indicates the effect of the separation of coarse pyroclasts may be significant.
To understand the implication of Figure 2c, we must somehow evaluate the effects of the above mechanisms. Particularly, to clarify the effect of the first mechanism, it is essentially important to describe the heterogeneity of the collapsing down-flow in the simulation results. For example, according to the definition in line 112 of Supplementary information, T c is calculated as an averaged temperature; however, it is not clear from the averaged value whether the collapsing down-flow is composed of a uniform mixture with a uniform temperature of T c or it is composed of heterogeneous mixture with a wide range of T . It is also helpful if we could know the variations in the temperature T , the density ρ and the mass fraction of ejected material ξ at each grid point before averaging; those values of T and ρ can be plotted as a function of ξ in a diagram like Figure 1 in this review report.
I suggest that the geological implication using Figure 3 is possible only after the above discussions on the mixture with low (T c − T a )/(T 0 − T a ). Judging from Figure  2c, I have a feeling that the above two mechanisms may explain a dense mixture with (T c − T a )/(T 0 − T a ) as low as 0.5; however, to explain the occurrence of PDCs with very low temperature (those of T = 100-300 • C in Figure 3), additional effects such as the presence of ground water would be necessary.

Specific comments
Some other specific comments follow.
(1) Lines 62-67: From the description of Supplementary Methods I could infer that the 3-D numerical simulations in this study are actually those for pressure-balanced jets with variable U 0 (i.e., so-called correctly expanded jets); however, the main text (as well as Table 1) gives an impression that the simulations are done for variable U v and P v , which is misleading. In this study, the physical quantities at P = P 0 (atmospheric pressure) are assumed to be uniform (i.e., top-hat distributions). On the other hand, recent numerical studies (Carcano et al., 2013;Koyaguchi et al., 2018) indicate that, when the jet is not pressure-balanced at the vent, the decompression and compression processes just above the vent strongly modify the distributions of the physical quantities at P = P 0 . As is discussed in Major comment, the distributions of the physical quantities in the column (e.g., the heterogeneity of T ) is considered to be critical to generate a dense mixture with a relatively low (T c − T a )/(T 0 − T a ). To avoid readers' misunderstanding, the author should explicitly state in the main text that the 3-D numerical simulations in this study are essentially those for pressure-balanced jet with variable U 0 , and, if possible, discuss how the assumption of top-hat distributions of physical quantities at P = P 0 affects the main conclusions of this paper.
(2) Lines 95-97: The definition of partial collapse based on Q c is rather indirect and it requires more explanations. I agree that Q c represent something important for understanding of partial collapse, but I am not fully convinced that Q c represents the fraction of collapsing mass as a PDC. There are at least two ambiguities. The separation of particles results in mass loss from an ascending column. Therefore, strictly speaking, the difference between the mass flow rate at the black dot and the total mass eruption rate does not represent the collapsing mass flow rate as a PDC. Considering the presence of large scale eddy around the top of fountain (H c ), a part of erupted mass may also be trapped by the large scale circulation even though the whole down-flow is not denser than the ambient air. These effects of the particle separations and circulations should be quantitatively evaluated to understand the meaning of Q c . In my opinion, to clarify whether Q c definitely represents the fraction of collapsing mass as a PDC, it is essentially important to show that the collapsing down flow is denser than the ambient air.
(3) Line 103: What is the definition of H jet ? Does it the level of the minimum point of horizontally integrated velocity in Figure 1? If so the level does not necessarily represent the level where the vertical momentum of the whole eruptive column is exhausted. The profile of horizontally integrated velocity is determined by competing effects of decreasing momentum due to negative buoyancy and increasing momentum of positive buoyancy. If the effect of increasing buoyancy due to positive buoyancy is large, the upward velocity may have a minimum point at a lower level before the initial momentum is exhausted.
(4) The result in Figure 2b (and related text from line 141 to 145) is interesting, although it may not be the main subject of this paper. In fact, similar tendencies have been pointed out in Costa et al. (2018) and Koyaguchi et al. (2018). Although the present result in Figure 2b is more comprehensive and looks more systematic, those previous works can be cited here.
(5) The caption of Supplementary Figure 3 needs to be corrected. There seem to be two important typos. Judging from the value of H c (i.e., the positions of the black dots), the black curves are the results of partial collapse and the pink curves are the results of the near total collapse. In the pink curves, the dotted curve shows nearly zero mass flow rate at high altitudes. According to the caption, the upper part of the this column (with pink color) is composed of coarse particles shown by dashed curve alone and almost fully depleted in fine particles, which does not make sense to me. I suspect that dashed and dotted curves represent the mass flow rate of fine and coarse particles, respectively.
Other minor points are annotated in the pdf file.

Appendix: Consistency and discrepancy with the column collapse condition proposed by Koyaguchi and Suzuki (2018) and Koyaguchi et al. (2018)
In the course of review process, I have confirmed how the present numerical results are consistent with the theoretical and numerical column collapse condition by Koyaguchi and Suzuki (2018) and Koyaguchi et al. (2018). Koyaguchi and Suzuki (2018) have shown that the column collapse condition of a pressure-balanced jet (e.g., Bursik and Woods, 1991) is expressed by a simple relationship as where M a is the Mach number after decompression (U 0 / √ y w R w T 0 using the present notation) andQ is the dimensionless mass eruption rate (see Koyaguchi and Suzuki (2018) for the normalization factor). An eruption column is stable to form a buoyant plume when M a /Q 1/5 > 1, whereas it collapses to form PDCs as M a /Q 1/5 < 1.
Using the values in Table 1, we can plot Q c (%) as a function of M a /Q 1/5 for the present numerical results (Figure 2). The present results are consistent with the prediction by Koyaguchi and Suzuki (2018) in the sense that Q c = 0 (i.e., fully buoyant) for M a /Q 1/5 > 1. The present results also show that Q c gradually increases with decreasing M a /Q 1/5 in the region of M a /Q 1/5 < 1. It is interesting to see the variation of Q c with M a /Q 1/5 collapses to a single relationship for Q 0 = 10 7 -10 9 kg/s; this may suggest that the present numerical study has captured a novel universal feature of the transition from a buoyant plume to total collapse.
It should be noted that the present numerical results and those of Koyaguchi et al. (2018) do not perfectly coincide with each other in a quantitative sense. The transitional state (i.e., partial collapse etc) is observed in the range of M a /Q 1/5 < 1 in the present study. On the other hand, it is typically observed in the range 1 < M a /Q 1/5 < 1.2 in the 3-D simulations of Koyaguchi et al. (2018); the range of the transitional zone varies with Q 0 and P v to a considerable extent in their simulations. This quantitative discrepancy may be caused by the fact that the present simulations take the effects of particle settling into account, whereas those of Koyaguchi et al. (2018) do not. The discrepancy may partly result from the difference in the condition at the level of P = P 0 . In Koyaguchi et al. (2018), the boundary condition of P v ̸ = P 0 is set at the crater top such that the physical quantities at the level of P = P 0 is heterogeneous because of the decompression/compression processes; this heterogeneity causes the diverse features in the transitional state in their simulations. In this study, on the other hand, a top hat distribution is imposed at the level of P = P 0 . It is suggested that further study would be necessary to fully understand the diverse features in the transitional state. The color of symbols are the same as the manuscript. Green: Q 0 = 10 7 kg/s, red: 10 8 kg/s, and balck: 10 9 kg/s.

Reviewer #1 -Erika Rader
This study uses a 3-d numerical model to simulate heights and convection processes within volcanic eruption columns with the intention of constraining the conditions that lead to column collapse and PDC formation. They found that reduced entrainment of ambient air leads to a hot, unstable column which can travel very far when it collapses. This variable is more likely to produce a long run-out from a PDC than the mass flux, which was previously thought to be the most significant factor for eruptions like Taupo. The conclusions are interesting and valuable to the volcanologic community and their work is convincing, although more discussion of what the modeling terms mean in real world scenarios (eg. What does a density of 19.3 kg/m3 mean and is it found in nature? See comment for Line 159 also) would make the results more impactful. There is not enough information in this manuscript for a non-expert to repeat the study, however, they list other studies which have more detailed methods sections, suggesting one could work it out if need be.
We sincerely thank the reviewer, Dr Erika Rader, for her positive comments about our manuscript. As we now specify in the revised manuscript (lines 77-78/373-376, manuscript with highlighted changes), the input parameters we have applied to our simulations were selected to be as realistic as Concerning the lack of information, we prefer to not provide an accurate description of the fluid dynamical model, as the numerical code has been already extensively described (and tested against numerical and experimental benchmarks) in previous papers, which we included as references in our manuscript. We rather choose to focus on describing the adopted initial and boundary conditions, and on defining the variables used to investigate the collapse dynamics, so as readers can easily assess our experimental procedure. We apologize for not being clear enough at this point. The intention was to say that a small amount of material rises buoyantly through the atmosphere as a vertical plume up to about 19 km, whereas 79% of the erupted material collapses to form a pyroclastic flow of long runout (>20 km). We have updated this information in the manuscript (lines 161-163, manuscript with highlighted changes).

Line 119 -fluctuations of what?
We simply meant density and temperature fluctuations, we have reworded accordingly (lines 169-174, manuscript with highlighted changes).
Line 127 -might be good to state what the "whole set of numerical simulations" includes here.
We are referring to all the simulations that are listed in Supplementary We revised this sentence so it now reads, "A corollary of our results is that progressively higher column collapse heights do not generate faster PDCs with longer runouts, in contrast with some previous views 38-40 ". See lines 260-262, manuscript with highlighted changes.
Line 145 -perhaps this isn't something to be added to the manuscript, but how did you test how atmospheric conditions might affect the likelihood of air entrainment? What conditions support producing these hotter PDC run outs?
We thank Erika Rader for her comment. We agree that different meteorological conditions may in part affect the efficiency of air entrainment during the jet phase, although it has recently been shown (e.g. Costa et al., 2016; Suzuki and Koyaguchi 2015) that, for example, atmospheric winds are mostly important in determining the total plume height for an assigned mass eruption rate, as they can shift downwind the buoyant and umbrella regions of the plume. This is particularly true in the case of weak plumes (mass eruption rate lower than 10 7 -10 8 kg/s), where the wind speed is comparable to the plume velocity. Given that our manuscript is the first study showing the thermal patterns leading to different column collapse regimes, we think that it is better and clearer to not consider the effects of wind, which will be investigated in future trials, as it is an interesting area for further extension. We decided to focus our study on the effect of the eruptions source parameters on the collapse efficiency, as we expect that those have a primary effect on the jet dynamics. For these reasons, our simulations are carried out in a stratified windless atmospheric condition (we have provided a new figure in Supplementary Materials to show the atmospheric profiles used for the simulations). Nonetheless, we expect the wind to have an effect too, even if weaker. Qualitatively, we expect an enhanced entrainment rate in strong wind conditions. It is worth noting that we also consider ash kinematic decoupling and preferential concentration of particles (e.g. Cerminara et al., 2016) to take into account the effect of coarse ash on the collapsing dynamics, and hence on the mixing efficiency. In the revised version, we highlight all the previous points (please refer to lines 70-72/366-369 in the manuscript with highlighted changes). Line 157 -I don't believe the citation you provided for this line states that the glass transition is going to be 74-82% of the eruption temperature. You are assuming that eruption temperature is the same as the experimental temperature Giordano used in that paper. Perhaps clarify this point a bit so someone doesn't assume volcanos only ever erupt the temperature where glass of that composition is a crystal-free bubble-free liquid, as they were in the experiments.  105-118 (2005) Line 159 -Your model provides no data in Figure 3. That appears to be literature search with no analysis using the model. Why not include estimates of air entrainment for each, or column densities for each or some output which your model gives by matching the temperature and runout length of the deposit? We explain this observation as a function of column collapse regimes. Since the original figure could not be directly fed with data from our numerical experiments, as suggested by Reviewer #1, we have edited this figure, by representing the temperature of the deposit proportionally to the eruptive temperature in order to include the percentage of collapsing mass derived by our models, which is linearly related to the temperature drop of the collapsing mass. We have therefore made changes to the explanatory text of the new version of Figure 4 at lines 293-328 in the manuscript with highlighted changes, reinforcing our findings, and explaining that: 1) the low mixing associated with near total collapse regimes fully justify welded pyroclastic flow deposits; 2) the highest mixing associated with partial collapse regimes cannot alone justify the occurrence of very low temperature pyroclastic flow deposits, and so it is necessary to consider at least an additional mechanism such as the involvement of external water during an eruption; 3) mixing associated with partial collapse -with percentages of collapse between 50% and 10%can likely be the principal factor to justify unwelded pyroclastic flow deposits with intermediate emplacement temperatures.
Line 301 -Include 'Qc' in the caption.

This has been added in the caption.
Line 322 -add 'symbols correspond to each eruptive center' to caption, or something similar.
This has been added in the caption.

Reviewer #2
This paper provides interesting and thoughtprovoking numerical results focusing on collapsing eruption columns, and discusses their applications to field data concerning temperatures of pyroclastic density currents (PDCs). The subject is important and some of the present results are novel. The present numerical results support the column collapse condition which has been recently proposed by Koyaguchi and Suzuki (2018) and Koyaguchi et al., (2018) (see Appendix of this review report). In addition, the results have quantitatively shown how the partial collapse develops in the transitional state of column collapse, which is certainly novel. I agree with the authors' conclusion that the present study provides a new insight into the mechanism of partial collapse of eruption column. On the other hand, there seem to be several points to be improved in this paper; some of the analyses of the simulation results look still preliminary and do not fully explain the observations in nature (e.g., those in Figure 3). Accordingly, I recommend the publication of this paper after some revisions. The points that can be improved are listed below.
We sincerely thank Reviewer #2 for the very helpful and encouraging comments. We appreciate his/her suggestions for improving the clarity and impact of our study. Overall, we recognize based on the Reviewer's comments that we could have described more the thermodynamic properties of the collapse region, which are important to understand the collapse mechanisms and, in turn, the variation in PDC's emplacement temperatures observed in nature. We have strived to address all the Reviewer's concerns and summarize them below. Figure 2c shows that the proportion of the collapsed flow, Q c , is correlated with the temperature of collapsing downflow, (T c -T a )=(T 0 -T a ). This is an interesting observation, and it is the key idea of this paper to explain the variation in temperature of PDCs in nature (i.e., Figure 3). However, the physics behind the relationship in Figure 2c is not clearly presented in the paper. As a consequence, a part of the geological implication using Figure 3 is not convincing as it stands.

Major comment
In fact, there is one important unresolved problem in Figure 2c. That is, how can the mixture with low (T c -T a )=(T 0 -T a ) be denser than air to form a PDC? The variations in Q c and (T c -T a )=(T 0 -T a ) observed in Figure 2c reflect the variation in the mass fraction of ejected magma (i.e., mixing ratio of ejected magma and air) in the collapsed flow (let us define this mass fraction as ξ); Q c and (T c -T a )=(T 0 -T a ) decrease as ξ decreases. The problem is that the mixture with low ξ is inevitably buoyant and it cannot form a PDC.
Let us consider a homogeneous mixture with ξ as a reference state. For simplicity, it is assumed that the effect of particle separation is negligible. In this reference case, the temperature of the eruption cloud, T, is calculated from where C Pair and C Pm are the specific heats of the air and the magma at constant pressures, respectively, and T a and T 0 are the temperatures of the ambient air and the magma, respectively. The specific heat of the magma is calculated from those of pyroclasts (C s ) and volcanic gas (C Pg ) as C Pm =y w C Pg +(1-y w )C s . This equation means that the vertical axis of Figure 2c is closely correlated with ξ (see the black curve in Figure 1 in this review report). On the other hand, substituting T in Eq.(1) into the equation of state, we can calculate the density of the mixture ρ as where P is the pressure, R air is the gas constant of the air. Here, we used the approximation that the volume fraction of pyroclasts is negligibly small in the eruption cloud. From Eq.(1) and Eq. (2) we obtain the relationship between the mass fraction of the ejected material ξ and the density ρ as the red curve in Figure 1. My preliminary calculation in Figure 1 suggests that the mixture with (T c -T a )=(T 0 -T a ) < 0.7 is expected to be buoyant ( Figure 1). (1) and (2) This simple consideration based on the reference state leads to a conclusion to explain the occurrence of low temperature PDC in Figure 3. The values of T and ξ for ρ=ρ air (the star in Figure  1) primarily depend on T 0 ; they decrease as T 0 becomes low, as is the case of phreatomagmatic eruptions (see Koyaguchi and Suzuki (2018) for more comprehensive discussions). Thus, it is inferred that very low temperature PDCs in Figure 3 resulted from phreatomagmatic eruptions with low T 0 . Figure 2c in this paper, on the other hand, indicates that the variation in PDCs' temperature can result from variation in U 0 (i.e., initial momentum) for given mass eruption rate Q 0 . As U 0 increases for given Q 0 , Q c decreases (see Appendix below), which generates a PDC with a low (T c -T a )=(T 0 -T a ). To apply this result to the observations in nature (e.g., the observations in Figure 3), we must specify the minimum value of (T c -T a )=(T 0 -T a ) for the mixture to be denser than air to form a PDC (i.e., whether it is as low as 0.5 like the case of Figure 2c or whether it can be even lower). For this purpose, we must understand the mechanism how the mixtures with (T c -T a )=(T 0 -T a ) < 0.7 can be denser than air. There may be some possible mechanism to form a dense mixture with low (T c -T a )=(T 0 -T a ) as follows.

Figure 1: The normalized temperature (black curve) and the density (ρ; red curve) as a function of the mass fraction of ejected materials in the eruption cloud (ξ) based on Eqs.
• The mixture of the erupted material and air is not necessarily a consequence of simple binary mixing, which is assumed in the calculations of Eqs. (1) and (2), but may be heterogeneous. If the thermal equilibrium in the mixture is incomplete, the density of such heterogeneous mixture may deviate from the red curve in Figure 1 to considerable extent (see the light blue region in Figure 1 in this report). • The above calculations of Eqs. (1) and (2) is based on the assumption that the effect of the particle separation is negligible, whereas Supplementary Figure 3 indicates the effect of the separation of coarse pyroclasts may be significant.
To understand the implication of Figure 2c, we must somehow evaluate the effects of the above mechanisms. Particularly, to clarify the effect of the first mechanism, it is essentially important to describe the heterogeneity of the collapsing down-flow in the simulation results. For example, according to the definition in line 112 of Supplementary information, T c is calculated as an averaged temperature; however, it is not clear from the averaged value whether the collapsing down-flow is composed of a uniform mixture with a uniform temperature of T c or it is composed of heterogeneous mixture with a wide range of T. It is also helpful if we could know the variations in the temperature T, the density ρ and the mass fraction of ejected material ξ at each grid point before averaging; those values of T and ρ can be plotted as a function of ξ in a diagram like Figure  1 in this review report.
I suggest that the geological implication using Figure 3 is possible only after the above discussions on the mixture with low (T c -T a )=(T 0 -T a ). Judging from Figure 2c, I have a feeling that the above two mechanisms may explain a dense mixture with (T c -T a )=(T 0 -T a ) as low as 0:5; however, to explain the occurrence of PDCs with very low temperature (those of T = 100-300 °C in Figure 3), additional effects such as the presence of ground water would be necessary.
We really appreciate this constructive and well-grounded critique. In response to your comment, we presented a new Figure 2, along with two new Supplementary Figures  (Supplementary Figures 2, and 3) and a paragraph to the results section (named: "Thermodynamic constraints for collapsing gas-particle mixtures") in the revised manuscript that give an explanation about the variations in temperature T, density ρ and mass fraction ξ of the collapsing down-flow mixture, necessary to explain the mechanisms by which a mixture with ( − )/( − ) < . can be denser that the surrounding atmosphere. It is revised as follow:

Thermodynamic constraints for collapsing gas-particle mixtures
To better understand the mechanism of temperature drop due to mixing with atmospheric air, we here derive, from thermodynamic considerations, the thermal conditions that would characterize the collapsing mixture. We assume an isobaric transformation between the two following states: initial) the erupted material and air are perfectly separated; final) the two phases are perfectly mixed and the resulting mixture is homogeneous. The final mixture temperature ( ) can be expressed, from enthalpy conservation, as a function of the mass fraction of the erupted material 13 . A dimensionless temperature ( ) can be written as: where is the initial temperature of the erupted material, is the initial air temperature, and is the specific heat at constant pressure of the air (thermodynamic properties used are listed in Supplementary Table 2). The specific heat at constant pressure of the erupted material is where is the initial gas fraction in the erupted material, and are the specific heats of the gas and solids, respectively. The ratio between the mixture density ( ) and the air density is: where and are the gas constant of the erupted gas and air, respectively. Equations (1) and (2) are derived under the following assumptions: 1) the enthalpy is conserved, thus the transformation is at constant pressure and the effects of dissipation and gravity can be disregarded; 2) the mixture in the final state is homogeneous; 3) particle separation can be disregarded; 4) air temperature and density are constant, i.e., atmospheric stratification is neglected; 5) the volume fraction of solid is negligibly small. In the following, we compare theoretical results with the 3D simulation, to evaluate to which extent these assumptions are correct cell by cell and in the whole collapsing region. From equations (1) and (2), the minimum temperature and mass fraction for which the mixture is non-buoyant ( ( ) > ) can be derived: In our case, ≈ ≈ 1 and equations (1)(2)(3)(4) simplify into: As pointed out by ref. 13 , and primarily depend on and . Using = 1100 K and = 290 K, we get ≈ ≈ 0.64. Below this value, a homogeneous mixture that cooled down by the sole effect of air entrainment would necessarily be positively buoyant.
In the volcanic case, stratification would be non-negligible: the air density drop along the collapsing region can be as high as 40% and the temperature drop is about 10% (see Supplementary Fig. 1). In Fig. 2, we show the comparison between the theoretical relationships (1) and (2) (solid curves), and the simulated values of non-dimensional temperature and density in each cell of the collapse region (points). The simulated dimensionless density is evaluated with equation (2), considering atmospheric temperature stratification ( ) for each cell in the collapsing region, while its theoretical value is obtained using the average air temperature between the crater and the collapsing height. The two end-member cases, with and without particle decoupling, are compared. In the equilibrium case (dusty-gas; Fig. 2a, c), the thermodynamics of the numerical solution almost perfectly follows the simplified theoretical model in every cell of the collapsing region. This result is two-fold: on one side, it shows that all the hypothesis used to obtain equations (1) and (2) are locally valid; on the other side, it confirms the ability of the numerical solver to accurately solve the enthalpy equation. In the case of particle decoupling (Fig. 2b, d), both temperature and density values deviate from equations (1) and (2), although following the expected theoretical trend. Since in this case we cannot use equations (1) and (2) to define and , these parameters are evaluated by searching the last point with negative buoyancy (Fig. 2b, d). Particle separation makes and smaller than those predicted by equations (3) and (4). Consequently, particle clustering and settling enhance and the temperature drop. For the two end-members, we find = 15% and = 79% with particle decoupling, while = 4% and = 70% without. We quantify this effect in Supplementary Fig. 3, where a larger temperature drop is observed because of particle decoupling (up to 20% for the partial collapse end-member).
The distribution of the points along the whole theoretical curves of Fig. 2 clearly indicates that the collapsing region is highly heterogeneous, both in terms of particle concentration and temperature. Thus, even if equations (1) and (2) hold locally, they cannot be used to estimate the global temperature and buoyancy of the collapsing region. To estimate the average temperature in the heterogeneous collapse region, we used instead a weighted average (described in the Methods section), based on the mass distribution in every cell. Its non-dimensional value is reported in Figure 2 (red line) as . In the partial collapse case (Fig. 2a, b), is smaller than , due to the heterogeneity and enhanced entrainment in the collapsing stream. In the fountaining regime (Fig. 2c, d) > , clearly indicating the establishment of full collapsing conditions with a smaller entrainment efficiency. Supplementary Fig. 4, for the partial collapse case, = 790 K and the average PDC temperature = 809 K, while for the total collapse end-member = 1090 K and = 1062 K. The good agreement between these estimates confirms the validity of our approach.

Figure 2.
Thermodynamic properties of the collapse region for the reference models. Dimensionless temperature and density / as a function of erupted material mass fraction in the collapsing region. Blue and red points are extracted from the numerical simulation data (temperature, mass fractions, density) and the atmospheric temperature in each cell of the collapse region (see Supplementary Fig. 1). Yellow and green solid lines are obtained from equations (1) and (2), using an average atmospheric temperature between the crater and the collapse height. The red solid line highlights the non-dimensional value of the average temperature of the collapsing region (computed as described in the Methods section), while the black solid line is the level of neutral buoyancy. Black dotted and dashed lines indicate and , respectively. All the blue points below and all the red points left of are less dense than the atmosphere at the same height. Panels a and b are for the partial collapse end-member, while c and d are for the total collapse. Panels a and c are for the dusty-gas model, while b and d take into account particle decoupling.

Supplementary Figure 3 | Comparison between dusty-gas and kinematic decoupling simulations.
Relationship between the percentage of collapsing mass and its corresponding average temperature value normalized to the initial magmatic temperature with respect to the atmospheric air temperature , for dusty-gas (black symbols) and kinematic decoupling (red symbols) simulations with the same initial eruptive conditions (mass eruption rate = 10 8 kg/s; temperature = 1123 K; water content 0.5< <2.5 wt%; pressure 0.1< <10 MPa). The black and red solid lines show the linear least square regression for dusty-gas and kinematic decoupling simulations, respectively. The red regression line is drawn using all data from kinematic decoupling simulations (it is the same shown in Figure 3c).

Supplementary Figure 4 | Time series of PDC temperature data for the reference models.
Temporal evolution of the mixture temperature (black solid line) in the bottom layer of PDC generated by the partial (a) and total (b) column collapse end-members. Temperature is measured at a distance from the injection point of 10 inlet radii ( = 10 ). Green solid lines indicate the time-averaged temperature of each pulse. The PDC temperatures (T pdc ) averaged over the reported time-window are 809 K and 1062 K for the partial and total collapse end-members, respectively. Such temperatures are very similar (within less than 5% in this case) to the mean temperatures at collapse.
The second mechanism contributes to the local (i.e. cell-by-cell) deviation from the theoretical law, and we were able to quantify such an effect (Figs. 2b,d). In the manuscript, we also listed: gravity effects, dissipation, atmospheric stratification, volume of the solid, as other assumption needed to get the Eqs. (1) and (2) suggested by the Reviewer. Moreover, we modify the theoretical equation for the density writing directly the density ratio ( ) to reduce the effects due to atmospheric stratification.
However, the main source of discrepancy between the predicted minimum temperature that a homogeneous collapsing column can reach, and the calculated averaged temperature is due to the heterogeneity of the collapsing region. Since the collapsing region is highly heterogeneous one (as shown by the new Figure 2), the collapsing temperature needs to be defined by some averaging procedure. We decided to weight-average the temperature with the mass contained in each cell. As a further check of the robustness of the procedure, we compared the mean temperature T c in the collapsing region with the temperature at the base of the forming PDCs T pdc (Supplementary Fig. 4). It turns out that ~ for both the end-member cases. The good agreement between these estimates confirms the validity of our approach. In Fig. 2 we also quantify the effect of particle decoupling on the conservation of the enthalpy. Decoupling enhances column collapse.
We believe that the new analytical and numerical analysis proposed by Reviewer #2 strengthens the understanding of the emplacement temperature variability of pyroclastic flow deposits shown in Figure 4 (which, as suggested by Reviewer #1, we amended to include the percentage of collapsing mass derived by our models, hoping that this addition makes comparison between the observed temperature data and our model results more quantitative).
We would like to thank you again for allow us to verify that the ASHEE code properly solve the energy balance equation.

Specific comments:
(1) Lines 62-67: From the description of Supplementary Methods I could infer that the 3-D numerical simulations in this study are actually those for pressure-balanced jets with variable U 0 (i.e., so-called correctly expanded jets); however, the main text (as well as Table  1) gives an impression that the simulations are done for variable U v and P v , which is misleading. In this study, the physical quantities at P = P 0 (atmospheric pressure) are assumed to be uniform (i.e., top-hat distributions). On the other hand, recent numerical studies (Carcano et al., 2013;Koyaguchi et al., 2018) indicate that, when the jet is not pressure-balanced at the vent, the decompression and compression processes just above the vent strongly modify the distributions of the physical quantities at P = P 0 . As is discussed in Major comment, the distributions of the physical quantities in the column (e.g., the heterogeneity of T) is considered to be critical to generate a dense mixture with a relatively low (T c -T a )=(T 0 -T a ). To avoid readers' misunderstanding, the author should explicitly state in the main text that the 3-D numerical simulations in this study are essentially those for pressure-balanced jet with variable U 0 , and, if possible, discuss how the assumption of top-hat distributions of physical quantities at P = P 0 affects the main conclusions of this paper.
We confirm (and we have better stated in the text at Lines 78−80, manuscript with highlighted changes) that source conditions for simulations are set for a pressure-balanced (correctlyexpanded) jet with a top-hat profile (i.e., uniform inlet conditions). The choice of limiting to these conditions was motivated: 1) to avoid simulating the decompression above the conduit exit, requiring very high-resolution grids at the inlet (that would have needed an enormous computational effort) and 2) to simplify the source parameters when focusing on the collapse dynamics. We are aware, especially after reading the more recent Koyaguchi and Suzuki 2018 paper, of the potential impact that a more accurate description of the compression/decompression patterns might have on the collapse dynamics. We have therefore welcomed Reviewer #2 suggestion to comment about this issue (Lines 81−86, manuscript with highlighted changes), which will be the topic for a future work.
We also agree with Reviewer #2, that we might have imposed directly the mixture parameters at the inlet, instead of imposing them at the conduit exit and computing analytically the conditions after the correct expansion, and that this can be somehow misleading for the reader. However, not only the two choices are completely equivalent, but we believe that ours allows a much direct linking to vent conditions deriving from conduit flow models. This is the approach that we followed also in our previous works (e.g., Esposti Ongaro et al., 2008).
We finally remark that the term vent in our paper refers to the conduit exit, not the crater exit. There is unfortunately not a general consensus in the literature about such nomenclature, and we preferred to be consistent with our previous works. (2) Lines 95-97: The definition of partial collapse based on Q c is rather indirect and it requires more explanations. I agree that Q c represent something important for understanding of partial collapse, but I am not fully convinced that Q c represents the fraction of collapsing mass as a PDC. There are at least two ambiguities. The separation of particles results in mass loss from an ascending column. Therefore, strictly speaking, the difference between the mass flow rate at the black dot and the total mass eruption rate does not represent the collapsing mass flow rate as a PDC. Considering the presence of large scale eddy around the top of fountain (H c ), a part of erupted mass may also be trapped by the large scale circulation even though the whole down-flow is not denser than the ambient air. These effects of the particle separations and circulations should be quantitatively evaluated to understand the meaning of Q c . In my opinion, to clarify whether Q c definitely represents the fraction of collapsing mass as a PDC, it is essentially important to show that the collapsing down flow is denser than the ambient air.
We thank Reviewer #2 for this comment. As suggested, we now explain in more detail our motivations to choose as a measure of the collapsing intensity. In the Methods section, we also quantify the effect of the recirculation eddy (lines 430-451 in the manuscript with highlighted changes). It is instead very difficult to directly quantify the amount of mass going into the PDC, because while the mixture spreads horizontally it continuously loses the mass that becomes positively buoyant, as a result of the entrainment and mixing with the atmosphere. In other words, where is the edge between the collapsing region and the beginning of the PDC?
For this reason, we decided to use a clearly defined method where the collapsing mass is obtained by comparing the mass flow rate of the ascending column above the collapsing region with the mass eruption rate. Thanks to conservation of mass, all the mass not present in the ascending column should have been lost by the collapsing region and the recirculation eddy.

As we describe in our previous reply to the Reviewer's major comment, the collapsing region is heterogeneous, so it contains both positively and negatively buoyant parcels (cells).
(3) Line 103: What is the definition of H jet ? Does it the level of the minimum point of horizontally integrated velocity in Figure 1? If so the level does not necessarily represent the level where the vertical momentum of the whole eruptive column is exhausted. The profile of horizontally integrated velocity is determined by competing effects of decreasing momentum due to negative buoyancy and increasing momentum of positive buoyancy. If the effect of increasing buoyancy due to positive buoyancy is large, the upward velocity may have a minimum point at a lower level before the initial momentum is exhausted.
Thank you for highlighting this point. We have clarified the explanation of our H jet measurement strategy further, as we realized how this can be misleading. As we now specify in the manuscript (lines 143-146/419-429 in the manuscript with highlighted changes), we decided to evaluate H jet by using a method that is not based on any assumption on the entrainment profile, but it is only based on the velocity profile. The other two methods that we now refer to in the paper, the Morton length scale generalized to the non-Boussinesq case and the method defined by Koyaguchi and Suzuki (2018), require the definition of the entrainment coefficient (to be either 0 or a constant value). In any case, we show that our methodology to evaluate H jet gives broadly comparable results to those obtained with the other two methods. (4) The result in Figure 2b (and related text from line 141 to 145) is interesting, although it may not be the main subject of this paper. In fact, similar tendencies have been pointed out in Costa et al. (2018) and Koyaguchi et al. (2018). Although the present result in Figure  2b is more comprehensive and looks more systematic, those previous works can be cited here.
We have now expanded the discussion around Figure 2b (please refer to lines 264-276/332-335 in the manuscript with highlighted changes), referring to the works suggested by the Reviewer. Based on Figure 2b, we also emphasize the importance of not only considering the effects of wind but also the collapse regime for determining the total plume height, which, for a given mass eruption rate, can be substantially reduced in the case of a partial collapse.
(5) The caption of Supplementary Figure 3 needs to be corrected. There seem to be two important typos. Judging from the value of Hc (i.e., the positions of the black dots), the black curves are the results of partial collapse and the pink curves are the results of the near total collapse. In the pink curves, the dotted curve shows nearly zero mass flow rate at high altitudes. According to the caption, the upper part of the this column (with pink color) is composed of coarse particles shown by dashed curve alone and almost fully depleted in fine particles, which does not make sense to me. I suspect that dashed and dotted curves represent the mass flow rate of fine and coarse particles, respectively.

Many apologies -Yes, the Reviewer is right. Thanks so much for catching this misstatement about the dashed and dotted curves that have now been corrected.
Minor points annotated in the manuscript: Line 39 -The initial momentum (U 0 ) is the other critical factor (e.g., Koyaguchi and Suzuki, 2018). See also the appendix of review report.
We have now added the reference suggested by the Reviewer (Lines 42-44, manuscript with highlighted changes).
Line 103 -How did you determined H jet ?

See comment (3).
Line 104 -I could not follow the logic of this part.
We rephrased this sentence in the revised manuscript (Lines 147-151 in the manuscript with highlighted changes) to remove the ambiguity. What we are saying is that when the jet becomes gradually more dilute due to turbulent shearing at its margins, then H un >> H jet , since the jet core remains unmixed. This is in accordance with what was previously suggested by Suzuki and Koyaguchi (2012) and Koyaguchi and Suzuki (2018).