## Abstract

Most models of volcanic ash flows assume that the flow is either dilute or dense, with dynamics dominated by fluid turbulence or particle collisions, respectively. However, most naturally occurring flows feature both of these end members. To this end, a two-layer model for the formation of dense pyroclastic basal flows from dilute, collapsing volcanic eruption columns is presented. Depth-averaged, constant temperature, continuum conservation equations to describe the collapsing dilute current are derived. A dense basal flow is then considered to form at the base of this current owing to sedimentation of particles and is modelled as a granular avalanche of constant density. We present results which show that the two-layer model can predict much larger maximum runouts than would be expected from single-layer models, based on either dilute or dense conditions, as the dilute surge can outrun the dense granular flow, or vice versa, depending on conditions.

## 1. Introduction

Pyroclastic currents are highly hazardous flows that consist of a suspension of particles in hot gas travelling at speeds up to hundreds of metres per second. They are commonly observed to have two main components (figure 1); a dense basal part, herein called the ‘pyroclastic flow’ with volumetric particle concentrations of tens of per cent, overlain by a dilute turbulent component, called the ‘pyroclastic surge’ with concentrations of the order of 1 per cent (Sparks *et al.* 1973; Druitt 1998).

Pyroclastic currents are frequently produced by the collapse of an eruption column (Sparks & Wilson 1976; Wilson 1976) that initially produces a dilute current from which the dense basal flow develops by mass transfer (Druitt & Sparks 1982). Exchange of mass and heat between the two layers allows for a spectrum of flows, being dominated by either the dense or dilute part (Fisher 1979; Druitt 1998). For example, in the 1997 explosions of Soufriére Hills Volcano, Montserrat, columns collapsed from heights of 400–750 m and generated dilute pyroclastic currents. These rapidly decelerated as particles sedimented into the basal regions, resulting in dense basal flows emerging from beneath at tens of metres per second (Druitt *et al.* 2002).

Most models of pyroclastic density currents consider either a dilute turbulent suspension (e.g. Dartevelle *et al.* 2004; Neri *et al.* 2007) or a dense pyroclastic flow (e.g. Druitt 1998; Kelfoun *et al.* 2009). The dilute models take either a multi-phase approach (e.g. Neri *et al.* 2007) or consider a ‘bulk-continuum’ pseudo-fluid mixture (e.g. Denlinger 1987). The simplest models use a depth-averaged assumption, similar to the classic ‘shallow-water’ equations, by assuming the horizontal extent of the current greatly exceeds its vertical height (e.g. Bursik & Woods 1996). The suspending turbulent current is considered well mixed with no fluid phase loss, and the concentration decreases owing to particles settling through a basal viscous sub-layer (Martin & Nokes 1989).

Dense pyroclastic flows have been modelled as Bingham fluids with a constant fluid viscosity and yield strength (e.g. Sparks 1976), as a non-deforming, frictional sliding block (e.g. Hayashi & Self 1992; Druitt 1998), as a rapid granular flow via discrete element simulations (e.g. Staron & Hinch 2005) or as a pseudo-fluid (e.g. Huppert & Dade 1998; Mangeney-Castelnau *et al.* 2003). In the latter, simplifications come from depth-averaging, and by treating the material as a constant density, incompressible, deformable Coulomb continuum with a stress-free surface.

These single-layer models do not capture the dynamics of many volcanic ash flows, which feature both dilute and concentrated regions within which turbulent suspension mechanics and granular interactions are important. Dartevelle *et al.* (2004) and Dufek & Bergantz (2007) have demonstrated that multi-phase models need to consider the whole spectrum of possible granular concentrations, to simulate mass exchanges between the dense and dilute regions of pyroclastic currents. Denlinger (1987) and Takahashi & Tsujimoto (2000) have developed models for the formation of upper dilute clouds from dense block-and-ash flows, and Wadge *et al.* (1998) and Widiwijayanti *et al.* (2009) have introduced extensions to map ash-cloud surge impacts beyond the dense basal flow. Two-layer models have been developed to model turbidity currents in the ocean (e.g. Drago 2002), and powder snow-dense flow avalanches (e.g. Eglit 1998; Issler 1998; Zwinger *et al.* 2003). These latter models consider the generation of an upper dilute powder-snow avalanche from the top of a basal dense snow flow avalanche through processes of re-suspension within a small interface layer, incorporating shear between the layers, elutriation and entrainment.

Here, we expand upon the study of Doyle *et al.* (2008) and present additional developments of the two-layer model for the formation of these dense basal flows from the upper dilute column collapse currents, including an analysis of the steady-flow problem. Observations of column collapse currents (e.g. Nairn & Self 1978; Cole *et al.* 1998; Druitt *et al.* 2002) and deposit interpretations (e.g. Sparks *et al.* 1973; Fisher 1979) indicate that a distinct interface exists between these dense and dilute layers. The bulk depositional flow unit, which is composed of light pumice and dense lithic clasts, is commonly overlain by fine, well-sorted and cross-bedded deposits from the upper turbulent cloud. Surface-wave-derived features on pumice flow deposit surfaces further illustrate that a sharp density contrast must exist between the layers (E. S. Calder 2009, personal communication).

The model presented herein has been specifically designed to elucidate the principles that govern the formation of these dense basal flows from the overlying dilute layer. There is assumed to be a distinct interface between the two layers, thus permitting their independent motion, and we develop separate models for their dynamics. We present the full derivation of this model in §2, and the numerical method in §3. Example calculations are presented in §4, and the steady-flow situation is examined in §5. In §6, we discuss the importance of adopting this model idealization to capture the maximum potential runout of these hazardous flows.

## 2. Development of the two-layer model

The two-layer model is depicted in figure 1. Within the dilute layer (*c*<*z*<*f*), the volume fraction of particles, *ϕ*, is small but varying temporally and spatially and the particles sediment relative to the suspending, turbulent air. In contrast, within the dense layer (*b*<*z*<*c*), the volume fraction of particles is high and is assumed to be constant (*ϕ*=*ϕ*_{b}); the relatively heavy particles are supported by the particle collisions and both phases have the same velocity. Therefore, the layers exhibit different balances of forces and hence dominant dynamics. However, for both layers, we retain the shallow-flow assumption and therefore negligible accelerations normal to the base, and the velocity is predominantly parallel to the underlying boundary.

Calculations and observations suggest that pyroclastic material cools by a maximum of 350 K in the very proximal region of a collapsing eruption column (Sparks *et al.* 1978), owing to a combination of air entrainment, degassing of juvenile clasts and incorporation of cold lithics (Calder *et al.* 2000). The temperature along the pyroclastic flow path is found to be generally constant (Calder *et al.* 2000) as the solid particles have a high thermal energy and large heat capacity, and so the hot pyroclastic current can heat and expand entrained air with little cooling of the flow (see discussions in Kieffer 1981; Valentine 1987; Ishimine 2005).

We thus model the dilute layer as a turbulent suspension of particles with a constant average temperature of 850 K based upon a vent temperature of 1100 K (Sparks *et al.* 1978, 1997). It will be further assumed to be incompressible (Huppert & Dade 1998). In our model idealization, we do not consider entrainment of the surrounding gas during the subsequent motion of the current. This has the potential to alter the density difference between the flowing current and the surrounding ambient. However, Hallworth *et al.* (1996) show that the buoyancy of the current is conserved under this mixing owing to a balance between reducing density difference and increasing flow height (i.e. thickness). The entrained ambient fluid must be accelerated up to the velocity of the current, and this exerts an effective drag on the motion. We neglect such effects in this investigation.

The density of the solid and interstitial gas phases are denoted by *ρ*_{s} and *ρ*_{g}, respectively, while the ambient is of density *ρ*_{a}. Thus, the bulk density of the current is given by *β*=*ρ*_{s}*ϕ*+*ρ*_{g}(1−*ϕ*). It is important to note that even in the dilute layer, while the particle concentration may be very dilute (*ϕ*≤0.01), the high density of the particles (up to 2400 kg m^{−3}) significantly enhances the bulk density and thus unlike many previous studies of dilute, particle-driven flows, these flows are non-Boussinesq. Models for Boussinesq flows commonly impose a finite height at the flow front, often described by a Froude condition (e.g. Huppert & Simpson 1980). However, this is not appropriate for non-Boussinesq models and so, in the manner of Ancey *et al.* (2007), we impose a negligible front height (§3*c*).

The key feature of the dilute layer is the potential for the relatively heavy particles to settle under gravity through the interstitial gas. On the assumption that the inertia of the individual particles is sufficiently small, the velocities of the solid and gas phases, denoted by **v**_{s} and **v**_{g}, respectively, are given by , where and are unit vectors aligned with the coordinate axes, is the gradient of the underlying boundary and the settling speed *ω*=|*ω*| may be determined from the expression (Sparks *et al.* 1997)
2.1
Here, *d* represents particle diameter and *C*_{d} is a drag coefficient, which is assumed close to unity (Woods & Bursik 1991). In what follows, we assume that the density-induced component of the gas flow far exceeds the component of the settling velocity downslope , and so will be neglected from the expression for the solid phase velocity above.

The flow in the dilute layer is turbulent, which suspends the relatively heavy particles. Thus, the governing equations model the properties of the flow averaged over fluctuations. This process introduces turbulent fluxes that are non-vanishing and contribute both to the effective shear stress generated by the flow and to the suspension of particles. When the flow is sufficiently dilute (*ϕ*≪1), and the velocity and solids concentration perturbations within the flow are incorporated via Reynolds averaging, the equations that express mass conservation of each phase are given by
2.2
where *K* denotes the turbulence-induced sediment diffusivity (Dyer & Soulsby 1988). The equation of motion for the mixture is given by
2.3
where **g** denotes gravitational acceleration, *P* the fluid pressure and *τ* the shear stress generated by the dilute mixture, for which the stresses generated by direct interactions of particles are neglected.

Turning to the dense basal layer, which is generated by the accumulation of particles sedimenting from the base of this dilute current, we assume a well-defined interface with the upper dilute current, as inferred from observational evidence of flows and deposits (§1). The key dynamical feature of the basal layer is that particle collisions generate significant stresses that may no longer be neglected. In particular, the flow generates sufficient normal stress to support the weight of the flowing material, and thus it is assumed that the volume fraction of particulate remains constant and that the gaseous and solid phase flow at the same velocity (**v**_{s}=**v**_{g}≡**v**_{b}). In the same manner as the dilute layer, all effects of heat transfer are neglected and the motion is assumed to be incompressible. Mass conservation and momentum balance are given by
2.4
where *β*_{b}=*ρ*_{g}(1−*ϕ*_{b})+*ρ*_{s}*ϕ*_{b} and *σ* denotes the stresses generated by the interactions between particles, and the fluid-generated stress is ignored.

We must now define suitable boundary conditions. At the upper surface of the dilute current (*z*=*f*(*x*,*t*), figure 1), it is assumed that, in the absence of entrainment, the interface is advected with the interstitial gas velocity and that there is no flux of particles across it. This implies that
2.5
where is a unit normal to the interface and the velocity fields have components **v**_{g}=(*u*_{g},*w*_{g}) and **v**_{s}=(*u*_{s},*w*_{s}). Additionally, as will be discussed below, we impose a continuous pressure field and vanishing shear stress at this interface. The boundary underlying the dense basal current is rigid and impermeable to the flow, and additionally we do not allow mass transfer between the flow and the underlying bed (i.e. there is assumed to be no erosion or deposition). Thus, denoting the basal velocity by **v**_{s}=(*u*_{s},*w*_{s}), the kinematic boundary condition at *z*=*b*(*x*) is given by
2.6
Further, we allow the solid phase to slip at the base and generate a Coulomb shear stress, which will be formulated below. At the interface between the dilute and dense layers, *z*=*c*(*x*,*t*), we impose conditions that conserve the mass flux of each species and the combined momentum flux. Denoting the unit normal to the interface by and the velocity of the interface by **v**_{i}, we impose
2.7
2.8and
2.9
where […]^{+}_{−} denotes the evaluation of the terms within the brackets either side of the interface.

The above presentation provides a description of the two-layer flow. However, we now simplify the models for the dynamics in the two layers by first imposing the shallow-layer approximation and then by depth integration. Under the shallow-layer approximation, the motion is taken to be predominantly parallel to the underlying boundary, and the vertical accelerations of each phase are assumed negligible. Then to leading order, the dominant balances in the vertical component of the combined momentum equation are
2.10
in the dilute current (*c*<*z*<*f*) and dense current (*b*<*z*<*c*), respectively. In addition, the assumption of shallowness implies that the interfaces are locally parallel with the underlying boundary and so . Thus, the continuity of the mass fluxes of each phase at the interface *z*=*c* implies that
2.11
where *ϕ*_{c} denotes the concentration at *z*=*c*^{+}. These may be further simplified given that the upper layer is dilute (*ϕ*_{c}≪1,*ϕ*_{c}≪*ϕ*_{b}) to yield
2.12
We may now integrate the expression for the conservation of fluid mass to *O*(*ϕ*_{c}),
2.13
where the height of the layer, *H*=*f*−*c*. Next, we integrate the equation governing the evolution of the volume fraction of particles, after first neglecting the streamwise diffusive flux (*K*∂*ϕ*/∂*x*) in comparison to the vertical diffusive flux (*K*∂*ϕ*/∂*z*), which follows owing to the assumption of shallowness. Further, as described above, we assume that *u*_{s}=*u*_{g} and and impose the boundary conditions of no flux of particles through the top of the layer (*z*=*f*) and no erosion of particles at *z*=*c* (where *K*∂*ϕ*/∂*z*=0). Thus, we find that to leading order,
2.14
We now analyse the streamwise component of the combined momentum equation for the dilute layer. First, we construct the pressure field, which is in hydrostatic balance to leading order. Imposing *P*=*P*_{0} at some distant elevated, horizontal surface , we find
2.15
Then, on the assumption that the streamwise velocity, *u*_{g}, is independent of depth, we find that the leading-order expression for momentum balance is given by
2.16
where ∂*τ*_{xx}/∂*x* is negligible relative to ∂*τ*_{xz}/∂*z* and the shear stress on the upper surface is neglected. The terms on the right-hand side of this momentum equation represent the downslope acceleration, the momentum lost from the dilute current owing to the particulate settling and the interfacial drag, which we model as proportional to the square of the velocity difference between the layers
2.17
This parametrization is analogous to a basal drag (e.g. Parker *et al.* 1986; Hogg *et al.* 2005), and we adopt a small drag coefficient *C*_{D}=0.001.

Now we depth-integrate the equations that model the dense basal layer, which has a height of *h*=*c*−*b*. First, mass conservation yields
2.18
The interfacial conditions at *z*=*c* imply, to leading order in terms of the volume fraction, *ϕ*, that both the shear and normal stress are continuous. Thus, *σ*_{zz}(*c*)=−*P*(*c*), and so integrating the expression for hydrostatic balance, we find that
2.19

Next, we assume that the stress tensor is isotropic and so *σ*_{xx}=*σ*_{zz} (see Pouliquen & Forterre 2002; Gray *et al.* 2003; Mangeney-Castelnau *et al.* 2003). Then, integrating the combined expression of momentum balance for both phases yields
2.20
The terms on the right-hand side of this equation represent the downslope acceleration, the pressure gradient associated with variations in the height of the overlying dilute layer, the flux of streamwise momentum from the dilute to the dense layer associated with the mass flux across the interface and the interfacial and basal drags. The latter is modelled using a Coulomb-like description so that
2.21
where *δ* is the dynamic angle of friction (Savage & Hutter 1989). The basal drag *σ*_{xz}(*z*=*b*) far exceeds the interfacial drag *σ*_{xz}(*z*=*c*)=*τ*_{c}, and so the latter may be neglected, unless the height of the basal flow is small (§4).

We now replace the integrals with averaged quantities. First, we define
2.22
This implies an averaged bulk density of . We then assume that
2.23
These expressions are exact identities only when the volume fraction and velocity fields are vertically uniform. Finally, we evaluate the volume fraction and velocity fields at the interface *z*=*c*, and to this end we assume that they adopt the average values and . Some studies have included factors to account for non-uniformities in the vertical structure of these fields (Parker *et al.* 1986; Hogg & Pritchard 2004). However, we neglect such factors in this study.

Finally, we present the system of governing equations for this two-layer model of pyroclastic flows. For the dilute layer, we have the following expressions for the mass conservation of the gas phase (2.24) and the combined mixture (2.25), and the corresponding combined balance of momentum (2.26):
2.24
2.25and
2.26
Note that in the derivation of the pressure terms, we have neglected curvature terms of *O*(∂_{x}*θ*(*x*)). However, changes in the average slope angle and thus weak curvature effects are still incorporated implicitly, via *θ*(*x*) (Savage & Nohguchi 1988).

For the dense basal layer, we may further simplify the governing equations by noting that *β*_{b}≫*ρ*_{a} and thus obtain the following equations that govern the evolution of the mass and momentum of the flow:
2.27
and
2.28
The second term on the right-hand side of equation (2.28) represents the pressure gradient on the basal layer exerted by variations in the height of the overlying dilute layer. This effect is negligible provided , a condition that is usually satisfied.

## 3. Numerical method

We use the finite-volume method (Leveque 2002) to solve these governing equations, adopting a first-order upwind Godunov approach (for full details, see Doyle 2007). This conservative approach is commonly adopted for shallow-water-type flows (e.g. Denlinger & Iverson 2004), and we use it owing to its shock-capturing capabilities, and its ability to apply the boundary conditions non-invasively. This numerical method solves a system of hyperbolic equations of the nonlinear form *q*_{t}+*f*(*q*)_{x}=*ψ*(*q*,*x*,*t*), where *q*=*q*(*x*,*t*)∈*R*^{m}. The terms *ψ*(*q*,*x*,*t*) are referred to as the source terms. We adopt the Godunov fractional step approach (Leveque 2002) to include these terms, where the full equation is split into two problems: (i) *q*_{t}+*f*(*q*)_{x}=0 and (ii) *q*_{t}=*ψ*(*q*), and solved using the procedure described in appendix A.

### (a) Basal-flow numerical method

To use the entire structure of the Riemann solver in the Godunov method is expensive (appendix A). Thus, for the basal-flow equations (2.27) and (2.28), we use the standard shallow-water Roe averages (Leveque 2002; Denlinger & Iverson 2004; Larrieu *et al.* 2006). The source terms are included via the fractional step approach, where they are further split into linear (basal friction, sedimentation, gravitational driving) and spatial (topographic) contributions. We use the flux-difference splitting method of Hubbard & Garcia-Navarro (2000) to solve the spatial source contributions, which are decomposed in a manner similar to the main finite-volume step (Hubbard & Dodd 2002; Denlinger & Iverson 2004). This method ensures that the appropriate equilibria of the underlying mathematical model are maintained by the numerical scheme. The linear source terms are then solved via the TRBDF2 (second-order trapezoidal rule and backward difference formula) method (Leveque 2002), which is composed of a two-step Runge–Kutta and implicit trapezoid BDF (backward difference formula) method.^{1}

### (b) Dilute current numerical method

The first-order, upwind, Godunov method is also used for the dilute current (2.24)–(2.26), again adopting a fractional-step approach. However, it is not possible to find a suitable conservative Roe average state that obeys the required conditions for an approximate Riemann solver (Doyle 2007). Thus, we use the alternative ‘f-wave’ approach (Leveque 2002; Bale *et al.* 2003), which guarantees numerical conservation when any linearization of the governing equations is made. This approach, described in appendix B, still splits the problem into two steps; however, any spatial source terms are included in the fluctuation calculations of the initial Godunov step and not treated separately. Linear source terms are handled in the same manner as the basal flow.

### (c) Code development and validation

For both layers, the numerical code is developed in C++ (Doyle 2007). Following (Leveque 2002), we impose a Courant–Friedrich–Levy condition to preserve stability, ensure convergence and control the time step. In the basal-flow layer, a Harten–Hyman entropy fix prevents the generation of entropy-violating shock waves, and a modified version of this fix is used for the f-wave method in the upper dilute current (appendix C). For both layers, we impose a stationary negligible prelayer of height *ϵ*=10^{−8} m, to prevent non-physical numerical solutions in both the Roe and f-wave solvers, owing to division by zero height (*h* or *H*=0) (e.g. Larrieu *et al.* 2006). Solutions do not depend upon minor variations to this prelayer value. For the dilute current, this prelayer is composed of a very dilute concentration of particles (*ϕ*_{0}=10^{−8}) suspended in the ambient air *ρ*_{a}.

A reflective boundary condition is imposed at the left origin (*x*=0 m), and an outflow condition is assumed for the right boundary. The dense basal front *x*_{gf} is determined by the location where the height *h* reduces to the minimum grain size of 100 μm, and the dilute current front *x*_{βf} where the bulk density equals that in the prelayer *β*_{a}(*ϕ*_{a}) (§2). Dilute current calculations continue until the average bulk density , along its length, is within 10 per cent of the ambient density *β*_{a}=*f*(*ϕ*_{a}).

The numerical solver for the dense basal layer is validated against the classic shallow-water dry dam-break solution (Ritter 1892), and the analytical results of Mangeney *et al.* (2000) for flows down a slope with friction (figure 2*a*–*c*). The numerical solver for the dilute current is also validated against the analytical dam-break solution (Ritter 1892; Leveque 2002), demonstrating that the modified entropy fix produces accurate results (figure 2*d*). In addition, the correct solution is found for a dam-break solution with a tailwater (Stoker 1957). Finally, we validate the numerical method for the inclusion of the basal topography terms against the example test cases of Leveque (1998), Hubbard & Garcia-Navarro (2000) and Caleffi *et al.* (2003), with a difference of less than 1.5 per cent between the height profiles produced by both the dense basal and dilute current flows (Doyle 2007). For both numerical layers, runouts differ by less than 0.5 per cent for grids of 2 or 5 m. Thus, calculations with initial collapsing column heights of *H*_{0}≤1100 m use 2 m cells, and calculations with taller columns use 5 m cells.

## 4. Examples of discrete eruption-column collapses

We now present two contrasting model calculations, to illustrate dense and dilute end-member flow behaviours. For these, we assume and an aspect ratio *a*=*H*_{0}/*x*_{0}=3, where *x*_{0} is the half-width of the instantaneous, constant volume, column collapse. Although this initial aspect ratio is relatively large, it is reasonable to assume that the shallow-water model provides an accurate description of the motion after the earliest stages of propagation. The derived basal flow has a bulk concentration of *ϕ*_{b}=0.5, consistent with field estimates (Sparks 1976; Druitt 1998), with a basal friction angle of *δ*=10^{°} (e.g. Freundt 1999). The basal drag *σ*_{xz}(*z*=*b*) far exceeds the interfacial drag *σ*_{xz}(*z*=*c*)=*τ*_{c}, and so the latter is neglected (§2), and we investigate the effect of this assumption in §4. We additionally neglect the pressure gradient imposed on the basal layer owing to the overlying dilute-layer height (2.28). All columns are released from rest over a simplified terrain, encompassing a slope inclined at *θ*=13^{°} for 10 km flanked by a 1^{°} plateau to 25 km, where there is a gradual, even, decrease in slope from 13^{°} to 1^{°} over a distance of 1 km.

For a short (*H*_{0}=550 m), coarse-grained (*d*=8 mm) column, the current transfers all its mass to the basal flow while still on the inclined slope (figure 3*a*). At this time, the derived basal flow has a wedge-shaped morphology, thinning towards the basal-flow front (figure 3*b*). The slope angle exceeds the basal friction angle *δ*, and thus independent motion of the basal flow occurs (figure 3*a*). This results in a thickening of the basal flow towards its front (figure 3*c*). The basal flow does not finally come to rest until it reaches the low angled plateau, where frictional deceleration occurs (figure 3*a*,*d*). Thus, the column collapse is dominated by a dense basal flow for the majority of its total propagation time. For a tall (*H*_{0}=1600 m), fine-grained (*d*=550 μm) column, the basal flow exhibits very different behaviour (figure 4*a*). The dilute current retains significant suspended mass when it reaches the low angled plateau, and frictional forces dominate the basal flow, preventing it from outrunning its parent flow. The final deposit thins towards its front (figure 4*b*), with a localized thickening at the back owing to material settling from the steeper slopes in that region. For these initial conditions, the propagation of the collapse is dominated by a dilute current for the majority of its total propagation time.

For these two simulations, the basal flow’s interfacial drag with the upper dilute current *σ*_{xz}(*z*=*c*)=*τ*_{c} is neglected, on the assumption that the basal drag far exceeds it (§2). However, as discussed by Pudasaini & Hutter (2007), if mass is added to the top of the basal flow, this mass must also be accelerated to reach the motion of the particles at the top, and the traction may no longer vanish at the upper surface. If instead we assume a continuous stress across the interface, the imposed traction on the surface of the basal flow results in an increase in the basal velocity. For the short, coarse-grained column, this decreases the final runout time by 8 per cent, and runout distance by less than 1 per cent (figure 5). For the tall, fine-grained column, this causes the basal-flow front to propagate independently beyond the termination distance of the dilute current, decelerating on the plateau until it comes to rest (figure 5). This is not observed for the stress-free calculations (figure 4*a*). However, the continuous-stress assumption has a minimal effect upon the upper dilute current, with only a minor decrease in its basal drag. For the tall, fine-grained column, this results in a slight increase in the propagation velocity of the dilute current towards the end of its propagation, with a decrease in stoppage time of less than 1 per cent. The maximum propagation distance of the dilute current, for both initial column heights, is negligibly different when either stress assumption is adopted.

## 5. Steady flow

Pyroclastic currents can establish a steady state, if fed continuously for the duration of a sustained column-collapse eruption, which can last several hours (Sparks *et al.* 1978; Branney & Kokelaar 1992). It is thus insightful to examine the steady, spatially evolving, two-layer motion predicted by this model, as these flows clearly illustrate the important coupling between the two layers and the possibility for independent dynamics. To this end, we seek steady solutions for the heights and velocities of each layer as well as the average bulk density of the upper dilute layer . Under these steady conditions, the expressions for the conservation of mass decouple from the momentum equations and may be integrated to give
5.1
5.2and
5.3
where *q* and *ϕ*_{0} denote the volume flux of fluid per unit width and the initial concentration of particles fed into the dilute ash cloud at the source. It is further assumed that there is initially no flux of material in the underlying dense layer. The key feature of the flow is that the volume fraction of particles in the dilute current decays exponentially (5.2), and thus delivers a progressively diminishing flux into the dense layer.

To complete the solutions, we must integrate the momentum equations, noting that, very close to the source, the height of the dense layer *h*(*x*) is sufficiently small so that the interfacial drag with the overlying dilute layer cannot be neglected. Further, we deduce that if the drag is written as equation (2.17), then at the source, the basal velocity is equal to the ash-cloud velocity . Finally, we specify this source velocity provided the conditions at the source are supercritical ; if they are subcritical, then the problem is already fully determined.

For flows along a planar surface (∂*b*/∂*x*≡0), we find that the momentum equations are given by
5.4
and
5.5
We integrate these equations until the bulk density of the dilute layer falls to the density of the surrounding atmosphere, at which point any residual particles in the ash cloud are lofted into the atmosphere. This occurs at
5.6
The spatial evolution of the flow, for a given choice of parameters, is illustrated in figure 6. We note that over this distance the dilute ash cloud accelerates and thins, and that the dense basal flow also accelerates but thickens owing to the source of particulate from the overlying layer. The acceleration of the dilute current is due to the downslope gravitational acceleration exceeding the deceleration associated with the sedimentation-related reduction of bulk density. For other parameter values, we can find flows that may initially accelerate, but then rapidly slow as the loss of suspended particles strongly affects the overall dynamics. It is noteworthy that for the calculated fields plotted in figure 6*a*–*c*, the effects of interfacial drag are negligible, apart from being very close to the source. In addition, the height of the dense layer remains much smaller than that of the dilute layer because the density of the dilute layer is so much smaller than that of the dense layer (*β*≪*β*_{b}). This implies that, to leading order, the velocity of the dense layer is given by
5.7
and thus, for some parameter values, it is possible for the velocity of the dense layer to exceed that of the dilute layer (figure 6*d*,*e*), illustrating the independent layer dynamics.

After the dilute ash-cloud lofts, the volume flux of particles per unit width in the dense layer then remains constant and is given by 5.8 Thereafter, the dynamics evolve according to equation (5.5) and the effects of interfacial drag are now negligible. When , the flow continues to accelerate down the plane. However, when the surrounding plateau is reached, the basal drag now exceeds the downslope acceleration and the flow begins to decelerate. The ensuing dynamics are complicated because the motion becomes time dependent and may feature bores that adjust the fast-moving granular layer to a more slowly moving, or arrested, pile of grains.

## 6. Discussion and conclusions

Dominant flow behaviour in the pyroclastic current is expected to be very different for dense flows and dilute surges, and to change owing to topography and blocking of the dense underflow (e.g. Fisher 1990; Browne & Gardner 2005). Evidence of spatial changes in deposits can be used to infer changes in the character of the flow with time and distance. For example, proximal lag breccia and cross-stratified surge deposits have been interpreted as being formed when the flow is in an initial expanded and dilute state (Sparks 1976; Druitt & Sparks 1982; Calder *et al.* 2000). Distal changes to massive deposits, commonly characterized by light pumice at the top and dense lithics at the base, suggest transformation to predominantly dense flows with distance, as occurs during eruption-column collapse (§1). These spatial variations in deposit and flow properties require consideration within the framework of models that incorporate the dual character of the flows, capturing the fundamentally different physics of both the dilute and dense regions. Our theory provides an approach to coupling the two components based on mass transfer from the dilute part to the dense flow. This two-layer model has the capability to capture spatial changes in flow dominance, where the independent propagation of each flow layer can be demonstrated for both steady and discrete flows (§§4 and 5).

This model idealization focuses on the formation of dense pyroclastic flows by concentrating on sedimentation processes. There are thus several caveats implicit to these results. The model is only in the depositional regime, the upward entrainment or erosion of particles from the basal flow back into the upper dilute current, and the effect that this has on the stress coupling at the interface has not been treated. The formation of further upper ash clouds from the surface of the dense basal flow can be incorporated into the model via shear, elutriation or an upward gas flux (Denlinger 1987; Takahashi & Tsujimoto 2000). This upward migration of particles could be incorporated via a negligibly thin re-suspension interface layer (§1, Zwinger *et al.* 2003). However, observations indicate that the interface is more probably a sharp transition than a transition zone (§1), and that sedimentation dominates flow propagation.

In our model idealization, we do not consider air entrainment (§2). Entrainment does not alter the buoyancy of the current (Hallworth *et al.* 1996), as any thickening of the current is balanced by a reduction in density owing to a dilution of the suspended particles, but it does induce an effective drag as entrained fluid is accelerated to the velocity of the dilute current. Entrainment-related dilution of the dilute current is expected to reduce the flux of particulate settling into the dense basal layer, and thus decrease the rate at which the basal layer forms. In turn, the upper dilute current may actually travel farther than that predicted by our model. Future development of this model could incorporate entrainment effects in the manner of Ellison & Turner (1959) and Parker *et al.* (1986), and compressibility in the manner of Timmermans *et al.* (2001).

Further development of the upper turbulent dilute-current model could also include a thermal-energy equation (e.g. Denlinger 1987; Takahashi & Tsujimoto 2000), a hindered settling velocity to represent how particle settling can be dependent upon the concentration of fines present (Druitt 1995), and an erosional term. Future models will thus need to develop criteria to discriminate between erosional basal flow or deposit formation. Further sophistication of the numerical method at the flow front could include the use of the kinetic scheme of Mangeney-Castelnau *et al.* (2003) that can handle discontinuous solutions, instead of applying a negligible prelayer throughout the domain (§3).

The primary goal of most pyroclastic current models is to be used as a tool for hazard analysis and assessment. It is thus vital to model these pyroclastic currents correctly. Only by including both the initial dilute cloud and the formation of derived dense basal flows can the maximum potential runout distance be calculated. Our focus here has been the development of such a two-layer model. The results demonstrate that if only dilute-current physics or single-layer models are adopted, then they are likely to underestimate the maximum runout of the current, as they do not model the generation of dense basal flows. This may have severe consequences for cases where the basal flow out-runs the dilute surge by large distances.

## Acknowledgements

We gratefully thank R. Stephen J. Sparks for providing invaluable advice during project inception and development. E.E.D. was supported by a UK Natural Environment Research Council PhD grant no: NER/S/A/2003/11201, and acknowledges code-development advice from Michael Hudson-Doyle. This work used the computational facilities of the Advanced Computing Research Centre and the Department of Mathematics, both at Bristol University.

## Footnotes

↵1 The TRBDF2 formulae is erroneously reported in Leveque (2002), and should be and for the stiff ordinary differential equation

*q*_{t}=*ψ*(*q*), where*Q*is the solution sought at time step*n*+1 and grid node*i*.

- Received July 29, 2010.
- Accepted October 20, 2010.

- This journal is © 2010 The Royal Society

## Appendix A. The Godunov numerical method

The first hyperbolic problem *q*_{t}+*f*(*q*)_{x}=0 is approximated to its linear form by finding the Jacobian *A*(*x*,*t*) for the system of equations and setting *q*_{t}+*A*(*x*,*t*)*q*_{x}=0 (Leveque 2002). A Riemann solver is then specified for any two cell-average states in space *Q*_{i−1} and *Q*_{i}, returning a set of *M*_{w} waves with speeds that satisfy . These waves are defined by , where is the *p*th eigenvector of *A*_{i−1/2}, and is a vector of coefficients defined by *α*^{p}=*l*^{p}(*q*_{r}−*q*_{l})=*l*^{p}(*Q*_{i}−*Q*_{i−1}), and are the corresponding left eigenvectors. This exact Riemann solver thus provides left- and right-going fluctuations to split the flux into and out of each cell via . The solution at each cell centre for each time step is then found via the first-order, upwind Godunov method (Leveque 2002): . After completion of this step, the source terms *q*_{t}=*ψ*(*q*) are usually included via standard ODE solvers, which are discussed further in Leveque (2002).

## Appendix B. The f-wave approach

The flux difference between cells is split into waves , travelling at speeds of , and directly decomposed as a linear combination of the eigenvectors of the Jacobian *A*_{i−1/2}, via (Bale *et al.* 2003)
B 1
where *Ψ*_{si−1/2} is the discretized spatial source term, the coefficients are defined by and *R* is the matrix of right eigenvectors of the Jacobian. The dilute current’s (2.24)–(2.26) Jacobian is
B 2
with the eigenvalues , *λ*_{2}=*u*, simplified to *λ*_{1,3}=*u*∓*c* and *λ*_{2}=*u*, and (Doyle 2007). The right-going eigenvectors are thus
B 3
and the left-going eigenvectors are (where *d*=2*ρ*_{a}−*β*)
B 4

Adopting a linearization of the governing equations, such that , leads to a set of left- and right-going fluctuations across the interface between the cells at *i*−1 and *i*,
B 5
These are then used in the upwind Godunov algorithm (Leveque 2002). The chosen average state is not critical to this f-wave method, thus we use the averages
B 6
and to define the eigenvalues and eigenvectors, using the values in the cell to the left (*l*=*i*−1) and right (*r*=*i*) of the cell interface. The spatial source terms (2.24)–(2.26) are then defined by
B 7
which is discretized at the cell edges to (Lowe 2005) and included in the f-wave decomposition (equation (B1)). The averages at the cell interface , and are used instead of the exact values at the cell centre to ensure that the appropriate equilibria of the underlying mathematical model are maintained by the numerical scheme (Hubbard & Garcia-Navarro 2000). Finally, the linear source terms *ψ*_{l}(*q*,*x*) (sedimentation, gravitational driving, basal drag) are incorporated via an additional fractional step, which we solve using the fourth-order Runge–Kutta method.

## Appendix C. Incorporating an entropy fix in the dilute-current numerical method

When the exact eigenvalues, describing the numerical-solution wave speeds, calculated in the left and right cells of the interface obey *λ*_{l}<0<*λ*_{r}, then a transonic rarefaction exists (Leveque 2002). This spreads partly to the left and right of the cell interface and, if incorrectly modelled, can lead to entropy-violating shocks in the numerical solutions. For the basal flow, we employ a Harten–Hyman entropy fix, which modifies the eigenvalues by partitioning them into the correct left- and right-going directions. The f-wave version of the Godunov method is defined directly from the waves themselves (equation (B5)). Thus, for the upper dilute current, a modification of this Harten–Hyman entropy fix is required. We employ a partitioning of the f-wave fluctuation directly, such that
C 1
and
C 2
using the Jacobians, eigenvalues and source terms defined in appendix B. This is an extension on the fix proposed by R. J. Leveque, available at http://www.amath.washington.edu/~claw/extensions/clawman/2d/examples/sphere/ (last viewed 11 June 2007). We have introduced the modifiers and to produce stable and accurate results. To prevent division by zero, and correctly partition the waves, the small number *η* has been introduced into expressions (C1) and (C2) to control the summation. If the average eigenvalue at the cell interface is , all waves are left-going, assuming that . If , all waves are right-going and thus if , the waves are split equally between left- and right-going fluctuations. However, if the exact eigenvalue calculated for this wave in the cell to the left of this interface also obeys *λ*_{l}<−*η* and to the right obeys *λ*_{r}>*η*, then the wave is a transonic rarefaction and the above entropy fix (C1) and (C2) must be applied.