EnvironmentalGeology(1993) 22:141-152

9 Springer-Verlag 1993

P. Caruso - M. T. Pareschi

Estimation of lahar and lahar-runout flow hydrograph on natural beds

Received: 17 September 1992 / Accepted:22 February 1993 Abstract A numerical model is used to simulate channelized lahars, flowing with a near-constant sediment concentration and volume. Water and debris are usually mobilized in short times and the peak discharge of lahars may reach many thousands of cubic meters per second in a few minutes. A relation for the energy dissipation term must be provided in the model, which in turn depends on debris flow rheology and shape and status of channel bed. A discussion of the form of this term is performed through the simulation of some historical (among the welldocumented) lahars and lahar-runout flows with concentrations ranging in a wide spectrum up to 70 percent by volume and irregularly shaped sand and coarser particles dispersed in a fluid matrix of water and fine material. As concentration increases and turbulence decreases, the dissipative term, which, in the first case, is proportional to the square average flow velocity (Manning or Chezy formulation) is well described by a linear dependence on flow velocity, as expected in the laminar case. The numerical reproduction of the examined historical cases suggests that the model can be used for hazard assessment, if some hypotheses are made about lahar hydrograph at the source, its volume, and the shape of the dissipative term. Key words Debris flow-- L a h a r - - Mudflow-- Numerical simulation of channelized debris flows

Introduction Recently, lahars have been defined by Pierson (1985) as "volcanic mudflows or debris flows," comprised of"dense, viscous slurries of poorly sorted gravel, sand, mud, and P. Caruso CNR-GNV M. T. Pareschi CNR-CSGDSA, Via Santa Maria 53, 1-56100 Pisa, Italy Offprint requests to: M. T. Pareschi

water". The term "lahar" as used in this paper refers loosely to rapidly flowing mixtures of sediments and muddy water that originate on the slope of a volcano, directly or indirectly triggered by volcanic eruptions. They have received particular scientific attention since the Mr. St. Helens (Washington, USA) eruption of 1980. In the first minutes of that eruption, the erosion of summit glacier tongues by pyroclastic flows triggered the Pine Creek and Muddy River lahars southeastward and the South Fork lahar westward. Five hours later, the North Fork lahar occurred, triggered by the seismic shaking of a debris avalanche on the north side of the volcano (Cummans 1981a, b, Janda and others 1981; Pierson 1985; Scott 1988). Although the South Fork, Pine Creek, and Muddy River lahars did little damage, the North Fork lahar, characterized by a volume of over 100 million cubic meters and by a long duration of the peak discharge, destroyed two large equipment storage facilities, six large and three small bridges, and more than 100 houses~ The interest in lahars was renewed by the occurrence, on 13 November 1985, of the catastrophic lahars of Nevado del Ruiz volcano (Colombia) which killed more than 22,000 people. Pyroclastic flows and surges, occurring during the eruption, melted tens of millions of cubic meters of glacier ice near the summit of the volcano, triggering catastrophic lahars along Azufrado, Lagunillas, and Guali valleys on the east side of the volcano and Molinos/Nereidas valley on the west side (Pierson and others 1990). Table 1 reports fatalities due to lahars triggered by some notable volcanic eruptions since the year 1000 AD (Yokoyama and others 1984). Sediment content, which is among the main factors affecting the rheological properties of a debris flow, ranges from 75-90 weight percent for lahars to 40-70 weight percent for hyperconcentrated flows (lahar-runout flows), characterized by an oily, plany smooth appearance, suppression of secondary turbulence, and a smaller average grain dimension (< 1 mm) (Pierson and Scott 1985; Pierson 1986). Flow transformations from lahar to lahar-runout flow and vice versa occur for many reasons, including (1) water floods progressively incorporate more sediment and

142

-ff

Table 1. Fatalities due to lahars triggered by some volcanic eruptions since 1000 AD (Yokoyama and others 1984) Volcano

Country

Year

0 tt) tw kLI

Fatalities

O ~ )...

0 H

Kelut Awu Cotopaxi Galunggung Nevado del Ruiz Awu Cotopaxi Awu Kelut Nevado del Ruiz

Indonesia Indonesia Ecuado r Indonesia Colombia Indonesia Ecuador Indonesia Indonesia Colombia

1586 1711 1741 1822 1845 1856 1877 1892 1919 1985

10,000 3,200 1,000 4,000 1,000 3,000 1,000 1,530 5,110 > 22,000

E o o9

>.e-~ t'~t~O

Io

i

,

,

.

i

9

E E o ~

,:1

(Slg LU ) e 6 ~ e q o s ! O

-g Total

E

v

o

O

o

r

53,900

E

o o v 113

,,e>-o . UJ 03

become lahars, and (2) lahars are diluted by the deposition of sediments (in concomitance, for example, with strong slope variations) or by the addition of water (coming from tributary rivers or preexisting water along the channel path). The South Fork lahar of 1980 shows, for example, a transition to hyperconcentrated flow about 45 km away from the Mt. St. Helens crater, owing to the confluence with North Fork River (Scott 1988). Lahars are generated in many ways (Neall 1976). Volcanic eruptions can catastrophically release water contained in crater lakes. Avalanches of water-satured debris can transform themselves into lahars. Pyroclastic flows can melt ice and snow and form debris flows. Torrential rainfalls, in some cases originated by rainstorms related to convecting eruption columns, on recently deposited tephra or other unconsolidated material, can trigger lahars. Generally, owing to the triggering mechanisms, great volumes of water and debris can be mobilized in short times (up to tens of thousands of cubic meters per second in a few minutes). According to the triggering mechanisms, lahar hydrographs have different shapes, which modify downvalley according to channel morphology, lahar volume, and sediment content. Hydrographs reconstructed near the source for lahars triggered by hot pyroclastic flows show a roughly triangular shape (Fig. 1b-g). In contrast, the initial hydrograph of the North Fork lahar, triggered, as already mentioned, by a harmonic tremor, had a trapezoidal shape with a rapid increase, a persistence of the maximum discharge for about 289 h (the duration is related to that of the triggering tremor), and a slower decrease (Fig. la) (Fairchild 1987; Pierson 1985; Pierson and others 1990). The velocities of historical lahars have varied greatly owing to differences in volume, channel dimensions, and grain-size distribution. Figure 2 shows the velocity range for some historical lahars: the 1980 lahars at Mt. St. Helens (Pierson 1985; Pierson and Scott 1985; Laenen and Hansen 1988); the Mt. Ruapehu (New Zealand) lahars of 1975 along the Wangaehu, Whakapapa, and Mangaturuturu valleys, triggered by the drainage of the Crater Lake (Paterson 1976); and the 1985 lahars at Nevado del Ruiz (Pierson and others 1990).

oVco I-O

o I

o O O O ~

O O O O '~"

00

9

0 0 0 C~

~

1

9

i

0 0 0 0

9

0 0 0 Q

i

,

0 0 0 0

i

.

0 0 0 0

(S/~m)aSJeq

w

0

~Sl(]

(Slgm)ef~eq~s!o E E

Er

o tt)

,,o

L,. c

o

E ~

o t=

O

,o| 9"r)...

o

O

tl3 ~ O

o03i= f

o

o

I-u

i

9

,

i

9

o

,

O

o oo

CO

C,J

~-

co

i

g oo

o v--

i

o g

o

(Slgm}egJ eqos![l

~e4

(s/~oa)o~Jeqos!o

co

E

tm o

o

tt)

mE~ n," 0

o

/

o

V0

o o o

o 0 o

o

u

~o 03

(Sl gWl oS~ e4o~!O

9

,

i

o

~g

,

!

9

oo

o

,

o

o co

(Sly: tu )o6Jetl os!o

Fig. la-g. Hydrographs near the source for past lahars, a-d refer to 1980 lahars at Mt. St. Helens: a North Fork lahar (after Fairchild 1987); b South Fork lahar (after Fairchild 1987); c Pine Creek lahar (after Pierson 1985); d Muddy River lahar (after Pierson 1985). Peak discharge of North Fork lahar is one order of magnitude smaller that the others, e-g lahars of 1985 at Nevado del Ruiz, Colombia (after Pierson and others 1990): e Molinos/Nereidas lahar; f Guali lahar; g Azufrado lahar

143 50

20 AZUFRADO ~

~" 40

~.

MUDDY RIVER

--

PINE CREEK

----o----

.~ 30

E

15

--

10

GUALI

NORTH FORK SoLFrH FORK

20 a.

lo

&

I

i

i

i

20

40

60

80

Distance

from

source

(kin)

0

b

~ 20

4=0

~ 60

8~0

~ 100

Distance

from

source

(km)

"1 120

12~

----a---

10

E 8 ~

WHAKAPAPA

=

MANGATURUTURU

9

WHANGAEHU

6

0.

c

Distance

i

i

100

200

from

source

(kin)

Lahar volumes vary greatly in magnitude (Fig. 3). In the 1985 Azufrado, Guali, and Molinos/Nereidas lahars of Nevado del Ruiz, triggered by pyroclastic flows, volumes first increased and then decreased (Pierson and others 1990) due to the strong variation in channel slope (Fig. 4). One of the largest eruptions in the last few decades at Mt. Ruapehu, that of 1975, mobilized, 1.8, 0.6, and 0.9 million

Fig. 2a-c. Peak velocities vs distance from source for some past lahars, a 1980 Mt. St. H e l e n s l a h a r s (after Pierson 1985; Pierson a n d Scott 1985; L a e n e n and H a n s e n 1988); b 1985 N e v a d o del Ruiz lahars (after Pierson and others 1990); c 1975 Mt. Ruapehu lahars (after Paterson 1976)

cubic meters of water/sediments along the Whangaeu, Mangaturuturu and Whakapapa valleys, respectively, against a total lake volume of 7 x 1 0 6 m 3 (only a fraction of this water overflowed from the crater, however) (Weir 1982; Vignaux and Weir 1990). In 1980, Mt. St. Helens lahars had volumes ranging from 13.5 and 14 million cubic meters, for Pine Creek-Muddy River and South Fork

40 5000

g" E

3O

WHANGAEHU

WHANGAEHU

........

9

NORTH FOR<

-----c:---

9

SOUTH

9

AZUFRADO

9

GUALI

20

FOF~K

GUALI

--

LAGUNILLAS 3000

MOL./NER.

LAGUNILLAS

==

,,

NORTH FORK

AZUFP&I~

4000

MOLINOS/NEREIDAS

2ooo "6 ;=

10 I

1000 D i

0

1 00 Distance

from

200 source

0

a

100

{kin) Distance

Fig.3. Volumes ofsome recent lahars (after Paterson 1976;Fairchild 1987; Pierson 1985; Pierson and others 1990). T h e v a l u e s o f the 1980 North Fork event must be multiplied by ten

from

200 crater

(kin)

Fig. 4. Longitudinal profiles of some l a h a r c h a n n e l s (after Paterson 1976; Fairchild 1987; Scott 1988; Pierson and others 1990)

144 lahars, respectively, triggered by pyroclastic flows and surges, to 120-140 million cubic meters for the North Fork lahar, via an avalanche deposit of sediments and ice (Fairchild 1987; Pierson 1985; Scott 1988). According to a reconstruction of the 1877 lahar along Rio Cutuchi, triggered by pyroclastic flows on Cotopaxi volcano (Ecuador), a value of about 150 million cubic meters has been estimated (Barberi and others 1992). This paper addresses the possibility of simulating channelized lahars in terms of velocity, discharge, flow depth, arrival times, etc., by an open channel numerical model, through the presentation of some historical cases. Ascertaining this possibility, even with some restrictions (such as constant volume and therefore the absence of significant erosion and deposition effects), is of primary importance, not only because information would be gained on phenomena for which direct measurements are scarce (consider their suddenness and the impassability of certain volcanic zones) but also because it would enable a better evaluation of the hazard in inhabited areas than is possible, following the traditional approach, by the merely historical reconstruction of the lahars of a certain volcano. Changes in the state of the volcano can be considered and thus lead to neither overestimated nor underestimated hazard evaluations. Think, for example, of the case of Cotopaxi volcano (Ecuador) where, over the last 100 years, a considerable reduction of the summit glacier, the main water supply reservoir of the voleano's lahars, has occurred. The lahar expected along the Rio Cutuchi--which runs through densely populated cities before disappearing, after over 150 km, into the Amazon forest--will therefore reasonably have a smaller volume than that of 1877, in the case of an eruption similar to that of 120 years ago. It is on this volume that an estimate of the hazard must be made (Barberi and others 1992).

correction factor, Q = discharge (= UA), h = flow depth, 9 = gravity acceleration, S= = channel slope (S= = tan 0 = bed slope). The method of solving equations 1 and 2 might be similar to that used to solve the dam-break flood wave problem if the shock equations are incorporated with the one-dimensional unsteady debris flow equations (Macedonio and Pareschi 1992). The differences between clear-water flows and debris flows are in the momentum correlation factor, even if the assumption of fl = 1 seems justified, unless there is a rapid change in A, such as in an alluvial fan, where twodimensional depth-averaged flow equations would be more accurate (Chen 1987). Another fundamental difference is in the dissipation term Sr which, for steady water flows, can be expressed, in terms of the average flow velocity U and the hydraulic radius Rn (equal for large channel to flow depth h), by the Manning formula (n is the Manning coefficient): n2U 2

SI = R~ 3

(3)

or by the Chezy formula: 02 SI

-

(4)

C2RH

where C is the Chezy coefficient, which can be expressed as (Smart 1984; Takahashi 1987):

where = a + 5.75 log S~I + 2%)_]

for z. > 0.2

(5)

Unsteady lahars and lahar-runout flows Lahars flowing down a narrow valley, when characterized by a constant sediment concentration and no significant differences in velocity between water and particles, can be described by one-dimensional dynamic equations, analogous to those for clear water (Chen 1987; Macedonio and Pareschi 1992):

OA

+

~Q

=0

(1)

and 8Q aflUQ ~h ~ + ~ + gA Ux = gA(S= - S~)

(2)

in which A = wet cross-section area, t = time, x -- distance along channel, U = average velocity (over A) along x direction, SI = energy dissipation term, fl = momentum

~-=A+

5.751og[0"2(~.~-S~-1) 1

A = 0.04Sx2; A=6.0;

for z. < 0.2

S x > 0.08

S=<0.08

In the formulae, z. = uZ./[(a/p - 1)gd] is the nondimensional tractive force, a and p are the densities of the particles and the fluid, u. is the shear velocity = x/gS=h and f the Darcy-Weisbach friction factor. The factors that affect the theological properties of a debris flow are: concentration of solids, clay content; type of clay; absolute size of the solid material; size distribution of clay, silt, sand, and gravel fraction; characteristics of clastic materials, such as shape, size and density; packing arrangement; proportion of fine-grained material to coarse clastic material such as rocks and boulders; and electrochemical characteristics of the liquid phase. These factors

145

illustrate that it is difficult to develop theoretical equations to predict the rheological properties of flows with natural grain materials. However, ifa simple rheological law z = Zo + yl(dv/dz)~ between shear rate (dv/dz) and shear stress ~ is assumed (Chen 1987), where Zo is the yield strength, Ss can be expressed as (Macedonio and Pareschi 1992):

U2 (7)

for dilatant flows. In the expression for shear stress, T, the flow behavior index, t/, is equal to 1 for a viscous flow and 2 for a dilatant behavior. If t/-- 1, #1 is the plastic viscosity. Laboratory experiments on slurries consisting of < 2 mm sediments from the 1980 North Fork deposit exhibit, for shear rates above 5 s -1, a viscoplastic behavior (linear variation of shear rate with shear stress) or a weakly dilatant one (Major and Pierson 1992). As sand volumetric concentration exceeds 0.2, for shear rates lower than 5 s -1, transient fluctuations and hysteresis effects are observed, probably due to grain rubbing, interlocking, and collision, and a simple rheological model cannot be used. Figure 12 (below) shows viscosity for different fines/sand ratios as a function of volumetric sediment concentration (Major and Pierson 1992); plastic viscosity shows a dependence on sediment concentration of the shape A e Bc, where A and B depend on silt-and-clay to sand ratio and it exhibits order-of-magnitude variation when sediment concentration changes as little as 2-4 percent (Major and Pierson 1992). For a turbulent muddy debris flow (hyperconcentrated flow), where the shear stress is due to particle collision and to the large scale turbulence mixing of the interparticle fluid, Takahashi (1991) proposes: U2

(8)

S I - (ghpZ)

where )_x/l+z+x/Zl;

+ Z --= .~2( a

sin or

\W dfto

"7 \ 101 o"

'c--0.3 1

10 0 t0 0

4 Pg h3 - n 2 ~ -

p=lEln(H+X/l+z

woter flow

(6)

for viscous flows and: 25 #1 U2

o 0.1 o 0.2 9

/

2v

Ss-

10 2

2

(9)

in which H is equal to a o / R , for a hydraulically smooth bed and to bk/h for a rough bed, R . = u . h / v , a o = 1/9.025, b = 1/30, k is the roughness height, v the kinematic viscosity of water, c is the sediment concentration in volume, c. is the maximum packed sediment concentration, a sin ~ = 0.02, 7,. = (a - p)c + p, d is the particle diameter, 2 is

........

t ........ 1 101 10 2 U/U* versus h/d

,

, , , ,~. h/d 103

Fig. 5. Average flow velocity over shear velocity vs flow depth over particle diameter, using Manning formula (dotted line) and formula 8 for a turbulent muddy debris flow (continuous lines). In formula 3, the Manning coefficienthas been computed by formulae 4 and 5. Dots are experimental values measured by Takahashi (1987)

the linear concentration = [ ( c . / c ) 1/3 - 1] -1, and • is the Karman constant; its value is reduced due to the effects of suspended sediments (to = 0.25 for c = 0.3 and ~c = 0.28 for c = 0.1) (Takahashi 1991). Figure 5 shows the resistance to flow in turbulent water flow (Manning law, with n = R~6/C, where C is computed according to formula 5), and turbulent mudflows (equations 8 and 9) compared with experimental values (Takahashi 1987). It can be seen that, as the ratio between flow depth h and particle diameter d increases, at low sediment contents, formulae 3-5 and 8 and 9 give the same behavior for turbulent water and turbulent mudflows.

Simulation of past lahars

In spite of a rich and extensive bibliography on historical lahar occurrences, only for the most recent ones (for example, the 1980-1982 Mt. St. Helens lahars, the 1985 Nevado del Ruiz lahars, the Mt. Ruapehu lahars in recent decades) quantitative information is available on lahar velocities, discharges, density, volumes, and flow depth. Unfortunately, some of these values are often derived from inferences, not from direct measurements, and so estimates could be affected by large errors. Here we discuss the applicability of equations 1 and 2 to the simulation of lahars and lahar-runout flows. In particular, we consider the Muddy River and Pine Creek lahars of 1980, triggered, as already mentioned, by pyroclastic surges and flows during the 18 May 1980 Mt. St. Helens eruption, and the 1975 Mt. Ruapehu lahar along the Whangaehu valley, triggered by the overflowing of the Crater Lake at the top of Mt. Ruapehu, as a consequence of an underwater explosion. For these lahars, the constraints imposed by the model, constant volume and homogeneous flow, were roughly verified along the considered paths. In addition,

146

PINE CREEK 30

MUDDY

Max. hydraulic depth (m) Peak discharge (103 m3/s)

[]

25

200

e .c

RIVER

E

2O

100, t--

d 15

B

B[]r~

[]

[] []

10

[]

9

50

!

8

0

[]

10

,

12

i

9

14

Distance

i

i

16

18

from

9

source

[]

e-

!

i

20

22

-100

(kin)

-200 13.9

MUDDY RIVER [] 9

25-

Max. hydraulic depth (m) I Peakdischarge(10 3 m3/s)

J

23.9 (kin)

from

source

PINE C R E E K

= 20d

1 8.9 Distance

30-

200

15 E

10

[] e

O'

t-

[]

5

[] []

[] 9

9

i

i

i

10

20

30

Distance

from

source

{km)

Fig. 6. Maximum hydraulic depth and peak discharge of 1980 Pine Creek lahar and 1980 Muddy River lahar (after Pierson 1985)

-$ el,CO J~

o

-200 9.9

for those lahars, we were able to collect, from the literature, a sufficient minimum amount of data to try a simulation. Figure 6 shows some values of maximum hydraulic depth and peak discharge along the channel for 1980 Pine Creek and Muddy River lahars (Pierson 1985). Figures 7 and 8 show, respectively, channel widths and channel slopes for Pine Creek and Muddy River. The first sections in the simulation are located, respectively, 9.9 and 13.9 km downstream from Mt. St. Helens crater, where the flow is already channeled and no longer erosive. Figures 9 and 10 show numerical peak discharge and maximum flow depth for Pine Creek and Muddy River. The volumes of the two lahars are 8 and 5.5 million cubic meters (Macedonio and Pareschi 1992). In the figures, field data [values of peak discharges are inferred from indirect observations and can be affected by errors of 15 percent or greater (Pierson 1985)] are also reported. Initial values are Q = 0, h = 0 (dry channel) at all the sections. Two, one, or zero boundary conditions are used at the first and last section, according to subcritical or supercritical flow. The hydrographs in the first section, an input to the simulation, are triangular with peak values from Fig. 6; their durations, Tin,x, have been deduced from known peak discharge and lahar volumes. When an additional boundary condition is needed, flow depth is assigned as derived assuming a linear increasing and then decreasing trend in the interval 0-Tmax, with a maximum value deduced from field observation (15.2 m for

, , , , , 11.9 13.9 15.9 17.9 19.9 D i s t a n c e f r o m s o u r c e (kin)

Fig. 7. Channel width along Pine Creek valley and Muddy River valley (after Pierson 1985)

1400 1200

A

4a----

MUDDYRIVER

9

PINE CREEK

1000

E 800 "13

600 < 400 200 0 0

I

i

I

!

I

i

10

20

30

40

50

60

source

(kin)

Distance

from

Fig. 8. Longitudinal profiles of Pine Creek and Muddy River channels (after Pierson 1985)

147

MUDDY RIVER LAHAR

PINE CREEK LAHAR 30000

E GI

0

m

0

I

I

I

I

I

I

I

I

6" A

E

16 t

e,.

5

"o 4

]: O

#-

lo 1

"- 3 ~4

81

:

8

10

1

:E

1

2

1

12 14 16 18 D i s t a n c e from s o u r c e

20 (kin)

22

12

A

.... -It-.-. m

,IF .... -~-'-'

Field data variable nv

I

14

I

I

16 18 20 22 Distance f r o m s o u r c e

"1'"

24 (kin)

' '

26

Field data variable nv nv=0.1 ml/2 sl/2

nv=0.2 ml/2 sl/2

nv=0.2ml/2 sl/2

nv=0.8 ml/2 sl/2

Fig. 9. Field data (big dots; Pierson 1985) and numerical values of peak discharges vs distance (top) and maximum flow depths vs distance (bottom) from source for 1980 Pine Creek lahar. The simulations were performed with variable n~ as deduced from Table 2 and constant n~ equal to 0.2 and 0.8 ml/2s 1/2, respectively

Table 2. n v calculated from field data for 1980 Pine Creek lahar (Pierson 1985; Macedonio and Pareschi 1992)

Location

Distance from crater (km)

n v (r/= 1) m V2s1/2

P2 P2.1 P3 P4 P5 P6 P7 P8 P9 P11

9.9 10.1 10.8 11.5 12.2 14.1 14.8 16.8 19.5 20.9

1.09 0.70 0.81 0.87 0.93 0.70 0.41 0.4 0.5 0.3

Fig. 10. Field data (big dots; Pierson 1985) and numerical values of peak discharges vs distance (top) and maximum flow depths vs distance (bottom) from source for 1980 Muddy River lahar. The simulations were performed with variable nv as deduced from Table 3 and constant nv equal to 0.1 and 0.2 ml/2s 1/2, respectively

Pine Creek and 4.7 m for Muddy River). Equation 6 is used for the dissipative term, with the values of nv reported in Tables 2 and 3. The values of n~ in the tables are deduced from field information (at peak), by imposing S~ = SI. In Figs. 9 and 10, dotted lines pertain to runs obtained using the values of nv reported in Tables 2 and 3. Extreme cases Table 3. n v calculated from field data for 1980 Muddy River lahar (Pierson 1985; Macedonio and Pareschi 1992)

Location

Distance from crater (km)

nv (~ = 1) m 1/2s1/2

M1 M2 M3 M4 M5 M6

13.9 18.1 21.4 22.7 24.5 28.1

0.17 0.14 0.15 0.09 0.1 0.34

148

with constant values of n~ = 0.1 and 0.2 m 1/2S1/2for Muddy River lahar and n~ = 0.2 and 0.8 ml/2s 1/2 for Pine Creek lahar are also reported. Figure 11 reports the envelopes of cumulative grain-size curves, derived from peak flow deposit (Pierson 1985), and solid concentrations by weight range from 0.84 to 0.86 for Pine Creek and from 0.83 to 0.86 for Muddy River lahar (Pierson 1985). No information on shear rate is available, even if observations on natural mudflows show rates of grain

100

99.5 99 98

P Q. J~

>

E :3 0

10

I

(ram)

diameter

1

I

0.1

i

0.01

ii

~\

:))

shear rarely exceeding 20 s -t, and perhaps are more commonly _<10 s -~ (O'Brien and Julien 1988). For the measured concentrations and fines to sand ratios, values for n~, deduced by extrapolating data from Fig. 12, are in the range 0.05-0.25 mms ~/2, comparable with those reported in Table 3. The greater values of n~ of Table 2 cannot be easily explained, since no macroscopic differences between grain size distributions and sediment concentration exist between the 1980 Pine Creek and Muddy River lahars. A possible explanation is the presence, in Pine Creek, of sharp bends, decreasing downvalley, which could have significantly slowed the lahar (Pierson 1985), that is, in this case, n, is not simply related to a rheological behavior, but also includes the slowing effect of sharp bends (Macedonio and Pareschi 1992). In principle, the

95 90

0.4

80 70 60 50 40 30 20

//~\\~.~ l ~ ~ , Z

l" /.:

10 5

'7,

"--.~

-8

I

I

6

-4

"1

Sx=.l

- .........

sx=2

Muddy River

I

"::.':/ .'.;.i

h

". .,'.'.~ .':.';.'1

2 1 0.5

',

..........

E r

Pine creek

.7

Sx=.01 Sx=.05

Iti/ '~

.7

/...;....;. 9

. , ,, :

(m)

8

Fig. 13. Manning coefficient deduced from formula 5 for different values of flow depth and bed slope I

I

I

I

I

I

"2

0

2

4

6

8

10

grain diameter (r scale)

6000 ~"

b

5000

a

4000 &

3000 2000

Fig. 11. Range of grain size distribution for Pine Creek and Muddy River lahars (after Pierson 1985)

1000 0

100

9

o th [3_

,

9

,

'

I

I

I

i 3.5:1 .5:1 2:1

9

400

1:1

600

Time

~IBSl~m

800

1000

(s)

12000-

10 9 1:4.5 t/I 0 t) t~

200

0

a

10000' / 0 / ~

>

9c

8000

-C

6000

O

U1

4000

n

2000 0.1 0.4 0.5 0.6 0.7 Volumetric sediment concentration Fig. 12. Viscosity at different fines to sand concentrations (after Major and Pierson 1992). Data refer to laboratory experiments on slurries consisting of < 2 - m m sediments from the 1980 North Fork deposit

0

....

I ....

50 Distance

I ....

100 from

I ....

150 source (km)

I

200

Fig. 14a, b. a Four different initial hydrographs used to simulate the 1975 Whangaehu lahar. The source is located 16 km from the Crater Lake. b Lag between front and peak arrival time for 1975 Whangaehu lahar (after Paterson 1976)

149 5000

-

n

=

0.030

4OOO

II

3000

r

L

2000

Q.

&

s

E

12

n

= 0.030

m-1/3

a

e'l

E

m-1/3

1000 0

........

b

..........

c

....................

d

9

s

a

=-

10 8

........

b

..........

c

....................

d

o

field data

E

6

4 i

i

5000 n

=

0.035

4000

m-1/3

....

b

.....

c

..........

d

2000

9

'

'

'

'

i

i

n

=

0.035

v

s

"e

~_

field data

m-1/3

s

a

E

a

3000

a 2 ~- 12-

E

10.

....

b

.....

c

..........

d

8-

6-

I,

field data

4

1000 2'

,'-'5

i

5000-

n = 0,040

4000 " 3000

-

0

i

A 12E

m-1/3 s

a

E

b

i

n

=

0.040 - -

m-1/3

....

b

....

b

.....

c

.....

c

..........

d

..........

d

9

field data

,o

8

o

6"

A

s

a

field data

2000 E 4"

Q

E

1000

=~

0

.

.

.

.

0 Distance

~i

i

50000 from f i r s t

100000 section

20-

I

Fig. 15a-c. Field data (dots; Paterson 1976)and computed peak discharges (continuous lines) for 1975 Whangaehu lahar. Simulations were performed with a n = 0.030 m-1/3s, b n = 0.035 m-~/3s, c n = 0.040 m-~/3s. In each figure, different lines refer to hydrographs a-d of Fig. 14

influence of these sharp bends should be considered by introducing a correction term in the equation of momentum (for example, assuming that a certain fraction of momentum is lost at a sharp bend) rather than by modifying n v. The 1975 lahar along Whangaehu valley (New Zealand) was triggered, as already mentioned, by the overflowing of the Crater Lake at the top of Ruapehu volcano. The volume changed from 1.8 million cubic meters at Karioi, 57 km from the Crater Lake, to 1.6 million cubic meters at K a u a n g a r o a Bridge, 169 km from the Crater Lake (Paterson 1976) and so can be considered approximately constant. The flow was turbulent, black, and oily in appearance (Paterson 1976). Weir (1982) reports that Ruapehu lahars are nonstructural or, equivalently, have a large water content (that is, they are more exactly hyperconcentrated flows or lahar-runout flows). Turbulence and the low sediment content indicate that formula 8 could be

,--1

5O00O

(m) Distance

|tom

first

1OO0O0 section

(m)

Fig. 16a-e. Field data (dots; Paterson 1976) and computed maximum flow depth (continuous lines) for 1975 Whangaehu lahar. Simulations were performed with a n = 0.030 m-1/3s,b n = 0.035 m-U3s, e n = 0.040 m-1/3s. In each figure, different lines refer to hydrographs a d of Fig. 14

more appropriate to estimate Ss than formulae 6 or 7. As seen in Fig. 5, for large values of h/d (flow depth over particle diameter), formula 8 has the same asymptotic values as formula 3, that is the behavior of turbulent water and turbulent mud flows is similar. In the figure, the ratio of flow velocity over shear velocity is computed in the first case by Manning formula and in the second case by equation 8 at different sediment contents (10 and 30 percent by volume). The simulations were thus performed using formula 3, with Manning coefficients, respectively, equal to 0.03, 0.035, and 0.04 m-t/3s. These values can be derived from Fig. 13, where the Manning coefficient, n, is deduced from formula 5, for different values of flow depth and bed slope. The first section in the simulation is located 16 km from the Crater Lake, because upstream the Whangaehu Valley is a wide flood plain (Paterson 1976). In the simulation, due to lack of more precise information, the Whangaehu

150 12-

3oooo 1

n = 0.030 m-1/3, s

........

to

8

~= 20000 ~

.~

4

a.

2-

a

0

d

a

t

O,03Dam_l/3 &

a

"E

lOOOO

[] ,

a

,

__T

n = 0.035

Or-,',

9

,

,

.

.

.

.

field data

,

30000

s I#1

E '~

.

9

2-

b

0

. 35 m-113 s a

field data

b C

4-

g.

=

20000

"E

10000'

d field data

b ,

9

~

'

I

0

I

I

30000-

10

n = 0.040 m-1/3 a ....

E

!i

13

20000 .......... d

?,

6'

9

field dala 10000

4-

d field data

2-

0

,

,

Distance

'

'"'

'

I

50000 from first

I

section

100000 (m)

0 Distance

50000 from first

section

100000 (m)

Fig. 17a-e. Field data (dots; Paterson 1976) and computed peak velocity (continuous lines) for 1975 Whangaehu tartar Simulations were performed with a n = 0.030 m-~/as, b n = 0.035 m-~/~s, e n = 0.040 m-~i3s. In each figure, different lines refer to hydrographs a-d of Fig. 14

Fig. 18a-e. Field data (dots; Paterson 1976) and computed arrival times (continuous lines) for 1975 Whangaehu lahar, Simulations were performed with a n = 0.030 m-1/3s, b n = 0,035 m-t/3s, e n = 0.040 m-l/~s. In each figure, different lines refer to hydrographs a-d of Fig. 14. In the figure, time equal to zero corresponds to the arrival time of the lahar to lhe first section of the simulation

channel has been arbitrarily assumed to be prismatic, with a triangular section 60 m wide at a depth of 5 m, according to an average estimate by Vignaux and Weir (1990), and sparse information by Paterson (1976). Twenty-nine sections have been used, spread along about 125 k m Figure 4 shows Whangaehu bed slopes. No direct information is available on hydrograph shape and maximum peak d~scharge in the first section of the simulation. To deduce these data, the following considerations hold. The estimated total volume was, as already mentioned, 1.8 million cubic meters, while the peak discharge at 6 km from the crater, where flow was not channeled, was 5000 m3/s, and at 44 km (where flow was channeled) it was 580 m3/s (Paterson 1976). Figure 14a curve c shows the initial hydrograph: a value of 3700 ma/s is obtained by a linear interpolation between 5000 and 580 ma/s. In the figure, the duration of the rising phase of about 5 - 6 rain has been obtained by extrapolating information

on the lag between front and peak arrival time (Fig. 14b, after Paterson 1976). Figures 15-18 show measured peak discharges, maximum flow depths, velocities, and arrival times of peak discharges (Paterson t976) compared with numerical results. Some simulations with a different shape of the initial hydrograph have been performed to analyze the sensitivity of the results to this input datum. Runs with a peak discharge of 5000 m 3/s of the hydrograph in the initial section (curve a in Fig. 14a) and peak discharge of 5000 and 3700 m3/s and rising time of 1 rain (curves b - d in Fig. 14a) instead of 350 s show differences in peak velocity, maximum flow depth, and travel times below 3 percent 10 km away from the source (Figs. 15-18). The agreement between field and computed data, here poorly sensitive to changes in the shape of the initial hydrograph, confirms the statement of Weir (1982) that Ruapehu lahar-runout flows may behave similarly to water flows.

151 Discussion and conclusions

Under the assumptions of constant sediment content and volume, an open channel model has been used to simulate both lahars and lahar-runout flows. Due to the variety of factors that affect the rheological properties of a lahar, mainly sediment/water content, average particle dimension, and sand-to-fines ratio, but also shear rate, suitable forms of the energy dissipation term must be adopted, depending on the case. Available data (on velocity, discharge, flow depth, etc.) from historical lahars, briefly reviewed in the introduction, are too poor to try to infer rheological laws. These must be derived from theoretical considerations and validated by laboratory experiments, for which it is possible to realize controlled conditions (stationary regime, regular and constant particle dimensions, constant bed slope, wide channels, etc.). The approach of Chen (1988) and Takahashi (1991) goes in this direction. However, it is important to verify, in natural cases, the laws derived from experiments and to control that effects, such as lobate snout at the flow front, presence of boulders above all at the front, lateral levees of coarse particles, etc., not considered in the model, play a minor role in the dynamics of these phenomena. The simulations presented for lahars in natural beds in two extreme, typical cases--that of a hyperconcentrated, turbulent flow (Ruapehau event of 1975) and that of viscous, much less turbulent flows (Pine Creek and Muddy River lahars of 1980)--seem to confirm that the model used is able to simulate the main features of these events. Additional roles may be played by channel meanderings, whose effects can be evaluated a posteriori by simulating, for a given channel, past events. The simulation of the Ruapehu case seems to suggest that, for lahar-runout flows, the Manning formula (based on a dependence of the friction term on square velocity) for turbulent water gives acceptable results, and the value of the Manning coefficient is obtainable with current

MUDDY RIVER R2=0.94

16-

methods. In these kinds of flow, in fact, as observed by Pierson and Scott (1985), the suspended load suppresses secondary turbulence (small waves, bubbles, froth) and the flow surface has an oily, glassy, smooth appearance, but bed-form activity produces significant large-scale turbulence. Measurements in flow resistances for 1982 hyperconcentrated flow along North Fork River (Mt. St. Helens) were probably not significantly different from those for normal streamflow values (Pierson and Scott 1985). Similar conclusions have also been derived by some authors (Weir 1982, Vignaux and Weir 1990) who developed and applied a kinematic wave model to Ruapehu lahars (of hyperconcentrated type) using a Manning-type law of friction. As the solids content increases, turbulence is further reduced. Visual observations by Li and Luo (1981) and Li and others (1983) on mudflows establish that a dense turbulence mixture occurs only at the head and that all other parts have laminar flow. As deduced from Tables 2 and 3, for Muddy River and some sections of Pine Creek lahar, if we used as viscosity that derived from laboratory experiments on natural lahar deposits of North Fork labar of 1980 at Mt. St. Helens (Major and Pierson 1992), owing to the high solids content, the Reynolds number is probably lower than the critical value for the onset of turbulence. In this case, a linear dependence of the dissipative term Ss on velocity is expected, and used, with acceptable results, to simulate those natural events. The regression of the peak values U vs h2Sx, deduced for the 1980 Muddy River lahar from field information at locations of Table 3 (Pierson 1985), shows a good R 2 (= 0.94), as reported in Fig. 19. It is important to emphasize that for large sediment content, also a dilatant behavior [i.e., a dependence of shear stress on square shear rate (du/dz)] could be, in principle, expected. According to Takahashi (1991) a modified Bagnold number N m = [ a 2 2 d 2 ( d u / d z ) ] / l l r , where/~r is the total viscosity, addresses the relative importance of the inertial stress (dilatant behavior) to the viscous stress (viscous behavior). For lahars where small size particles are a large fraction of the solids content and therefore the average particle dimension is of the order of a few millimeters or less, N,, could be smaller than 450, the experimental threshold value for the dilatant behavior.

14-

Acknowledgments We thank two anonymous reviewers provided helpful criticisms of earlier versions of this paper.

~-, 1 2 -

'=ElOj 86

References

4 2

0.0

i

i

0.1

0.2 Sx h 2

i

i

i

0.3

0.4

0.5

(rrl 2 )

Fig. 19. Peak velocities plotted with respect t o h2Sx (here h is the maximum flow depth, equivalent, for wide channels, to the hydraulic radius) for Muddy River lahar of 1980 (Pierson 1985)9The regression equation used, characteristics of laminar uniform flows (equation 2, with Sx = Sy, t?/t?x = 0, 0/& = 0) displays a good R2 (= 0.94)

Barberi F, Caruso P, Macedonio G, Pareschi MT and Rosi M (1992) Reconstruction and numerical simulation of the lahar of the 1877 eruption of Cotopaxi volcano (Ecuador). Acta Vulcanol 2:35-44 Chen C (1987) Comprehensive review of debris flow modelling concepts in Japan. Rev Eng Geol 7:13-29 Chen C (1988) General solution of the viscoplastic debris flow. J Hydrual Eng 114:259-282 Cummans J (1981a) Mudflows resulting from the May 18, 1980, eruption of Mt. St. Helens, Washington. US Geological Survey Circular 850-B. 16 pp

152 Cummans J (198 lb) Chronology ofmudflows on the South Fork and North Fork Toutle River following the May 18 eruption. In: Lipmann PW and Mullineaux DR (Eds), The 1980 Eruption of Mount St. Helens, Washington. US Geological Survey Professional Paper 1250. pp 479-486 Fairchild LH (1987) The importance of lahar initiation process. In Costa JE and Wieczorek GF (Eds), Debris flows/avalanches: Process recognition and mitigation. Rev Eng Geol 7:51-61 Janda RJ, Scott KM, Nolan KM and Martinson HA (1981) Lahar movement, effects, and deposition. In Lipmann PW and Mullineaux DR (Eds), The 1980 Eruption of Mount St. Helens, Washington. US Geological Survey Professional Paper 1250. pp 461 478 Laenen A and Hansen RP (1988) Simulation of three lahars in the Mount St. Helens area, Washington using a one-dimensional, unsteady-state streamflow model. US Geological Survey, WaterResource Investigations Report 88-4004. 20 pp Li J and Luo D (1981) The formation and characteristics of mudflow and flood in the mountain area of the Dachao River and its prevention. Z Geomorphol 25 :470-484 Li J, Yuan J, Bi C and Luo D (1983) The main features of the mudflow in Jiang-Jia Ravine. Z Geomorphol 27:325-341 Macedonio G and Pareschi MT (1992) Numerical simulation of some lahars from Mt. St. Helens. J Volcanol Geotherm Res (in press) Major JJ and Pierson TC (1992) Debris flow rheology: Experimental analysis of fine-grained slurries. Water Resour Res 28(3): 841-857 Neall VE (1976) Lahars as major geological hazards. Bull Int Assoc Eng Geol 14:233-240 Paterson BR (1976) The effects of lahar from the 1975 April Mt Ruapehu eruption and the threat of future eruption on Tongariro power development. NZ Geological Survey Unpublished Engineering Geology Report EG 230

Pierson TC (1985) Initiation and flow behavior of the Pine Creek and Muddy River lahars, Mt. St. Helens, Washington. Geol Soc Am Bull 96:1056 1069 Pierson TC (1986) Flow behavior of channelized debris flows, Mount St. Helens, Washington. In Abrahams AD (Ed), Hillslope processes. Boston: Allen & Unwin. pp 269-296 Pierson TC and Scott KM (1985) Downstream dilution of a lahar: Transition from debris flow to hyperconcentrated streamflow. Water Resour Res 21(10): 1511-1524 Pierson TC, Janda RJ, Thouret JC and Borrero CA (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. J Volcanol Geotherm Res 41:17 66 Scott KM (1988) Origin, behaviour, and sedimentology oflahars and lahar-runout flows in the Toutle-Cowlitz River system. US Geological Survey Professional Paper 1447-A. 74 pp Smart GM (1984) Sediment transport formula for steep channels. J Hydraul Div ASCE 110(3):267-276 Takahashi T (1987) High velocity flow in steep erodible channels. In Fluvial Hydraulics, Proceedings of the 22nd IAHR Congress, Lausanne, Switzerland. pp 42-53 Takahashi T (1991) Debris flow. Balkema AA (ed.), Rotterdam: IAHR Monograph Series. 165 pp Vignaux M and Weir GJ (1990) A general model for Mt. Ruapehu lahars. Bull Volcanol 52:381-390 Weir JG (1982) Kinematic wave theory for Ruapehu lahars. NZ J Sci 25:197 203 Yokoyama I, Tilling RI and Scarpa R (1984) International mobile early-warning systems(s) for volcanic eruptions and related seismic activities. Paris: UNESCO FP/2106-82-01(2286). 102 pp