The run-out distance of large-scale pyroclastic density currents: A two-layer depth-averaged model

Accepted Manuscript The run-out distance of large-scale pyroclastic density currents: A two-layer depth-averaged model

Hiroyuki A. Shimizu, Takehiro Koyaguchi, Yujiro J. Suzuki PII: DOI: Reference:

S0377-0273(18)30536-5 https://doi.org/10.1016/j.jvolgeores.2019.03.013 VOLGEO 6579

To appear in:

Journal of Volcanology and Geothermal Research

Received date: Revised date: Accepted date:

21 November 2018 18 March 2019 22 March 2019

Please cite this article as: H.A. Shimizu, T. Koyaguchi and Y.J. Suzuki, The run-out distance of large-scale pyroclastic density currents: A two-layer depth-averaged model, Journal of Volcanology and Geothermal Research, https://doi.org/10.1016/ j.jvolgeores.2019.03.013

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT Journal of Volcanology and Geothermal Research

The run-out distance of large-scale pyroclastic density currents: A two-layer depth-averaged model Hiroyuki A. Shimizu∗, Takehiro Koyaguchi, Yujiro J. Suzuki

RI PT

Earthquake Research Institute, The University of Tokyo, Japan

Abstract

CE

PT E

D

MA

NU

SC

We have developed a new two-layer pyroclastic density current (PDC) model that considers the effects of the entrainment and thermal expansion of ambient air in the upper dilute layer in order to investigate the dominant factors controlling the run-out distance of large-scale PDCs. Numerical simulations show a dilute current spreading radially from the collapsing eruption column edge and forming a thin, dense, basal current through particle settling; i.e., a two-layer PDC forms. Forward propagation of the two-layer PDC stops owing to lift-off of the dilute current and/or deposition of the dense current, and eventually the two-layer PDC converges to a steady state. Our parametric study classifies the flow patterns of steady-state two-layer PDCs into one of two flow regimes according to the relative magnitude of the run-out distances of the dilute and dense currents. This relative magnitude critically depends on the ratio of the deposition speed at the base of the dense current (D) to the speed of particles settling from the dilute current to the dense current (Ws ). The run-out distance of the whole two-layer PDC is determined by that of the dilute current when D/Ws is large (≳ 4 × 10−3 ). Otherwise, when D/Ws is small (≲ 4 × 10−3 ), the run-out distance of the dense basal current exceeds that of the dilute current; the former increases as D/Ws decreases. The run-out distance of the dilute current is too short to explain some of the long run-out distances of PDCs observed in the field (e.g., the 1991 Pinatubo and 2014 Kelud eruptions). It is suggested that the dense current traveling beyond the lift-off point of the parent dilute current plays a significant role in the emplacement of PDCs with such long run-out distances.

AC

Keywords: Explosive volcanism; Pyroclastic density current; Two-layer model; Run-out distance Highlights:

• We developed a new two-layer pyroclastic density current (PDC) model. • The model comprises a dilute current over a dense basal current. • Run-out distance of a PDC is determined from that of the dilute or dense current. • Run-out distances of dilute currents are greatly reduced by the entrainment of air. • Dense currents play a role in the emplacement of PDCs with long run-out distances.

∗ Corresponding author at: 1-1-1 Yayoi, Bunkyo-ku, Tokyo 113-0032, Japan. E-mail address: [email protected] (H.A. Shimizu).

1

ACCEPTED MANUSCRIPT Journal of Volcanology and Geothermal Research

1

Introduction

SC

RI PT

During volcanic eruptions, a hot mixture of volcanic particles and gas is ejected from the volcanic vent and can flow along the ground surface as a pyroclastic density current (PDC) (see the reviews by Druitt, 1998 [22]; Branney and Kokelaar, 2002 [5]; Sulpizio et al., 2014 [61]; Dufek, 2016 [24]). A PDC can devastate urban areas due to its high dynamic pressures and high temperatures (e.g., Cole et al., 2015 [14]), with the perimeter of the hazardous area determined by the run-out distance. The run-out distance depends on a number of geological factors and can extend over 100 km from the source vent during a large-scale eruption (Wilson, 1991 [75]; Streck and Grunder, 1995 [60]; Wilson et al., 1995 [76]; Henry et al., 2012 [32]; Roche et al., 2016 [52]). According to previous numerical and experimental works, it is affected by the eruption conditions such as magma discharge rate (e.g., Bursik and Woods, 1996 [9]; Dade and Huppert, 1996 [17]; Dufek and Bergantz, 2007 [25]), the physical processes of PDCs such as entrainment of ambient air and sedimentation (e.g., Lube et al., 2007 [45]; Roche et al., 2008 [53]; Andrews and Manga, 2012 [2]), and the topography (e.g., Doyle et al., 2008 [20]; Andrews and Manga, 2011 [1]). Because of the interplay between the different factors, it is difficult to predict the run-out distances of PDCs.

AC

CE

PT E

D

MA

NU

The complex interplay between the various factors comes from the fact that a PDC has an extremely wide range of particle concentrations. A PDC generally becomes stratified as the particle concentration and density show a remarked increase toward the base of the current (Valentine, 1987 [70]; Branney and Kokelaar, 2002 [5]; Burgisser and Bergantz, 2002 [8]; Breard et al., 2016 [7]). It is composed of a voluminous dilute turbulent suspension with a low particle concentration (≲ 1 vol.%) over a thin dense fluidized granular flow with a high particle concentration (∼ 40–50 vol.%). The dilute-part flow is driven by the density contrast with the ambient air, and its density decreases as particles settle and ambient air is entrained, heated, and expanded during transport. Finally, the frontal region of the dilute part can become lighter than the ambient air, so that its buoyancy lifts it off the ground (e.g., Sparks et al., 1993 [59]; Andrews and Manga, 2012 [2]). The dynamics of the dense part, on the other hand, are controlled by the different physical processes from those for the overlying dilute part; i.e., it is affected mainly by particle-particle collision, frictional interaction between the current and the ground, and deposition/erosion processes at the base (e.g., Lube et al., 2007 [45]; Roche et al., 2008 [53]; Girolami et al., 2010 [29]; Roche et al., 2013 [54]). The deposition reduces progressively the mass of particles moving above the static-mobile interface that migrates upward, and the dense-part flow stops when the mobile region is consumed (cf. Roche, 2012 [51]). The run-out distance of a PDC with strong density stratification depends on the dynamics of both the dilute and dense parts (e.g., lift-off of the dilute part and deposition of the dense part), and the interactions between them (e.g., particle transfer from one to the other). The essential mechanisms controlling the run-out distance of a PDC with strong density stratification can be examined using a two-layer dynamical model (Doyle et al., 2008 [20], 2010 [21], 2011 [19]; Kelfoun, 2017 [38]). The two layers, an upper dilute current over a dense basal current, evolve separately but are coupled through mass and momentum exchanges such as inter-layer particle transfer. Therefore, the above factors in both currents (e.g., lift-off of the dilute current and deposition of the dense current), as well as the interactions between the currents, can be evaluated in the two-layer model. The previous two-layer models showed that the strong density stratification causes a wide variation in the run-out distance of small-scale PDCs generated by collapses of Vulcanian eruption column and lava dome. This paper is concerned with the run-out distance of large-scale PDCs generated by collapse

2

ACCEPTED MANUSCRIPT Journal of Volcanology and Geothermal Research

RI PT

of Plinian eruption column. Bursik and Woods (1996) [9] developed a dynamical model that calculates the run-out distance of a large-scale PDC flowing on a (nearly) flat surface as a function of eruption conditions. Assuming the PDC to be a fully dilute current, they proposed that the distance at which the current lifts off the ground owing to particle settling and thermal expansion of entrained air is regarded as the run-out distance, and estimated how the runout distance of the large-scale PDC depends on magma discharge rate at the source. The relationship between run-out distance and magma discharge rate is useful to estimate eruption conditions (e.g., the intensity of explosive eruption) from field observations (e.g., the extent of PDC deposits). However, the model requires further assessment, as it does not consider the effects of strong density stratification (i.e., the presence of the dense basal current).

2.1

Method Basic equations

PT E

2

D

MA

NU

SC

Here, we aim to discuss how the previously obtained relationship between run-out distance and magma discharge rate is modified by the strong density stratification, based on a twolayer model. The previous two-layer models neglect the effects of the entrainment and thermal expansion of ambient air in the dilute current, so that they cannot be meaningfully compared with the model of Bursik and Woods (1996) [9]. We therefore develop a new two-layer PDC model that considers the effects of the entrainment and thermal expansion of ambient air in the upper dilute current, in addition to the deposition in the lower dense current and the mass and momentum transfer between the two currents. We carry out an extensive parametric study to clarify the dominant factors controlling the run-out distance of the two-layer PDC (e.g., to identify the parameters controlling which of lift-off of the dilute current and deposition of the dense current determines the run-out distance of the whole current). Although the present model deals with PDC dynamics on an idealized flat surface for comparison with the model of Bursik and Woods (1996) [9], we briefly discuss possible topographic effects on the run-out distance for large-scale PDCs.

AC

CE

Large-scale PDCs commonly form near-symmetrical, circular distributions of deposits around the source vent (e.g., Branney and Kokelaar, 2002 [5]). We therefore design a model to describe an axisymmetric PDC spreading radially from a central source vent (Fig. 1). In the model, an eruption column consisting of volcanic particles, volcanic gas, and entrained air collapses to generate a dilute current from the column edge at a constant mass flow rate on a horizontal ground surface. Although the processes associated with column collapse are complex (e.g., Sweeney and Valentine, 2017 [66]; Valentine and Sweeney, 2018 [71]), the collapsing column is tentatively assumed to be composed of a homogeneous dilute mixture for comparison with the model of Bursik and Woods (1996) [9]. The two layers used in the present model to describe the dynamics of PDC with strong density stratification are separately evolving depth-averaged dilute and dense currents that are coupled through mass and momentum exchanges as suspended particles in the dilute current settle into the dense current (Doyle et al., 2008 [20], 2010 [21], 2011 [19]). A deposit progressively aggrades upward from the base of the dense and/or dilute currents (e.g., Branney and Kokelaar, 2002 [5]; Lube et al., 2007 [45]; Girolami et al., 2010 [29]; Roche, 2012 [51]). The basic equations of the dilute and dense currents and the deposit are shown below; some of their details and the numerical procedures for solving them are given in Appendix A. Summaries of the constant and varied parameters used in this paper are listed in Tables 1 and 2, respectively.

3

ACCEPTED MANUSCRIPT

Edge of collapsing eruption column

(b) z zf

(c) . M0

zc zb 0

Dilute PDC

h

Dense PDC Deposit

hH r

0

r

PT E

r0

D

MA

NU

t=0

SC

(a)

RI PT

Journal of Volcanology and Geothermal Research

AC

CE

Figure 1: Illustration of the two-layer pyroclastic density current (PDC) model. (a) A PDC is generated by a steady supply of a dilute current from the edge of the eruption column collapsing above the volcanic vent (i.e., r = r0 ), and spreads radially from the source over a horizontal ground surface. (b) The dilute current with thickness h (red region) forms a dense basal current with thickness hH (blue region). A deposit (gray region) progressively aggrades upward from the base of the dense current (and/or the dilute current). (c) The supply of a homogeneous dilute current with a constant mass flow rate M˙ 0 starts suddenly at t = 0.

4

ACCEPTED MANUSCRIPT

RI PT

Journal of Volcanology and Geothermal Research

Table 1: Common input parameters and constants for simulations.

na0 a0 Ri0 ε

CE

ns0

AC

MA

NU

SC

Meaning Gravitational body force per unit mass Gas constant of volcanic gas Gas constant of air Specific heat of solid particles Specific heat at constant pressure of volcanic gas Specific heat at constant pressure of air Temperature of ambient air Density of solid particles (In Appendix F) Density of ambient air Imposed frontal Froude number Volume fraction of solid particles in dense current Ch´ezy drag coefficient of dilute current Ch´ezy drag coefficient of dense current Temperature of dilute current at r0 (In Appendix F) Mass fraction of solid particles in dilute current at r0 (In Appendix F) Mass fraction of entrained air in dilute current at r0 (In Appendix F) Aspect ratio of h0 to r0 (In Appendix F) Richardson number of dilute current at r0 (In Appendix D) Non-dimensional thickness of artificial bed

D

ρa F rN ϕsH Cdc Cdb T0

Value [unit] 9.81 [m/s2 ] 462 [J/(kg·K)] 287 [J/(kg·K)] 1100 [J/(kg·K)] 1810 [J/(kg·K)] 1004 [J/(kg·K)] 273 [K] 1000 [kg/m3 ] (1000 or 2500) 1.23 [kg/m3 ] 1.19 0.4 10−4 10−4 950 [K] (800 or 950) 0.9 (0.9 or 0.97) 0.07 (0 or 0.07) 0.2 (0.1 or 0.2) 1 (0.99 or 1.01) 10−10

PT E

Variable g Rv Ra Cs Cpv Cpa Ta ρs

5

ACCEPTED MANUSCRIPT Journal of Volcanology and Geothermal Research

CE

D

Ws /(u0 a0 ) 1.85 1.16 0.735 0.464 0.293 0.185 0.116 0.0735 0.0735 0.0464 0.0293 0.0185 0.0116 0.00735 1.85 1.16 0.735 0.464 0.293 0.185 0.116 0.0735 0.0464 0.0293 0.0185 0.0116 0.00735

NU

SC

D [m/s] 3.0 × 10−3 3.0 × 10−3 3.0 × 10−3 3.0 × 10−3 3.0 × 10−3 3.0 × 10−3 3.0 × 10−3 3.0 × 10−3 3.0 × 10−4 3.0 × 10−4 3.0 × 10−4 3.0 × 10−4 3.0 × 10−4 3.0 × 10−4 9.0 × 10−3 9.0 × 10−3 9.0 × 10−3 9.0 × 10−3 9.0 × 10−3 9.0 × 10−3 9.0 × 10−3 9.0 × 10−3 9.0 × 10−4 9.0 × 10−4 9.0 × 10−4 9.0 × 10−4 9.0 × 10−4

MA

Ws [m/s] 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 0.3 0.3 0.3 0.3 0.3 0.3 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 0.3 0.3 0.3 0.3 0.3

PT E

M˙ 0 [kg/s] 105 106 107 108 109 1010 1011 1012 107 108 109 1010 1011 1012 105 106 107 108 109 1010 1011 1012 108 109 1010 1011 1012

AC

Simulation Run 1 Run 2 Run 3 Run 4 Run 5 Run 6 Run 7 Run 8 Run 9 Run 10 Run 11 Run 12 Run 13 Run 14 Run 15 Run 16 Run 17 Run 18 Run 19 Run 20 Run 21 Run 22 Run 23 Run 24 Run 25 Run 26 Run 27

RI PT

Table 2: Input parameters and flow regimes.

6

D/Ws 1 × 10−3 1 × 10−3 1 × 10−3 1 × 10−3 1 × 10−3 1 × 10−3 1 × 10−3 1 × 10−3 1 × 10−3 1 × 10−3 1 × 10−3 1 × 10−3 1 × 10−3 1 × 10−3 3 × 10−3 3 × 10−3 3 × 10−3 3 × 10−3 3 × 10−3 3 × 10−3 3 × 10−3 3 × 10−3 3 × 10−3 3 × 10−3 3 × 10−3 3 × 10−3 3 × 10−3

Regime 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

ACCEPTED MANUSCRIPT

Table 2: (continued).

CE

D

Ws /(u0 a0 ) 1.85 1.16 0.735 0.464 0.293 0.185 0.116 0.0735 0.0464 0.0293 0.0185 0.0116 0.00735 1.85 1.16 0.735 0.464 0.293 0.185 0.116 0.0735 0.0735 0.0464 0.0293 0.0185 0.0116 0.00735

NU

SC

D [m/s] 1.5 × 10−2 1.5 × 10−2 1.5 × 10−2 1.5 × 10−2 1.5 × 10−2 1.5 × 10−2 1.5 × 10−2 1.5 × 10−2 1.5 × 10−3 1.5 × 10−3 1.5 × 10−3 1.5 × 10−3 1.5 × 10−3 3.0 × 10−2 3.0 × 10−2 3.0 × 10−2 3.0 × 10−2 3.0 × 10−2 3.0 × 10−2 3.0 × 10−2 3.0 × 10−2 3.0 × 10−3 3.0 × 10−3 3.0 × 10−3 3.0 × 10−3 3.0 × 10−3 3.0 × 10−3

MA

Ws [m/s] 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 0.3 0.3 0.3 0.3 0.3 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 0.3 0.3 0.3 0.3 0.3 0.3

PT E

M˙ 0 [kg/s] 105 106 107 108 109 1010 1011 1012 108 109 1010 1011 1012 105 106 107 108 109 1010 1011 1012 107 108 109 1010 1011 1012

AC

Simulation Run 28 Run 29 Run 30 Run 31 Run 32 Run 33 Run 34 Run 35 Run 36 Run 37 Run 38 Run 39 Run 40 Run 41 Run 42 Run 43 Run 44 Run 45 Run 46 Run 47 Run 48 Run 49 Run 50 Run 51 Run 52 Run 53 Run 54

RI PT

Journal of Volcanology and Geothermal Research

7

D/Ws 5 × 10−3 5 × 10−3 5 × 10−3 5 × 10−3 5 × 10−3 5 × 10−3 5 × 10−3 5 × 10−3 5 × 10−3 5 × 10−3 5 × 10−3 5 × 10−3 5 × 10−3 1 × 10−2 1 × 10−2 1 × 10−2 1 × 10−2 1 × 10−2 1 × 10−2 1 × 10−2 1 × 10−2 1 × 10−3 1 × 10−2 1 × 10−2 1 × 10−2 1 × 10−2 1 × 10−2

Regime 2a 2a 2a 2a 2a 2a 2a 2a 2a 2a 2a 2a 2a 2b 2b 2b 2b 2b 2b 2b 2b 2b 2b 2b 2b 2b 2b

ACCEPTED MANUSCRIPT Journal of Volcanology and Geothermal Research

2.1.1

Dilute PDC

We model the dilute current as a highly turbulent suspension flow consisting of solid particles, volcanic gas, and entrained ambient air (e.g., Bursik and Woods, 1996 [9]). The depth-averaged conservation equations of mass and momentum at radius r and time t in a dilute current with thickness h(r, t), velocity u(r, t), and density ρ(r, t) are (1) (2)

RI PT

∂ 1 ∂ ns (ρh) + (ρuhr) = ρa E|u| − ρWs , and ∂t r ∂r nsH ( ) ∂ 1 ∂ ( 2 ) ∂ ρ − ρa 2 ∂zc ns (ρuh) + ρu hr + gh = −(ρ − ρa )gh − ρuWs − τc , ∂t r ∂r ∂r 2 ∂r nsH

D

MA

NU

SC

respectively, where ρa is the density of ambient air, E is the entrainment coefficient, ns is the mass fraction of solid particles in the dilute current, nsH is the mass fraction of solid particles in the dense basal current, Ws is the averaged settling speed of solid particles at the base of the dilute current, g is the gravitational body force, zc is the height of the basal contact, and τc is the basal drag. The effect of air entrainment is taken into account in the mass equation (i.e., the first term on the right-hand side (RHS) of Eq. (1)). The effect of particle settling is considered in both the mass and momentum equations (i.e., the second terms on the RHSs of Eqs. (1) s ρWs is the sum of the settling rate of particles (ns ρs Ws ) and that of the and (2)); here nnsH 1−nsH interstitial gas ( nsH ns ρWs ). The effect of basal drag is represented in the momentum equation (i.e., the third term on the RHS of Eq. (2)). The first term on the RHS of Eq. (2) represents the acceleration over the basal slope generated by the presence of the dense current. Note that the present model is time-dependent, unlike the steady-state model of Bursik and Woods (1996) [9]. Depth-averaged equations (i.e., Eqs. (1) and (2)) provide a good approximation of the dynamics of currents where turbulent mixing is sufficiently intense to maintain vertically uniform concentration (e.g., Bonnecaze et al. 1993 [3]). See Appendix A.1 for details of E, Ws , and τc .

PT E

The bulk density ρ(r, t) is given by the thermal equation of state (e.g., Bursik and Woods, 1996 [9]): ns T 1 = + (na Ra + (1 − ns − na )Rv ) , ρ ρs p

(3)

AC

CE

where ρs is the density of solid particles, Ra is the gas constant of air, Rv is the gas constant of volcanic gas, and T is the temperature of the dilute current. Here, p is the hydrostatic pressure of the dilute current, which is written as p = ρa Ra Ta −ρa gzf + 12 ρgh, where Ta is the temperature of ambient air and zf is the height of the free surface of the dilute current. In this study, the value of p is approximated as the atmospheric pressure ρa Ra Ta in Eq. (3) for simplicity (e.g., Bursik and Woods, 1996 [9]). The mass fraction of entrained air, na (r, t), and that of solid particles, ns (r, t), evolve temporally and spatially on the basis of their mass conservation equations: ∂ 1 ∂ (na ρh) + (na ρuhr) = ρa E|u|, and ∂t r ∂r ∂ 1 ∂ (ns ρh) + (ns ρuhr) = −ns ρWs , ∂t r ∂r

(4) (5)

respectively. To consider the change in temperature T (r, t) due to air entrainment, we solve the depth-averaged enthalpy conservation equation of the dilute current: ∂ 1 ∂ 1 − nsH (ρCp T h) + (ρCp T uhr) = ρa E|u|Cpa Ta − ns ρWs Cs T − ns ρWs Cpv T ∂t r ∂r nsH u2 gh gh ns + ρa E|u| + ρa E|u| + ρWs , (6) 2 2 2 nsH 8

ACCEPTED MANUSCRIPT Journal of Volcanology and Geothermal Research

where Cs is the specific heat of solid particles, Cpa is the specific heat at constant pressure of air, Cpv is the specific heat at constant pressure of volcanic gas, and Cp (r, t) is the specific heat at constant pressure of the dilute current, which is given by Cp = ns Cs + na Cpa + (1 − ns − na )Cpv .

(7)

RI PT

The first term on the RHS of Eq. (6) is the thermal energy added by air entrainment; the second and third terms are the losses of thermal energies of settling particles and volcanic gas, respectively. The fourth term on the RHS represents the rate at which kinetic energy is transferred to the entrained air. The fifth and sixth terms on the RHS represent the changes in thermal energy by conversion from potential energy. Eq. (6) is derived from the depth-averaged total energy conservation equation (cf. Bursik and Woods, 1996 [9]; see Appendix B for details).

NU

SC

To describe realistic dilute PDC dynamics, a balance between the buoyancy pressure driving the flow front and the resistance pressure caused by the acceleration of the ambient air at the front of the dilute current should be considered (Ungarish, 2009 [69]; Shimizu et al., 2017 [56]). This dynamical balance is written as the front condition: √ ρN − ρa drN = uN = F rN ghN at r = rN (t), (8) dt ρa

MA

where the subscript N denotes the front, and F rN is the imposed frontal Froude number and is set to 1.19 (Huppert and Simpson, 1980 [36]).

Dense PDC

PT E

2.1.2

D

When the dilute current becomes lighter than the ambient air, it lifts off the ground to produce a co-ignimbrite ash plume (e.g., Andrews and Manga, 2012 [2]). The present model does not calculate the dynamics of the dilute current after lift-off, as the set of Eqs. (1)–(8) is valid only where the dilute current is denser than the ambient air.

CE

The dense basal current, generated by particles settling from the dilute current, is modeled as a homogeneous fluidized granular flow consisting of solid particles and volcanic gas. The depthaveraged conservation equations of mass and momentum in a dense current with thickness hH (r, t) and velocity uH (r, t) are

AC

∂hH 1 ∂ ns ρ + (uH hH r) = Ws − D, and (9) ∂t r ∂r nsH ρH ( ) ) ∂ 1 ∂ ( 2 ∂ 1 ρH − ρa 2 (uH hH ) + uH hH r + ghH ∂t r ∂r ∂r 2 ρH ∂zb hH ∂ ns ρ τc τb ρH − ρa ghH − ((ρ − ρa )gh) + uWs − uH D + − , (10) =− ρH ∂r ρH ∂r nsH ρH ρH ρH respectively, where D is the deposition speed at the base, zb is the height of the contact between the dense current and the deposit, and τb is the basal drag. The dense current is assumed to have a constant bulk density ρH = ρs ϕsH + ρgH (1 − ϕsH ), where ϕsH (≡ nsH ρH /ρs ) is the volume fraction of solid particles, ρgH (≡ p/(Rv T0 )) is the density of the gas phase, and T0 is the initial temperature of the overlying dilute current (see Subsection 2.2 below for details). The particle supply from the overlying dilute current and the deposition at the base of the dense current are accounted for in the mass and momentum equations (i.e., the RHS of Eq. (9), and the third and fourth terms on the RHS of Eq. (10)). On the RHS of Eq. (10) the first term represents the acceleration over the basal slope generated by the deposition, and the second term represents 9

ACCEPTED MANUSCRIPT Journal of Volcanology and Geothermal Research

the pressure gradient on the dense current exerted by variations in the load of the dilute current (Doyle et al., 2011 [19]); the fifth and sixth terms are the interfacial drag with the overlying dilute current and the basal drag, respectively.

2.1.3

RI PT

For the dense current, the resistance of ambient air at the flow front is not explicitly considered, because the momentum lost due to frontal resistance is negligibly small with respect to the current’s momentum (Shimizu et al., 2017 [56]). The particle transfer from the dense current to the overlying dilute current is tentatively ignored in the present model, but it can generate a dilute ash cloud (cf. Fisher, 1979 [28]; Nakada and Fujii, 1993 [49]) and may play an important role in PDC dynamics (e.g., Kelfoun, 2017 [38]); we will assess its effect in later sections (Sections 3 and 5). Deposit

PT E

D

MA

NU

SC

A deposit progressively aggrades upward from the base of the dense or dilute current (e.g., Branney and Kokelaar, 2002 [5]; Lube et al., 2007 [45]; Girolami et al., 2010 [29]; Roche, 2012 [51]). The mass conservation equation of the dense current (i.e., Eq. (9)) indicates that the deposit is formed either by aggradation at the base of the dense current or directly by aggradation at the base of the dilute current. The aggradation of the deposit from the dilute current occurs when the particle-settling rate at the base of the dilute current is smaller than the deposition rate of the dense current at the position where the dense current is absent (i.e., the two conditions that the RHS of Eq. (9) < 0 and hH = 0 are simultaneously satisfied). Assuming that the deposit is not erodible and has the same bulk density as the dense current, we can write the aggradation rate of material in the deposit as  D (Aggradation from the dense current) ∂zb  = (11) ns ρ  ∂t Ws (Aggradation from the dilute current).  nsH ρH Given the assumption of similar properties to the dense current, the mass fraction of solid particles and the density of the deposit are the same terms as used above (in Eqs. (9) and (10)) for the dense current (nsH and ρH , respectively).

AC

CE

In our model (i.e., Eq. (11)), the deposition at the base of the dense current reduces progressively the mass of particles moving above the static-mobile interface that migrates upward (i.e., progressive aggradation; Branney and Kokelaar, 1992 [4]), and the dense current stops when the mobile region is consumed. Generally, the processes of basal deposition also depend on the rheology of dense current and the steepness of slope of the ground surface. For example, in some previous models where the basal deposition is not taken into account (e.g., Doyle et al., 2008 [20], 2011 [19]; Kelfoun, 2017 [38]), the dense current stops to form the deposits when its momentum is consumed due to the basal drag modeled by a Coulomb or plastic rheology (i.e., en masse deposition; Sparks, 1976 [58]). The en masse deposition appears to be qualitatively different from the progressive aggradation; however, it can be interpreted as a special case of the progressive aggradation (Roche, 2012 [51]). The diverse depositional processes between these two cases can be expressed by a combination of the basal drag τb modeled by a Ch´ezy-type drag (see Eq. (A.3)) and the basal deposition (i.e., the deposition speed D described in Eq. (11)) on the basis of the experimental studies by Roche et al. (2008) [53]. The possible ranges of D and the limitations of this model are discussed in Appendix A.1.

10

ACCEPTED MANUSCRIPT Journal of Volcanology and Geothermal Research

2.2

Boundary and initial conditions

ns0 Cs Tm + (1 − ns0 − na0 )Cpv Tm + na0 Cpa Ta , Cp0

(12)

SC

T0 =

RI PT

We assume that a homogeneous dilute current is steadily supplied during time t ≥ 0 from the edge of the collapsing eruption column (i.e., r = r0 ) as a boundary condition (Fig. 1; cf. Bursik and Woods, 1996 [9]). We are concerned with the PDCs generated during magmatic eruptions; the dilute current at the column edge consists of magma (i.e., the volcanic particles and gas erupted from the volcanic vent) and entrained ambient air. Given the mass fractions of entrained air and solid particles at the column edge (i.e., na0 and ns0 , respectively), the specific heat at constant pressure at the column edge, Cp0 , is determined by Eq. (7). Note that the mass fraction of volcanic gas at the column edge is given by 1 − ns0 − na0 ; i.e., the mass fraction of volcanic gas in the magma is written as (1 − ns0 − na0 )/(1 − na0 ). The temperature of the dilute current at the column edge, T0 , is estimated by conservation of heat between the magma of temperature Tm and the entrained air of temperature Ta (e.g., Bursik and Woods, 1996 [9]):

D

MA

NU

where ns0 Cs Tm is the enthalpy of solid particles, (1 − ns0 − na0 )Cpv Tm is the enthalpy of volcanic gas, and na0 Cpa Ta is the enthalpy of entrained air. Given the density of solid particles, ρs , the bulk density of the dilute current at the column edge, ρ0 , is determined by the thermal equation of state (i.e., Eq. (3)); the volume fraction of solid particles at the column edge, ϕs0 , is given by ns0 ρ0 /ρs . We specify na0 , ns0 , T0 , and ρs , and use them to calculate the other material properties of the dilute current at the collapsing column edge (i.e., Cp0 , ρ0 , and ϕs0 ). In this paper, we fix na0 , ns0 , T0 , and ρs at typical values used in Bursik and Woods (1996) [9] (na0 = 0.07, ns0 = 0.9, T0 = 950 K, and ρs = 1000 kg/m3 ; see Table 1) roughly corresponding to ρ0 = 2.98 kg/m3 and ϕs0 = 2.68 × 10−3 . These values also roughly correspond to the mass fraction of volcanic gas in magma, (1 − ns0 − na0 )/(1 − na0 ) = 0.032, and the temperature of the magma, Tm = 1000 K.

AC

CE

PT E

The other physical values of the dilute current at the collapsing column edge (i.e., the thickness h0 , velocity u0 , and radius of the collapsing column r0 ) are obtained from the mass flow rate at the column edge, M˙ 0 (≡ 2πr0 ρ0 u0 h0 ), the Richardson number at the column edge, Ri0 (≡ (ρ0 − ρa )gh0 /(ρ0 u20 )), and the aspect ratio a0 (≡ h0 /r0 ):  ( )2 / ( )1/5  a M˙ ρ0 − ρa 0 0 , (13) g/Ri0 h0 =   2πρ0 ρ0 u0 =

r0 =

{

( )2 }1/5 a0 M˙ 0 ρ0 − ρa g/Ri0 , 2πρ0 ρ0 (  )2 / ( )1/5 1  a0 M˙ 0 ρ0 − ρa g/Ri0 .  a0  2πρ0 ρ0

(14)

(15)

In this paper, we fix Ri0 = 1 and a0 = 0.2 (as discussed in Appendix C; see Table 1), such that M˙ 0 becomes the only free parameter at the source boundary. Because we set na0 = 0.07, M˙ 0 can be directly related to magma discharge rate, M˙ m , as M˙ 0 = M˙ m /0.93. As mentioned above, we assume that the collapsing column generates a homogeneous dilute mixture in this paper for comparison with the model of Bursik and Woods (1996) [9]. On the other hand, some recent studies have shown that the behavior of particles in a dilute collapsing gas-particle mixture is complex so that the conditions of PDCs at the column edge can deviate from what is assumed here; for example, particles could rebound and/or concentrate at the 11

ACCEPTED MANUSCRIPT Journal of Volcanology and Geothermal Research

base of the collapsing column to generate a dense basal current at the column edge (Sweeney and Valentine, 2017 [66]; Valentine and Sweeney, 2018 [71]). The effects of such processes in the collapsing column on the PDC dynamics can be investigated using our model by changing conditions at the column edge; a further parametric study for this purpose is in progress.

2.3

Non-dimensional parameters

SC

RI PT

The above basic equations are non-dimensionalized with the quantities at the collapsing column edge (see Appendix C). The dimensionless basic equations have a number of non-dimensional parameters. In this paper, we fix all the non-dimensional parameters determined by the material properties (i.e., na0 , ns0 , nsH , ρa /ρ0 , ρH /ρ0 , Cpa Ta /(Cp0 T0 ), Cs /Cp0 , and Cpv /Cp0 ) at their typical values. We also fix the basal drag coefficients of the dilute and dense currents (i.e., Cdc and Cdb ) and the imposed frontal Froude number F rN on the basis of laboratory experiments (see Appendix C for details). When the values of the Richardson number Ri0 and the aspect ratio a0 are fixed as mentioned above, Ws /(u0 a0 ), D/Ws , Eu0 /Ws , and gh0 /(2Cp0 T0 ) become variable non-dimensional parameters.

PT E

D

MA

NU

The above four non-dimensional parameters have the following physical meanings. The parameter Ws /(u0 a0 ) is a time-scale ratio that compares the propagation time of the current for a distance r0 (i.e., r0 /u0 ) with the particle settling time for a height h0 (i.e., h0 /Ws ). The parameter D/Ws represents the relative magnitude of the effect of the deposition at the base of the dense current to that of the particle supply from the overlying dilute current to the dense current. The parameter Eu0 /Ws represents the relative magnitude of the effect of air entrainment at the top of the dilute current to that of particle settling at the base of the dilute current. The parameter gh0 /(2Cp0 T0 ) represents the relative magnitude of mechanical energy to thermal energy (i.e., enthalpy), which is negligibly small under the condition of the mass flow rate M˙ 0 ≲ 1012 kg/s. This means that the fourth to sixth terms on the RHS of the enthalpy conservation equation of the dilute current (Eq. (6)) are negligible. Thus, Ws /(u0 a0 ), D/Ws , and Eu0 /Ws are the key parameters controlling the dynamics of the two-layer PDC. It is important to note that Ws /(u0 a0 ) is closely related to the eruption condition. From Eq. (14) it is given as a function of M˙ 0 as

CE

Ws Ws = u0 a0 a0

{

a0 M˙ 0 2πρ0

(

ρ0 − ρa g/Ri0 ρ0

)2 }−1/5 .

(16)

AC

In order to investigate how the present two-layer model modifies the run-out distance of the one-layer dilute PDC model of Bursik and Woods (1996) [9], we perform a parametric study for different values of Ws /(u0 a0 ) and D/Ws ; the value of Eu0 /Ws is calculated from the flow dynamics (see Eq. (A.1)). We set the range of the particle settling speed Ws as 0.3–3.0 m/s on the basis of terminal velocities of variously sized solid particles (see Appendix A.1). Considering the above range of Ws and a wide range of the mass flow rate at source, M˙ 0 = 105 –1012 kg/s, we set Ws /(u0 a0 ) as ranging from 7 × 10−3 to 2 × 100 (see Eq. (16)). It should be noted that the range of Ws /(u0 a0 ) reflects mainly that of M˙ 0 , although it also depends on the variation in Ws . The deposition speed D may be estimated from the hindered settling speed (i.e., Eq. (A.5); e.g., Cao et al., 2004 [10]; Li et al., 2018 [43]). We set various values of D/Ws ranging from 10−3 to 10−2 on the basis of Eq. (A.5) (see Appendix A.1). Summaries of the constant and varied parameters are listed in Tables 1 and 2, respectively.

12

ACCEPTED MANUSCRIPT Journal of Volcanology and Geothermal Research

3 3.1

Results General features

Our numerical simulations show that a two-layer PDC spreads radially from the collapsing column edge (i.e., r = r0 ), stops forward propagation because of lift-off of the dilute current and/or deposition of the dense current, and eventually converges to a steady state. A typical result for a two-layer PDC is shown in Fig. 2 (Run 5 with Ws /(u0 a0 ) = 0.292 and D/Ws = 1.0 × 10−3 ; see Supplementary Movie S1 for details).

MA

NU

SC

RI PT

In Fig. 2, after the dilute current is generated from the column edge, it forms a thick flow “head” followed by the flow “body” owing to frontal resistance of ambient air (see Fig. 2a; cf. Ungarish, 2009 [69]; Shimizu et al., 2017 [56]). As the dilute current spreads out from the source, the current density decreases through particle settling, air entrainment, and thermal expansion of the entrained air. The rate of the density decrease due to particle settling is smaller in the head than in the body, because the former is thicker than the latter (cf. Bonnecaze et al., 1993 [3]). Furthermore, the rate of entrainment is smaller in the head than in the body, because the former is slower and therefore has larger Ri than the latter (see Eq. (A.1)). For these reasons, the current density is smallest in the body just behind the head. As this minimum density becomes smaller than the ambient density, the current begins to lift off at this location, leaving an isolated head (see Fig. 2b). The density of the isolated head decreases with time and distance; as a result, it fully lifts off and disappears (see Fig. 2c and d). Finally, the body converges to a steady state (see Fig. 2b–f), the profile of which agrees with the supercritical steady solution proposed by Bursik and Woods (1996) [9] (see Appendix D).

AC

CE

PT E

D

The solid particles settling from the dilute current form a thin, dense, basal current. The dense current initially forms a thick flow head followed by the flow body (see Fig. 2a–d) because of the complex combinations of the mass and momentum supplied from the overlying dilute current as well as the radial fluxes of mass and momentum in the dense current. Subsequently, the dense current evolves downstream until its radial mass flux becomes zero. The radial mass flux of the dense current (i.e., the left-hand side (LHS) of Eq. (9)) is determined by two competing effects: it increases through mass supply from the dilute current (i.e., the first term on the RHS of Eq. (9)) and decreases through basal deposition (i.e., the second term on the RHS of Eq. (9)). In the distal area where the dilute current is absent, the radial mass flux decreases with distance due to deposition in the body, and is smallest just behind the head. As this minimum mass flux balances the deposition rate, the flow body stops forward propagation and converges to a steady state, leaving an isolated head (see Fig. 2e and f). Finally, the isolated head fully disappears (see Fig. 2f; see also Supplementary Movie S1). The results show that the run-out distance of each of the dilute and dense currents converges to a steady-state value (referred to as the steady run-out distance) after the transient “overshoot” (see Fig. 3). The maximum run-out distance of each current is determined by the farthest point reached by the isolated head. The difference between the maximum and steady run-out distances is referred to as the overshoot distance. The overshoot occurs because of the source boundary condition that a steady supply of a dilute current starts suddenly at t = 0 (see Fig. 1c). When the eruption duration is much longer than the time required for the two-layer PDC to reach a steady state, the distributions of the major deposits are determined by the steady-state current. Accordingly, although the evaluation of the overshoots is important from the viewpoint of hazard mitigation, herein we focus on the steady-state flow properties of two-layer PDCs.

13

ACCEPTED MANUSCRIPT Journal of Volcanology and Geothermal Research

(a) Time, t /T = 1.0 1.5

(b) Time, t /T = 1.3 1.5

Dilute current

Isolated head

0 0.003

0 0.003

Body Head Dense current

0 0.006

zb / h0

Deposit

0.003 2 3 Distance, r / r 0

0

4

h / h0

NU MA

0 0.003

D

0 0.006

1

PT E

0.003 2 3 Distance, r / r 0

1

0 0.003 0 0.006 0.003 0

4

1

(e) Time, t /T = 4.5

AC

0 0.003

h / h0 Isolated head

0 0.006 0.003 1

1 0.5

zb / h0

hH / h0

0.5

4

(f) Time, t /T = 8.0

hH / h0

1

2 3 Distance, r / r 0

1.5

CE

h / h0

1.5

zb / h0

4

0.5

zb / h0

h / h0

1

0

2 3 Distance, r / r 0

1.5

0.5

hH / h0

1

(d) Time, t /T = 2.5

(c) Time, t /T = 2.0 1.5

0

0.003

SC

1

0 0.006

RI PT

0.5

hH / h0

0.5

0

zb / h0

Lift-off

1

h / h0

1

hH / h0

zb / h0

hH / h0

h / h0

Body Head

2 3 Distance, r / r 0

0 0.003 0 0.006 0.003 0

4

1

2 3 Distance, r / r 0

4

Figure 2: Representative numerical results of a two-layer PDC at scaled times t/T = (a) 1.0, (b) 1.3, (c) 2.0, (d) 2.5, (e) 4.5, and (f) 8.0 from the beginning of propagation in Run 5. Scaled thicknesses of the dilute (h/h0 ; red) and dense (hH /h0 ; blue) currents and the deposit (zb /h0 ; black) are shown. 14

ACCEPTED MANUSCRIPT

Overshoot distance

4 3 Overshoot distance

Steady run-out distance

2 1

Steady run-out distance

Dense current Dilute current

0 1

2

3

4

5

7

8

NU

Time, t / T

6

SC

0

RI PT

Front position, r / r0

Journal of Volcanology and Geothermal Research

3.2

MA

Figure 3: Front positions of the dilute (red) and dense (blue) currents scaled with r0 as a function of scaled time t/T for Run 5.

Steady-state flow properties

PT E

D

The flow patterns of steady-state two-layer PDCs are classified depending on whether the steady run-out distance of the dense current is longer or shorter than that of the dilute current (referred to as “Regime 1” or “Regime 2”, respectively). Regime 2 includes two sub-regimes that either have (“sub-Regime 2a”) or lack (“sub-Regime 2b”) the dense current.

AC

CE

Fig. 4 shows representative steady-state profiles for Regime 1 (Run 5 with Ws /(u0 a0 ) = 0.292 and D/Ws = 1.0 × 10−3 ; see Fig. 2) and Regime 2 (Run 32 with Ws /(u0 a0 ) = 0.292 and D/Ws = 5.0 × 10−3 ; see also Supplementary Movie S2). Given that the parameters are identical for these two runs, other than D/Ws , they have identical profiles for the dilute current. On the other hand, the differences in the dynamics of their dense currents are explained by the difference in D/Ws ; a small value results in Regime 1, which has a thicker dense current than Regime 2 (see Fig. 4a). The velocities of the dense currents of both regimes are similar (see Fig. 4b), reflecting that they are determined mainly by the momentum of particles supplied from the dilute current. The velocities are slower than that of the dilute current because of basal drag at the dense current. In the distal area where the dilute current is absent, the velocity of the dense current decreases with distance owing to the basal drag. The above variation in the thickness and velocity of the dense current can be explained by the changes with distance of the radial mass flux of particles in the dense current (i.e., the LHS of Eq. (9)) (Fig. 4c). In Regime 1, the radial mass flux of the dense current increases with distance in the proximal area where the dilute current exists, because the rate of mass supply from the dilute current (i.e., the first term on the RHS of Eq. (9)) exceeds that of basal deposition (i.e., the second term on the RHS of Eq. (9)). In the distal area without the dilute current, the radial mass flux of the dense current decreases with distance owing to basal deposition. On the other hand, in Regime 2 the radial mass flux of the dense current decreases with distance, even in the 15

ACCEPTED MANUSCRIPT Journal of Volcanology and Geothermal Research

proximal area where the dilute current exists, as D/Ws is large. In Regime 2, the deposit is formed directly from the dilute current in the distal area.

RI PT

Note that the majority of particle mass is transported by the dilute current in the proximal area for both regimes (Fig. 4c). The major difference between them is that the particles are transported farther by the dense current in the distal area in Regime 1. In reality, elutriation from a part of this dense current in the distal area generates a new dilute ash cloud. However, the radial mass flux of particles in such an ash cloud is considered to be much smaller than that in the parent dense current. Thus, our classification based on the radial mass flux of particles (i.e., Fig. 4c) is robust.

Discussion

NU

4

SC

On the basis of an extensive parametric study, we constructed a flow-regime map in the parameter space of the two given non-dimensional parameters Ws /(u0 a0 ) and D/Ws (Fig. 5). The figure shows that the transition between the two regimes is determined by D/Ws , and it is independent of Ws /(u0 a0 ). As mentioned above, Ws /(u0 a0 ) is closely related to the mass flow rate at the source M˙ 0 (see Eq. (16)); therefore, it is suggested that the regime transition is determined by D/Ws regardless of M˙ 0 (i.e., the intensity of an explosive eruption).

CE

PT E

D

MA

This section discusses the mechanisms that determine the regime transition observed in Fig. 5 to clarify how the present two-layer model modifies the run-out distance of the one-layer dilute PDC model of Bursik and Woods (1996) [9]. The transition can be explained by considering the relationship between steady run-out distance and Ws /(u0 a0 ) (i.e., Fig. 6). This figure shows that the curves for the dilute and dense currents are sub-parallel, and that changes in D/Ws shift the curve for the dense current almost in parallel. As a result, the transition between the two regimes occurs at a constant value of D/Ws regardless of Ws /(u0 a0 ). We discuss this feature through the following two steps. First, we derive the analytical solution of the steady run-out distances from the mass conservation equations (i.e., Eqs. (5) and (9)) when neglecting the effect of air entrainment. The analytical solution provides a framework to explain the dependence of the steady run-out distances on Ws /(u0 a0 ) and D/Ws . Second, we compare the analytical solution with numerical results for the two-layer model considering the effect of air entrainment to evaluate the additional effect of air entrainment on the steady run-out distances.

AC

When the effect of air entrainment is neglected (i.e., E = 0), we can derive the analytical solution of the steady run-out distance of the two-layer PDC from the mass conservation equations. The analytical solution for the dilute current was derived by Bursik and Woods (1996) [9]. We extend it to that for the two-layer PDC by adding the solution for the dense current as √ [ ]/ ρ0 − ρg0 Ws r∞ = 1 + 2 ln [No air entrainment], (17) Dilute current: r0 ρa − ρg0 u 0 a0 √ /( ) ρ0 − ρa r∞H D Ws Dense current: = 1+2 [No air entrainment], (18) r0 ρH W s u 0 a0 where r∞ and r∞H are the steady run-out distances of the dilute and dense currents, respectively, and ρg0 is the density of the gas phase in the dilute current at the collapsing column edge (see Appendix E for derivations). Eqs. (17) and (18) explain the feature that the relationships between r∞ /r0 and Ws /(u0 a0 ) and between r∞H /r0 and Ws /(u0 a0 ) are parallel for a given D/Ws when r∞ /r0 ≫ 1 and 16

ACCEPTED MANUSCRIPT

10 1

(a)

10 0 10

Dilute current Dense current (D/Ws = 1×10-3 ) Dense current (D/Ws = 5×10-3 )

RI PT

Flow thickness

Journal of Volcanology and Geothermal Research

-1

10 -2 10 -3 10 -4

SC

10 -5

Dilute current Dense current (D/Ws = 1×10-3 ) Dense current (D/Ws = 5×10-3 )

(b)

NU

1.5 1 0.5

MA

Flow velocity

2

(c)

1

D

0.8 0.6 0.4

PT E

Radial mass flux of solid particles

0 1.2

0.2 0

2

3

4

Distance, r / r0

CE

1

Dilute current Dense current (D/Ws = 1×10-3 ) Dense current (D/Ws = 5×10-3 )

AC

Figure 4: Representative numerical results for steady-state two-layer PDCs in Regime 1 (Run 5 with D/Ws = 1 × 10−3 ; solid curves) and Regime 2 (Run 32 with D/Ws = 5 × 10−3 ; dashed curves). Red and blue curves indicate dilute and dense currents, respectively. Both regimes have identical profiles for the dilute current. (a) Scaled thicknesses of the dilute (h/h0 ) and dense (hH /h0 ) currents. (b) Scaled velocities of the dilute (u/u0 ) and dense (uH /u0 ) currents. (c) Scaled radial mass fluxes of solid particles in the dilute (ns ρuhr/(ρ0 u0 h0 r0 )) and dense (nsH ρH uH hH r/(ρ0 u0 h0 r0 )) currents.

17

ACCEPTED MANUSCRIPT

RI PT

Journal of Volcanology and Geothermal Research

-2

NU

D / Ws

SC

10

-3

MA

10

10-1

Eq. (20) sub-Regime 2a Eq. (19)

Regime 1

100

D

10-2

sub-Regime 2b

(High

PT E

Ws / (u0 a0) . M0

Low)

For a given W s

AC

CE

Figure 5: Flow-regime map of the two-layer PDCs in the parameter space of the two given nondimensional parameters Ws /(u0 a0 ) and D/Ws . Blue symbols indicate Regime 1. Red symbols indicate Regime 2; the circles and triangles represent sub-Regimes 2a and 2b, respectively. The solid line is the approximate analytical solution of the boundary between the two regimes (i.e., Eq. (19)), and the dashed line is the analytical solution of the boundary between the sub-regimes (i.e., Eq. (20)).

18

ACCEPTED MANUSCRIPT Journal of Volcanology and Geothermal Research

Eq. (18) (D / Ws = 1×10 -3 )

Eq. (17)

RI PT

10 1

(D /Ws = 1×10 -3) (D /Ws = 5×10 -3)

NU MA

10 -2

PT E

D

10 0

SC

Scaled steady run-out distance, r∞ /r0 and r∞H / r0

Dilute current r∞ /r0 Dense current r∞H /r0

10 -1

10 0

10 1

Ws / (u0 a0) . M0

(High

Low)

CE

For a given W s

AC

Figure 6: Steady run-out distances of the two-layer PDCs scaled with r0 as a function of the non-dimensional parameter Ws /(u0 a0 ). Symbols represent the numerical results of the two-layer PDC with air entrainment, and the dashed curves are the analytical results of the two-layer PDC without air entrainment (i.e., Eqs. (17) and (18)). Red represents the dilute current (r∞ /r0 ). Blue and cyan represent the dense currents (r∞H /r0 ) with D/Ws = 1 × 10−3 and 5 × 10−3 , respectively.

19

ACCEPTED MANUSCRIPT Journal of Volcanology and Geothermal Research

RI PT

r∞H /r0 ≫ 1 (dashed curves in Fig. 6). They also explain the negative correlations of r∞ /r0 and r∞H /r0 with Ws /(u0 a0 ). As the steady run-out distance of the dilute current, r∞ /r0 , is determined by the point of lift-off where the current density becomes smaller than the ambient density due to particle settling, it depends primarily on the ratio of the timescale of current propagation to that of particle settling (Ws /(u0 a0 )); r∞ /r0 decreases as Ws /(u0 a0 ) increases (i.e., the effect of particle settling becomes more significant; dashed red curve in Fig. 6). Similarly, as the steady run-out distance of the dense current, r∞H /r0 , is determined by the point where the current stops forward propagation due to basal deposition, it depends primarily on the ratio of the timescale of current propagation to that of deposition (D/(u0 a0 )); r∞H /r0 decreases as Ws /(u0 a0 ) increases for a given D/Ws (dashed blue curve in Fig. 6).

D

MA

NU

SC

In the second step, we discuss the additional effect of air entrainment on the steady runout distance. As mentioned in Subsection 2.3, the relative magnitude of the effect of air entrainment to that of particle settling is expressed by the non-dimensional parameter Eu0 /Ws . Although the value of the entrainment coefficient E locally changes (see Eq. (A.1)), our numerical analysis suggests that the profile of E does not vary for a wide range of initial conditions with different Ws /(u0 a0 ); therefore, the horizontally averaged entrainment coefficient ∫r Eave (≡ r2 2−r2 r0∞ (Er)dr) remains almost constant for different Ws /(u0 a0 ). The constant Eave ∞ 0 means that the relative magnitude of the effect of air entrainment (i.e., Eave u0 /Ws ) is inversely proportional to Ws /(u0 a0 ). Accordingly, when Ws /(u0 a0 ) is large (≳ 100 ), Eu0 /Ws is small, so that the steady run-out distance of the dilute current with air entrainment (red symbols in Fig. 6) approaches that without air entrainment (dashed red curve in Fig. 6). On the other hand, when Ws /(u0 a0 ) is small (≲ 100 ), Eu0 /Ws is large and thermal expansion of air entrained into the hot dilute current enhances the lift-off of the current, so that the steady run-out distance of the dilute current with air entrainment becomes much shorter than that without air entrainment.

CE

PT E

Of note, air entrainment also leads to a decrease in the steady run-out distance of the dense current in the two-layer PDC (blue symbols in Fig. 6). This is qualitatively explained by the fact that air entrainment into the dilute current reduces the particle supply from the dilute current to the dense current. Fig. 6 shows that the steady run-out distance of the dense current decreases with increasing Ws /(u0 a0 ) at almost the same rate as that of the dilute current within the range of the present parametric study.

AC

When the curves showing r∞ /r0 and r∞H /r0 as a function of Ws /(u0 a0 ) in Fig. 6 are parallel, the transition condition between Regimes 1 and 2 can be obtained from the analytical solution of the steady run-out distance of the two-layer PDC without air entrainment (i.e., from the relative magnitude of Eqs. (17) and (18)) as follows:  [ ] ρ0 −ρg0 ρ0 −ρa  ≲ / ln (Regime 1) D ρH [ ρa −ρg0 ] (19) Ws ≳ ρ0 −ρa / ln ρ0 −ρg0 (Regime 2). ρH

ρa −ρg0

This criterion explains the regime transition of the numerical results considering air entrainment (black solid line in Fig. 5). The boundary between sub-Regimes 2a and 2b is derived as follows from the mass conservation of the dense current (i.e., Eq. (9)) and the source boundary condition that the dense current is not directly generated from the collapsing column edge: { (sub-Regime 2a) D < ϕs0 /ϕsH (20) Ws ≥ ϕs0 /ϕsH (sub-Regime 2b),

20

ACCEPTED MANUSCRIPT Journal of Volcanology and Geothermal Research

where ϕs0 (≡ ns0 ρ0 /ρs ) is the volume fraction of solid particles in the dilute current at the collapsing column edge (black dashed line in Fig. 5).

Geological implications

SC

5

RI PT

The above results indicate that the presence of the density stratification affects the steady run-out distance of the whole PDC; if the dense current outruns the dilute current, the total distance becomes longer than the prediction by the one-layer dilute PDC model. How dense currents outrun dilute currents primarily depends on a non-dimensional parameter, D/Ws . This dependence on D/Ws reflects the fact that the run-out distance of a two-layer PDC is determined by the mass balance: that is, the global balance between the deposition rate at the base of the dense current (D) and the particle-supply rate from the dilute current to the dense current (Ws ); see Eq. (E.6). It should be noted that the results in Fig. 6 are based on the basic equations where the effect of slope of the ground surface (i.e., topography) is not considered. Nevertheless, the above results are robust in a qualitative sense, and have important geological implications, particularly for the run-out distance of large-scale PDCs.

D

MA

NU

This section discusses the relationship between the steady run-out distance of the two-layer PDC and the magma discharge rate M˙ m (≡ 0.93M˙ 0 ), which can be obtained from Fig. 6 and the definition of Ws /(u0 a0 ) in Eq. (16) (Fig. 7). Fig. 7 shows that the steady run-out distance of the dilute current (red region in the figure) increases with increasing M˙ m . The steady runout distance of the whole two-layer PDC is determined by that of the dilute current when D/Ws ≳ 4 × 10−3 . Otherwise, when D/Ws ≲ 4 × 10−3 the steady run-out distance of the dense basal current (blue region) exceeds that of the dilute current. The steady run-out distance of the dense current increases as D/Ws decreases; it is more than twice as long as that of the dilute current when D/Ws = 1 × 10−3 . We discuss below how the results in Fig. 7 modify previously obtained relationships.

AC

CE

PT E

Some previous relationships between the run-out distance and M˙ m have been proposed on the basis of the subcritical and supercritical steady solutions of the dilute PDC model on a flat surface (Bursik and Woods, 1996 [9]). In particular, the relationship based on the subcritical solution has been commonly used to estimate the magma discharge rates of explosive eruptions (i.e., M˙ m ) from field observations of the distributions of PDC deposits (e.g., Bursik and Woods, 1996 [9]; Carazzo et al., 2008 [11]). The present result for the run-out distances of dilute currents (red region in Fig. 7) coincides with that of the supercritical steady solution (see Fig. D.1). The run-out distance of the supercritical solution is much shorter than that of the subcritical solution (see Fig. D.1) because of the lift-off due to the thermal expansion of entrained air. A key implication of the present results is that the subcritical steady solution is unstable. The absence of the subcritical solution means that previous estimations based on this solution have overestimated the run-out distance of the steady-state dilute PDC. It is suggested that other previous dilute PDC models in which air entrainment was not considered (e.g., Dade and Huppert, 1996 [17]; Dade, 2003 [16]) also overestimate the run-out distance of the dilute PDC. The absence of the subcritical solution implies that the run-out distance of the dilute current cannot quantitatively explain some of the long run-out distances of PDCs observed in the field (e.g., the June 15, 1991 eruption of Pinatubo, and the February 13, 2014 eruption of Kelud; see Fig. 7). The PDC generated during the Pinatubo eruption reached up to ∼ 13 km from the source (Scott et al., 1996 [55]), and the magma discharge rate M˙ m was estimated to be 3 × 108 –1 × 109 kg/s from observed data such as column height and the total amount of the deposit divided by the duration of the eruption (Holasek et al., 1996 [34]; Koyaguchi and Ohno, 21

ACCEPTED MANUSCRIPT Journal of Volcanology and Geothermal Research

D/

Pinatubo

10

D

NU

Kelud

1

×

>~ 4 s W /

-3

0 ×1

Dilute current (Ws =3) (Ws =0.3)

Dense current (Ws =3) (Ws =0.3)

108 109 1010 Magma discharge rate, . . Mm (= 0.93 M0) [kg/s]

1011

PT E

D

107

Ws

=1

-3

10

RI PT

ic

SC

aph r g o Top ffects e

MA

Run-out distance of PDC [km]

100

AC

CE

Figure 7: Run-out distance of PDC as a function of magma discharge rate M˙ m (≡ 0.93M˙ 0 ). Circles represent the numerical results of steady run-out distances of the two-layer PDC model. Red and blue represent the dilute and dense currents, respectively. We set Ws = 0.3 and 3.0 m/s. The steady run-out distance of the whole two-layer PDC is determined by that of the dilute current when D/Ws ≳ 4 × 10−4 and that of the dense current when D/Ws ≲ 4 × 10−3 (here, D/Ws = 1 × 10−3 ). The source conditions other than the aspect ratio a0 are the same as those of Bursik and Woods (1996) [9] (see Table 1). Squares represent estimations for PDCs observed in the field (June 15, 1991 eruption of Pinatubo, and February 13, 2014 eruption of Kelud).

22

ACCEPTED MANUSCRIPT Journal of Volcanology and Geothermal Research

2001 [39]; Wiesner et al., 2004 [72], 2005 [73]; Suzuki and Koyaguchi, 2009 [63]). Similarly, the run-out distance for the Kelud PDC was ∼ 6 km (Maeno et al., in press [47]), and its M˙ m was estimated to be 3 × 107 –1 × 108 kg/s (Maeno et al., in press [47]; Suzuki and Iguchi, in press [62]). These observed run-out distances are much longer than that of the supercritical dilute current (see Fig. 7); this conclusion has been confirmed for a wide range of potential source conditions (ρs = 1000–2500 kg/m3 , ns0 = 0.9–0.97, na0 = 0–0.07, T0 = 800–950 K, and a0 = 0.1–0.2; see Appendix F). To explain the observed long run-out distances, we propose an alternative mechanism: that is, the dense current traveling beyond the lift-off point of the dilute current (i.e., the two-layer PDC of Regime 1).

MA

NU

SC

RI PT

Some geological data, such as those of sedimentary structures, suggest that dense (fluidized) granular flows can form long extents of PDC deposits (e.g., Wilson, 1985 [74]; Cas et al., 2011 [12]; Roche et al., 2016 [52]). Generally, the wide variety of sedimentary structures seen in PDC deposits result from differences in the flow-particle interactions within the bottom boundary layer between dilute and dense currents (Branney and Kokelaar, 2002 [5]). It is considered that deposition from the base of the dilute current produces stratified and/or cross-stratified facies due to tractional sorting and that deposition from the dense current produces massive, poorly sorted facies due to inhibited traction. Many of PDC deposits with long extents are predominantly massive from proximal to distal areas, suggesting that the dense current can transport the mass of particles to the distal area. Incidentally, Roche et al. (2016) [52] estimated that the dense current forming a widespread ignimbrite had low flow speeds ∼ 20 m/s in the distal area; this estimation is consistent with our result with M˙ 0 = 1011 kg/s (Run 7).

CE

PT E

D

It should be emphasized that the dense current in the present two-layer PDC model is different from previous dense-current models, in which the dense currents generated at the collapsing column are driven mainly by their own buoyancy force and/or the momentum obtained at the collapsing column. In such previous models, the majority of particle mass is transported by the dense current from the proximal to distal areas. Our results for the two-layer PDC model show that both dilute and dense currents coexist, and it is the dilute current that transports the majority of the particle mass in the proximal area (see Fig. 4c). Furthermore, Fig. 4b shows that the dense current is driven mainly by the momentum of particles supplied from the overlying dilute current. This means that the dense current can gain its driving force through particle supply from the overlying dilute current even when it does not have a sufficient buoyancy force (associated with its own thickness) and/or a sufficient initial momentum.

AC

Although this study uses the two-layer PDC model to examine the dynamics of currents flowing on a flat surface, the topography can be another important factor to determine the run-out distance of the PDCs; we qualitatively assess possible effects of topography on the run-out distance for large-scale PDCs. Generally, the influence of topography on the dynamics of the dilute and dense currents strongly depends on their flow thicknesses (e.g., Branney and Kokelaar, 2002 [5]; Andrews and Manga, 2011 [1]). Our results for large-scale PDCs show that the thicknesses of the dilute and dense currents are the orders of 102 –103 m and 10−1 –100 m, respectively. In this case, the major topographic effect is that the dilute current can surmount topographic obstacles (with heights of 101 –102 m) and produces stratified veneers on topographic highs, whereas the dense current’s travel is affected by local topography and produces massive valley-confined deposits with a “net-like” distribution (see Fig. 8). Such a topographic effect does not change the fact that the run-out distance of the dense current is determined by the global balance of the deposition rate and the particle-supply rate from the dilute current (i.e., D/Ws ). On the other hand, the topographic effect might reduce the effective value of D/Ws and thus increase the run-out distance of the dense current (Figs. 7 and 8). The particle-supply

23

ACCEPTED MANUSCRIPT Journal of Volcanology and Geothermal Research

rate is considered to be unaffected by topography, whereas the averaged value of D would be strongly reduced because of the low local deposition rate on steep slopes (e.g., Lusso et al., 2017 [46]) and the local erosion process promoted on steep topography (e.g., Mangeney et al., 2010 [48]; Farin et al., 2014 [27]).

SC

RI PT

The effects of topography on the run-out distance for small-scale PDCs are much more complex. There are some previous works that have explored how dense currents can outrun dilute currents on the basis of the two-layer models without the effects of entrainment and thermal expansion of ambient air (Doyle et al., 2008 [20], 2010 [21], 2011 [19]; Kelfoun, 2017 [38]). The results of these studies indicate that the presence of slope of the ground surface directly affects the momentums of both the dilute and dense currents in the case of relatively small-scale PDCs. Our results suggest that the effect of thermal expansion reduces the run-out distances of both the dilute and dense currents; however, a quantitative assessment of this additional effect is not possible at present. We need further studies on the basis of multi-dimensional simulations (e.g., Kelfoun, 2017 [38]) to evaluate quantitatively these topographic effects while comparing with geological data such as distributions and sedimentation patterns related to topography.

AC

(a)

CE

PT E

D

MA

NU

Finally, we briefly discuss the effect of the particle transfer from the dense current to the overlying dilute current, which is ignored in the present model. As mentioned above, a minor dilute ash cloud can be generated by elutriation from the dense current. This mass transfer has two influences in the context of the present analyses. First, it can reduce the effective value of the particle-supply rate from the dilute current to the dense current (i.e., Ws ). This effect increases the run-out distance of the dilute current and decreases that of the dense current. The other is the generation of a secondary dilute ash cloud from a dense current of Regime 1. Such a secondary ash cloud from a Regime 1 current would account for a wedge-shaped front observed in laboratory experiments (Lube et al., 2015 [44]; Breard and Lube, 2017 [6]). The secondary dilute ash cloud can even travel beyond the run-out distance of the parent dense current, which explains the fact that stratified facies are often observed in distal areas. The above two effects may be evaluated by a recent two-layer model of Kelfoun (2017) [38] where the particle transfer from the dense current to the dilute current is taken into account. Quantitative evaluation of these effects, particularly that of the propagation of a secondary dilute ash current, would be important from the viewpoint of hazard mitigation.

(b)

Figure 8: Illustrations of the effects of topography on a large-scale two-layer PDC: (a) flat topography, (b) valley and ridge topography.

24

ACCEPTED MANUSCRIPT Journal of Volcanology and Geothermal Research

6

Concluding remarks

SC

RI PT

This paper has presented a two-layer pyroclastic density current (PDC) model that considers the effects of the entrainment and thermal expansion of ambient air. On the basis of the new twolayer PDC model, we have discussed how the run-out distance of the one-layer dilute PDC model of Bursik and Woods (1996) [9] is modified by the strong density stratification. Our parametric study shows that the flow patterns of the steady-state two-layer PDCs critically depend on the ratio of the deposition speed at the base of the lower dense current (D) to the particle settling speed at the base of the upper dilute current (Ws ). The steady run-out distance of the whole two-layer PDC is determined by that of the dilute current when D/Ws ≳ 4 × 10−3 . Otherwise, when D/Ws ≲ 4 × 10−3 the steady run-out distance of the dense basal current exceeds that of the dilute current. The steady run-out distance of the dilute current is too short to explain some long run-out distances of PDCs observed in the field. It is suggested that the dense current traveling beyond the parent dilute current plays a significant role in the emplacement of PDCs with such long run-out distances.

AC

CE

PT E

D

MA

NU

One of our main conclusions is that the run-out distance of the dense current is determined by the global balance of the basal deposition rate and the particle-supply rate from the overlying dilute current (D/Ws ). Although the present two-layer model deals with PDC dynamics on an idealized flat surface, this conclusion allows us to infer a possible effect of topography on the run-out distance of large-scale PDCs; the run-out distance of the dense current along valleys may increase owing to the small effective values of D/Ws . If we are to verify the topographic effects, the present two-layer PDC model should be applied to a pseudo three-dimensional surface (i.e., a digital elevation model of volcanic topography; cf. Kelfoun, 2017 [38]). For an accurate estimation of the run-out distance of a dense current, increasingly sophisticated models of the deposition speed D would also be necessary. Although the range of D was tentatively estimated here from the hindered settling speed (Appendix A.1), the interaction of D with the flow dynamics (e.g., the effects of basal drag and pore pressure; cf. Iverson and Ouyang, 2015 [37]; Gueugneau et al., 2017 [30]) should be considered. The evaluation of these effects awaits further study.

25

ACCEPTED MANUSCRIPT Journal of Volcanology and Geothermal Research

Notation

PT E

D

MA

NU

SC

RI PT

Aspect ratio of h0 to r0 Ch´ezy drag coefficient Heat capacity at constant pressure (of dilute current), J/(kg·K) Heat capacity of solid particle, J/(kg·K) (Mean) diameter of solid particle, m Deposition speed at base of dense current, m/s Specific internal energy of dilute current, J/kg Entrainment coefficient Imposed frontal Froude number Gravitational body force per unit mass, m/s2 Thickness of (dilute) current, m Characteristic thickness scale of dense current, m Empirical exponent Mass flux, kg/s Mass fraction Thermodynamic pressure, Pa Distance (i.e., radius) from volcanic vent, m Gas constant, J/(kg·K) Particle Reynolds number Richardson number of dilute current Time, s Temperature (of dilute current), K Characteristic time scale, s Velocity component of (dilute) current in r direction, m/s Volume flux of dilute current (≡ 2πuhr), m3 /s Particle settling speed at base of dilute current, m/s Coordinate in vertical direction, m Non-dimensional thickness of artificial bed Dynamic viscosity of gas phase, Pa·s Mass density (of dilute current), kg/m3 Drag at base of current, kg/(m2 ·s) Volume fraction

CE

a0 Cd Cp Cs d D e E F rN g h HH m M˙ n p r R Res Ri t T T u V˙ Ws z ε ηg ρ τ ϕ

AC

Subscript a Air ave Used to emphasize a horizontally averaged variable b Upper surface of ground or deposit (i.e., base of (dense) current) c Upper surface of dense current (i.e., contact surface between dilute and dense currents) f Upper surface of dilute current g Gas phase (i.e., volcanic gas and entrained air) H Dense (i.e., high particle concentration) current m Material erupted from volcanic vent (i.e., magma) N Front (i.e., nose) of dilute current s Solid particle v Volcanic gas 0 Collapsing eruption column edge ∞ Run-out distance of steady-state (dilute) current Superscript ∗ Used to emphasize a scaled variable 26

ACCEPTED MANUSCRIPT Journal of Volcanology and Geothermal Research

Acknowledgements The manuscript was improved by comments from Olivier Roche, Greg A. Valentine, and an anonymous reviewer. Part of this study was supported by funding from the Ministry of Education Science and Culture of Japan (No. 17H02949).

A

Details of the two-layer PDC model

A.1

SC

RI PT

This appendix provides details of the two-layer PDC model. In the mass, momentum, and enthalpy equations of the dilute current (i.e., Eqs. (1), (2), and (4)–(6)), the effects of air entrainment, particle settling, and basal drag are taken into account through E, Ws , and τc , respectively. For the dense current, the effects of basal deposition and basal drag are considered in the mass and momentum equations (i.e., Eqs. (9) and (10)) through D and τb . The formulas for E, Ws , τc , D, and τb used in the basic equations are described below. The numerical procedures for solving the basic equations are also given.

Closure equations

0.075 , (1 + 718Ri2.4 )0.5

(A.1)

MA

E=

NU

The formulas for E, τc , and τb are given as follows. We adopt the entrainment coefficient based on the laboratory experiments of Parker et al. (1987) [50]:

D

where Ri (≡ (ρ − ρa )gh/(ρu2 )) is the local Richardson number. The basal drags of the dilute and dense currents (i.e., τc and τb ) are given by turbulent Ch´ezy-type drags (e.g., Doyle et al., 2011 [19]; Shimizu et al., 2017 [56]): (A.2)

τb = ρH Cdb uH |uH |,

(A.3)

PT E

τc = ρCdc (u − uH )|u − uH |, and

AC

CE

respectively, where Cdc and Cdb are the Ch´ezy drag coefficients. Note that τc becomes ρCdc u|u| when the dense current is absent at the base of the dilute current. The values of Cdc and Cdb (∼ 10−4 ) can be estimated on the basis of the empirical formula proposed by Hager (1988) [31] (cf. Hogg and Pritchard, 2004 [33]). According to experimental studies of an initially fluidized granular flow (Roche et al., 2008 [53]), the fluid-inertial behavior of the dense current is correctly reproduced by the Ch´ezy-type drag. The dense-current model using the Ch´ezy-type drag implies that the dense current is kept in a fluidized state; this assumption is justified by autofluidization due to air escape from the interstices of a rough substrate (Ch´edeville and Roche, 2014 [13]). In the present dense-current model, diverse depositional processes (i.e., progressive aggradation and en masse deposition) are expressed by a combination of the Ch´ezy-type drag and the mass balance between the basal deposition rate and the particle-supply rate from the overlying dilute current; this mass balance is represented by the parameter D/Ws , which is described below. Our results in Section 3 suggest that the behavior of the two-layer PDC strongly depends on the ratio D/Ws . The parameters Ws and D have various uncertainties associated with multiphase physics and polydispersity; the values of these parameters have wide ranges (Ws = 0.3–3.0 m/s and D = 3 × 10−4 –3 × 10−2 m/s). The value of Ws represents the averaged settling speed of particles in the dilute current; as the first approximation, it is estimated from the terminal

27

ACCEPTED MANUSCRIPT Journal of Volcanology and Geothermal Research

velocity of solid particles of diameter d (e.g., Kunii and Levenspiel, 1969 [41]):  (ρ −ρ )gd2 s g   18η  )  ( g 4(ρs −ρg )2 g 2 1/3 d Ws = 225ηg ρg  ( )    3.1(ρs −ρg )gd 1/2

(Res ≤ 6) (6 < Res < 500)

(A.4)

(Res ≥ 500),

ρg

SC

D = (1 − ϕsH )m Ws ,

RI PT

where ρg is the density of the gas phase composed of volcanic gas and entrained air, ηg is the dynamic viscosity of the gas phase, and Res (≡ ρg dWs /ηg ) is the particle Reynolds number. The range of Ws = 0.3–3.0 m/s roughly corresponds to mean particle diameters of d = 10−2 – 100 mm. The deposition speed D is determined by combinations of the deposition and erosion processes (e.g., Dufek et al., 2015 [26]; Iverson and Ouyang, 2015 [37]; Delannay et al., 2017 [18]). Although these processes are not fully understood, the value of D may be estimated from the hindered settling speed (e.g., Cao et al., 2004 [10]; Li et al., 2018 [43]), which is written as (A.5)

D

MA

NU

where ϕsH is the particle volume fraction (0.4–0.5; e.g., Breard et al., 2016 [7]), and m is an empirical exponent (typically 2 to 6 for very fine and well sorted materials with volcanic particles up to 0.25 mm, and 7 to 12 for poorly sorted materials with particles up to 4 mm; Druitt et al., 2007 [23]). Eq. (A.5) with m = 2–6 covers the values given in laboratory experiments for initially fluidized granular flows consisting of particles up to 0.25 mm (i.e., D = 4×10−2 –8×10−2 m/s; Girolami et al., 2010 [29]). On the other hand, the present study applies Eq. (A.5) with m = 7–12 to the estimation of the values of D, because PDC deposits are generally poorly sorted. Thus, we set various values of D/Ws ranging from 10−3 to 10−2 (see Fig. A.1); i.e., D = 3 × 10−4 –3 × 10−2 m/s.

A.2

AC

CE

PT E

We emphasize here that the possible range of D expressed by Eq. (A.5) is tentative one. As mentioned above, the parameter D provides a kinematic condition at the interface between static and mobile particles, and hence, its effective value of D may also depend on the rheology of dense currents and other mechanical processes of particles at the base of currents. For example, some laboratory experiments of granular flows show that, when the flows move on steep slopes, effective values of deposition speed are significantly reduced due to physical processes such as erosion of bed material (e.g., Mangeney et al., 2010 [48]; Farin et al., 2014 [27]; Lusso et al., 2017 [46]). Thus, to explain diverse depositional processes (i.e., progressive aggradation and en masse deposition) in a unified manner, sophisticated models of the basal deposition speed D would be necessary.

Numerical procedure

The present time-dependent two-layer model solves two sets of partial differential equations: the conservation of mass, momentum, and enthalpy in the dilute current (i.e., Eqs. (1), (2), and (4)–(6)), and the conservation of mass and momentum in the dense current (i.e., Eqs. (9) and (10)). Each of these sets can be written as a hyperbolic system with source terms (i.e., the terms on the RHSs of the conservation equations and the additional terms for computation with axisymmetric geometries). As each of these systems is nonlinear, shocks may develop in each current. Accordingly, the hyperbolic part of each system is numerically solved by a finite volume method with shock-capturing capability; the fluxes at the cell faces are calculated by the HLL approximate Riemann solver (e.g., Toro 2001 [67]; LeVeque, 2002 [42]). The source terms are numerically solved by a fractional-step method (e.g., LeVeque, 2002 [42]).

28

ACCEPTED MANUSCRIPT

RI PT

101 100

SC

m=7

-1

NU

10

-2

10

-3

D

10-4 10-5 10-2

m=12

MA

10

10-1

PT E

Ratio of deposition speed to particle settling speed, D/Ws

Journal of Volcanology and Geothermal Research

100

Particle volume fraction, fsH

AC

CE

Figure A.1: The ratio of the deposition speed at the base of the dense current, D, to the particle settling speed at the base of the dilute current, Ws , as a function of the volume fraction of solid particles in the dense current, ϕsH . The value of D is assumed to be the hindered settling speed (i.e., Eq. (A.5)). We set the empirical exponent m = 7–12 (Druitt et al., 2007 [23]). Dashed lines represent ϕsH = 0.4 and 0.5.

29

ACCEPTED MANUSCRIPT Journal of Volcanology and Geothermal Research

B

RI PT

The numerical treatment of the front propagation of the dilute current, where ρ/ρa = 100 – 101 , calculates the front condition (i.e., Eq. (8)) as a boundary condition at each time step (i.e., the Boundary Condition model of Shimizu et al., 2017 [56]). For the dense current, where ρH /ρa = 102 –103 , we can set a thin artificial bed ahead of the front and solve the mass and momentum conservation Eqs. (9) and (10) across the boundary between the flow front and the artificial bed to reproduce approximately the front propagation of the dense current (i.e., the Artificial Bed model of Shimizu et al., 2017 [56]). The thickness of the artificial bed is given by r0 s0 ρ0 εHH , where ε is the non-dimensional thickness of the artificial bed, and HH (≡ nnsH ρH Ws u0 ) is the characteristic thickness scale of the dense current. We set ε = 10−10 (Shimizu et al., 2017 [56]).

Enthalpy conservation equation (6) and total energy conservation equation

NU

SC

This appendix gives the derivation of the depth-averaged enthalpy conservation equation (6) from the total energy conservation equation. It is also shown that the total energy conservation equation of Bursik and Woods (1996) [9] is an approximation of our total energy conservation equation.

PT E

D

MA

Integrating the total energy conservation for the Euler equations (cf. Toro, 2009 [68]) from z = zc to zf (cf. Ungarish, 2009 [69]), we obtain the depth-averaged total energy conservation equation: ( ) ) ( ) ( ( ) h h ∂ 1 2 1 ∂ 1 3 + zc gh + + zc ughr + puhr ρeh + ρu h + ρ ρeuhr + ρu hr + ρ ∂t 2 2 r ∂r 2 2 ) ( ∂h ∂zc ns 1 ns 1 − ρgh = ρa E|u|gzf − ρWs gzc − ρu2 Ws − τc u + p − ρgh 2 ∂t ∂t nsH 2 nsH ( ) 1 − nsH +ρa E|u|Cpa Ta − ns ρWs Cs T + ns ρWs Cpv T , (B.1) nsH

AC

CE

where e is the specific internal energy of the dilute current. The first and second terms on the RHS of Eq. (B.1) represent the potential energies added and lost by air entrainment and particle settling, respectively. The third and fourth terms represent the kinetic energies lost by particle settling and basal drag, respectively. The fifth and sixth terms represent the thermal energies added and lost by air entrainment and particle settling, respectively. Substituting the depth-averaged mass and momentum equations (i.e., Eqs. (1) and (2)) and a caloric equation of state (i.e., e + p/ρ = Cp T ) into Eq. (B.1) yields the depth-averaged enthalpy conservation equation (6). The present total energy conservation (i.e., Eq. (B.1)) and hence the enthalpy conservation (i.e., Eq. (6)) include all the terms for mechanical energy. However, some of these terms are neglected in the energy equation of Bursik and Woods (1996) [9]. When the mass flow rate M˙ 0 is smaller than 1012 kg/s, the amount of change in thermal energy by conversion from mechanical energy (i.e., the fourth to sixth terms on the RHS of Eq. (6)) is negligibly small, so that the energy equation of Bursik and Woods (1996) [9] is approximately equivalent to Eqs. (6) and (B.1).

30

ACCEPTED MANUSCRIPT Journal of Volcanology and Geothermal Research

C

Non-dimensionalization

This appendix gives the non-dimensionalization of the basic equations for the present two-layer PDC model by the quantities at the collapsing column edge (i.e., the source boundary). We use r0 as the horizontal length scale, h0 as the thickness and height scales, u0 as the velocity scale, and T (≡ r0 /u0 ) as the time scale. We also use ρ0 as the reference density, T0 as the reference temperature, and Cp0 as the reference specific heat at constant pressure. Accordingly, we can switch all dimensional variables to scaled variables (denoted below by an asterisk) as follows: r = r0 r∗ , z = h0 z ∗ , h = h0 h∗ , hH = h0 h∗H , u = u0 u∗ , uH = u0 u∗H , t = T t∗ ,

RI PT

ρ = ρ0 ρ∗ , T = T0 T ∗ , Cp = Cp0 Cp∗ .

(C.1)

The mass fractions of entrained air and solid particles are also scaled by na0 and ns0 , respectively (cf. Ungarish, 2009 [69]): (C.2)

SC

na = na0 n∗a , and ns = ns0 n∗s .

NU

Applying Eqs. (C.1) and (C.2) to Eqs. (1), (2), (4)–(6), and (8) yields the non-dimensional basic equations of the dilute current:

(C.3)

(C.4)

(C.5)

(C.6)

AC

CE

PT E

D

MA

Conservation of entrained air mass: ∂ 1 ∂ ρa /ρ0 E ∗ |u |, (n∗ ρ∗ h∗ ) + ∗ ∗ (n∗a ρ∗ u∗ h∗ r∗ ) = ∂t∗ a r ∂r na0 a0 Conservation of solid particle mass: 1 ∂ Ws /u0 ∗ ∗ ∂ (n∗s ρ∗ h∗ ) + ∗ ∗ (n∗s ρ∗ u∗ h∗ r∗ ) = − ns ρ , ∗ ∂t r ∂r a0 Conservation of bulk mass: ( ) ∂ ∗ ∗ 1 ∂ ns0 Ws /u0 ∗ ∗ ρa /ρ0 Eu0 |u∗ | ∗ ∗ ∗ ∗ (ρ h ) + ∗ ∗ (ρ u h r ) = ns ρ −1 , ∂t∗ r ∂r nsH a0 ns0 /nsH Ws n∗s ρ∗ Momentum conservation: ( ∗ ) ∂ ∂ ∗ ∗ ∗ 1 ∂ ( ∗ ∗2 ∗ ∗ ) 1 ρ − ρa /ρ0 ∗2 (ρ u h ) + ∗ ∗ ρ u h r + Ri0 ∗ h ∂t∗ r ∂r ∂r 2 1 − ρa /ρ0 ρ∗ − ρa /ρ0 ∗ ∂zc∗ ns0 Ws /u0 ∗ ∗ ∗ Cdc ∗ ∗ = −Ri0 h − ns ρ u − ρ (u − u∗H )|u∗ − u∗H |, ∗ 1 − ρa /ρ0 ∂r nsH a0 a0 Enthalpy conservation: ) ∂ ( ∗ ∗ ∗ ∗) 1 ∂ ( ρ Cp T h + ∗ ∗ ρ∗ Cp∗ T ∗ u∗ h∗ r∗ ∗ ∂t r ∂r ρa Cpa Ta E ∗ Cs Ws /u0 ∗ ∗ ∗ 1 − nsH Cpv Ws /u0 ∗ ∗ ∗ = |u | − ns0 ns ρ T − ns ρ T ρ0 Cp0 T0 a0 Cp0 a0 nsH /ns0 Cp0 a0 (( ( ) ) ) ρa gh0 /2 ∗ ρa E ∗ 1 u∗2 ns0 Ws /u0 ∗ ∗ + h |u | 1− + 1 + n ρ , s Cp0 T0 ρ0 a 0 ρ0 Ri0 h∗ nsH a0 Front condition: √ √ ∗ ∗ √ drN ρ0 ρN − ρa /ρ0 ∗ ∗ ∗ = u∗N = F rN Ri0 h at r∗ = rN (t ). ∗ dt ρa 1 − ρa /ρ0 N

31

(C.7)

(C.8)

ACCEPTED MANUSCRIPT Journal of Volcanology and Geothermal Research

Similarly, applying Eqs. (C.1) and (C.2) to Eqs. (9) and (10) yields the non-dimensional basic equations of the dense current:

SC

RI PT

Mass conservation: ( /( )) ∂h∗H D 1 ∂ ns0 ρ0 Ws /u0 ∗ ∗ ns0 ρ0 ∗ ∗ ∗ ∗ ∗ + ∗ ∗ (uH hH r ) = ns ρ 1 − n ρ , (C.9) ∂t∗ r ∂r nsH ρH a0 Ws nsH ρH s Momentum conservation: ( ) ∂ 1 ∂ ( ∗2 ∗ ∗ ) ∂ 1 1 − ρa /ρH ∗2 ∗ ∗ h (u h ) + ∗ ∗ uH hH r + Ri0 ∗ ∂t∗ H H r ∂r ∂r 2 1 − ρa /ρ0 H ( ∗ ) 1 − ρa /ρH ∗ ∂zb∗ ρ0 ∗ ∂ ρ − ρa /ρ0 ∗ = −Ri0 h − Ri0 hH ∗ h 1 − ρa /ρ0 H ∂r∗ ρH ∂r 1 − ρa /ρ0 ( /( )) ns0 ρ0 Ws /u0 ∗ ∗ ∗ D ns0 ρ0 ∗ ∗ u∗ + ns ρ u 1 − n ρ nsH ρH a0 Ws nsH ρH s u∗H Cdb ∗ ∗ ρ0 Cdc ∗ ∗ ρ (u − u∗H )|u∗ − u∗H | − u |u |. (C.10) + ρH a 0 a0 H H

MA

NU

The non-dimensional parameters in Eqs. (C.3)–(C.10) are classified into three groups: those determined from material properties (i.e., na0 , ns0 , nsH , ρa /ρ0 , ρH /ρ0 , Cpa Ta /(Cp0 T0 ), Cs /Cp0 , and Cpv /Cp0 ), those related to drags (i.e., Cdc , Cdb , and F rN ), and the others (i.e., a0 , Ri0 , Ws /(u0 a0 ), Eu0 /Ws , D/Ws , and gh0 /(2Cp0 T0 )). In this paper, we fix the non-dimensional parameters of the first and second groups and some of the parameters related to source conditions (a0 and Ri0 ) in the third group for the following reasons.

CE

PT E

D

Given that we are concerned with the PDCs generated by total column collapse during a magmatic eruption, the material properties of the dilute current at the collapsing column edge are set as typical values used in Bursik and Woods (1996) [9] (i.e., na0 = 0.07, ns0 = 0.9, and T0 = 950 K). Accordingly, all the non-dimensional parameters determined by material properties are fixed. We also fix all the non-dimensional parameters related to drags (i.e., Cdc , Cdb , and F rN ) as typical values on the basis of laboratory experiments: Cdc and Cdb are set as 10−4 by the empirical formula proposed by Hager (1988) [31] (cf. Hogg and Pritchard, 2004 [33]), and F rN is set to 1.19 (Huppert and Simpson, 1980 [36]).

AC

In the third group, a0 , Ri0 , and Ws /(u0 a0 ) are related to the source conditions. Among these three parameters, the variation in Ws /(u0 a0 ) directly reflects the wide variation in M˙ 0 (see Eq. (16)). The values of a0 and Ri0 also depend on the eruption column dynamics (cf. Suzuki and Koyaguchi, 2012 [64]; Hunt and Burridge, 2015 [35]; Koyaguchi et al., 2018 [40]); however, their physical processes are not fully understood. In this study, we tentatively set the aspect ratio a0 = 0.2 on the basis of three-dimensional numerical simulations (e.g., Suzuki et al., 2005 [65]; Costa et al., 2018 [15]; Koyaguchi et al., 2018 [40]). We also assume that the Richardson number at the source Ri0 is on the order of 100 when a PDC near the source area is driven by gravitational collapse of the eruption column (cf. Bursik and Woods, 1996 [9]; Slim and Huppert, 2011 [57]). Although most of the results in this paper are for Ri0 = 1, we also investigated the effects of variations in Ri0 from the viewpoint of the stability of the steady subcritical solution for the dilute current (see Appendix D for details). Further experimental and numerical studies are necessary to understand fully the relationships between the eruption column dynamics and the source conditions of PDCs.

32

ACCEPTED MANUSCRIPT Journal of Volcanology and Geothermal Research

D

Stability of the steady subcritical solution for dilute PDCs

Bursik and Woods (1996) [9] proposed on the basis of their steady-state model that the dilute current is supercritical (i.e., the local Richardson number Ri < 1 throughout the flow) when the source Richardson number Ri0 < 1, whereas it remains subcritical (i.e., Ri > 1) for Ri0 > 1 (black and gray dashed curves in Fig. D.1, respectively). This appendix shows that the subcritical solution of Bursik and Woods (1996) [9] is unstable for any value of Ri0 .

SC

RI PT

Given a steady source condition, the present time-dependent PDC model shows that the current eventually converges to a steady state. The steady-state profile of the present result agrees with the steady solution of Bursik and Woods (1996) [9] (i.e., the supercritical solution) for Ri0 < 1 (see the results for Ri0 = 0.99 in Fig. D.1). For Ri0 > 1, on the other hand, the present model shows that the initially subcritical flow undergoes a transition to a supercritical flow near the source; consequently, the present steady-state profile for Ri0 > 1 also approaches the supercritical solution of Bursik and Woods (1996) [9] (see the results for Ri0 = 1.01 in Fig. D.1). We confirmed that even when the steady subcritical solution of Bursik and Woods (1996) [9] is given as the initial condition of the present model, the subcritical solution disappears and eventually approaches the supercritical solution (see Supplementary Movie S3).

Analytical solution of steady run-out distance of two-layer PDC without air entrainment

PT E

E

D

MA

NU

The discrepancy regarding the solution for Ri0 > 1 arises from the fact that there is no physical mechanism sustaining the boundary condition (i.e., thickness and velocity) at the front of the subcritical solution of Bursik and Woods (1996) [9]. A small perturbation can cause a collapse of the flow front to form a radially spreading supercritical flow. The perturbation propagates backward; consequently, the whole current except for near the source becomes supercritical. This means that the steady subcritical solution of Bursik and Woods (1996) [9] is unstable.

CE

This appendix derives the analytical solution of the steady run-out distance of the two-layer PDC from the mass conservation equations while neglecting the effect of air entrainment. The analytical solution for the dilute current derived by Bursik and Woods (1996) [9] is extended here to that for the two-layer PDC by adding the solution for the dense current.

AC

Bursik and Woods (1996) [9] derived the analytical solution for the dilute current from the steady-state form of the mass conservation of solid particles in the dilute current (i.e., Eq. (5)): d (2πrns ρuh) = −2πrns ρWs . dr

(E.1)

When air entrainment is negligible, the radial volume flux of the current, V˙ (≡ 2πuhr), remains approximately constant along the length of the current (i.e., V˙ ≈ V˙ 0 ≡ 2πu0 h0 r0 ). Accordingly, Eq. (E.1) is reduced to d 2πWs r. (ns ρ) = −ns ρ dr V˙ 0

(E.2)

Eq. (E.2) is integrated from the collapsing column edge r = r0 to the steady run-out distance r∞ to yield [ ] ns0 ρ0 πWs 2 ln = (r∞ − r02 ), (E.3) ns∞ ρ∞ V˙ 0 33

ACCEPTED MANUSCRIPT Journal of Volcanology and Geothermal Research

(a)

This study [Ri 0 = 0.99 and 1.01]

10

Supercritical solution [Ri 0 = 0.99] Subcritical solution [Ri 0 = 1.01]

SC

10 0

-1

10

2

(b)

MA

10

3

D

10 1 10

PT E

10 0 -1

10 -3 1

CE

10 -2

AC

Richardson number, Ri

10 10 4

RI PT

Bursik and Woods (1996) 1

NU

Flow thickness, h / h0

10 2

2 Distance, r / r0

3

Figure D.1: Steady-state profiles of (a) the scaled thickness, h/h0 , and (b) the Richardson number, Ri, of the dilute current. Solid red curves represent the numerical results of the present time-dependent model using Ri0 = 0.99 and 1.01. (The profiles are identical.) The dashed black and gray curves represent the numerical results of the steady model of Bursik and Woods (1996) [9] using Ri0 = 0.99 and 1.01, respectively. The parameters other than Ri0 in these runs are the same as those in Run 45.

34

ACCEPTED MANUSCRIPT Journal of Volcanology and Geothermal Research

where the subscript ∞ denotes values at the steady run-out distance of the dilute current. Because ns ρ = ρs (ρ − ρg )/(ρs − ρg ) (i.e., Eq. (3)), ρg ≈ ρg0 , and ρ∞ = ρa , the steady run-out distance of the dilute current can be obtained by √ [ ] ρ0 − ρg0 M˙ 0 2 r∞ = r0 + ln . (E.4) πρ0 Ws ρa − ρg0

RI PT

We derive the analytical solution for the dense current from the steady-state form of its mass conservation (i.e., Eq. (9)): d ns (2πrρH uH hH ) + 2πrρH D = 2πr ρWs . dr nsH

SC

Eq. (E.5) is integrated from the collapsing column edge r = r0 to r∞H to yield ) ∫ r∞H ( ns 2 2 πρH D(r∞H − r0 ) = − ρWs dr, 2πr nsH r0

(E.5)

(E.6)

D

MA

NU

where r∞H is the steady run-out distance of the dense current. Eq. (E.6) means that the steady run-out distance of the dense current, r∞H , is determined by the global balance of the particlesupply rate from the dilute current (i.e., the RHS of Eq. (E.6)) and the deposition rate (i.e., the LHS of Eq. (E.6)) regardless of the presence or absence of air entrainment effect. When air entrainment into the overlying dilute current is negligible, the particle-supply rate on the RHS of Eq. (E.6) is given by the change in the radial mass flux of the overlying dilute current (see Eq. (1)), so that we obtain ∫ r∞H d 2 2 (E.7) (2πrρuh) dr. πρH D(r∞H − r0 ) = − dr r0

PT E

The integral of Eq. (E.7) can be obtained only numerically when r∞H < r∞ (i.e., in Regime 2). Otherwise, when r∞H ≥ r∞ (i.e., in Regime 1), Eq. (E.7) is rewritten as 2 πρH D(r∞H − r02 ) = M˙ 0 − M˙ ∞ ,

(E.8)

AC

CE

where M˙ ∞ is the radial mass flux of the dilute current at r∞ . As mentioned above, the radial volume flux of the dilute current remains approximately constant (i.e., V˙ ≈ V˙ 0 ), and the density of the dilute current at r∞ is the ambient density (i.e., ρ∞ = ρa ), so that M˙ ∞ ≈ ρa V˙ 0 = M˙ 0 ρa /ρ0 . Substituting M˙ ∞ = M˙ 0 ρa /ρ0 into Eq. (E.8) yields the steady run-out distance of the dense current: √ M˙ 0 ρ0 − ρa r∞H = r02 + . (E.9) πρH D ρ0

F

Sensitivity of supercritical dilute PDC to source conditions

In this appendix we investigate the sensitivity of our model to source conditions (i.e., the density of solid particles ρs , the mass fraction of solid particles ns0 , the mass fraction of entrained air na0 , the temperature T0 , and the aspect ratio a0 ; see Table F.1) to confirm the conclusion derived from the results in Fig. 7. Fig. F.1 shows that the steady run-out distance of the supercritical dilute current is substantially shorter than the run-out distances observed in the field (the June 15, 1991 eruption of Pinatubo, and the February 13, 2014 eruption of Kelud) for a wide range of ρs (1000–2500 kg/m3 ), ns0 (0.9–0.97), na0 (0–0.07), T0 (800–950 K), and a0 (0.1–0.2). 35

ACCEPTED MANUSCRIPT

10

RI PT

Pinatubo

SC

Kelud

Reference (see Table 1)

a0 = 0.1

108 109 Magma discharge rate, . . Mm (= (1- na0 ) M0 ) [kg/s]

MA

107

NU

1

rs = 2500 kg/m3 ns0 = 0.97; n = 0 a0 T0 = 800 K

1010

PT E

D

Steady run-out distance of supercritical dilute PDC [km]

Journal of Volcanology and Geothermal Research

AC

CE

Figure F.1: Sensitivity of steady run-out distance of two-layer PDC of Regime 2 (i.e., supercritical dilute current) to source conditions (i.e., magma discharge rate M˙ m (≡ (1 − na0 )M˙ 0 ), density of solid particles ρs , mass fraction of solid particles ns0 , mass fraction of entrained air na0 , temperature T0 , and aspect ratio a0 (≡ h0 /r0 )). Black dashed line (Runs 49–52) represents numerical results in the main text (ρs = 1000 kg/m3 , ns0 = 0.9, na0 = 0.07, T0 = 950 K, and a0 = 0.2; see Tables 1 and 2 for the other parameters). Pink solid line (Runs 55–58) represents the results with the same parameters as for Runs 49–52 except that ρs = 2500 kg/m3 . Green solid line (Runs 59–62) represents the results with the same parameters as for Runs 49–52 except that ns0 = 0.97 and na0 = 0. Blue solid line (Runs 63–66) represents the results with the same parameters as for Runs 49–52 except that T0 = 800 K. Purple solid line (Runs 67–70) represents the results with the same parameters as for Runs 49–52 except that a0 = 0.1. Note that the purple solid line almost overlaps the blue solid line. A summary of the parameters used in this sensitivity analysis is listed in Table F.1. Squares represent estimations for PDCs observed in the field (June 15, 1991 eruption of Pinatubo, and February 13, 2014 eruption of Kelud).

36

ACCEPTED MANUSCRIPT Journal of Volcanology and Geothermal Research

Table F.1: Input parameters for sensitivity analysis to source conditions.

D PT E CE AC

37

na0 0.07 0.07 0.07 0.07 0.00 0.00 0.00 0.00 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07

T0 [K] 950 950 950 950 950 950 950 950 800 800 800 800 950 950 950 950

RI PT

ns0 0.9 0.9 0.9 0.9 0.97 0.97 0.97 0.97 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9

SC

ρs [kg/m3 ] 2500 2500 2500 2500 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000

NU

M˙ 0 [kg/s] 107 108 109 1010 107 108 109 1010 107 108 109 1010 107 108 109 1010

MA

Simulation Run 55 Run 56 Run 57 Run 58 Run 59 Run 60 Run 61 Run 62 Run 63 Run 64 Run 65 Run 66 Run 67 Run 68 Run 69 Run 70

a0 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.1 0.1 0.1 0.1

ACCEPTED MANUSCRIPT Journal of Volcanology and Geothermal Research

References [1] Andrews, B.J., Manga, M., 2011. Effects of topography on pyroclastic density current runout and formation of coignimbrites. Geology 39 (12), 1099–1102. https://doi.org/10.1130/G32226.1. [2] Andrews, B.J., Manga, M., 2012. Experimental study of turbulence, sedimentation, and coignimbrite mass partitioning in dilute pyroclastic density currents. J. Volcanol. Geotherm. Res. 225–226, 30–44. https://doi.org/10.1016/j.jvolgeores.2012.02.011.

RI PT

[3] Bonnecaze, R.T., Huppert, H.E., Lister, J.R., 1993. Particle-driven gravity currents. J. Fluid Mech. 250, 339–369. https://doi.org/10.1017/S002211209300148X.

SC

[4] Branney, M.J., Kokelaar, P., 1992. A reappraisal of ignimbrite emplacement: Progressive aggradation and changes from particulate to non-particulate flow during emplacement of high-grade ignimbrite. Bull. Volcanol. 54 (6), 504–520. https://doi.org/10.1007/BF00301396. [5] Branney, M.J., Kokelaar, P., 2002. Pyroclastic density currents and the sedimentation of ignimbrites. Memoirs of the Geological Society of London, 27.

MA

NU

[6] Breard, E.C.P., Lube, G., 2017. Inside pyroclastic density currents – uncovering the enigmatic flow structure and transport behaviour in large-scale experiments. Earth Planet. Sci. Lett. 458, 22–36. https://doi.org/10.1016/j.epsl.2016.10.016. [7] Breard, E.C.P., Lube, G., Jones, J.R., Dufek, J., Cronin, S.J., Valentine, G.A., Moebis, A., 2016. Coupling of turbulent and non-turbulent flow regimes within pyroclastic density currents. Nat. Geosci. 9, 767–771. https://doi.org/10.1038/NGEO2794.

PT E

D

[8] Burgisser, A., Bergantz, G.W., 2002. Reconciling pyroclastic flow and surge: the multiphase physics of pyroclastic density currents. Earth Planet. Sci. Lett. 202 (2), 405–418. https://doi.org/10.1016/S0012-821X(02)00789-6. [9] Bursik, M.I., Woods, A.W., 1996. The dynamics and thermodynamics of large ash flows. Bull. Volcanol. 58 (2–3), 175–193. https://doi.org/10.1007/s004450050134.

CE

[10] Cao, Z., Pender, G., Wallis, S., Carling, P., 2004. Computational dam-break hydraulics over erodible sediment bed. J. Hydraul. Eng. 130 (7), 689–703. https://doi.org/10.1061/(ASCE)0733-9429(2004)130:7(689).

AC

[11] Carazzo, G., Kaminski, E., Tait, S., 2008. On the dynamics of volcanic columns: A comparison of field data with a new model of negatively buoyant jets. J. Volcanol. Geotherm. Res. 178 (1), 94–103. https://doi.org/10.1016/j.jvolgeores.2008.01.002. [12] Cas, R.A.F., Wright, H.M.N., Folkes, C.B., Lesti, C., Porreca, M., Giordano, G., Viramonte, J.G., 2011. The flow dynamics of an extremely large volume pyroclastic flow, the 2.08-Ma Cerro Gal´an Ignimbrite, NW Argentina, and comparison with other flow types. Bull. Volcanol. 73 (10), 1583–1609. https://doi.org/10.1007/s00445-011-0564-y. [13] Ch´edeville, C., Roche, O., 2014. Autofluidization of pyroclastic flows propagating on rough substrates as shown by laboratory experiments. J. Geophys. Res. 119 (3), 1764–1776. https://doi.org/10.1002/2013JB010554. [14] Cole, P.D., Neri, A., Baxter, P.J., 2015. Hazards from pyroclastic density currents. In: Sigurdsson, H. (Ed.), The Encyclopedia of Volcanoes (Second Edition), Academic Press, Amsterdam, 943–956. https://doi.org/10.1016/B978-0-12-385938-9.00054-7. 38

ACCEPTED MANUSCRIPT Journal of Volcanology and Geothermal Research

[15] Costa, A., Suzuki, Y.J., Koyaguchi, T., 2018. Understanding the plume dynamics of explosive super-eruptions. Nat. Comm. 9, 654. https://doi.org/10.1038/s41467-018-02901-0. [16] Dade, W.B., 2003. The emplacement of low-aspect ratio ignimbrites by turbulent parent flows. J. Geophys. Res. 108 (B4), 2211. https://doi.org/10.1029/2001JB001010. [17] Dade, W.B., Huppert, H.E., 1996. Emplacement of the Taupo ignimbrite by a dilute turbulent flow. Nature 381, 509–512. https://doi.org/10.1038/381509a0.

RI PT

[18] Delannay, R., Valance, A., Mangeney, A., Roche, O., Richard, P., 2017. Granular and particle-laden flows: from laboratory experiments to field observations. J. Phys. D: Appl. Phys. 50 (5), 053001. https://doi.org/10.1088/1361-6463/50/5/053001. [19] Doyle, E.E., Hogg, A.J., Mader, H.M., 2011. A two-layer approach to modelling the transformation of dilute pyroclastic currents into dense pyroclastic flows. Proc. R. Soc. A 467, 1348–1371. https://doi.org/10.1098/rspa.2010.0402.

SC

[20] Doyle, E.E., Hogg, A.J., Mader, H.M., Sparks, R.S.J., 2008. Modeling dense pyroclastic basal flows from collapsing columns. Geophys. Res. Lett. 35, L04305. https://doi.org/10.1029/2007GL032585.

NU

[21] Doyle, E.E., Hogg, A.J., Mader, H.M., Sparks, R.S.J., 2010. A two-layer model for the evolution and propagation of dense and dilute regions of pyroclastic currents. J. Volcanol. Geotherm. Res. 190 (3–4), 365–378. https://doi.org/10.1016/j.jvolgeores.2009.12.004.

MA

[22] Druitt, T.H., 1998. Pyroclastic density currents. In: Gilbert, J.S., Sparks, R.S.J. (Eds.), The Physics of Explosive Volcanic Eruptions. Geol. Soc. London Spec. Pub. 145, 145–182. https://doi.org/10.1144/GSL.SP.1996.145.01.08.

PT E

D

[23] Druitt, T.H., Avard, G., Bruni, G., Lettieri, P., Maez, F., 2007. Gas retention in finegrained pyroclastic flow materials at high temperatures. Bull. Volcanol. 69 (8), 881–901. https://doi.org/10.1007/s00445-007-0116-7. [24] Dufek, J., 2016. The Fluid Mechanics of Pyroclastic Density Currents. Annu. Rev. Fluid Mech. 48, 459–485. https://doi.org/10.1146/annurev-fluid-122414-034252.

CE

[25] Dufek, J., Bergantz, G.W., 2007. Dynamics and deposits generated by the Kos Plateau Tuff eruption: Controls of basal particle loss on pyroclastic flow transport. Geochem. Geophys. Geosyst. 8 (12), Q12007. https://doi.org/10.1029/2007GC001741.

AC

[26] Dufek, J., Esposti Ongaro, T., Roche, O., 2015. Pyroclastic Density Currents: Processes and Models. In: Sigurdsson, H. (Ed.), The Encyclopedia of Volcanoes (Second edition), Academic Press, Amsterdam, 617–629. https://doi.org/10.1016/B978-0-12-3859389.00035-3. [27] Farin, M., Mangeney, A., Roche, O., 2014. Fundamental changes of granular flow dynamics, deposition, and erosion processes at high slope angles: Insights from laboratory experiments. J. Geophys. Res. 119 (3), 504–532. https://doi.org/10.1002/2013JF002750. [28] Fisher, R.V., 1979. Models for pyroclastic surges and pyroclastic flows. J. Volcanol. Geotherm. Res. 6 (3–4), 305–318. https://doi.org/10.1016/0377-0273(79)90008-8. [29] Girolami, L., Roche, O., Druitt, T.H., Corpetti, T., 2010. Particle velocity fields and depositional processes in laboratory ash flows, with implications for the sedimentation of dense pyroclastic flows. Bull. Volcanol. 72 (6), 747–759. https://doi.org/10.1007/s00445010-0356-9. 39

ACCEPTED MANUSCRIPT Journal of Volcanology and Geothermal Research

[30] Gueugneau, V., Kelfoun, K., Roche, O., Chupin, L., 2017. Effects of pore pressure in pyroclastic flows: Numerical simulation and experimental validation. Geophys. Res. Lett. 44 (5), 2194–2202. https://doi.org/10.1002/2017GL072591. [31] Hager, W.H., 1988. Abflussformeln f¨ ur turbulente Str¨omungen. Wasserwirtschaft 78, 79– 84.

RI PT

[32] Henry, C.D., Hinz, N.H., Faulds, J.E., Colgan, J.P., John, D.A., Brooks, E.R., Cassel, E.J., Garside, L.J., Davis, D.A., Castor, S.B., 2012. Eocene-Early Miocene paleotopography of the Sierra Nevada-Great Basin-Nevadaplano based on widespread ash-flow tuffs and paleovalleys. Geosphere 8 (1), 1–27. https://doi.org/10.1130/GES00727.1. [33] Hogg, A.J., Pritchard, D., 2004. The effects of hydraulic resistance on dam-break and other shallow inertial flows. J. Fluid Mech. 501, 179–212. https://doi.org/10.1017/S0022112003007468.

SC

[34] Holasek, R.E., Self, S., Woods, A.W., 1996. Satellite observations and interpretation of the 1991 Mount Pinatubo eruption plumes. J. Geophys. Res. 101 (B12), 27635–27655. https://doi.org/10.1029/96JB01179.

NU

[35] Hunt, G.R., Burridge, H.C., 2015. Fountains in industry and nature. Annu. Rev. Fluid Mech. 47, 195–220. https://doi.org/10.1146/annurev-fluid-010313-141311.

MA

[36] Huppert, H.E., Simpson, J.E., 1980. The slumping of gravity currents. J. Fluid Mech. 99 (4), 785–799. https://doi.org/10.1017/S0022112080000894.

D

[37] Iverson, R.M., Ouyang, C., 2015. Entrainment of bed material by Earth-surface mass flows: Review and reformulation of depth-integrated theory. Rev. Geophys. 53 (1), 27–58. https://doi.org/10.1002/2013RG000447.

PT E

[38] Kelfoun, K., 2017. A two-layer depth-averaged model for both the dilute and the concentrated parts of pyroclastic currents. J. Geophys. Res. 122 (6), 4293–4311. https://doi.org/10.1002/2017JB014013.

CE

[39] Koyaguchi, T., Ohno, M., 2001. Reconstruction of eruption column dynamics on the basis of grain size of tephra fall deposits: 2. Application to the Pinatubo 1991 eruption. J. Geophys. Res. 106 (B4), 6513–6533. https://doi.org/10.1029/2000JB900427.

AC

[40] Koyaguchi, T., Suzuki, Y.J., Takeda, K., Inagawa, S., 2018. The condition of eruption column collapse: 2. Three-dimensional numerical simulations of eruption column dynamics. J. Geophys. Res. 123. https://doi.org/10.1029/2017JB015259. [41] Kunii, D., Levenspiel, O., 1969. Fluidization Engineering, John Wiley, New York. [42] LeVeque, R.J., 2002. Finite Volume Methods for Hyperbolic Problems. Cambridge University Press, Cambridge. [43] Li, J., Cao, Z., Hu, K., Pender, G., Liu, Q., 2018. A depth-averaged two-phase model for debris flows over erodible beds. Earth Surf. Process. Landforms 43 (4), 817–839. https://doi.org/10.1002/esp.4283. [44] Lube, G., Breard, E.C.P., Cronin, S.J., Jones, J., 2015. Synthesizing large-scale pyroclastic flows: Experimental design, scaling, and first results from PELE. J. Geophys. Res. 120 (3), 1487–1502. https://doi.org/10.1002/2014JB011666.

40

ACCEPTED MANUSCRIPT Journal of Volcanology and Geothermal Research

[45] Lube, G., Huppert, H.E., Sparks, R.S.J., Freundt, A., 2007. Static and flowing regions in granular collapses down channels. Phys. Fluids 19, 043301. https://doi.org/10.1063/1.2712431. [46] Lusso, C., Bouchut, F., Ern, A., Mangeney, A., 2017. A free interface model for static/flowing dynamics in thin-layer flows of granular materials with yield: Simple shear simulations and comparison with experiments. Appl. Sci. 7 (4), 386. https://doi.org/10.3390/app7040386.

RI PT

[47] Maeno, F., Nakada, S., Yoshimoto, M., Shimano, T., Hokanishi, N., Zaennudin, A., Iguchi, M., in press. A sequence of a plinian eruption preceded by dome destruction at Kelud volcano, Indonesia, on February 13, 2014, revealed from tephra fallout and pyroclastic density current deposits. J. Volcanol. Geotherm. Res. https://doi.org/10.1016/j.jvolgeores.2017.03.002.

SC

[48] Mangeney, A., Roche, O., Hungr, O., Mangold, N., Faccanoni, G., Lucas, A., 2010. Erosion and mobility in granular collapse over sloping beds. J. Geophys. Res. 115, F03040. https://doi.org/10.1029/2009JF001462.

NU

[49] Nakada, S., Fujii, T., 1993. Preliminary report on the activity at Unzen Volcano (Japan), November 1990-November 1991: Dacite lava domes and pyroclastic flows. J. Volcanol. Geotherm. Res. 54 (3–4), 319–333. https://doi.org/10.1016/0377-0273(93)90070-8.

MA

[50] Parker, G., Garcia, M., Fukushima, Y., Yu, W., 1987. Experiments on turbidity currents over an erodible bed. J. Hydraul. Res. 25 (1), 123–147. https://doi.org/10.1080/00221688709499292.

PT E

D

[51] Roche, O., 2012. Depositional processes and gas pore pressure in pyroclastic flows: an experimental perspective. Bull. Volcanol. 74 (8), 1807–1820. https://doi.org/10.1007/s00445012-0639-4. [52] Roche, O., Buesch, D.C., Valentine, G.A., 2016. Slow-moving and far-travelled dense pyroclastic flows during the Peach Spring super-eruption. Nat. Comm. 7, 10890. https://doi.org/10.1038/ncomms10890.

AC

CE

[53] Roche, O., Montserrat, S., Ni˜ no, Y., Tamburrino, A., 2008. Experimental observations of water-like behavior of initially fluidized, dam break granular flows and their relevance for the propagation of ash-rich pyroclastic flows. J. Geophys. Res. 113 (B12), B12203. https://doi.org/10.1029/2008JB005664. [54] Roche, O., Ni˜ no, Y., Mangeney, A., Brand, B., Pollock, N., Valentine, G.A., 2013. Dynamic pore-pressure variations induce substrate erosion by pyroclastic flows. Geology 41 (10), 1107–1110. https://doi.org/10.1130/G34668.1. [55] Scott, W.E., Hoblitt, R.P., Torres, R.C., Self, S., Martinez, M.M.L., Nillos, T., 1996. Pyroclastic flows of the June 15, 1991, climactic eruption of Mount Pinatubo. In: Newhall, C.G., Punongbayan, R.S. (Eds.), Fire and Mud: Eruptions and Lahars of Mount Pinatubo, Philippines, University of Washington Press, Seattle, 545–570. [56] Shimizu, H.A., Koyaguchi, T., Suzuki, Y.J., 2017. A numerical shallow-water model for gravity currents for a wide range of density differences. Prog. Earth Planet. Sci. 4 (1), 8. https://doi.org/10.1186/s40645-017-0120-2. [57] Slim, A.C., Huppert, H.E., 2011. Axisymmetric, constantly supplied gravity currents at high Reynolds number. J. Fluid Mech. 675, 540–551. https://doi.org/10.1017/jfm.2011.71. 41

ACCEPTED MANUSCRIPT Journal of Volcanology and Geothermal Research

[58] Sparks, R.S.J., 1976. Grain size variations in ignimbrites and implications for the transport of pyroclastic flows. Sedimentology 23 (2), 147–188. https://doi.org/10.1111/j.13653091.1976.tb00045.x. [59] Sparks, R.S.J., Bonnecaze, R.T., Huppert, H.E., Lister, J.R., Hallworth, M.A., Mader, H., Phillips, J., 1993. Sediment-laden gravity currents with reversing buoyancy. Earth Planet. Sci. Lett. 114 (2–3), 243–257. https://doi.org/10.1016/0012-821X(93)90028-8.

RI PT

[60] Streck, M.J., Grunder, A.L., 1995. Crystallization and welding variations in a widespread ignimbrite sheet; the Rattlesnake Tuff, eastern Oregon, USA. Bull. Volcanol. 57 (3), 151– 169. https://doi.org/10.1007/BF00265035. [61] Sulpizio, R., Dellino, P., Doronzo, D.M., Sarocchi, D., 2014. Pyroclastic density currents: state of the art and perspectives. J. Volcanol. Geotherm. Res. 283, 36–65. https://doi.org/10.1016/j.jvolgeores.2014.06.014.

SC

[62] Suzuki, Y.J., Iguchi, M., in press. Determination of the mass eruption rate for the 2014 Mount Kelud eruption using three-dimensional numerical simulations of volcanic plumes. J. Volcanol. Geotherm. Res. https://doi.org/10.1016/j.jvolgeores.2017.06.011.

NU

[63] Suzuki, Y.J., Koyaguchi, T., 2009. A three-dimensional numerical simulation of spreading umbrella clouds. J. Geophys. Res. 114 (B3), B03209. https://doi.org/10.1029/2007JB005369.

MA

[64] Suzuki, Y.J., Koyaguchi, T., 2012. 3-D numerical simulations of eruption column collapse: Effects of vent size on pressure-balanced jet/plumes. J. Volcanol. Geotherm. Res. 221–222, 1–13. https://doi.org/10.1016/j.jvolgeores.2012.01.013.

PT E

D

[65] Suzuki, Y.J., Koyaguchi, T., Ogawa, M., Hachisu, I., 2005. A numerical study of turbulent mixing in eruption clouds using a three-dimensional fluid dynamics model. J. Geophys. Res. 110 (B8), B08201. https://doi.org/10.1029/2004JB003460. [66] Sweeney, M.R., Valentine, G.A., 2017. Impact zone dynamics of dilute mono- and polydisperse jets and their implications for the initial conditions of pyroclastic density currents. Phys. Fluids 29 (9), 093304. https://doi.org/10.1063/1.5004197.

CE

[67] Toro, E.F., 2001. Shock-Capturing Methods for Free-Surface Shallow Flows. Wiley, Chichester.

AC

[68] Toro, E.F., 2009. Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction. Springer, Berlin. [69] Ungarish, M., 2009. An Introduction to Gravity Currents and Intrusions. CRC Press, Boca Raton. [70] Valentine, G.A., 1987. Stratified flow in pyroclastic surges. Bull. Volcanol. 49 (4), 616–630. https://doi.org/10.1007/BF01079967. [71] Valentine, G.A., Sweeney, M.R., 2018. Compressible flow phenomena at inception of lateral density currents fed by collapsing gas-particle mixtures. J. Geophys. Res. 123 (2), 1286– 1302. https://doi.org/10.1002/2017JB015129. [72] Wiesner, M.G., Wetzel, A., Catane, S.G., Listanco, E.L., Mirabueno, H.T., 2004. Grain size, areal thickness distribution and controls on sedimentation of the 1991 Mount Pinatubo tephra layer in the South China Sea. Bull. Volcanol. 66 (3), 226–242. https://doi.org/10.1007/s00445-003-0306-x. 42

ACCEPTED MANUSCRIPT Journal of Volcanology and Geothermal Research

[73] Wiesner, M.G., Wetzel, A., Catane, S.G., Listanco, E.L., Mirabueno, H.T., 2005. Grain size, areal thickness distribution and controls on sedimentation of the 1991 Mount Pinatubo tephra layer in the South China Sea. Bull. Volcanol. 67 (5), 490–495. https://doi.org/10.1007/s00445-005-0421-y. [74] Wilson, C.J.N., 1985. The Taupo eruption, New Zealand. II. The Taupo Ignimbrite. Phil. Trans. R. Soc. Lond. A 314 (1529), 229–310. https://doi.org/10.1098/rsta.1985.0020. [75] Wilson, C.J.N., 1991. Ignimbrite morphology and the effects of erosion: a New Zealand case study. Bull. Volcanol. 53 (8), 635–644. https://doi.org/10.1007/BF00493690.

AC

CE

PT E

D

MA

NU

SC

RI PT

[76] Wilson, C.J.N., Houghton, B.F., Kampt, P.J.J., McWilliamst, M.O., 1995. An exceptionally widespread ignimbrite with implications for pyroclastic flow emplacement. Nature 378 (6557), 605–607. https://doi.org/10.1038/378605a0.

43