Pulse propagation in gravity currents

P. A. Allen, a) R. M. Dorrell, O. G. Harlen, R. E. Thomas, and W. D. McCaffrey School of Computing, Faculty of Engineering and Physical Sciences, University of Leeds, Leeds LS2 9JT, UK Energy and Environment Institute, University of Hull, Hull HU6 7RX, UK School of Mathematics, Faculty of Engineering and Physical Sciences, University of Leeds, Leeds LS2 9JT, UK School of Earth and Environment, Faculty of Environment, University of Leeds, Leeds LS2 9JT, UK


I. INTRODUCTION
Gravity currents are flows driven by pressure gradients resulting from density differences. These may be the result of temperature, suspended sediment or salinity differences. Gravity currents form significant geophysical flows in atmospheric, terrestrial and subaqueous environments. Examples include, landslides, avalanches, turbidity currents, pyroclastic flows and lahars 1 . Pulses are a common feature in gravity currents and may result from flow instabilities, variable supply of dense material 2 , combining of flows from different sources, flow splitting and recombining 3 , or flow interactions with topography 4 . For example, failure mechanisms for landslides and similar events are varied 5 and can result in pulsed flows: an initial failure of an embankment or dam can create a steep main scarp as the supporting material slides away, which, in turn, can lead to a further 'retrogressive' failure. The process may repeat creating a significantly larger event comprised of many smaller pulses. Surges and pulses internal to gravity currents can have a significant impact on the hazards associated and flow properties when compared to a single release of the same volume. This is particularly significant for compositional gravity currents where the variations in velocity affect the deposition or erosion that can occur 6 .
A pyroclastic flow generated by the 1997 eruption of the Soufriére hills volcano contain three distinct major flow surges over its 25 minute duration 7 . Deposits from the first two major surges partially filled the main drainage channel which the pyroclastic flow flowed along. This caused the third a) scpaa@leeds.ac.uk to overspill and travel into a region considered to be at low risk. The release dynamics of these gravity currents and potential evolution downstream in multiple surges impacts the dynamics of the flow. The deposits left by previous surges, and information they contain, can be destroyed by subsequent surges making flow dynamics difficult to identify. Further, run-out length, inundation zones and hazards are affected by the internal dynamics.
Turbidity currents sub-marine gravity currents in which the density difference is caused by suspended sediment and have significant economic impact; they can travel with head speeds as high as 19 ms −18 and cause significant damage to sea floor equipment, including pipelines, oil rig moorings 9 , and seafloor telecommunication cables 10,11 . Turbidity currents are a key mechanism of sediment distribution from shallow to deep marine environments through the oceans and lead to create of sedimentary rocks (turbidites) which are linked to hydrocarbon reservoirs 12 . Deposits of seismically generated turbidity currents at the Cascadia margin, Washington USA record multiple currents that combined downstream at as many as seven confluences 13 . The separation time between the flows can be negligible or of the order a few hours. Further, experimental modeling has demonstrated that the signature of individual turbidity currents can be destroyed after the different events interact 14 .
Gravity currents have been extensively studied by theoretical and experimental approaches, often based around the idealized lock-release or lock-exchange problem, where dense fluid is released by the rapid removal of a gate, providing a well-controlled initial condition. This method provides a suitable means to create repeatable, spatially and temporally varying, fixed-volume currents that allow meaningful comparison to theoretical models [15][16][17][18][19] . Recently, surge effects have This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI:10.1063/1.5130576 2 been studied in the laboratory for turbulent saline flows using multiple lock-releases, in which a series of lock-boxes positioned behind the first were released at set times after the initial release 14,20,21 . The second current intruded into the first release and propagated towards the head of the current.
Dimensional arguments and simple integral models are able to predict the position of the head x N as a function of the initial conditions by assuming a dominant balance between buoyancy and either inertial or viscous forces 17,22,23 . By exploiting the large length L and shallow depth h of these flows, a first order asymptotic approximation in the aspect ratio δ = h/L followed by depth-averaging produces the hyperbolic shallowwater equations 24 . After an initial transience the flows generally enter a self-similar phase where the flow structure is independent of time and length scales. These similarity solutions capture some of the horizontal features of the flow. However, similarity solutions do not always exist. For example, in the axisymmetric spreading of a compositional gravity current, the tail of the current is self-similar but the head is time-dependent 25 .
The neglected vertical gradients in the shallow water model become significant at the head of the flow where turbulent drag and three-dimensional flow structures dissipate momentum. To capture this dissipation, shallow-water models are supplemented with an imposed dimensionless flow velocity at the head via a densimetric Froude number condition 26 . The densimetric Froude number is a dimensionless velocity of the current Fr = u/ √ g ′ h, where g ′ is the buoyancy adjusted gravity, u is the streamwise velocity and h is the flow depth, and is ratio between flow speed and infinitesimal long surface waves on the gravity current 27 . The buoyancy adjusted gravity expressed in terms of the current and ambient densities, ρ c and ρ a , respectively, and gravitational acceleration g is g ′ = g(ρ c − ρ a )/ρ c . Theoretical values of the Froude number can be determined through application of Bernoulli's theorem and a momentum balance of far upstream and far downstream of the head of the current in rectangular channels 28,29 or by a vorticity balance without a dissipation assumption 30 for Boussinesq flows. These differing approaches have been demonstrated to only differ by an assumption about the dissipation and can be reconciled within the same framework 31 . In addition, densimetric Froude number conditions have been calculated for a stratified ambient 32 and non-rectangular cross-sections 33 . The densimetric Froude number condition can also be determined experimentally to be 1.2 16 . For non-Boussinesq flows, where the density difference between the current and the ambient is large, the Froude number can be large and even tend to infinity. This corresponds to zero depth at the front of the flow 24 .
Self-similar solutions of the lock-release problem provide a valuable tool for numerical model validation and enable development of these models to explore more complicated scenarios. In the majority of applications, numerical integration is required to obtain solutions. Direct numerical simulation (DNS) of the Navier-Stokes equations has also been conducted and demonstrates that depth-averaged models are able to capture the principal physics of gravity current dynamics with significantly reduced computational time 9 . Shallow-water models can capture shocks (pulses) which appear as discontinuities in the solution. In reality, the viscosity of the flow smooths the discontinuity, but the shocks are still able to accurately capture the velocity of such disturbances 34 . Bonnecaze, Huppert, and Lister 34 employed a Lax-Wendroff finite-difference scheme to numerically integrate their shallow-water model that included sediment transport for a lock-release problem. The scheme accurately captured internal shocks which occurred behind the head downstream from the lock.
Hogg 27 employed the method of characteristics to solve the problem of a single-lock release flow. The shallow-water equations yield two families of characteristic curves along which the Riemann invariants α = u + √ g ′ h and β = u − √ g ′ h are conserved. Hogg 27 showed that the structure of the solution depended qualitatively on the densimetric Froude number, Fr. For Fr ≥ 2, characteristics leaving the back of the head never reached the back of the lock-box and an internal shock formed for Fr > 2. For Fr < 2, a structured solution exists in which the characteristic (x,t)-space is split into regions where both, one or neither of the two characteristics variables are constant. Hogg 27 analytically determined the boundary between these regions and the solution when at least one of the characteristic variables is constant. These models have not been extended to the multiple release case.
The goal of this paper is to extend the lock-release problem for the shallow-water equations to a double-release case, where a second equally sized lock-box is released subsequently. The second release will create a shock that will propagate through towards the head current. A Lax-Wendroff finite-difference scheme based on the implementation of Bonnecaze, Huppert, and Lister 34 is employed to solve the governing shallow-water equations. The characteristics are then computed from this solution afterwards, in order to describe the form of the solution in (x,t) space. For a double-release problem with identical lock-boxes, there is an additional parameter as well as the densimetric Froude number: the dimensionless release time t * re = t re √ g ′ h lock /l, where t re is the gate separation time, h lock is the lock depth and l is the lock length. The work presented here explores the (Fr,t re ) parameter space for Fr < 2. Simulations of the single-release case are compared to the analytical solution of Hogg 27 for validation. The double-release simulations reveal a variety of distinct regions in the (Fr,t re ) parameter space with qualitatively different behavior in the shock velocity. For t re → ∞, the two releases behave as two non-interacting events, whereas for t re ≤ 1 the flow is effectively a single discrete event of twice the volume. However, the two events interact, affecting pulse propagation, for intermediate t re -values and a range of qualitatively different solutions are obtained. Regions are separated by whether or not they are affected by the amount of material in the second lock-box and the path through the single release solution structure. The velocity of the shock has implications for the dynamics of the pulsed gravity currents flows discussed earlier and the assessment of their hazards.
The paper is structured as follows: The shallow water model is presented for both the single-and double-release configurations in §II; Results from our numerical model and This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI:10.1063/1.5130576
3 the parametric study are presented in §III; Implications are discussed in section §IV; Finally, we conclude in §V.

II. THEORETICAL MODELING
Consider a constant volume gravity current propagating over a fixed horizontal rigid surface in two spatial dimensions (x, z), where x and z are the horizontal and vertical coordinates, respectively, with the time from the first release defined as t, Fig. 1. The current has density ρ c and the deep ambient has density ρ a . Thus, the buoyancy-adjusted gravity for the current may be expressed as g ′ = (ρ c − ρ a )g/ρ c , where g = 9.81 ms −2 . In the lock-gate configuration, the flow quickly reaches a state where the height h(x,t) is much smaller than the length of the current x N (t). Therefore, in considering the horizontal momentum of the flow we can average over the depth and assume purely hydrostatic pressure. Further, the flows are assumed inviscid (with no basal drag), entrainment is negligible, and the ambient is quiescent and infinitely deep. These assumptions allow us to apply the simplified depth-averaged shallow-water equations 35 , where (3) are respectively the depth-integrated momentum per unit mass and the depth of the flow expressed in terms of the horizontal velocity v(x, z,t) 36 . The shallow water equations are first order approximations in terms of the aspect ratio between the depth and length of the current δ and contain no source terms, i.e. drag and entrainment. However, at the head of the current x = x N (t) the dissipation is accounted for via a densimetric Froude number condition 29 . Further, a dynamic boundary condition is imposed at the head where Fr is a constant and subscript N denotes a value at the head. Both lock-boxes are assumed to be of the same length l and filled to a depth of h lock . Initially, no flux boundary conditions are imposed at the back of both lock-boxes x = 0, l. After the second gate is released at t = t re , the no flux condition at x = l is removed. From the momentum equation (2), no flux is equivalent to ∂ h/∂ x = 0 and so Similarly, the initial conditions, Fig. 1, are defined as The mass and momentum conservation equations (1) where * denotes a dimensionless variable. Similarly, the dimensionless boundary conditions (5) and (4) are: and initial conditions (6) become A. Analysis of characteristics: single release case This section will discuss the behavior of the flow for t < t re , which is equivalent to the single release solution of Hogg 27 . From this point the * are neglected from the dimensionless variables for brevity, unless stated otherwise. The system of equations (7)(8) can be transformed into its characteristic form 24 by changing to characteristic variables Thus, α and β , the Riemann invariants, are constant along characteristics curves with gradients in (x,t)-space of u + c and u − c, respectively. The gradients u ± c are the eigenvalues of the system of equations (7) & (8) and, provided the flow depth is non-zero, are real and distinct everywhere. Thus, the system is hyperbolic and the method of characteristics may be applied. The values of the characteristic variables can be determined from boundary or initial conditions that the characteristics pass through. This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. Initially, α = 2 and β = −2 on all positive and negative characteristics that originate from 0 < x < 2, t = 0, see Fig.  2. Thus, with the imposed densimetric Froude number condition u N = Frc N , positive characteristics collide with the head at a finite time provided Fr → ∞ (in the limit Fr → ∞, the head corresponds to the leading characteristic). While positive characteristics arrive at the head with α = 2, the velocity and wave speed at the head are constant and take value

PLEASE CITE THIS ARTICLE AS
Thus, negative characteristics emanating from the back of head have constant An expansion fan of negative characteristics emanates from (2,0) and connects the two regions where β is constant, see Fig. 2b. These negative characteristics correspond to straight lines through the origin and satisfy The curves x ref and x fan collide at t = (2 + Fr) 3/2 / √ 8. Beyond this point, x fan enters a region of varying α and therefore has a non-constant gradient. In contrast, x ref enters a region of constant β and thus has constant gradient until reaching the head. Hogg 27 gives the x ref characteristic between the back wall and the head as because of the no flux boundary condition (9), while in regions U 2n+1 , S 2n+2 and U 2n+2 because of the densimetric Froude number condition (10). As n → ∞, α → 0, and β → 0, and thus u → 0 and c → 0. When Fr ∼ 0 , λ ∼ 1 and therefore the flows will interact for a large number of reflections and hence longer release times. For Fr ∼ 2 there is minimal interaction between the events. Critically, x ref and x fan partition the single-release solution into three distinct regions for any fixed t where the behavior of α This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI:10.1063/1.5130576 6 and β are qualitatively different.

B. Extension to double release
From the structure described in section II A we can determine the nature of the solution when the second gate is released at t = t re . If t re ≤ 1 then trivially the solution behaves identically to a single-release of lock-box length 2, because the backwards traveling disturbance has insufficient time to reach x = 1. For t re > 1, a shock is created where positive characteristics in x < 1, having α = 2, collide with positive characteristics from x ≥ 1 having α < 2. Depending on the release time, t re , the shock is released into a region of constant depth (uniform region) or varying depth and ve-locity (complex region) and so the relative position of the two curves, x ref and x fan , determine the initial motion of the shock. Imposing that mass and momentum fluxes are conserved across the shock, the shock velocityẋ s can be obtained from the Rankine-Hugoniot conditions of the shallow water equations, because they are in a conservative form. In terms of the characteristic variables ahead of and behind the shock, α + , β + , α − , β − , respectively, the shock velocityẋ s iṡ from mass conservation (7) anḋ from momentum conservation (8). Initially α − = 2 as the positive characteristics come from unperturbed fluid. Shocks that propagate into a uniform region have both α + and β + constant and thus, by the Rankine-Hugoniot conditions (21) & (22), bothẋ s and β − are also constant. In simple and complex regions, the shock will accelerate or decelerate and values of β − will vary. In uniform regions adjacent to x = 1, the boundary condition u(1) = 0 (9) implies that u(x) = 0 throughout the uniform region. Thus, and whilst α − = 2, the problem replicates the wet dam break 24 and the shock velocity can be calculated explicitly throughout the uniform region. The Rankine-Hugoniot conditions (21) & (22) provide an implicit relation for the constant shock velocity,ẋ s , for uniform regions adjacent to the head. For a shock of positive velocity,ẋ s , causality implies that positive characteristics cannot be emitted by the shock. Thus, the shock represents the furthest point in the domain that has been affected by the release of the second gate.
Similar to the line x ref (t), an additional line x fin (t) is introduced for the second release. The first branch tracks the backwards propagating disturbance of the second release, i.e. the fastest negative characteristic from (1,t re ). On this characteristic β = −2, and positive characteristics intersecting it arrive from unperturbed fluid and therefore α = 2. Thus and so At t = t re + 1, the fluid at the back of the second lock starts to be affected by the gate release and beyond this time α < 2 at x = 0. The last α = 2 characteristic leaves x = 0 at t = t re + 1, which is denoted as the continuation of the line x fin . The second branch of the curve x fin defines the part of the solution affected by the finite length of the second lock-box. If this characteristic intersects the shock, then α − < 2 thereafter. Both β − and the shock velocityẋ s are constant when the shock propagates through a uniform region. Thus, for a shock propagating in a uniform region, a region of constant β is created, which in turn creates another region of constant α < 2 upon reaching the back of the lock-box. The structure of the characteristic space is displayed for a shock released into a uniform region, Fig. 4, and a complex region, Fig. 5 . For shocks released into either uniform or complex regions negative characteristics will have gradients greater than −1 (14) and another expansion fan of negative characteristics must exist at (1,t re ). A shock released into a uniform region will initially travel at a constant speed. Further, the flow depth and velocity will be constant either side of the shock. This will hold until the shock intersects one of the three curves x ref , x fan , or x fin , with each possibility leading to a different structure behind the shock. The example drawn in Fig. 4  The shallow-water equations (7) & (8) coupled with boundary (9) & (10) and initial (11) conditions are solved using the This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. At the point of release, an additional expansion fan is bounded between x fin (· · · ) and the slowest backwards traveling characteristic from the second release (· · −). Online version in color.

PLEASE CITE THIS ARTICLE AS
method of Bonnecaze, Huppert, and Lister 34 . For stability the initial depth, h N , and velocity, u N , at the head were set to slumping phase values; Further, the shock initiates at (x,t) = (1,t re ), where u = 0, and thus the positive and negative characteristics ahead of the shock takes values α + = 2 h(1,t re ) and β + = −2 h(1,t re ).
Together with α − = 2, the Rankine-Hugoniot conditions (21) & (22) were used to determine the initial shock velocity and depth,ẋ s (t re ), and these were imposed at the node coinciding with the shock for the first time step after release. The shallow water equations are remapped to the unit interval using the change of variables (ζ , τ) = x/x N (t),t). This removes the moving boundary condition simplifying the application of the Froude number condition. For full details see appendix A of Bonnecaze, Huppert, and Lister 34 . The transformed equations are solved using a Lax-Wendroff finite different scheme. Numerical integration with a upwind finitedifference scheme is used to determine the positive characteristic that arrives at the head, ζ = 1, at the next time step and provide, together with the Froude number condition, a second equation for u and c at the head. This enabled us to determine suitable values used as boundary conditions at the head for the next time step. Model validation is performed by comparing the numerical solution at t = 1 and the curve x ref (t) with the exact solution given by Hogg 27 . This was chosen, because of the significance of the three lines x ref (t), x fan (t) and x fin (t) to assess the ability of the simulations to capture the distinct regions where behavior changes. This validation is presented in the Appendix This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. and x fan (· · · −), the head of the current x N (-) and the shock x s (--). At the point of release, an additional expansion fan is bounded between x fin (· · · ) and the slowest backwards traveling characteristic from the second release (· · −). Online version in color.

III. RESULTS AND DISCUSSION
In this section, the range of possible solutions from the shallow-water model are classified in terms of qualitative differences in the shock velocity throughout its motion towards the head of the current. The generic structure of the (Fr,t re ) phase space is presented from a parametric study, and the solution types are presented.

A. Shock Evolution
The numerical solutions are first distinguished by the region the shock is released within, C 1 , U 2 , C 3 etc. If the shock at x s (t) is released into a uniform region, the shock velocity remains constant until it enters the simple wave region. If the shock is released into a complex region, its velocity varies from the outset. The characteristics that bound the region of varying or constant α and β , x fan (t) and If the shock is released into a complex region, three distinct paths through (x,t)-space may occur: i) The shock first intersects x fan and enters a simple wave region in which β is constant. The shock then intersects x ref after it has reflected off the front, entering another complex region; ii) The shock intersects x fan first, but intersects x ref before it has reflected off the head. Thus, the shock enters a uniform region until colliding with the head; iii) The shock intersects x ref , entering a simple region within which α is constant. It then intersects This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI:10.1063/1.5130576
x fan , entering a uniform region. The three possibilities for C 1 are displayed in Fig. 6a. A similar distinction is drawn for the uniform cases, for example U 2 in Fig. 6b: i) β varies within the simple wave region before becoming constant again; ii) β starts varying first followed by α; iii) α starts varying followed by β . Further, if the shock is affected by the finite length of the domain lock-box, i.e. x fin intersects with the shock before it reaches the head of the flow, then this is equivalent to α − < 2 beyond this point. The first time the shock and x fin coincide is defined as t fin , such that x fin (t fin ) = x s (t fin ). If the shock is affected by the finite length of the lock-box, the case is labeled 'F', otherwise 'N'. For example U 2 Fi is a shock that is released into the U 2 region of type i that is affected by the finite length of the lock-box.
Numerical solutions for the U 2 N shocks are displayed in Fig. 7. For case i, Fig. 7a&b, β varies upon entering the S 3 region and the shock accelerates at approximately a constant rate until entering U 3 , where it returns to a constant but higher velocity until colliding with the head of the current. In contrast the shock velocity increases throughout in cases ii, Fig.  7c&d, and iii, Fig. 7e&f, as first one characteristic starts varying and then the other (β first for case ii and α first for case iii). The shock velocity at t fin is largest for case iii, with case i being the slowest.
For larger densimetric Froude numbers, the shock has further to travel before reaching the front and therefore is more likely to be affected by the finite length of the lock-box. The U 2 F shocks are qualitatively similar to the U 2 N cases, Fig.  8, until x fin collides with the shock at t = t fin , after which the velocity decreases and the maximum shock velocity occurs before the head. Case U 2 Fi does not exist. i.e., the shock can only be affected by the finite length of the lock-box if it enters a region where α is varying. Further, our simulations reveal that this case U 4 Fi does not exist. Shocks released into complex regions immediately increase in velocity as both α + is decreasing and β + is increasing. Example C 3 N cases are presented in Fig. 9. As expected, the acceleration of the shock decreases after t fan for case i, because β takes a constant value until t ref , after which β starts decreasing again and the acceleration of the shock increases. For cases ii and iii, the acceleration of the shock decreases when it enters the S 3 region before becoming zero when entering U 3 until reaching the head of the current. For C 3 and subsequent complex and uniform cases, the shock may feel the affect of the finite length of the lock-box before colliding with x ref or x fan . The acceleration of the shock still decreases after it is affected by the finite length of the domain, but the range of possibilities becomes more complex, Fig. 10. This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. run at these parameter values in order to resolve the parameter space accurately. The distinct regions of the parameter space are presented in Fig. 12, with the i, ii, and iii distinctions only shown for cases C 1 , U 2 and C 3 . As Fr → 0, the boundaries collapse onto odd integers and the regions of complex shocks tend to zero. This is to be expected, as x ref and x fan coincide for Fr = 0, and, in fact, u = 0 everywhere so positive and negative characteristics have gradient 1 or -1, respectively. Cases C 1 ii and C 1 iii are the only two possible cases where the shock collides with x ref with α = 2. Our simulations reveal that neither of these cases exist.

PLEASE CITE THIS ARTICLE AS
The classification of each case determines exactly which three regions the shock travels through. the Rankine-Hugoniot conditions (21) & (22) ensure that the shock velocity is constant in uniform regions where the shock is not affected by the finite length of the lock-box α − = −2. For shocks in simple wave regions, one of α + or β + will take constant value given in equation (19) or (20), respectively, while the other will vary monotonically between two constant values dependent on Fr (from the neighboring uniform regions). In a complex region, the values of α + and β + both vary between the constant values given in equations (19) & (20).
Solving the Rankine-Hugoniot conditions (21) & (22) numerically for a fixed α − = 2 and varying Fr, α + and β + enables us to explore the regions of the flow where the shock x s may accelerating or decelerating. Studying the regions up to S 7 and C 6 revealed that the shock accelerates for all possible parameter values in regions of fixed α + (S 3 , S 5 and S 7 ) and for all values of α + above a critical value of the densimetric Froude number Fr c in regions of fixed β + (S 2 , S 4 and S 6 ). Similarly one variable was fixed and the other varied to determine Fr c -values for complex regions. This analysis re-This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. Online version in color. This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI:10.1063/1.5130576
vealed that apart from case C 3 Ni in the C 4 region and C 3 Niii in the C 3 region, all the Fr c -values lie in the F region of the parameter space and therefore, all shocks in the N region up to case U 6 have non-decreasing velocities,ẋ s . Sampling the C 3 Niii case did not reveal any decelerating shocks. Decelerating shocks were found in C 3 Ni near the N-F boundary and just before the shock reaches the head. However, these were also observed for other cases near the N-F boundary. As the numerical method is not exact, the maximum shock velocity may arise just before the reflection of the backwards traveling disturbance x fin reaches the shock, Fig. 10b,d&f. The velocity maximum occurs just before t fin . This is believed to be a numerical artifact from the dissipative nature of the Lax-Wendroff scheme and that N shock do no necessarily decelerate.
The velocity of the flow behind the shock u − = (α − − β − )/2 is related to the velocity of the shock through the Rankine-Hugoniot conditions (21) & (22). The simulations reveal that the acceleration of the flow behind the shock has the same sign as the shock acceleration, except when the shock is propagating through S 2n regions of the flow, i.e. regions where β + takes a constant value. In these regions the acceleration has the opposite polarity. Critically, this means that an increasing shock speed is equivalent to an increasing fluid velocity behind the shock for all other regions in the single release flow.
Depth h and velocity u = m/h are compared for cases U 2 Ni and C 3 Ni in Fig. 11. The position of the curves x ref ( ) and x fan ( ) are displayed along the x-axis to highlight the boundaries between uniform, simple, and complex regions. After the second release the position of the shock ( ) and x ref (•) are also shown. Both cases are for the same Fr, so the initial dynamics are identical. After the second release, the flow depth is deeper, but slower, behind the shock for the U 2 Na case when compared with the C 3 Ni case. When the shock reaches the head of the flow the expected constant depth and velocity is observed for the U 2 Na case. Whereas for the C 3 Ni case, the velocity increases to a maximum at the head of the flow. This maximum is larger than the maximum velocity for the U 2 Na case.

C. Momentum of the Head
Both the destructive potential, i.e. the amount of damage it can cause, of the gravity current and its run-out length can be affected by the fluid momentum at the head of the current m(x N ,t). In Fig. 13 we compare the momentum at the time when the shock collides with the head at the same time, m(x N ,t col ), for four different Froude numbers Fr against pulse separation times t re . Also shown on this plot is the corresponding momentum for a single release of twice the size at (x,t) = (x N ,t col ). For a fixed Fr, the momentum is significantly lower for a single release than the corresponding double release. Although both cases contain the same amount of material, the depth at the head of the flow starts decreasing later for the single release and therefore the dissipation is larger. Further, the dissipation is significantly larger at the head of the flow than at the shock.
In contrast to the single release, where increasing the densimetric Froude number decreases the dissipation at the head, plotted against separation time the momentum in the corresponding double release is lower for higher densimetric Froude numbers. This is a consequence of the head moving faster for higher Froude numbers. If instead the momentum at t = t col is plotted against the position of the head x N (t col ), Fig.  13b the expected trend is observed with higher Froude numbers being less dissipative. Critically, although the distinction is less with larger Froude numbers Fr, a single release event is more dissipative than a double release.

IV. IMPLICATIONS
We can now reflect on the implications of our model findings on the geophysical events discussed in the introduction.For compositional flows, the particle size distribution is the dominant control on when erosion and deposition rates are in balance 6 . For dilute flows, this equilibrium point defines a particle concentration above or below which the flow is erosional or depositional, respectively. The models discussed in Dorrell et al. 6 are increasing functions of the shear velocity at the bed u * , which can be empirically modeled as proportion to, or related to a positive power of the depth-averaged velocity u = m/h [37][38][39] . Thus, in regions where the flow speed is increasing it is more net erosional and conversely for regions where the flow speed is decreasing it is more net depositional. When the particle concentration is large, for example lahars and landslides, the particles are in constant contact and suspended by matrix strength rather than the fluid. Therefore there is no simple relationship between erosion and deposition rates and the acceleration of the flow for these flows. Instead erosion and deposition are controlled by a stress balance at the bed 40 .
The hazard assessments conducted for the Soufriére Hills Volcano, Montserrat assume a continuous steady release of material from the lava dome during the 1997 eruption creating a pyroclastic surge flow 7 . However, at the 1997 eruption the dome failed retrogressively producing a continuous release with three distinct peaks in volume flux identified as distinct pulses. Sedimentation from pyroclastic surges can create a complex flow structure with a dense, basal pyroclastic flow component with high-particle concentration, which is overlain by and an upper, less dense pyroclastic surge. In the 1997 event, 1.4 km downstream from source the third pyroclastic flow associated with the third pulse overspilled the drainage channel and went on to hit the villages of Streatham and Windy Hill. Despite the bend coinciding with a constriction in the channel cross-sectional area these villages were not considered at risk. However, the previous two pulses had left significant deposits thus increasing the risk of surge detachment and overspilling.
The pulsed nature of this flow affected its run-out and the inundation zone of this flow when compared to the a continuous release of the same size 7 . The shallow-water model discussed in this paper suggests that the separation time between This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. the pulses may have a crucial affect on whether the pulses were erosional or depositional at a particular location. Pulses that are close together (case N) are non-decelerating and may entrain more material downstream increasing the depth of the channel or increase the velocity/momentum of the flow. For large separation times the pulses may transition from erosional to depositional downstream filling in the channel and making it more susceptible to future flows overspilling.

V. CONCLUSIONS
We have explored the affect of pulses on gravity current propagation using an extension of the shallow-water model for the single lock-release case studied by Hogg 27 . The range of solutions are classified in terms of two parameters: the Froude number at the head of the current, Fr, and a dimensionless pulse separation time, t re . For t re ≤ 1 the problem is identical to a single release, whilst the limits t re → ∞ and t re → 1 + This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.  Variations in pulse velocity affect the rate of energy dissipa-tion, and thus of the energy transferred through to the head of the current, which may enable the flow to transition from laminar to turbulent behavior. For pulse-prone, Boussines compositional flows such as pyroclastic flows, the dynamics of the flow depend on dynamics of the release and the changes in flow velocity may have implications for hazard prediction models, which sometimes neglect the release dynamics and the subsequent pulses created.

PLEASE CITE THIS ARTICLE AS
Verification of the above scheme was conducted by comparing simulations of a single lock-release problem to the solution presented in Hogg 27 . Depth and velocity profiles at t = 1 are compared for fixed ∆t and varying ∆x and vice versa, Fig. 14. The error of the variable E f is calculated via the ℓ 1 -norm of the variable | f a − f n |, (A.1) where f = h or f = u ≡ m/h and N int is the number of nodes over the averaging interval. The numerical, subscript n, and analytical, subscript a, solutions are interpolated onto an equally spaced grid with N int = 10 4 . For fixed ∆x the error quickly converges to the spatial error and therefore a numerical solution ∆t = 10 −6 and ∆x = 0.005 is used instead of the analytical solution for comparison. After t = 1, complex regions start appearing in the characteristic space and a depth and velocity profile across the length of the current cannot be explicitly written down everywhere. However, the curve x ref has been expressed in closed form in equation (18) and this expression is compared to a positive characteristic emanating from (1,1) calculated from the numerical solution. An excellent agreement between the analytical and numerical expression for x ref is observed, Fig. 15. Although the Lax-Wendroff finite-difference scheme is formally second order accurate, the model verification suggests that it is first order accurate in both time and space. This is to be expected as the analytic solution of the dam break does not have a continuous derivative everywhere. Although not displayed here, for lower resolutions the contribution to the error is largest at x = 0 and x = 1.
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.