Geomorphology 53 (2003) 25 – 44 www.elsevier.com/locate/geomorph
A smoothed-particle hydrodynamic automaton of landform degradation by overland flow M. Bursik a,*, B. Martı´nez-Hackert a, H. Delgado b, A. Gonzalez-Huesca b b
a Department of Geology, University at Buffalo, SUNY, Buffalo, NY 14260, USA Instituto de Geofı´sica, Universidad Nacional Auto´noma de Mexico, Mexico D.F., Mexico
Received 13 February 2001; received in revised form 4 September 2001; accepted 4 October 2001
Abstract A novel model for calculating hillslope degradation by overland flow is based on a coupling of the cellular automaton and smoothed-particle hydrodynamics theories. Simulated debris-laden floods were initiated on digital elevation models (DEMs). The flows were created by solving equations for diffusion routing, and erosion and deposition. We separated advection and diffusion terms using operator splitting, then solved the advection term via automaton and the diffusion term via smoothed-particle hydrodynamics. The technique was tested on geometrically simple synthetic slopes, against data for erosion and deposition by debris-laden flows in a gully developing on a scoria cone (Arizona), and against data on a recent moderate-sized debris-laden flood generated by possible partial glacier melting at Popocate´petl volcano (Mexico). The flow and erosion/deposition patterns on the simple synthetic slopes appear intuitively reasonable. On the synthetic scoria-cone slope, the model reproduced some of the observed spatial patterns of event-scale erosion and deposition. Reasonable flow speeds, runup heights, and runout lengths were produced for the debris-laden flood at Popocate´petl. The coupled cellular automaton – smoothed-particle hydrodynamics technique represents some advance in the modeling of landform degradation and channel development by cellular automata (CA). Use of the technique may aid our understanding of the link between hillslope-scale and landform-scale degradation patterns. D 2002 Elsevier Science B.V. All rights reserved. Keywords: Landform degradation; Channel development; Overland flow; Debris flow; Simulation; Smoothed-particle hydrodynamics; Cellular automaton; Popocate´petl; San Francisco Volcanic Field; Arizona; Mexico
1. Introduction Cellular automata (CA) have proven capable of simulating some of the features of landform evolution. The ‘‘precipiton’’ model of Chase (1992) included an implementation of both degradation by linear-diffusive * Corresponding author. E-mail address:
[email protected] (M. Bursik).
processes, such as rain-splash, and degradation by overland flow. Overland flow was implemented by allowing sediment movement between adjacent, landform-scale, hillslope cells. Chase was able to show that some of the fractal features of entire landforms, developed over hundreds of thousands to millions of years, could be explained with simple cellular rules. Because the drainage length was restricted to extend only between adjacent cells, the model was not scalable.
0169-555X/02/$ - see front matter D 2002 Elsevier Science B.V. All rights reserved. doi:10.1016/S0169-555X(02)00346-X
26
M. Bursik et al. / Geomorphology 53 (2003) 25–44
Favis-Mortlock (1998) used a cellular automaton to characterize the initiation and evolution of rill systems at the hillslope scale. The model of Favis-Mortlock (to a degree scale independent) simulated some of the network and growth features of rill systems and intrarill geometry, avoiding some of the requirements of standard hillslope models—such as the Water Erosion Prediction Project (WEPP) and the European Soil Erosion Model (EUROSEM)—to specify rill geometry and flow partitioning. Between the landform and the rill scales, and long and short timespans, lies a ‘‘noman’s land.’’ Is it possible to use CA to bridge this gap and help us to understand how local, short-term processes relate to long-term landform modification? A cellular automaton consists of an array, the elements or cells of which can assume a finite set of values. These values evolve through time by a set of fixed algebraic rules applied to the old value and the values of nearby cells. In geomorphology, CA have been used to simulate phenomena as diverse as the development of braided stream patterns (Murray and Paola, 1994), the large-scale degradation of landforms (Chase, 1992), and the development of wind-ripple stratigraphy (Forrest and Haff, 1992). In replicating some of the features of natural systems, CA encapsulate the physics of a phenomenon within algebraic relationships rather than a system of differential equations. In classical CA as outlined by Wolfram (1983), cells were acted upon simultaneously by one set of rules. Chase (1992) applied the rules to cells sequentially and randomly to develop the complex topography seen in large-scale landforms. Murray and Paola (1994) used the concept of tracking two linked variables through time, in their case topography and river discharge, i.e., essentially running two linked CA simultaneously. The present contribution builds on these earlier works by applying the rules sequentially to cells containing a nonzero value of flow depth, and by tracking both topography and flow depth as they develop. Typically, finite difference or finite element techniques have been used in computational fluid dynamics to model open-channel flow of river floods and debris flows (Chaudhry, 1993; Molinaro and Natale, 1994; Chen, 1997). Both techniques rely on an array of cells, at each of which the values of the changing variables are calculated through time (particle-in-cell or PIC technique). In contrast to PIC techniques,
smoothed-particle hydrodynamics (SPH) relies on a group of fluid ‘‘particles’’ or ‘‘parcels’’ that are not located at fixed grid points, i.e., they flow from one position to any other in the framework (Monaghan, 1992). As a result, the values of parameters associated with the cells occupied by fluid parcels alone are calculated, and the rest of the grid is left unchanged. The location, velocity, momentum, and energy of each particle are evaluated through use of a smoothing kernel, which allows the properties of each particle to rely on those of other particles, with nearer particles having a greater influence than more distal particles. The goal of the present work is to develop an automaton model of overland flow, channel formation, and channel evolution on hillslopes, which might be able to bridge between the models of Chase (1992) and Favis-Mortlock (1998). With respect to the Chase model, relaxation of restrictions on the relationship between drainage length and landform size should make the model scalable. With respect to the FavisMortlock model, accommodating variable flow width, divergent flow, variable velocity, and deposition, and running in the time domain (water packets interact) should advance model capabilities. To accomplish these advances, the principles of smoothed-particle hydrodynamics are coupled with cellular automata in a formulation of the shallow water equations. The CA/ SPH model is tested first using simplified topographic systems, such as ramps and diverging channels. Next, we undertake a qualitative test of the model against field examples of a gully, the morphology of which has been monitored for 5 years, and of a debris-laden flood on Popocate´petl volcano in Mexico. Volcanoes represent ideal environments to test such models because of the steep slopes and the high frequency of degradational events on their flanks. The results represent the first application of smoothed-particle hydrodynamics to the shallow-water formulation, as well as the first coupling of automata and smoothedparticle hydrodynamics.
2. Analytical model 2.1. Equations of motion Flood routing is the calculation of depth and velocity as a fluid propagates downslope, usually within a
M. Bursik et al. / Geomorphology 53 (2003) 25–44
channel. Most routing is based on solution of ‘‘shallow-water’’ or St. Venant equations, depth-averaged forms of the continuity and Navier –Stokes equations. A variety of formulations of the shallow-water equations for different natural and engineered channel geometries and material properties are discussed in several recent texts and compilations (Chaudhry, 1993; Molinaro and Natale, 1994; Singh, 1996; Chen, 1997). For the current situation, the formulation captures some of the characteristics of a sediment laden flood that approaches pure-water at one extreme, and a grain flow at the other extreme, as occurs in debris flows. The equations of motion can be written as Bh Bhu þ ¼0 Bt Bs
ð1Þ
for continuity, and Bu Bu u þu ¼ gS0 g ð1 hw =hÞcos#tanub Bt Bs AuA u B 4vm 2 egcos# h Bs ðkðh hw Þ þ hw Þ ð2Þ for momentum (Savage and Hutter, 1989; Chaudhry, 1993; Iverson, 1997). In Eqs. (1) and (2), t is time, s is the local downstream spatial coordinate, h is flow depth, u is flow speed, S0 is surface slope (S0 = tan # = Bz / Bs c sin #), # is slope angle, hw is depth of water alone, ub is the internal friction angle of the mixture of clasts, v is the volume fraction of water, with m the kinematic viscosity of the fluid plus fine particles, e the flood-wave aspect ratio, and k the ‘‘earth pressure coefficient.’’ One simplification of the present work is that h = hs + hw, where hs is the depth of particles. This simplification is strictly true for mixtures in which there is no density contrast between phases and no support from the framework of particles. Thus, the simplification poorly characterizes debris flows, and an improvement of the present work should include a relaxation of this assumption. The equations are formulated for a channel—or in our case a water parcel— of constant width. Lateral spreading of water on open slopes is accomplished by transfer of water among parcels, not creation of new parcels. Two resistance terms occur in the momentum equation—one static (second term on the right-hand side) and the other viscous (third term on the right-hand side).
27
By rearrangement and dividing through by g, Eq. (2) can be put in a form analogous to the St. Venant equation: u u ð1 hw =hÞcos#tanub Sf ¼ 4vm 2 þ h AuA B 1 Bu u Bu ¼ S0 ecos# ðkðh hw Þþ hw Þ Bs g Bt g Bs ð3Þ where Sf is the friction slope. In the end-member case of dilute flow, a Manning fluid-resistance term might be more appropriate, and hw would everywhere be equivalent to h. The substitution yields the St. Venant equation of momentum: Sf ¼
n2m u2 Bh 1 Bu u Bu ¼ S0 ecos# Bs g Bt g Bs h4=3
ð4Þ
where nm is the Manning roughness parameter. The continuity and momentum equations can be combined into a single expression, in which the local (Bu/Bt) and convective (uBu/Bs) accelerations are assumed negligible, to solve for the flux per unit width, uh. The results of the simplification will be equations of the diffusion-wave type (Singh, 1996). Diffusion-wave equations adequately characterize unsteady flows in long reaches. They do not characterize conditions at channel transitions, such as hydraulic jumps. From Eq. (3): gh3 u uh ¼ ð1 hw =hÞcos#tanub S0 AuA 4vm Bhw Bh þ ecos# ð1 kÞ k Bs Bs
ð5Þ
Putting the expression for uh from Eq. (5) into Eq. (1), we get an expression of the form (ignoring higherorder terms): Bh Bh B 2 hw B2 h þ 3jS0 h2 jeh3 ð1 kÞ k Bt Bs Bs2 Bs2 u B tanub ½h2 ðh hw Þ ¼j ð6Þ AuA Bs where j = 4g/vm. The term on the right-hand side describes the resisting effects of static frictional forces within the flow.
28
M. Bursik et al. / Geomorphology 53 (2003) 25–44
2.1.1. Asymptotic behavior The characteristics of Eq. (6) can be explored by considering the end-member cases of dilute and highly concentrated flow. For the flow of dilute mixtures, Eq. (6) reduces to Bh Bh B2 h þ 3jS0 h2 jeh3 2 ¼ 0 Bt Bs Bs
stream floods, all showing a steep rise with a peak discharge (or depth) value near the beginning of the event and a gradually decreasing depth following the peak value (Pierson, 1995). A variety of discretized shapes satisfying these general rules were used as initial conditions for simulations.
ð7Þ
(The approximation used here is similar to that used to derive Eq. (4), that h f hw.) Eq. (7) is a form of diffusion routing equation. The diffusion wave described in Eq. (7) has a celerity, c, of 3jS0h2 and a diffusivity, D, of jeh3, where the celerity is the speed of the diffusion wave (Singh, 1996) and the diffusivity describes the rate at which the wave attenuates. As the flood elongates and e ! 0, the third term of Eq. (7) approaches zero and the equation collapses to that for a kinematic wave, as it should. Near the runout limit of a concentrated flow, the water depth, hw decreases to values below hs. In this case, a crust of coarse particles is thought to develop at the flow perimeter, hindering propagation (Major, 1997). For this situation we have hwbh or hw f 0, and Eq. (6) becomes Bh u Bh B2 h þ 3jh2 S0 tanub jeh3 k 2 ¼ 0 Bt AuA Bs Bs ð8Þ The celerity is c = 3jh2(S0 (u/AuA)tan ub) wherein (u/AuA)tan ub acts as a friction slope to slow the flow. Once the slope angle drops to this value, the flow no longer propagates and the sediments are deposited en masse.
3. Solution method The method used to solve the diffusion routing equation and to propagate synthetic hydrographs downstream couples a cellular automaton with smoothed-particle hydrodynamics to track the movement of discrete parcels of water plus sediment on a digital elevation model (DEM), consisting of an evenly spaced grid of elevation values. Despite tracking discrete flood parcels, smoothed-particle hydrodynamics provides a smoothing or averaging technique to account for the interaction between each fluid parcel and its neighbors and to approximate continuum behavior. The cellular method is based on principles similar to those described by Murray and Paola (1994) for the evolution of channel systems in streams heavily laden with sediment (braided streams). In the current formulation, we assume that water is generated at the source in a relatively instantaneous event and that there is sufficient water to maintain some sediment in motion as long as there is a sufficient downslope gradient. Eqs. (7) and (8) are of the form Bh ¼ L1 h þ L2 h Bt
ð11Þ
The boundary condition for a diffusion-wave equation consists of a hydrograph specified at both t = 0 and s = 0 as localized functions of both variables, which become vanishingly small at a finite distance from the origin:
where L1 and L2 are spatial linear operators. The technique of operator splitting is used to divide the spatial operators in Eqs. (7) and (8) from one another, allowing the advection and diffusion of water to be treated separately in arriving at the final updated values for position and depth (Press et al., 1986).
hðt; 0Þ ¼ HðtÞ
ð9Þ
3.1. Cellular automaton
hð0; sÞ ¼ GðsÞ
ð10Þ
2.2. Boundary conditions
Initial and near-source hydrographs from a number of debris-laden floods are similar to those for normal
The automaton handles the advection operator, moving each fluid particle from cell to cell in the framework. In the automaton, each water parcel is
M. Bursik et al. / Geomorphology 53 (2003) 25–44
moved in turn from its current cell into the adjoining cell with the lowest value of elevation plus water depth (Fig. 1). If more than one cell has the same lowest elevation plus water depth, then the parcel flows into one of these lower cells at random. This effectively prohibits two water parcels from occupying the same cell simultaneously. The cellular rule is
29
expressed in terms of the successive positions, !rt, of a water particle: !rtþ1 ¼!rt þ m!ct Dt i i
ð12Þ
where i is the index for each fluid parcel, m is an index of time steps to determine whether the displacement is
Fig. 1. Schematic diagram showing the nature of the computational grid on a landform (a) and the movement of a water parcel within the cellular grid (b). Bold arrow for rit shows a potential path that the fluid parcel could take; medium dashed line shows actual path taken, including step from rit to rit + 1.
30
M. Bursik et al. / Geomorphology 53 (2003) 25–44
zero or one grid cell, and ! c obtains its magnitude from Eq. (7) or (8) and its direction from the automaton. !ct Dt is greater than one-half the distance Once m between grid points, the water parcel is propagated into a new grid position. The reduced random space algorithm (Miyamoto and Sasaki, 1997) ensures that speeds in cornercentered cell propagation are the same as those in face-centered cell propagation. The inclusion of this algorithm is an improvement for the future. Currently, the variation in speed caused by cell geometry is 40%. 3.2. Smoothed-particle hydrodynamics For the water particle SPH, we assume that a flood of water plus sediments proceeds downslope as a diffusion wave following Eqs. (7) and (8). The water parcels change at time t as htþ1 ¼ hti þ Dhi i where Dhi is given by the SPH solution to B2 h Bh Dhi ¼ D 2 c Dt Bs Bs which is, for the current model, 2 2hk sk 2 2 c esk =d Dhi ¼ Rk pffiffiffi 2 2D d pd
ð13Þ
ð14Þ
velocity and net diffusivity arise from the nature of the topography encountered. With the second algorithm, the full diffusion-routing equations are used and the celerity is solved as c ¼ 3jh2 ðS0 ð1 vÞtanub Þ and the diffusivity as D ¼ jeh3 k
ð17Þ
both directly from Eqs. (7) and (8) (Appendix A). Thus, as v ! 1, c ! 3jh2S0 (Eq. (7)) and as v ! 0, c ! = 3jh2(S0 tan ub) (Eq. (8)). 3.3. Sediment load, hs A thorough treatment of the complex process of erosion is beyond the scope of this contribution. However, to estimate a reasonable value for erosion, we assume that the substrate is easily eroded so that the flood carries a load equal to capacity. For normal stream flow, the sediment load has been found to depend on slope angle and flow speed. The sediment load carried by parcel i at time t, hsit , is therefore (Dingman, 1984; Murray and Paola, 1994) htsi ¼ J ðcðS0 Sc ÞÞn
ð15Þ
pffiffiffi given a Gaussian weighting kernel, i.e., ð1= pdÞ expððsk =dÞ2 Þ: (The reader is referred to Monaghan, 1992 for a full explanation of the SPH technique, including the critical idea of the weighting kernel. Briefly though, the kernel determines the distance over which one fluid parcel influences others.) In Eq. (15), k are the different water particles (not the different cells in the neighborhood), sk is the distance from the ith particle (being tracked) to the kth particle, and d is the length scale over which particle smoothing is done. The celerity and diffusivity have been treated using two different algorithms to understand the effect of the simulated topography on flow propagation, and to explore how a simple model might capture essential flow characteristics. With the first algorithm, m!c and D are taken to be constants. A water parcel moves if a nearby cell has a lower elevation. The resultant net
ð16Þ
ð18Þ
where J and n are dimensionless parameters related to the size of particles, and Sc is a threshold slope angle below which there is deposition for water parcel i. Eq. (18) is based on empirical observations of the sediment capacity of normal stream flow. Debris flows display abnormal bulking relative to normal stream flow. Both the exponent and the multiplicand in Eq. (18) have been varied to simulate different types of bulking behavior (see Appendix A), such as might occur in debris flows. Besides the deposition that occurs when a parcel reaches its sediment capacity, deposition also occurs in a debris flow if the fluid pressure drops sufficiently that particle support vanishes (Eq. (8)). By these erosion and deposition mechanisms, sediment mass is conserved within a flow using the rule t htþ1 si ¼ hsi þ Bj
X
j
ð19Þ
where Bj and Rj are the amount of sediment eroded and deposited by parcel i in path grid cell j.
M. Bursik et al. / Geomorphology 53 (2003) 25–44
3.4. Boundary conditions The initial hydrograph is specified as the depth of water and sediment in each parcel at t = 0. In practice, an initial water depth is assigned to each of N parcels (h0 = hw0, N = 9 in the current simulations) over a restricted area within the DEM at time zero. Because N is fixed in the current implementation, the DEM must be scaled to the size of the initial inundation area occupied by the water parcels. The scaling factor relates the initial volume of water, V, and the hillslope length, L, so that ðV 1=3 =LÞmodel cðV 1=3 =LÞdata :
ð20Þ
4. Results and analyses Realizations of the model run on different synthetic landscapes, from simplified to natural, are compared on a qualitative basis in the following subsections with intuitive expectations and data. In the first subsection, we discuss general features of simulation results. Subsequent subsections treat comparisons of simulation output with field data. Details regarding parameter values and functions used in the simulations are given in Appendix A (bulking and rheological characteristics) and Appendix B (characteristic speed). 4.1. Simple test slopes As they should, the simulated flows propagate at speeds that are dependent on local slope angle. For simulations in which the celerity is fixed, pffiffiffi instantaneous speed is either 0, 1.0 or 1.4 ( 2 ) times the characteristic speed, c (Appendix B), depending on whether there is no lower cell into which to move or the lower cell is at a corner or face of the current cell (Fig. 2). The net velocity measured over any given number of time steps varies from these instantaneous values, depending mostly on the mean slope between the beginning and ending cell. The reason for the dependence of the net velocity on the slope angle can be rationalized as follows. Consider water flowing down a steep, intermediate and shallow slope (Fig. 2). On a steep slope, in which there is always only one
31
cell lower than the initial cell, water will flow consistently along one pathline. As long as this pathline is approximately rectilinear, net velocity will be near the maximum possible value for the given pathline azimuth (i.e., 1.4c if at a diagonal, and 1.0c if on the cell faces). As the slope angle decreases, the probability that a pathline will have start- and endpoints close to one another increases. For intermediate slopes, a relatively high probability still exists that net velocity magnitude will be high. For shallow slopes, however, the probability of a high net velocity magnitude is low. We can formalize the relationship between net velocity and the grid by imagining what would happen on an ideal horizontal plane (Fig. 3). What is the probability that the net velocity magnitude will take on a given value? A certain number of end cells exist in which a water parcel could terminate its flow, and the measured net velocity would take on a given value. A number of paths could be taken to arrive at each of these end cells. Multiplying the number of cells at the given distance or implicit velocity magnitude times the number of paths to each cell yields the number of possible ways in which a given velocity magnitude could arise. The sum of all these yields the total number of paths that the water parcel could take. The probability that a given velocity magnitude will be achieved is then the number of paths with that magnitude divided by the sum of the total number of paths (Fig. 4). The implicit velocity is more likely to be low rather than high, although there is a finite probability that even the maximum possible velocity magnitude will be achieved. The probability is 7% if averaging is over two time steps that the maximum possible velocity magnitude will be measured, whereas the probability is < 1% if averaging is over three time steps. A nonzero velocity magnitude will more likely be measured than a zero velocity magnitude. A low value for velocity magnitude will more likely be measured on a surface that dips at a shallow angle than on a more steeply dipping surface. The simulations with variable celerity show quite similar behavior to those with fixed celerity. Overall, as the slope angle decreases and the flow spreads, instantaneous celerity decreases. The net down-channel velocity (the true velocity as averaged over a number of time steps) also decreases with slope angle.
32
M. Bursik et al. / Geomorphology 53 (2003) 25–44
Fig. 2. Schematic diagram showing the dependence of velocity on (a) mean slope angle, and (b) mean slope angle and slope irregularity. Cells with horizontal stripes signify preceding time step; cells with diagonal stripes signify current time step; cells with stipple signify potential path; cells with vertical stripes signify final position at end of measurement interval. The maximum possible velocity is ustraight if flow is across cell face centers.
M. Bursik et al. / Geomorphology 53 (2003) 25–44
33
The flows on this test slope typically show a downstream branching pattern as they develop that is only weakly seen in the final runout pattern. The branching and resulting crosscutting lobate pattern are typical for braided streams and debris flow fronts observed in the field, where they often arise by levee breaching and overtopping. In the simulations, they arise because of the patchy deposition pattern, which creates local high spots that were not encountered by the water parcels that passed earlier. These local, often transient high spots must be circumnavigated by the water parcels. 4.2. Gully development on Walker Lake scoria cone
Fig. 3. Schematic diagram illustrating the probability that a given velocity magnitude will be attained on a horizontal plane. Same pattern scheme for cells as in Fig. 2.
Results are discussed more thoroughly in the section on Popocate´petl. Flow paths for both channelized and unconfined flow were investigated by performing simulations on a synthetic slope that had channelized and unconfined segments (Fig. 5). In these experiments, flow was initially confined in a narrow channel. Flow was bank-to-bank within this narrow reach. Where the channel doubled in width, the flow meandered between the channel banks. Where the channel width increased to a value much greater than the initial channel width, the flow swung in a single, broad loop across the slope, showing a complex braiding and occasional lobe-forming pattern as it developed. If the reach had been longer, it may have meandered between edges. Overall, the active flow pattern qualitatively resembled the meandering, bifurcating pattern common to flow on rough slopes, heavily laden stream flow, or debris flow (e.g., Murray and Paola, 1994). The flow pattern did not resemble that observable in laboratory experiments and on smooth engineered slopes in which monotonic thinning and broadening occurred at transitions in channel width (Chaudhry, 1993).
Erosion and deposition in a natural rill and gully system on a small volcanic scoria (cinder) cone can be used to validate hillslope-scale erosion and deposition in the model system. The San Francisco Volcanic Field, Arizona, is a neogene scoria cone field consisting of late Miocene to Holocene volcanic rocks. The scoria cones show different stages of erosion, partly due to great age diversity. In this semiarid region, heavy monsoon rainfall in July to September profoundly modifies the parched regolith. Walker Lake scoria cone (N 35j23V30U and W 111j44V50U), a tuff ring (Blauvelt, 1998) of age 2.01 F 0.22 ma (Wolfe et
Fig. 4. The graph illustrates the probability that a given velocity magnitude will be measured if averaging velocity over two or three time intervals. N(DT) is the length of the interval over which averaging is done, in multiples of time steps. The maximum possible velocity is ustraight if flow is across cell face centers. The maximum possible is 1.4 times this value (flow across cornercentered cells).
34
M. Bursik et al. / Geomorphology 53 (2003) 25–44
Fig. 5. Schematic view of an experiment to test the effects of channel widening and flow over an open slope on flow pattern. The walls of the channel have been removed to be able to visualize the flow. The water parcels flowed down a rough ramp of constant average slope angle. The scale of the flow deposit (areas of raised relief) has been exaggerated relative to that for the topography for illustrative purposes. No scale.
al., 1987) in the late Volcano to early Planeze stage of Kear (1957), is of particular interest for model evaluation. A forest fire destroyed much of its vegetation in 1996, shortly before the monsoon season began. Rills and gullies developed rapidly at that time on the 170-m high, 800-m wide volcano. 4.2.1. Description of the gully system Measurements of surface degradation and gullying were taken from the first monsoon season following the 1996 fire, to the most recent season in 2000. Gully systems on the WSW flank have been particularly active. Rills initiate on the upper midslope, which consists of a burned area on which unstable, infiltration-excess overland flow develops during rainfall events of 2 –4 cm h 1 and f 1.5 cm total. Gullyheads occur on the lower midslope, where midslope return flow meets converging rill flow. One of several gullies was chosen for detailed investigations. The gullyhead consists of two main channels of average depth 0.5 to 1.2 m. The headscarp region was 7 m wide in 1996, and has increased in width and depth, and migrated upslope with each subsequent monsoon season. At 6 m downslope from the gullyhead, the two channels converge into one channel of triangular
cross-section that conducts water to the footslope break, where the channel splits into a distributary system of rills on the alluvial footslope. Erosion and deposition processes have been recorded with erosion pins, which are 0.5 –0.7 m long iron bars set up along the length of the rill –gully– footslope system. The erosion pins were driven into the regolith, and the surface height change between rainfall events was measured as changes in surface elevation against the relatively stable erosion pins. 4.2.2. Simulation of gully development A simple, scaled DEM of an ungullied hillslope was made in TIF image format. The 10-m horizontal, 0.2-m vertical spacing of the DEM allowed us to resolve relatively small changes in elevation at a horizontal spacing comparable to that of the spacing between erosion pins on the natural hillslope. The slope angle of the DEM averaged f 6j. Although this is the same as the average slope angle within the natural test area, it is lower than the angle ( f 15j+) within that part of the test area with which we have been able to make comparisons. Water parcels were initiated on the upper part of the DEM in numerous realizations of rainfall events.
M. Bursik et al. / Geomorphology 53 (2003) 25–44
Although most parcels were initiated within a restricted area, some parcels were started at other positions on the DEM to attenuate the unnatural regularity of the DEM. A number of realizations were run to initiate and develop a gully. Realizations were run until the gully appeared to have approximately the same degree of development and scaled dimensions as the natural gully in planform. The model gully runs almost straight across contours at a consistent depth, then shallows and breaks into a distributary system of subgullies or rills. This breakup occurs at a distance of 90 m downslope from the gullyhead, which is similar to the breakup of the natural system 70 m downslope. In nature, a small change in slope angle is associated with the termination of the gully. No intergully slope angle change occurs in the model system. The model gully was developed by this means until it was 2 –4 m deep; the natural gully averages about 1 m deep. On a local scale, the bottom of the model gully is rougher than the bottom of the natural gully. Once the model gully had been developed to about the same morphology as the natural gully, a realization was run for comparison with the measurements obtained from the erosion pins. The results of this
35
realization compare favorably, in terms of erosion and deposition patterns, with natural erosion and deposition patterns generated during storms (Fig. 6). The simulated amount of water in the one test event was 28 m3. The real gullyhead begins 150 m downslope from the drainage divide. Gullies are spaced f 10 m laterally across the scoria cone hillslope. Thus, simulated runoff of 28 m3 is equivalent to the runoff of f 2 cm of rainfall from the gully catchment of f 1500 m2. Our measurements indicate that infiltration-excess overland flow begins during moderate to large rainstorms once f 1.5 cm of rain has fallen. Two cm of rainfall runoff is thus equivalent to f 3.5 cm total rainfall. In the small natural storm used for comparison, 1.3 cm of rain fell, and in the large storm, 8.7 cm of rain fell. The simulated erosional pattern, which is intermediate to the erosional pattern seen in the natural events in both amplitude and extent, is therefore qualitatively consistent with the simulated rainfall total. We know that erosion and deposition in the natural gully system are controlled by local slope angle. To test the effect of changes in slope angle along the runout path, one model run was realized on a slope having an upper section at 30j and lower section at 0j. The final erosion and deposition pattern is spotty
Fig. 6. Results for the gully on the SW side of Walker Lake cone, San Francisco volcanic field, Arizona. Hillslope surface-height changes were measured during episodes of drizzle to light rainfall (3 mm h 1, small storm) and moderate to heavy rainfall (15 mm h 1, large storm) (Laws, 1941). The data for two events are compared to one realization of the model. In previous realizations, the model gully was initiated and developed to a geometry similar to the natural geometry. In the realization illustrated, the model, event-scale erosion, and deposition pattern are shown to be of a scale intermediate to that for the two natural events.
36
M. Bursik et al. / Geomorphology 53 (2003) 25–44
or patchy (Fig. 7) as observed in the natural gully, and as mapped elsewhere in the field for braided streams and debris flows (e.g., Janda et al., 1981; Pierson et al., 1990). During flow, patches of erosion and deposition were transient; and the final erosion and deposition pattern arose only after the last water parcels had passed. The net degradation pattern showed the greatest amount of erosion in the headwaters and the least amount near the distal margin. Deposition was concentrated in a thick lobe on the lower horizontal segment. These results show that the model responds in the correct manner to changes in slope angle. The model does, however, produce a deposit that is more localized than was observed in the natural gully. This is because the criterion used for deposition by debris flow causes almost all material to be deposited in one grid cell. 4.3. Popocate´petl debris-laden flow On 21 December 1994, Popocate´petl volcano—50 km from downtown Mexico City—reawakened and has continued to erupt to the present (October, 2001). The largest eruption of the current episode was on 30 June 1997. Early on 1 July 1997, a debris-laden flow reached the village of Xalitzintla, within the Alseseca channel on the NE flank of the volcano (Fig. 8). We
used field studies of this event to test and validate the manner in which the current model handles landformscale phenomenology. 4.3.1. Description of the event Inspection of the erosion and deposition patterns resulting from the flood of July 1 were carried out in the months following the event. Tracing the deposits upstream from Xalitzintla, we found that the flood originated near the outlet tongue of the Ventorillo glacier, at an elevation of ca. 4800 m. The total runout length was 17 km. The maximum runup was 8 m in the proximal region (0 –5 km from the glacier tongue), 5 m in the medial region (5– 12 km), and 0.8 m in the distal region (12 – 17 km) (Fig. 9). The proximal flow speed was f 20 m/s. The speed was estimated from a superelevation of 7.3 m across a radial distance of 25.6 m between flow-levee margins on a bend of radius 140 m (Pierson, 1985). Grain size data showed that the flood bulked with clays as it propagated downstream so that the deposits have the appearance of typical lahar debris in the medial and distal reaches (Fig. 10). Thus, the flood evolved with distance into a true debris flow. The flood both eroded and deposited in the proximal reach, primarily deposited in levees in the medial reach, and deposited a terminal lobe in the distal reach. In the following section, we compare
Fig. 7. Detail of a fully scaled flow realization for an experiment on a ramp (slope angle 30j) that debouches onto a horizontal plane for testing of the response of the model erosion and deposition characteristics to changes in slope angle. Flow was unconfined on the ramp and on the plane.
M. Bursik et al. / Geomorphology 53 (2003) 25–44
37
Fig. 8. (a) Map of the Popocate´petl region showing separate debris-laden flow realizations (pixelated gray areas) given different initial hydrographs and starting points on the glacier. Most simulated flows go down the Alseseca channel toward Xalitzintla, although a fraction propagates down the Nexapa channel and down the broader channel toward the region of Atlautla, depending mostly upon starting position. Northing and easting are shown in UTM coordinates. The topography used in the current simulations is shown at a 500-m contour spacing. Major towns and villages around the northern half of the volcano are shown in the grid pattern. Some of the stream channels tributary to the Alseseca channel are shown in light gray. Two facies of the 1600-year BP San Nicolas lahar deposit (Gonzalez-Huesca et al., 1997) are shown in horizontal lined pattern for comparison. The flow path of the 1 July 1997 lahar is shown in bold dark gray, also for comparison. (b) The expanded view shows the upper region of the volcano with sites (Figs. 9 and 10) shown by filled squares.
model realizations with these data to determine whether the model yields physically reasonable results over spatially large and complex regions. We emphasize that the model is not being developed to ‘‘compete’’ in the field of analysis of lahar hazards and inundation zones; numerous numerical and quantita-
tive techniques already exist (e.g., Honda and Egashira, 1997; Iverson et al., 1998). An elevation model of the volcano was digitized from government topographic maps and saved as an ASCII text file of integer elevation values. Only the largest channels, such as the Alseseca, were resolved
38
M. Bursik et al. / Geomorphology 53 (2003) 25–44
Fig. 9. Channel cross-sections measured in January 1998, at sites shown in Fig. 8. Stipple pattern represents material deposited in the July, 1997, event. Most of the depositional levees are related to the July event; some of the erosion in the channel bottom is related to activity subsequent to July, 1997.
M. Bursik et al. / Geomorphology 53 (2003) 25–44
39
Fig. 10. Grain size analyses for samples from the Popocate´petl debris-laden flow deposit for different reaches, at localities given in Fig. 8. Note the fining of the deposit with downstream distance.
with the 100-m horizontal, f 4-m vertical spacing of the DEM. Slope angle transitions as low as f 4% and features as small as 100 m in lateral extent were distinguishable. 4.3.2. Comparison of model and data Several tests were run to determine general features of the model and DEM before simulations of the Popocate´petl debris-laden flow. To test the fidelity of the DEM in terms of the locations of channels,
we inspected whether simulated runout patterns did follow known channels—including the Alseseca— that lead downslope from the glacier. The simulated floods did follow the larger channels (Fig. 8). Once onto the plain at the base of the cone, where gradients are lower and water parcels propagate in a variety of directions, the simulated floods typically spread over a relatively wide region, as exhibited by a prehistoric lahar downstream from the village of Xalitzintla.
40
M. Bursik et al. / Geomorphology 53 (2003) 25–44
Fig. 11. Log – log plot of runout length as a function of total volume Shi for a number of debris-laden flood realizations. Runout length of the water flood is defined as the total distance traveled along channel by any of the water parcels. Runout length of the debrisladen flow is defined as the distance between the source region and the most distal downstream pixel showing deposition. As discharge increases, debris-laden flow runout length increases because of the increased ability of the water to carry debris past channel reaches that have a gradient lower than the critical slope angle Sc.
Simulated debris-laden flood runout lengths, defined as the distance from the initiation zone to the most distal deposition site, depended on total volume normalized by cell size, Rihi (Fig. 11), as is the case in nature (Pierson, 1995). In contrast, the runout distance of debris-free water floods was not dependent on discharge. The motion of simulated debris-free stream flow was arrested only by ponding in closed depressions. The difference between debris flows and floods is seen on volcanoes, where debris flows are often transformed by deposition into water floods that continue further downslope (Janda et al., 1981). 1 July 1997 event: We can compare some features of model simulations with data for the July 1997 event. A number of poorly constrained parameters are associated with debris flows. Of particular concern is the fraction of sediment, and how that fraction affects viscosity and other intrinsic parameters (Appendix A). A detailed evaluation of these parameters is beyond the scope of this work. Simulated paths and runout distances show good agreement with the data. Model and data are similar for flow realizations with about 105 – 106 m3 of water initially (Fig. 8). The amount of rainfall and glacial
ablation available for the flow are not known because no measurements are available, but the amount of water used in the realizations represents runoff from about 1– 10 cm of rainfall and ice-melt on the uppermost flanks of the volcano. Looking at one realization of the model as an example, we compare proximal, medial and distal runup heights. Unlike some of the other parameters, the true inundation heights differ markedly from the inundation heights in the realization (Fig. 12). This is partly because of the resolution of the DEM, but also because the simulated flood waves peak to unusually high values near the source region. Given that the model and data runup heights show the same downstream development, some adjustment of the parameters of the model could yield a better fit. Hydrographs that develop through time as the flood wave proceeds downslope both attenuate and disperse (Fig. 13), as has been observed for debris flow hydrographs from numerous volcanoes (Pierson, 1995). For comparison of the effects of the diffusion coefficient on hydrographs, we use simulations with fixed diffusion coefficient. Measured at the same distance from the source region, the flow with the greater diffusion coefficient, D, shows greater attenuation and dispersion, as expected. A nondispersive flood wave could be simulated by adjusting to D = 0 (kinematic routing).
Fig. 12. Runup or inundation heights as a function of distance, comparing model results with data. In the models, inundation height is taken to be the depth of flow for the water parcel at the peak of the hydrograph. Brackets indicate time averages over 21 time steps.
M. Bursik et al. / Geomorphology 53 (2003) 25–44
Fig. 13. Model initial hydrograph and hydrographs after given number of time steps for simulation with constant c and D, showing the diffusion of the hydrograph with distance and time. Lambda (k) is flood wave length, and hmax(0) is peak flood height at t = 0. The diffusion of the flood wave is similar to that seen in real events. Negative value of depth is the result of numerical error.
Flood velocity is consistent with the one estimate from the event (Fig. 14). We can check the model flow speed at the point where true speed was estimated to be f 20 m/s. Instantaneous flow velocities at this distance are of the same magnitude as was estimated from the field data. The field data indicate bulking with distance. For the realization used above, bulking is seen in the
41
Fig. 15. Fraction of water in flow, v, as a function of distance showing bulking behavior up to 4 km from source, typical for lahars that originate from breakout floods or precipitation (Scott et al., 1995).
model up to about 4 km from the source (Fig. 15). Beyond this distance, much of the sediment is lost from the flow and bulking values are low. Again, with some adjustment of the relatively poorly known values for the model parameters, a good fit to the data could probably be achieved.
5. Discussion and conclusions
Fig. 14. Variation in c and u with distance in the proximal region. Filled triangle shows the one estimate of velocity from the field at locality 2. The lines are magnitudes averaged over 21 time steps (indicated by brackets around the symbols); symbols are samples of instantaneous values. This is for a realization with variable c and D.
The present work includes advances in the modeling of hillslope flow, landform degradation, and channel development using cellular automata. A new, coupled cellular automaton – smoothed-particle hydrodynamic technique has been developed to simulate the propagation of sediment-charged floods and debris flows as components in an evolving hillslope degradation system. The model appears to yield reasonable flow speed dependencies, paths, and erosion/deposition patterns on channelized and open synthetic hillslopes. Although the numerous erosion relationships that were tested are heuristic, these relationships were readily adjustable to yield an erosion/deposition pattern on a simulated hillslope gully that was similar to the pattern observed in a real gully on a hillslope in Arizona. Some characteristics of a larger simulated flow were consistent with observations of and estimates from a recent debrisladen flood in Mexico. Inundation heights were greater
42
M. Bursik et al. / Geomorphology 53 (2003) 25–44
than those observed in the real event; but flow speed, runout length and hydrograph evolution with distance were consistent with observation. The model advances the automaton technique in several respects. (The enumeration below to some extent follows that of Favis-Mortlock (1998) to help clarify the contributions of this work.) (i) The model is scale independent, although the number and size of water parcels should be scaled to the topographic grid. The model in some respects bridges the gap between the landform model of Chase (1992) and the hillslope model of Favis-Mortlock (1998) and its progenitors. (ii) The water parcels (‘‘precipitons’’ of Chase, 1992; ‘‘raindrops’’ of Favis-Mortlock, 1998) can be of different sizes. (iii) Numerous water parcels are released onto the elevation grid at one time. (iv) The water parcels interact dynamically. (v) One algorithm implemented as a result of this work showed that it is possible to obtain a dependence of flow speed on slope angle without specifying this dependence explicitly in the equations of motion. This result may have ramifications for construction of models intended to propagate flows over DEMs displaying irregular topopgraphy. (vi) Flow velocity and dispersion follow the diffusion wave equation. (vii) Channels can develop and evolve fully in threedimensional form. (viii) Deposition occurs both when capacity is exceeded in normal stream flow, and when the inclination drops below the friction angle in debris flow. Shortcomings in the current effort are: (i) The stability and convergence characteristics of the technique are not treated in this work. (ii) The model does not treat rapidly varied flow conditions. (iii) Several parameters were left fixed or simplified, such as flow depth or pressure. (iv) No new water parcels can be added once a realization has been initiated.
(v) No infiltration is assumed. (vi) No spatial variations in substrate characteristics have been included. (vii) Although the model runs in the time domain, time-varying hillslope parameters have not been used. (viii) Only overland flow is modeled. The model does not treat an entire degrading hillslope system, which also includes transport by rainsplash, soilfall, etc. The model is currently being developed to address some of the shortcomings (Appendix C). One advanced treatment includes the effects of rainsplash transport. Because the model is based on algebraic relationships and is implemented as linked automata on image data, local properties can be readily defined and varied, spatially and temporally, to allow for substrate and flow path heterogeneities.
Acknowledgements This research was supported in part by grants from NASA (NAG53142) to study landscape evolution, NSF (ACI-0121254) to model geophysical flows, UNAM (DGAPA-IN103393), and Secretaria de Gobernacion (Ministry of the Interior) through the National Center for Disaster Prevention (CENAPRED). Two anonymous reviewers and the editors are acknowledged for useful suggestions that improved the text. We thank B. Pitman, J.L. Macı´as, T. Pierson, J. Major, and A. Patra for useful suggestions about the work. J.J. Monaghan kindly provided reprints and preprints of many of his works. The DEM was based on two separate digitizations of the standard topographic maps for the region, one by M.M. Brugman and the other with the aid of M.F. Sheridan. Occasional corrective measures to our thinking by A.D. Abrahams are greatly appreciated.
Appendix A A first-order effect on the rheological properties of a water and sediment mixture is exerted by the fraction of sediment, particularly clay (see references in Scott et al., 1995). We have not characterized the
M. Bursik et al. / Geomorphology 53 (2003) 25–44
clay content specifically in this work but only require that rheology depend upon sediment content, which is defined as for a dilute mixture: v¼
hw ðhw þ hs Þ
ðA:1Þ
Perhaps the most important variable that depends on sediment content is viscosity or for a turbulent flow, some variable such as friction factor, eddy viscosity or Manning’s constant. A smooth variation in viscosity with water content can be constructed with a relationship such as: mðvÞ ¼ k1 þ
k2 : 1 þ ðv=k3 Þn
ðA:2Þ
The parameter k1 is set to the minimum value that m can attain over the interval of possible v (i.e., f 0.2 < v < 1.0 for sorted mixtures of fluid and solids). As v ! 1,m ! k1. The value of k1 + k2 is near the maximum possible value that viscosity can attain as v ! 0. The values of n and k3 are adjusted to allow convergence to asymptotic values within some tolerance level. Because the manner of dissipation in the flow changes with changing water content, m is an effective viscosity. The particular function used in the current analysis causes a greater viscosity with the addition of sediment. However, an effective viscosity might actually decrease in value with increasing soil content as the flow transforms from a turbulent to a laminar mode of dissipation. As noted in the text, a standard relationship that has been found effective for sediment load in normal stream flow is htsi ¼ J ðcðS0 Sc ÞÞn
2hti 1 þ 21 2 c =v
ðA:5Þ
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi cðSl Sc Þ þ v cðSl Sc Þ ðA:6Þ
in which Sl was tried as either S0 or Sf in a given simulation. (Favis-Mortlock, 1998 discusses a relationship appropriate to rill flow.) All the above functions have the property of converging to a constant value of sediment load as slope, speed (celerity), or water volume increase. In general, Eq. (A.4) provided the most stable solutions, with Eq. (A.5) somewhat less stable. For the current simulations, the only other major parameters that were unconstrained were e and tan ub. A constant value of unity was used for e, while values between 0.0 and 0.3 were tried for tan ub. Perhaps the best results were obtained using a value of f 0.1. Thus, the models were in fact developed with relatively few free parameters. The values of some of these parameters and their dependence on more fundamental variables are quite poorly constrained at this time.
Appendix B For simulations in which the celerity is kept fixed with distance, the simulated speed is scaled to a characteristic flow speed. The characteristic speed is taken to be the normal speed for the flow at the source, un (Singh, 1996). The normal depth, hn, at the source satisfies Eq. (2) for the condition of steady, uniform flow. Eq. (2) in this situation reduces to, assuming a dilute flood, 0 ¼ gS0 4vm
ðA:3Þ
At low concentrations, this relationship has been found to be accurate. However, it cannot hold for the high concentrations and therefore different rheological and suspension behavior that occurs in debris flows. We have experimented, therefore, with three possible, heuristic relationships ! ðSl Sc Þ2 t t hsi ¼ 2hi exp ðA:4Þ cv2c htsi ¼
htsi ¼ 0:02ð1 vÞ
43
u h2
ðB:1Þ
or un ¼ jS0 h2n
ðB:2Þ
Because flow depth rather than discharge is specified at the source, the normal speed corresponds to that for the initial flow depth or un ¼ jS0 h20
ðB:3Þ
With typical values for the parameters in Eq. (B.1) of v f 1, m f 1 m2/s, S0 f 0.1, and h0 f 10 m, un f 25 m/s.
44
M. Bursik et al. / Geomorphology 53 (2003) 25–44
Appendix C The models discussed in the text were written as Pascal macro scripts for the free, open-source program, NIH Image. Digital elevation models consisted of either 8- or 16-bit integer image or text files. The latest version of the macros and images upon which they can operate are available for download at http:// www.volcano.buffalo.edu/pub/automaton. Image (Scion Image) is available for Windows platforms from http://www.scioncorp.com and for the Macintosh at http://rsb.info.nih.gov/nih-image. New versions of the script are currently being coded as Java macros for NIH ImageJ—the Javabased successor to Image, and in C, for model scaling. The Java macros will run through the web.
References Blauvelt, D., 1998. Evolution of processes during landform degradation: examples from the San Francisco Volcanic Field, AZ. MA thesis, Dept. of Geology, State University of New York at Buffalo, Buffalo, NY, USA. Chase, C.G., 1992. Fluvial landsculpting and the fractal dimension of topography. Geomorphology 5, 39 – 57. Chaudhry, M.H., 1993. Open-Channel Flow. Prentice-Hall, Englewood Cliffs, NJ. 483 pp. Chen, C.-L., 1997. Debris-Flow Hazards Mitigation: Mechanics, Prediction and Assessment. Amer. Soc. Civ. Eng., New York. 817 pp. Dingman, S., 1984. Fluvial Hydrology. Freeman, New York. Favis-Mortlock, D., 1998. A self-organizing dynamic systems approach to the simulation of rill initiation and development on hillslopes. Computers & Geosciences 24, 353 – 372. Forrest, S.B., Haff, P.K., 1992. Mechanics of wind ripple stratigraphy. Science 255, 1240 – 1243. Gonzalez-Huesca, A.E., Delgado, H., Urrutia-Fucugauchi, J., 1997. The San Nicolas lahar at Popocatepetl volcano (Mexico): a case study of a glacier-icemelt-related debris flow, triggered by a blast at the onset of a plinian eruption. IAVCEI General Assembly Abstracts, Gobierno de Jalisco, Guadalajara, Mexico. Honda, N., Egashira, S., 1997. Prediction of debris flow characteristics in mountain torrents. In: Chen, C.-L. (Ed.), Debris-Flow Hazards Mitigation: Mechanics, Prediction and Assessment. Amer. Soc. Civ. Eng., New York, pp. 707 – 716. Iverson, R.M., 1997. The physics of debris flows. Reviews of Geophysics 35, 245 – 296. Iverson, R.M., Schilling, S.P., Vallance, J.W., 1998. Objective delineation of lahar-inundation hazard zones. Geological Society of America Bulletin 110, 972 – 984.
Janda, R.J., Scott, K.M., Nolan, K.M., Martinson, H.A., 1981. Lahar movement, effects, and deposits. In: Lipman, P.W., Mullineaux, D.R. (Eds.), The 1980 Eruptions of Mount St. Helens, Washington. U. S. Geol. Surv. Prof. Pap., vol. 1250, pp. 461 – 478. Washington. Kear, D., 1957. Erosional stages of volcanic cones as indicators of age. New Zealand Journal of Science and Technology 38, 671 – 682. Laws, J.O., 1941. Measurement of fall-velocity of water drops and raindrops. Transactions-American Geophysical Union 22, 709 – 721. Major, J.J., 1997. Depositional processes in large-scale debris-flow experiments. Journal of Geology 105, 345 – 366. Miyamoto, H., Sasaki, S., 1997. Simulating lava flows by an improved cellular automata method. Computers & Geosciences 23, 283 – 292. Molinaro, P., Natale, L., 1994. Modelling of Flood Propagation Over Initially Dry Areas. Amer. Soc. Civ. Eng., New York. 373 pp. Monaghan, J.J., 1992. Smoothed particle hydrodynamics. Annual Review of Astronomy and Astrophysics 30, 543 – 574. Murray, A.B., Paola, C., 1994. A cellular model of braided rivers. Nature 371, 54 – 57. Pierson, T.C., 1985. Initiation and flow behavior of the 1980 Pine Creek and Muddy River lahars, Mount St. Helens, Washington. Geological Society of America Bulletin 96, 1056 – 1069. Pierson, T.C., 1995. Flow characteristics of large eruption-triggered debris flows at snow-clad volcanoes: constraints for debris-flow models. Journal of Volcanology and Geothermal Research 66, 283 – 294. Pierson, T.C., Janda, R.J., Thouret, J.-C., Borrero, C.A., 1990. Perturbation and melting of snow and ice by the 13 November 1985 eruption of Nevado del Ruiz, Colombia, and consequent mobilization, flow and deposition of lahars. Journal of Volcanology and Geothermal Research 41, 17 – 66. Press, W.H., Flannery, B.P., Teukolsky, S.A., Vetterling, W.T., 1986. Numerical Recipes. Cambridge Univ. Press, Cambridge, UK. 818 pp. Savage, S.B., Hutter, K., 1989. The motion of a finite mass of granular material down a rough incline. Journal of Fluid Mechanics 199, 177 – 215. Scott, K.M., Vallance, J.W., Pringle, P.T., 1995. Sedimentology, behavior, and hazards of debris flows at Mount Rainer, Washington. U.S. Geological Survey Professional Paper, vol. 1547. Washington. Singh, V.P., 1996. Kinematic Wave Modeling in Water Resources: Surface-Water Hydrology. Wiley-Interscience, New York. 1399 pp. Wolfe, E.W., Ulrich, G.E., Holm, R.F., Moore, R.B., Newhall, C.G., 1987. Geologic map of the central part of the San Francisco volcanic field, North Central Arizona. U.S. Geologic Survey Miscellaneous Field Studies Map MF-1957, Washington. Wolfram, S., 1983. Cellular automata. Los Alamos Science 9, 2 – 21.