A MODEL OF GRANULAR FLOWS OVER AN ERODIBLE SURFACE ...
DISCRETE AND CONTINUOUS DYNAMICAL SYSTEMS–SERIES B Volume 3, Number 4, November 2003
Website: http://AIMsciences.org pp. 589–599
A MODEL OF GRANULAR FLOWS OVER AN ERODIBLE SURFACE E.B. PITMAN, C.C. NICHITA A.K. PATRA, A.C. BAUER M. BURSIK AND A. WEBB
Department of Mathematics, Department of Mechanical and Aerospace Engineering, Department of Geology, University at Buffalo, Buffalo, NY 14260
To Dave Schaeffer, who shows how mathematics can impact science. Abstract. We present a framework for modeling a dry geophysical mass of granular material – a debris or volcanic avalanche or landslide – flowing over an erodible surface. We also describe a computing environment that incorporates topographical data into a parallel, adaptive mesh computational algorithm that solves the model equations.
1. Introduction. Slow moving geophysical mass flows – debris flows, block and ash flows, volcanic avalanches, or landslides – may be initiated by volcanic activity. In these flows, constituent particles are typically centimeter to meter sized, and the flows, sometimes as fast as tens of meters per second, propagate distances of tens of kilometers. As these flows slow, the particle mass sediments out, yielding deposits that can be a hundred meters deep and many kilometers in length [9]. These flows can contain O(106 − 1010 ) m3 or more of material, and many also include a significant amount of water. The range of scales and the complexity of the rheology for geological materials, coupled with the mathematical problem of describing and computing a free surface flow, is a significant challenge. We lack a full understanding of how these mass flows are initiated, but there is a growing understanding of processes governing flow, once that motion has been initiated. This paper describes one effort in modeling and simulating geophysical mass flows. Our efforts try to strike a balance between fidelity to physics on the one hand and mathematical and computational tractability on the other. The starting point for the modeling effort here is the pioneering work of Savage and Hutter [21]. Our contribution is to incorporate into the modeling the effect of erosion on the dynamics of dry, flowing granular materials. In a first approximation we assume that changes in elevation of the basal surface due to erosion are small, and we solve the equations of motion with the bed surface frozen. For the largescale flows that are our primary interest, this assumption seems justified, but the assumption may not be accurate when comparing model results against small tabletop experiments. 1991 Mathematics Subject Classification. 86A60, 86-08, 74H15. Key words and phrases. Geological problems, Computational methods, Numerical approximation of solutions.
589
590
PITMAN, NICHITA, PATRA, BAUER, BURSIK AND WEBB
The governing equations consist of a system of hyperbolic balance laws. Using a parallel, adaptive mesh computational framework, simulations were performed and results tested against data from small scale experiments of granular avalanches on inclined plane [15], and against data from a debris flow event at the Little Tahoma Peak on Mt. Rainier, Washington [24]. The next section contains a derivation of our model equations. We then describe the TITAN computing environment used to simulate flows. We close with computational results and comparisons with experimental data. 2. Modeling Issues. In this section we summarize one approach to modeling geophysical mass flows. Although interstitial fluid is an important ingredient in many real flows, we concentrate here on dry granular flows without fluid. The approach we adopt was originally introduced by Savage and Hutter [21]. Starting from equations of mass and momentum balance on an inclined plane, using a Coulomb constitutive relation and making scaling arguments, they developed a one-dimensional model system similar to the shallow water equations. They also performed a series of table-top experiments, to test model predictions. Savage and Hutter later extended their work to two dimensions, found a similarity solution to the governing equations, and examined flows over more general basal surfaces [6]. Iverson [7, 8, 4] extended the modeling to include the effects of pore fluid, and employed modern shock-capturing numerical techniques to solve the model system. Gray [5] re-derived the governing equations for general basal surfaces. Pitman [18] wrote a Godunov solver for one-dimensional flows that incorporates general topographical features into a simulation. A heuristic derivation of the model equations proceeds in a manner similar to the derivation of the shallow water equations. A thin layer of granular material at constant density is assumed to flow over a smooth basal surface; motion normal to the bed is neglected. By depth-averaging the momentum balance laws and the equation of incompressibility, one derives a system of PDEs for the flow depth h and the flow velocities in the down-slope and cross-slope directions. There are three significant differences distinguishing geophysical mass flow models from shallow water. First, the assumption of Mohr-Coulomb plasticity is far too complex to use in a depth-averaged model. A simplifying assumption is made, namely that all shear and longitudinal stresses are proportional to the normal (lithostatic) stress; similar assumptions appear in soil mechanics and have their roots in the classic work of Terzaghi [25] and Rankine [19]. Second, basal friction and basal surface curvature lead to drag terms that are more complicated than the frictional drag in shallow water models. Finally, the equations introduced below include a model of bed erosion. It is estimated that as much as 50% of a typical flow deposit is comprised of eroded material; this paper offers one approach to capturing this important physical effect. A complementary modeling effort can be found in [11], which uses the hydrodynamiclike constitutive theory for rapid granular flows, also in the context of a thin layer assumption. As an aside, it is interesting to note that a time-dependent model of an incompressible granular material with a Mohr-Coulomb constitutive relation – the starting point for the modeling – has been shown to be unconditionally ill-posed, similar to a mixed forward-backward heat equation qt = qxx − qyy [22]. The governing equations consist of a system of hyperbolic balance laws. To solve these equations we have developed a numerical simulation framework that imports
GRANULAR FLOWS OVER AN ERODIBLE SURFACE
591
digital elevation data, solves the governing equations, and assists in visualization. A brief description of the code will be presented below; the interested reader can find a more detailed description of this parallel, adaptive mesh, finite volume Godunov solver elsewhere [15, 16]. The model equations and computational environment is used to study the effect of erosion on the dynamics of flowing granular masses. By assuming that changes in the basal surface due to erosion are small, we ignore changes to the surface slopes and curvature that erosion might induce, and we do not alter the digital elevation map that defines the terrain. Thus the net effect of erosion consists only in changing the net mass and momentum of the flow. Because of the accuracy of widely available elevation data, this assumption seems justified. The erosion model has been tested against experimental and field data [2], which the reader can consult for details. 3. Governing Equations. In this section we summarize the system of equations governing our erosion model. For a detailed derivation of these equations the reader is referred to [2]. The fundamental balance laws for mass and momentum for an incompressible continuum may be written ∇·v = 0
ρ(∂t v + v · ∇v) + ∇ · T = ρg
(1)
Here v is the velocity, ρ is the (constant) density of the granular material, and T is the stress tensor. The granular material is assumed to be an incompressible continuum satisfying a Mohr-Coulomb constitutive relation. The Mohr-Coulomb assumption states that slip planes appear inside the bulk when the internal state of stress equals the yield criteria, σt /σn = tan ϕint , where σn , σt are the normal and shear stresses acting on a plane element inside the granular material and ϕint is the internal friction angle of the medium. We assume that the interface between flowing granular material and the basal surface z = b(x, t) can move as a result of an exchange of sediment particles. The normal speed of the interface ∂t b + vx ∂x b + vy ∂y b − vz is given by an erosion rate es . We postulate an empirical formula for this rate, assuming there is a threshold shear stress below which there is no erosion, and that the amount of eroded material increases monotonically with the bed shear stress above this threshold. Changes to the surface are assumed small, and are ignored when determining surface slope and curvature. A friction law is assumed at the interface between the flowing mass and the basal r surface: Tb nb − nb (nb · Tb nb ) = |vvr | (nb · Tb nb ) tan ϕbed . Here vr is the slip velocity (r denoting the velocity component), at the upper side of an infinitesimally thin boundary layer that forms at the basal contact surface, Tb is the stress at this basal surface, nb the surface normal, and ϕbed the bed friction angle. Although it might be argued that, if the surface is eroding, frictional stresses would be dramatically reduced, we assume that the sliding friction relationship is not changed during the erosion process. The upper free surface z = h(x, t), is assumed to be stress-free, and a kinematic boundary condition ∂t h + vx ∂x h + vy ∂y h − vz = 0 is imposed on the motion of the surface.
592
PITMAN, NICHITA, PATRA, BAUER, BURSIK AND WEBB
The three dimensional system as described gives a detailed description of the flow, but the computational complexity of such a model, especially in the treatment of the free surfaces, precludes its use in our applications. Before delving into the details of the model derivation, we mention that the Equations (1) could be scaled in both dependent and independent variables, giving insight into the relative magnitude of terms is the governing system. Such a scaling lies at the heart of the original derivation of Savage and Hutter [21, 6]. In particular, the thin layer model equations we solve are predicated on a long wave assumption. That is, changes in the downslope direction are assumed to occur on a lengthscale L, while changes in the direction normal to the basal surface on a scale H, and we assume H/L = ǫ ≪ 1. Our derivation below does not explicitly rely on this scaling, but the long wave assumption does influence the choices made during the modeling. The interested reader should also see [8]. Depth-averaging is applied to an arbitrary element in a local coordinate system that has the Oz axis directed normal to the bed surface, and Ox in the direction of maximum local slope. In this coordinate system the equations for the basal surface and for the free surface are z = 0, and z = h(x, y, t), respectively. We now briefly illustrate the depth-averaging procedure. Integrate the incompressibility equation through the layer to find Z h(x,y,t) (∂x vx + ∂y vy + ∂z vz ) dz = 0 . 0
The Leibnitz formula allows the interchange of differentiation and integration. Then apply the kinematic conditions on the upper and lower surfaces to find (2) ∂t h + ∂x (hv x ) + ∂y (hv y ) = es . Rh Rh The depth-averaged velocities are defined by hv x = 0 vx dz, and hv y = 0 vy dz. q 2 2 The erosion rate es is an empirical factor, es = α (hv x ) + (hv y ) /h if (T σ ≥ T 0 ), or es = 0 if (T σ < T 0 ). In these formulae, T σ is the total shear stress at the given point and T 0 is a threshold stress, and α is a proportionality constant to be fitted to experimental results. Turning now to the momentum equations, the x momentum equation first, integrate in the normal direction again using the Leibnitz formula to find, £ ¢ ¡ ¤ h ρ ∂t (hvx ) + ∂x hv 2x + ∂y (hv x v y ) − ρ [vx (∂t z + vx ∂x z + vy ∂y z − vz )]0 = h
−∂x (hT xx ) − ∂y (hT yx ) + [Txx ∂x z + Tyx ∂y z − Tzx ]0 + ρgx h
(3)
The overbar on T denotes a depth-averaged stress component, and a term like gx is the component of gravity in the x-direction. The kinematic and stress condition at the upper surfaces reads ¯ ¯ ¯ ¯ T · n ¯z=h(x,y,t) = Txx ¯z=h(x,y,t) ∂x h + Tyx ¯z=h(x,y,t) ∂y h − Tzx ¯z=h(x,y,t) = 0
Combining this with the basal stress condition in the depth-averaged momentum equation yields: £ ¢ ¡ ¤ ρ ∂t (hv x ) + ∂x hv 2x + ∂y (hv x v y ) = ρes vr − ∂x (hT xx ) − ∂y (hT yx ) + Tzx |z=0 + ρgx h
(4)
Experimental observations suggests that the downslope velocity profile is blunt, and all shear is confined to a very thin layer near the basal surface. Based on these
GRANULAR FLOWS OVER AN ERODIBLE SURFACE
593
notions, we assume the sliding velocity vr can be approximated by the average velocity v x . Because of the shallowness assumption, we assume a lithostatic relation for the normal stresses in the z direction. The Mohr-Coulomb then theory provides a highly nonlinear relation among the stresses, a theory is too complex for use here. Instead we employ a relation from soil mechanics that has its roots in the work of Terzaghi [25] and Rankine [19] – assume the longitudinal and shear stresses are proportional to the normal stress, with constant proportionality: Txy = ρhgz x Tyx = −sgn( ∂v ∂y )kap 2 sin ϕint ; similar relations hold for the xz- yz-, xx- and yystresses. In this formula, the stress coefficient kap differs in the active and passive stress states. The active or passive state of stress is developed if an element of material is elongated or compressed, and the formula for the corresponding states can be derived from the Mohr diagram [23] kap
£ ¤1/2 1 ± 1 − cos2 ϕint (1 + tan2 ϕbed ) −1 =2 cos2 ϕint
(5)
in which the negative square root corresponds to an active state (∂v x /∂x+∂v y /∂y > 0), and the positive root to the passive state (∂v x /∂x + ∂v y /∂y < 0). Combining all these equations, the depth averaged x-momentum equation can be written µ ¶ 1 2 2 ∂t (hv x ) + ∂x hv x + kap gz h + ∂y (hv x v y ) = es vx + gx h 2 · µ ¶¸ vx vx q gz h 1 + −hkap sgn(∂v x /∂y)∂y (gz h) sin ϕint − tan ϕbed (6) rx gz v 2x + v 2y Here rx is the radius of curvature in the x-direction of the basal surface. The relation for the y-momentum equation is similar, and can be obtained by interchanging x and y in the equation above. Writing the relations for mass and momentum balance, the governing equations form a system of hyperbolic balance laws that can be written in matrix form ∂U ∂F ∂G + + =S ∂t ∂x ∂y
(7)
where U = (h, hvx , hvy )t , F = (hvx , hvx2 +0.5kap gz h2 , hvx vy )t , G = (hvy , hvx vy , hvy2 + 0.5kap gz h2 )t , and S = (es , Sx , Sy )t and where ¶¸ · µ vx vx tan ϕbed Sx = es vx +gx h−hkap sgn(∂vx /∂y)∂y (gz h) sin ϕint − q gz h 1 + rx gz v 2x + v 2y ¶¸ · µ vy tan ϕbed gz h 1 + Sy = es vy +gy h−hkap sgn(∂vy /∂x)∂x (gz h) sin ϕint − q ry gz v 2x + v 2y vx
This system has a structure similar to the shallow water equations, except for the complexity of the dissipation terms on the right-hand side. The system is strictly hyperbolic and genuinely nonlinear away from the “vacuum state” h = 0, with √ characteristic speeds λ± = u ± βh, where β = kap .
594
PITMAN, NICHITA, PATRA, BAUER, BURSIK AND WEBB
4. TITAN2D code. TITAN2D is a parallel adaptive-mesh finite volume computational environment based on a Godunov solver, for solving the governing PDEs Eqns (7). A detailed description of the code can be found in [15]; here we give a brief review of important features of the code. TITAN2D begins by importing digital elevation data from the GRASS data (Geographic Resource Analysis Support System) geographic information systems (GIS) data source. This GRASS data is then geo-rectified and coded into a grid, in a manner similar to the GRID module of ARC/INFO [1]. From this data, the elevation slopes and curvatures, which are parameters in the governing equations, are computed and a regular Cartesian grid on which Eqns. (7) will be solved is defined. The equations are solved by a finite-volume method using a Godunov method in an explicit time-integration scheme. The basic approach employed here is due originally to Davis [3]. The dependent variables are considered as cell-averages, and their values are advanced by a predictor-corrector method; the source terms are included in these updates, and no splitting – neither for the source nor the multiple dimensions – is necessary. ¶ µ Gi,j+1/2 − Gi,j−1/2 Fi+1/2,j − Fi−1/2,j n+1 n (8) + + Si,j ui,j = ui,j − ∆t ∆x ∆y In this update formula, both the numerical fluxes F, G and the source terms S are evaluated at the mid-time n + 21 ∆t values, which are obtained in a predictor step: ∆t (9) (Ai,j ∆x u + Bi,j ∆y u + Si,j ) 2 Here, A and B are the Jacobians of the physical flux functions F and G, respectively, and ∆x u and ∆y u are slope-limited approximations to the x- and y-derivatives, respectively. The numerical fluxes are then evaluated at the cell edges using an approximate Riemann solver. We have experimented with Riemann solvers due to Davis, Roe, Russanov, and the HLL solver [3, 10, 27], and all give similar results. Testing suggests that a correct handling of the snout and tail of the pile is the most important of the several factors influencing the computations. In particular, robust results rely on sufficiently fine grids and a proper accommodation of waves propagating into the h = 0 vacuum at the front and back of the pile. Although numerical difficulties often accompany hyperbolic systems with source terms, the approach just described has provided acceptable results especially when used with fine enough grids [15] (see also [27] for results from shallow water computations, and [29] regarding so-called positivity preserving schemes). In addition, the simplicity and small operation count of the method recommend its use. This numerical method is integrated into a parallel implementation coupled with adaptive meshing, allowing for efficient processing of large data sets without degrading numerical accuracy. In particular, a special scheme based on a space-filling curve is used to index grid cells, an important feature when cells are refined and unrefined [14, 15]. n+ 21
ui,j
= uni,j −
5. Results and Experiments. We have tested our model equations (without erosion) and numerical solver against results from table-top experiments and against field measurements of historical flows over natural terrain; see, for example, [15, 16, 24]. Here we report on numerical results of the erosion model system, and compare these against flow on a simple inclined plane flow. We also show a comparison of flow with and without erosion at the Little Tahoma Peak on Mt. Ranier, Washington.
GRANULAR FLOWS OVER AN ERODIBLE SURFACE No Erosion 1.4
Flow front distance from origin(meters) Flow tail distance from origin(meters) Maximum pile height(scaling factor 35 mm) Maximum velocity (meters/s)
1.2
1
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0
(A)
Erosion(alpha=0.1,H0=0.01) 1.4
Flow front distance from origin(meters) Flow tail distance from origin(meters) Maximum pile height(scaling factor 35 mm) Maximum velocity (meters/s)
1.2
0
0.1
0.2
0.3
0.4
595
0.5
Time(s)
0.6
0.7
0.8
0.9
0
(B)
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Time(s)
Figure 1. Location of pile front and tail, and maximum pile height and velocity, for sand flowing down an inclined plane in numerical experiments. A) Erosion is not active, while B) erosion with a threshold of 0.1 is active. Consider then the dynamics of a pile of sand, about 400 g in mass, released on an inclined plane. The experimental plane is about 12 meter wide and 1 meter long, made of a composite particle board with coarse sandpaper glued on the top. At a distance of 0.8 meters from the top, the incline abruptly meets a horizontal plane. We have performed two sets of experiments, one with an erodible surface, one without erosion. For runs with erosion, we have covered the plane surface with a layer of sand about 1 cm. deep. Because the friction angle between the sand and the sandpaper is large, this layer remains stable for a range of inclination angles at which we run experiments. We investigate velocities at a series of time steps, as well as runout distance, sediment deposition, and final mass of sediment load. Physical experiments and detailed comparisons are reported in [2]. Figure (1) gives results from TITAN simulations; Figure(2) presents results from the experiment. In Figure (1) we show the position of the nose and tail of the flowing pile as it descends an incline, as well as velocity and pile height during the run. Plots A and B compare results of simulations without erosion and with erosion activated. In the erosion simulation, the threshold for the initiation of erosion is 0.1 - that is, 10% of the maximum initial pressure - and the erosion rate α = 0.1. In Figure(2), the location of the front of the pile is plotted, as is the maximum speed of the pile. For the physical experiments, the internal friction is 36.6 degrees, the basal friction 29 deg. The plane inclined at an angle of 44.3 degrees in the run with no erosion, and 42.5 deg. in the run with erosion. In the numerical runs, the initial pile was ellipsoidal, with a maximum height scaled to unity. The various friction and inclination angles were set to match the experimental values. Experimental and numerical flows show similar overall characteristics, and all stop rapidly soon after the incline meets the horizontal. The maximum pile height in the experiments was 1 cm in the case of no erosion, and 1.5 cm with erosion (a 50% difference). These values compare with final maximum pile heights of 0.7 cm in the case of no erosion, and 0.8 cm with erosion for the numerical results (a 15% difference). As intuition might suggest, maximum pile height in both laboratory and numerical runs is somewhat larger for the case of erosion, but not significantly. The most obvious difference and perhaps the most important difference between the
596
(A)
PITMAN, NICHITA, PATRA, BAUER, BURSIK AND WEBB
(B) Figure 2. Location of pile front, and maximum speed, for sand in laboratory experiments flowing down an inclined plane. (A) Erosion is not active, while (B) erosion is active.
cases of erosion and no erosion is in maximum downstream speed, which is 15−20% larger without erosion in both laboratory and numerical experiments. In comparing experimental and numerical results, perhaps the biggest discrepancy is the time scale on which flow occurs. In the numerical simulations with and without erosion, the front of the pile reaches the end of the plane - and effectively halts - after about 0.45 sec. On the other hand, in experiments this takes about 0.6 sec. Comparison between the numerical and experimental finding of the maximum flow speed are not well established. It does appear that, in the experiments, the maximum speed without erosion is perhaps as much as 20% larger than with erosion, consistent with the simulation results. But these experiments did not provide sufficient detail to determine the time course of the maximum velocity profile, particularly in the no erosion experiments. There are several shortcomings with this analysis, including differences in friction and slope angles among the comparison datasets, but the results suggest that the model mimics several characteristics of real flows, at least real, small-scale flows. We turn attention now to simulations of a 1963 rockfall avalanche at the Little Tahoma Peak. Field work suggests that a flank of this peak gave way during a series of avalanches over several weeks, with a total of about 107 m3 of material falling. Numerical results here simulate one large avalanche flowing over the terrain as measured in recent GIS data. At the top of the mountainside is Emmons Glacier, while lower down the basal surface is made of geologic material. This changing terrain presents a challenge to our simulations - we have not yet included a variable basal friction in the code. For simulations reported here, we use a small basal friction angle, intermediate between that for geologic material (∼ 30 degrees) and ice (perhaps ∼ 5−10 degrees). We show the peak velocities computed in simulations with and without erosion in Figure 3. The peak velocity is seen to be about 10% smaller when erosion is activated. In Figure 4 we show contours of flow simulations at several times down Little Tahoma, with and without erosion. The pile height during these simulations are differ, especially during the latter stages of the flow. The final deposit with erosion is about twice as large as when erosion is not active. During the simulations, the moving pile sweeps out an area that is in good agreement with field reconstruction of the deposit from the ’63 events. The simulations
GRANULAR FLOWS OVER AN ERODIBLE SURFACE
597
Figure 3. Maximum velocity of a simulated flow down Little Tahoma Peak, with and without erosion. Erosion is activated whenever the shear stress is greater than 1% of the initial gravitational pressure.
Figure 4. Contour of the flowing mass down Little Tahoma Peak at three different times. For each pair of figures, in the frame on the left erosion is absent, and erosion is active in the frame on the right, with a threshold of 1% of the peak initial stress.
show qualitative agreement with the flow history and run-out, as estimated from field studies, but the match is not quantitatively right. Nonetheless, the comparison is encouraging; see [24] for further details.
6. Summary. We provide a summary of a mathematical system of PDEs describing dry geophysical mass flows over realistic, erodible surfaces. We then briefly describe the TITAN2D computer code for solving this model system, and offer a comparison of simulations results with experiments. There are significant modeling issues yet to be addressed, such as the presence of pore fluid in the moving mass. Nonetheless, for dry flows, preliminary analysis suggests TITAN provides a computational environment with adequate fidelity to the physical flows it intends to simulate; experiments and computational validation is continuing. The code is designed to help scientists and civil protection authorities assess the risk of and mitigate hazards due to dry debris and avalanche flows. The the core numerical solver can be adapted as more complete rheological models are developed.
598
PITMAN, NICHITA, PATRA, BAUER, BURSIK AND WEBB
We now have in place a computational environment for simulating the system of equations modeling single-phase geophysical mass flows, including a model of erosion, and a compute environment that integrates accurate topographical data into a parallel, adaptive mesh Godunov solver. Simulation results match experimental findings qualitatively, but detailed quantitative comparison, especially largescale fieldwork against which to compare simulations, is still required. The MohrCoulomb theory which lies at the heart of the physical modeling presented here requires only a small number of physical parameters (internal and basal friction angles, and two erosion parameters), these parameters are only loosely defined, and the measurements of these quantities are crude. The thin layer modeling, and our numerical solution methodology is capable of extension to constitutive theories other than the Mohr-Coulomb relation used here. One may rightly ask “Does the specific constitutive theory employed in such a thin layer model qualitatively change simulation results?”. It is also important to recognize that, by necessity, any “thin layer” model hides important information through the averaging process. For instance, one open problem is to recover information about the internal flow dynamics. The numerical method described here, and similar approaches by others [4, 5], are (or can be) formally second order accurate. The sheer size of the regions on which we wish to simulate flow, however, necessitates large grid cells. Adaptive gridding and a parallel computing framework provides the capacity to perform computations on reasonably sized grid cells. Our experience shows that a grid with sufficiently small cells is necessary for reasonable simulations. At the same time, elevation data must be accurate enough to allow these fine grids to be faithful to the actual topography. In some cases, this topographical data is available down to 5m × 5m-scale with a ±1m vertical accuracy. A remaining challenge is to develop a strategy for computational mesh refinement coupled to a hierarchy of refined topographical mappings. Finally, the model used here accounts only for the deformation of solids material. Missing is the effect of interstitial fluid, especially water mixed into the granular material. Pore fluid may also be important in initiating mass flows. Iverson [7, 8] has taken a first step by in including fluid into a thin layer model, but more work is required before we can be confident of any approach. 7. Acknowledgments. This work has been supported by NSF grants ITR 0121254 and EHR 0087665. Computing support from the Center for Computational Research at the University at Buffalo is gratefully acknowledged. REFERENCES [1] ArcInfo http://www.esri.com [2] Bursik, M.I., Nichita, C.C., and Pitman, E.B. Geophysical Research Letters, in preparation. [3] Davis S.F., (1988) “Simplified second order Godunov type methods”, SIAM J. Sci. Stat. Comp. 9 445-473. [4] Denlinger, R.P and Iverson, R.M (2001) “Flow of variably fluidized granular material across three-dimensional terrain 2. Numerical predictions and experimental tests”, J. Geophys. Res. 106 553-566 [5] Gray, J.M.N.T., Wieland, M. and Hutter, K. (1999) “Gravity-driven free surface flow of granular avalanches over complex basal terrain” Proc. Roy. Soc. London Serial A 455 18411874 [6] Hutter, K., Siegel, M., Savage, S.B., and Nohguchi, Y. (1993) “Two dimensional spreading of a granular avalanche down an inclined plane; Part 1: Theory” Acta Mech. 100 37-68 [7] Iverson, R.M. (1997) “The physics of debris flows”, Rev. Geophys. 35 245-296
GRANULAR FLOWS OVER AN ERODIBLE SURFACE
599
[8] Iverson, R.M. and Denlinger, R.P, (2001) “Flow of variably fluidized granular material across three-dimensional terrain 1. Coulomb mixture theory”, J. Geophys. Res. 106 537-552 [9] Legros, F. (2002) “The mobility of long runout landslides”, Eng. Geology 63 301-331 [10] Leveque, R. J. (2002) Finite Volume Methods for Hyperbolic Problems Cambridge University Press [11] Louge, M.Y. and Keast, S. “On dense granular flows down flat frictional inclines” (2001)Phys. Fluids 13 1213-1233 [12] National Research Council (1991) “The eruption of Nevado del Ruiz Volcano, Colombia, South America, November 13, 1985”, The Committee on Natural Disasters [13] National Research Council (1994) “Mount Ranier: Active Cascade Volcano”, The US Geodynamics Committee [14] Patra, A. Laszloffy, A. and Long, J. “Data Structures and Load Balancing for Parallel Adaptive hp Finite Element Methods”, submitted Comp. Method. in App. Mech. and Engineering. [15] Patra, A.K., Bauer, A.C., Nichita, C.C., Pitman, E.B., Sheridan, M.F., Bursik, M., Rupp, B., Webber, A., Stinton, A., Namikawa, L., and Renschler, C. “Parallel Adaptive Numerical Simulation of Dry Avalanches over Natural Terrain” submitted J. Volcanology and Geothermal Research [16] Pitman, E.B. , Patra, A.K., Bauer, A.C., Nichita, C.C., Sheridan, M. and Bursik, M. “Computing Debris Flows” accepted, Physics of Fluids [17] Pitman, E.B. Patra, A.,. Bauer, A., Sheridan, M. and Bursik, M. “Computing Debris Flows and Landslides: Towards a Tool for Hazard Risk Evaluation”, to appear, Proc. Hyp ’02 T. Hou and E. Tadmor, eds. [18] Pitman, E. B. and Tulyakov, S. (1998) “Granular debris flows in 1 dimension” UB Math Department report [19] Rankine, W.J.M. (1857) “On the stability of loose earth”, Philos. Trans. Roy. Soc. London 147 9-27. [20] Rodriguez-Elizarraras, S., Siebe, C., Komorowski, J.-C., Espindola, J.M., and Saucedo, R. (1991) “Field observations of pristine block-and-ash flow deposits emplaced April 16-17, 1991 at Volcan de Colima, Mexico” J. Volcanology and Geothermal Research 48 399-412. [21] Savage, S.B. and Hutter, K. (1989) “The motion of a finite mass of granular material down a rough incline”, J. Fluid Mech. 199 177-215 [22] Schaeffer, D.G. (1987) “Instability in the evolution equations describing granular flow”, J. Diff. Eqn. 66 19-50. [23] Schofield, A.N. and Wroth, C.P. (1968), Critical State Soil Mechanics McGraw-Hill, London. [24] Sheridan, M.F., Stinton, A.J., Patra, A.K., Pitman, E.B., Bauer,A.C., and Nichita, C.C. “Evaluating TITAN2D mass-flow model using the 1963 Little Tahoma Peak avalanches, Mount Rainier, Washington”, submitted Journal of Volcanology and Geothermal Research. [25] Terzaghi, K. (1936) “The shearing resistance of saturated soils and the angle between planes of shear”, Proc. 1st Int. Conf. Soil Mech. 54-56 [26] Tiling, R. I., “Volcanic Hazards”, American Geophysical Union Short Course, Volume I, 1989. [27] Toro, E. F. (2001) Shock Capturing Methods for Free Surface Shallow Flows Wiley and Sons [28] Vallance, J. W. (1999) “Post-glacial lahars and potential hazards in the White Salmon River system on the Southwest flank of Mt. Adams, Washington” USGS Bulletin 2161 [29] Wang, C. and Liu, J.-G. (2003) “Positivity property of second order flux-splitting schemes for the compressible Euler equations”, Disc. Cont. Dynam. Sys. 3 201-228
Received March 2003; revised July 2003. E-mail address: [email protected], [email protected]; [email protected], [email protected], [email protected]; [email protected]
Website: http://AIMsciences.org pp. 589–599
A MODEL OF GRANULAR FLOWS OVER AN ERODIBLE SURFACE E.B. PITMAN, C.C. NICHITA A.K. PATRA, A.C. BAUER M. BURSIK AND A. WEBB
Department of Mathematics, Department of Mechanical and Aerospace Engineering, Department of Geology, University at Buffalo, Buffalo, NY 14260
To Dave Schaeffer, who shows how mathematics can impact science. Abstract. We present a framework for modeling a dry geophysical mass of granular material – a debris or volcanic avalanche or landslide – flowing over an erodible surface. We also describe a computing environment that incorporates topographical data into a parallel, adaptive mesh computational algorithm that solves the model equations.
1. Introduction. Slow moving geophysical mass flows – debris flows, block and ash flows, volcanic avalanches, or landslides – may be initiated by volcanic activity. In these flows, constituent particles are typically centimeter to meter sized, and the flows, sometimes as fast as tens of meters per second, propagate distances of tens of kilometers. As these flows slow, the particle mass sediments out, yielding deposits that can be a hundred meters deep and many kilometers in length [9]. These flows can contain O(106 − 1010 ) m3 or more of material, and many also include a significant amount of water. The range of scales and the complexity of the rheology for geological materials, coupled with the mathematical problem of describing and computing a free surface flow, is a significant challenge. We lack a full understanding of how these mass flows are initiated, but there is a growing understanding of processes governing flow, once that motion has been initiated. This paper describes one effort in modeling and simulating geophysical mass flows. Our efforts try to strike a balance between fidelity to physics on the one hand and mathematical and computational tractability on the other. The starting point for the modeling effort here is the pioneering work of Savage and Hutter [21]. Our contribution is to incorporate into the modeling the effect of erosion on the dynamics of dry, flowing granular materials. In a first approximation we assume that changes in elevation of the basal surface due to erosion are small, and we solve the equations of motion with the bed surface frozen. For the largescale flows that are our primary interest, this assumption seems justified, but the assumption may not be accurate when comparing model results against small tabletop experiments. 1991 Mathematics Subject Classification. 86A60, 86-08, 74H15. Key words and phrases. Geological problems, Computational methods, Numerical approximation of solutions.
589
590
PITMAN, NICHITA, PATRA, BAUER, BURSIK AND WEBB
The governing equations consist of a system of hyperbolic balance laws. Using a parallel, adaptive mesh computational framework, simulations were performed and results tested against data from small scale experiments of granular avalanches on inclined plane [15], and against data from a debris flow event at the Little Tahoma Peak on Mt. Rainier, Washington [24]. The next section contains a derivation of our model equations. We then describe the TITAN computing environment used to simulate flows. We close with computational results and comparisons with experimental data. 2. Modeling Issues. In this section we summarize one approach to modeling geophysical mass flows. Although interstitial fluid is an important ingredient in many real flows, we concentrate here on dry granular flows without fluid. The approach we adopt was originally introduced by Savage and Hutter [21]. Starting from equations of mass and momentum balance on an inclined plane, using a Coulomb constitutive relation and making scaling arguments, they developed a one-dimensional model system similar to the shallow water equations. They also performed a series of table-top experiments, to test model predictions. Savage and Hutter later extended their work to two dimensions, found a similarity solution to the governing equations, and examined flows over more general basal surfaces [6]. Iverson [7, 8, 4] extended the modeling to include the effects of pore fluid, and employed modern shock-capturing numerical techniques to solve the model system. Gray [5] re-derived the governing equations for general basal surfaces. Pitman [18] wrote a Godunov solver for one-dimensional flows that incorporates general topographical features into a simulation. A heuristic derivation of the model equations proceeds in a manner similar to the derivation of the shallow water equations. A thin layer of granular material at constant density is assumed to flow over a smooth basal surface; motion normal to the bed is neglected. By depth-averaging the momentum balance laws and the equation of incompressibility, one derives a system of PDEs for the flow depth h and the flow velocities in the down-slope and cross-slope directions. There are three significant differences distinguishing geophysical mass flow models from shallow water. First, the assumption of Mohr-Coulomb plasticity is far too complex to use in a depth-averaged model. A simplifying assumption is made, namely that all shear and longitudinal stresses are proportional to the normal (lithostatic) stress; similar assumptions appear in soil mechanics and have their roots in the classic work of Terzaghi [25] and Rankine [19]. Second, basal friction and basal surface curvature lead to drag terms that are more complicated than the frictional drag in shallow water models. Finally, the equations introduced below include a model of bed erosion. It is estimated that as much as 50% of a typical flow deposit is comprised of eroded material; this paper offers one approach to capturing this important physical effect. A complementary modeling effort can be found in [11], which uses the hydrodynamiclike constitutive theory for rapid granular flows, also in the context of a thin layer assumption. As an aside, it is interesting to note that a time-dependent model of an incompressible granular material with a Mohr-Coulomb constitutive relation – the starting point for the modeling – has been shown to be unconditionally ill-posed, similar to a mixed forward-backward heat equation qt = qxx − qyy [22]. The governing equations consist of a system of hyperbolic balance laws. To solve these equations we have developed a numerical simulation framework that imports
GRANULAR FLOWS OVER AN ERODIBLE SURFACE
591
digital elevation data, solves the governing equations, and assists in visualization. A brief description of the code will be presented below; the interested reader can find a more detailed description of this parallel, adaptive mesh, finite volume Godunov solver elsewhere [15, 16]. The model equations and computational environment is used to study the effect of erosion on the dynamics of flowing granular masses. By assuming that changes in the basal surface due to erosion are small, we ignore changes to the surface slopes and curvature that erosion might induce, and we do not alter the digital elevation map that defines the terrain. Thus the net effect of erosion consists only in changing the net mass and momentum of the flow. Because of the accuracy of widely available elevation data, this assumption seems justified. The erosion model has been tested against experimental and field data [2], which the reader can consult for details. 3. Governing Equations. In this section we summarize the system of equations governing our erosion model. For a detailed derivation of these equations the reader is referred to [2]. The fundamental balance laws for mass and momentum for an incompressible continuum may be written ∇·v = 0
ρ(∂t v + v · ∇v) + ∇ · T = ρg
(1)
Here v is the velocity, ρ is the (constant) density of the granular material, and T is the stress tensor. The granular material is assumed to be an incompressible continuum satisfying a Mohr-Coulomb constitutive relation. The Mohr-Coulomb assumption states that slip planes appear inside the bulk when the internal state of stress equals the yield criteria, σt /σn = tan ϕint , where σn , σt are the normal and shear stresses acting on a plane element inside the granular material and ϕint is the internal friction angle of the medium. We assume that the interface between flowing granular material and the basal surface z = b(x, t) can move as a result of an exchange of sediment particles. The normal speed of the interface ∂t b + vx ∂x b + vy ∂y b − vz is given by an erosion rate es . We postulate an empirical formula for this rate, assuming there is a threshold shear stress below which there is no erosion, and that the amount of eroded material increases monotonically with the bed shear stress above this threshold. Changes to the surface are assumed small, and are ignored when determining surface slope and curvature. A friction law is assumed at the interface between the flowing mass and the basal r surface: Tb nb − nb (nb · Tb nb ) = |vvr | (nb · Tb nb ) tan ϕbed . Here vr is the slip velocity (r denoting the velocity component), at the upper side of an infinitesimally thin boundary layer that forms at the basal contact surface, Tb is the stress at this basal surface, nb the surface normal, and ϕbed the bed friction angle. Although it might be argued that, if the surface is eroding, frictional stresses would be dramatically reduced, we assume that the sliding friction relationship is not changed during the erosion process. The upper free surface z = h(x, t), is assumed to be stress-free, and a kinematic boundary condition ∂t h + vx ∂x h + vy ∂y h − vz = 0 is imposed on the motion of the surface.
592
PITMAN, NICHITA, PATRA, BAUER, BURSIK AND WEBB
The three dimensional system as described gives a detailed description of the flow, but the computational complexity of such a model, especially in the treatment of the free surfaces, precludes its use in our applications. Before delving into the details of the model derivation, we mention that the Equations (1) could be scaled in both dependent and independent variables, giving insight into the relative magnitude of terms is the governing system. Such a scaling lies at the heart of the original derivation of Savage and Hutter [21, 6]. In particular, the thin layer model equations we solve are predicated on a long wave assumption. That is, changes in the downslope direction are assumed to occur on a lengthscale L, while changes in the direction normal to the basal surface on a scale H, and we assume H/L = ǫ ≪ 1. Our derivation below does not explicitly rely on this scaling, but the long wave assumption does influence the choices made during the modeling. The interested reader should also see [8]. Depth-averaging is applied to an arbitrary element in a local coordinate system that has the Oz axis directed normal to the bed surface, and Ox in the direction of maximum local slope. In this coordinate system the equations for the basal surface and for the free surface are z = 0, and z = h(x, y, t), respectively. We now briefly illustrate the depth-averaging procedure. Integrate the incompressibility equation through the layer to find Z h(x,y,t) (∂x vx + ∂y vy + ∂z vz ) dz = 0 . 0
The Leibnitz formula allows the interchange of differentiation and integration. Then apply the kinematic conditions on the upper and lower surfaces to find (2) ∂t h + ∂x (hv x ) + ∂y (hv y ) = es . Rh Rh The depth-averaged velocities are defined by hv x = 0 vx dz, and hv y = 0 vy dz. q 2 2 The erosion rate es is an empirical factor, es = α (hv x ) + (hv y ) /h if (T σ ≥ T 0 ), or es = 0 if (T σ < T 0 ). In these formulae, T σ is the total shear stress at the given point and T 0 is a threshold stress, and α is a proportionality constant to be fitted to experimental results. Turning now to the momentum equations, the x momentum equation first, integrate in the normal direction again using the Leibnitz formula to find, £ ¢ ¡ ¤ h ρ ∂t (hvx ) + ∂x hv 2x + ∂y (hv x v y ) − ρ [vx (∂t z + vx ∂x z + vy ∂y z − vz )]0 = h
−∂x (hT xx ) − ∂y (hT yx ) + [Txx ∂x z + Tyx ∂y z − Tzx ]0 + ρgx h
(3)
The overbar on T denotes a depth-averaged stress component, and a term like gx is the component of gravity in the x-direction. The kinematic and stress condition at the upper surfaces reads ¯ ¯ ¯ ¯ T · n ¯z=h(x,y,t) = Txx ¯z=h(x,y,t) ∂x h + Tyx ¯z=h(x,y,t) ∂y h − Tzx ¯z=h(x,y,t) = 0
Combining this with the basal stress condition in the depth-averaged momentum equation yields: £ ¢ ¡ ¤ ρ ∂t (hv x ) + ∂x hv 2x + ∂y (hv x v y ) = ρes vr − ∂x (hT xx ) − ∂y (hT yx ) + Tzx |z=0 + ρgx h
(4)
Experimental observations suggests that the downslope velocity profile is blunt, and all shear is confined to a very thin layer near the basal surface. Based on these
GRANULAR FLOWS OVER AN ERODIBLE SURFACE
593
notions, we assume the sliding velocity vr can be approximated by the average velocity v x . Because of the shallowness assumption, we assume a lithostatic relation for the normal stresses in the z direction. The Mohr-Coulomb then theory provides a highly nonlinear relation among the stresses, a theory is too complex for use here. Instead we employ a relation from soil mechanics that has its roots in the work of Terzaghi [25] and Rankine [19] – assume the longitudinal and shear stresses are proportional to the normal stress, with constant proportionality: Txy = ρhgz x Tyx = −sgn( ∂v ∂y )kap 2 sin ϕint ; similar relations hold for the xz- yz-, xx- and yystresses. In this formula, the stress coefficient kap differs in the active and passive stress states. The active or passive state of stress is developed if an element of material is elongated or compressed, and the formula for the corresponding states can be derived from the Mohr diagram [23] kap
£ ¤1/2 1 ± 1 − cos2 ϕint (1 + tan2 ϕbed ) −1 =2 cos2 ϕint
(5)
in which the negative square root corresponds to an active state (∂v x /∂x+∂v y /∂y > 0), and the positive root to the passive state (∂v x /∂x + ∂v y /∂y < 0). Combining all these equations, the depth averaged x-momentum equation can be written µ ¶ 1 2 2 ∂t (hv x ) + ∂x hv x + kap gz h + ∂y (hv x v y ) = es vx + gx h 2 · µ ¶¸ vx vx q gz h 1 + −hkap sgn(∂v x /∂y)∂y (gz h) sin ϕint − tan ϕbed (6) rx gz v 2x + v 2y Here rx is the radius of curvature in the x-direction of the basal surface. The relation for the y-momentum equation is similar, and can be obtained by interchanging x and y in the equation above. Writing the relations for mass and momentum balance, the governing equations form a system of hyperbolic balance laws that can be written in matrix form ∂U ∂F ∂G + + =S ∂t ∂x ∂y
(7)
where U = (h, hvx , hvy )t , F = (hvx , hvx2 +0.5kap gz h2 , hvx vy )t , G = (hvy , hvx vy , hvy2 + 0.5kap gz h2 )t , and S = (es , Sx , Sy )t and where ¶¸ · µ vx vx tan ϕbed Sx = es vx +gx h−hkap sgn(∂vx /∂y)∂y (gz h) sin ϕint − q gz h 1 + rx gz v 2x + v 2y ¶¸ · µ vy tan ϕbed gz h 1 + Sy = es vy +gy h−hkap sgn(∂vy /∂x)∂x (gz h) sin ϕint − q ry gz v 2x + v 2y vx
This system has a structure similar to the shallow water equations, except for the complexity of the dissipation terms on the right-hand side. The system is strictly hyperbolic and genuinely nonlinear away from the “vacuum state” h = 0, with √ characteristic speeds λ± = u ± βh, where β = kap .
594
PITMAN, NICHITA, PATRA, BAUER, BURSIK AND WEBB
4. TITAN2D code. TITAN2D is a parallel adaptive-mesh finite volume computational environment based on a Godunov solver, for solving the governing PDEs Eqns (7). A detailed description of the code can be found in [15]; here we give a brief review of important features of the code. TITAN2D begins by importing digital elevation data from the GRASS data (Geographic Resource Analysis Support System) geographic information systems (GIS) data source. This GRASS data is then geo-rectified and coded into a grid, in a manner similar to the GRID module of ARC/INFO [1]. From this data, the elevation slopes and curvatures, which are parameters in the governing equations, are computed and a regular Cartesian grid on which Eqns. (7) will be solved is defined. The equations are solved by a finite-volume method using a Godunov method in an explicit time-integration scheme. The basic approach employed here is due originally to Davis [3]. The dependent variables are considered as cell-averages, and their values are advanced by a predictor-corrector method; the source terms are included in these updates, and no splitting – neither for the source nor the multiple dimensions – is necessary. ¶ µ Gi,j+1/2 − Gi,j−1/2 Fi+1/2,j − Fi−1/2,j n+1 n (8) + + Si,j ui,j = ui,j − ∆t ∆x ∆y In this update formula, both the numerical fluxes F, G and the source terms S are evaluated at the mid-time n + 21 ∆t values, which are obtained in a predictor step: ∆t (9) (Ai,j ∆x u + Bi,j ∆y u + Si,j ) 2 Here, A and B are the Jacobians of the physical flux functions F and G, respectively, and ∆x u and ∆y u are slope-limited approximations to the x- and y-derivatives, respectively. The numerical fluxes are then evaluated at the cell edges using an approximate Riemann solver. We have experimented with Riemann solvers due to Davis, Roe, Russanov, and the HLL solver [3, 10, 27], and all give similar results. Testing suggests that a correct handling of the snout and tail of the pile is the most important of the several factors influencing the computations. In particular, robust results rely on sufficiently fine grids and a proper accommodation of waves propagating into the h = 0 vacuum at the front and back of the pile. Although numerical difficulties often accompany hyperbolic systems with source terms, the approach just described has provided acceptable results especially when used with fine enough grids [15] (see also [27] for results from shallow water computations, and [29] regarding so-called positivity preserving schemes). In addition, the simplicity and small operation count of the method recommend its use. This numerical method is integrated into a parallel implementation coupled with adaptive meshing, allowing for efficient processing of large data sets without degrading numerical accuracy. In particular, a special scheme based on a space-filling curve is used to index grid cells, an important feature when cells are refined and unrefined [14, 15]. n+ 21
ui,j
= uni,j −
5. Results and Experiments. We have tested our model equations (without erosion) and numerical solver against results from table-top experiments and against field measurements of historical flows over natural terrain; see, for example, [15, 16, 24]. Here we report on numerical results of the erosion model system, and compare these against flow on a simple inclined plane flow. We also show a comparison of flow with and without erosion at the Little Tahoma Peak on Mt. Ranier, Washington.
GRANULAR FLOWS OVER AN ERODIBLE SURFACE No Erosion 1.4
Flow front distance from origin(meters) Flow tail distance from origin(meters) Maximum pile height(scaling factor 35 mm) Maximum velocity (meters/s)
1.2
1
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0
(A)
Erosion(alpha=0.1,H0=0.01) 1.4
Flow front distance from origin(meters) Flow tail distance from origin(meters) Maximum pile height(scaling factor 35 mm) Maximum velocity (meters/s)
1.2
0
0.1
0.2
0.3
0.4
595
0.5
Time(s)
0.6
0.7
0.8
0.9
0
(B)
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Time(s)
Figure 1. Location of pile front and tail, and maximum pile height and velocity, for sand flowing down an inclined plane in numerical experiments. A) Erosion is not active, while B) erosion with a threshold of 0.1 is active. Consider then the dynamics of a pile of sand, about 400 g in mass, released on an inclined plane. The experimental plane is about 12 meter wide and 1 meter long, made of a composite particle board with coarse sandpaper glued on the top. At a distance of 0.8 meters from the top, the incline abruptly meets a horizontal plane. We have performed two sets of experiments, one with an erodible surface, one without erosion. For runs with erosion, we have covered the plane surface with a layer of sand about 1 cm. deep. Because the friction angle between the sand and the sandpaper is large, this layer remains stable for a range of inclination angles at which we run experiments. We investigate velocities at a series of time steps, as well as runout distance, sediment deposition, and final mass of sediment load. Physical experiments and detailed comparisons are reported in [2]. Figure (1) gives results from TITAN simulations; Figure(2) presents results from the experiment. In Figure (1) we show the position of the nose and tail of the flowing pile as it descends an incline, as well as velocity and pile height during the run. Plots A and B compare results of simulations without erosion and with erosion activated. In the erosion simulation, the threshold for the initiation of erosion is 0.1 - that is, 10% of the maximum initial pressure - and the erosion rate α = 0.1. In Figure(2), the location of the front of the pile is plotted, as is the maximum speed of the pile. For the physical experiments, the internal friction is 36.6 degrees, the basal friction 29 deg. The plane inclined at an angle of 44.3 degrees in the run with no erosion, and 42.5 deg. in the run with erosion. In the numerical runs, the initial pile was ellipsoidal, with a maximum height scaled to unity. The various friction and inclination angles were set to match the experimental values. Experimental and numerical flows show similar overall characteristics, and all stop rapidly soon after the incline meets the horizontal. The maximum pile height in the experiments was 1 cm in the case of no erosion, and 1.5 cm with erosion (a 50% difference). These values compare with final maximum pile heights of 0.7 cm in the case of no erosion, and 0.8 cm with erosion for the numerical results (a 15% difference). As intuition might suggest, maximum pile height in both laboratory and numerical runs is somewhat larger for the case of erosion, but not significantly. The most obvious difference and perhaps the most important difference between the
596
(A)
PITMAN, NICHITA, PATRA, BAUER, BURSIK AND WEBB
(B) Figure 2. Location of pile front, and maximum speed, for sand in laboratory experiments flowing down an inclined plane. (A) Erosion is not active, while (B) erosion is active.
cases of erosion and no erosion is in maximum downstream speed, which is 15−20% larger without erosion in both laboratory and numerical experiments. In comparing experimental and numerical results, perhaps the biggest discrepancy is the time scale on which flow occurs. In the numerical simulations with and without erosion, the front of the pile reaches the end of the plane - and effectively halts - after about 0.45 sec. On the other hand, in experiments this takes about 0.6 sec. Comparison between the numerical and experimental finding of the maximum flow speed are not well established. It does appear that, in the experiments, the maximum speed without erosion is perhaps as much as 20% larger than with erosion, consistent with the simulation results. But these experiments did not provide sufficient detail to determine the time course of the maximum velocity profile, particularly in the no erosion experiments. There are several shortcomings with this analysis, including differences in friction and slope angles among the comparison datasets, but the results suggest that the model mimics several characteristics of real flows, at least real, small-scale flows. We turn attention now to simulations of a 1963 rockfall avalanche at the Little Tahoma Peak. Field work suggests that a flank of this peak gave way during a series of avalanches over several weeks, with a total of about 107 m3 of material falling. Numerical results here simulate one large avalanche flowing over the terrain as measured in recent GIS data. At the top of the mountainside is Emmons Glacier, while lower down the basal surface is made of geologic material. This changing terrain presents a challenge to our simulations - we have not yet included a variable basal friction in the code. For simulations reported here, we use a small basal friction angle, intermediate between that for geologic material (∼ 30 degrees) and ice (perhaps ∼ 5−10 degrees). We show the peak velocities computed in simulations with and without erosion in Figure 3. The peak velocity is seen to be about 10% smaller when erosion is activated. In Figure 4 we show contours of flow simulations at several times down Little Tahoma, with and without erosion. The pile height during these simulations are differ, especially during the latter stages of the flow. The final deposit with erosion is about twice as large as when erosion is not active. During the simulations, the moving pile sweeps out an area that is in good agreement with field reconstruction of the deposit from the ’63 events. The simulations
GRANULAR FLOWS OVER AN ERODIBLE SURFACE
597
Figure 3. Maximum velocity of a simulated flow down Little Tahoma Peak, with and without erosion. Erosion is activated whenever the shear stress is greater than 1% of the initial gravitational pressure.
Figure 4. Contour of the flowing mass down Little Tahoma Peak at three different times. For each pair of figures, in the frame on the left erosion is absent, and erosion is active in the frame on the right, with a threshold of 1% of the peak initial stress.
show qualitative agreement with the flow history and run-out, as estimated from field studies, but the match is not quantitatively right. Nonetheless, the comparison is encouraging; see [24] for further details.
6. Summary. We provide a summary of a mathematical system of PDEs describing dry geophysical mass flows over realistic, erodible surfaces. We then briefly describe the TITAN2D computer code for solving this model system, and offer a comparison of simulations results with experiments. There are significant modeling issues yet to be addressed, such as the presence of pore fluid in the moving mass. Nonetheless, for dry flows, preliminary analysis suggests TITAN provides a computational environment with adequate fidelity to the physical flows it intends to simulate; experiments and computational validation is continuing. The code is designed to help scientists and civil protection authorities assess the risk of and mitigate hazards due to dry debris and avalanche flows. The the core numerical solver can be adapted as more complete rheological models are developed.
598
PITMAN, NICHITA, PATRA, BAUER, BURSIK AND WEBB
We now have in place a computational environment for simulating the system of equations modeling single-phase geophysical mass flows, including a model of erosion, and a compute environment that integrates accurate topographical data into a parallel, adaptive mesh Godunov solver. Simulation results match experimental findings qualitatively, but detailed quantitative comparison, especially largescale fieldwork against which to compare simulations, is still required. The MohrCoulomb theory which lies at the heart of the physical modeling presented here requires only a small number of physical parameters (internal and basal friction angles, and two erosion parameters), these parameters are only loosely defined, and the measurements of these quantities are crude. The thin layer modeling, and our numerical solution methodology is capable of extension to constitutive theories other than the Mohr-Coulomb relation used here. One may rightly ask “Does the specific constitutive theory employed in such a thin layer model qualitatively change simulation results?”. It is also important to recognize that, by necessity, any “thin layer” model hides important information through the averaging process. For instance, one open problem is to recover information about the internal flow dynamics. The numerical method described here, and similar approaches by others [4, 5], are (or can be) formally second order accurate. The sheer size of the regions on which we wish to simulate flow, however, necessitates large grid cells. Adaptive gridding and a parallel computing framework provides the capacity to perform computations on reasonably sized grid cells. Our experience shows that a grid with sufficiently small cells is necessary for reasonable simulations. At the same time, elevation data must be accurate enough to allow these fine grids to be faithful to the actual topography. In some cases, this topographical data is available down to 5m × 5m-scale with a ±1m vertical accuracy. A remaining challenge is to develop a strategy for computational mesh refinement coupled to a hierarchy of refined topographical mappings. Finally, the model used here accounts only for the deformation of solids material. Missing is the effect of interstitial fluid, especially water mixed into the granular material. Pore fluid may also be important in initiating mass flows. Iverson [7, 8] has taken a first step by in including fluid into a thin layer model, but more work is required before we can be confident of any approach. 7. Acknowledgments. This work has been supported by NSF grants ITR 0121254 and EHR 0087665. Computing support from the Center for Computational Research at the University at Buffalo is gratefully acknowledged. REFERENCES [1] ArcInfo http://www.esri.com [2] Bursik, M.I., Nichita, C.C., and Pitman, E.B. Geophysical Research Letters, in preparation. [3] Davis S.F., (1988) “Simplified second order Godunov type methods”, SIAM J. Sci. Stat. Comp. 9 445-473. [4] Denlinger, R.P and Iverson, R.M (2001) “Flow of variably fluidized granular material across three-dimensional terrain 2. Numerical predictions and experimental tests”, J. Geophys. Res. 106 553-566 [5] Gray, J.M.N.T., Wieland, M. and Hutter, K. (1999) “Gravity-driven free surface flow of granular avalanches over complex basal terrain” Proc. Roy. Soc. London Serial A 455 18411874 [6] Hutter, K., Siegel, M., Savage, S.B., and Nohguchi, Y. (1993) “Two dimensional spreading of a granular avalanche down an inclined plane; Part 1: Theory” Acta Mech. 100 37-68 [7] Iverson, R.M. (1997) “The physics of debris flows”, Rev. Geophys. 35 245-296
GRANULAR FLOWS OVER AN ERODIBLE SURFACE
599
[8] Iverson, R.M. and Denlinger, R.P, (2001) “Flow of variably fluidized granular material across three-dimensional terrain 1. Coulomb mixture theory”, J. Geophys. Res. 106 537-552 [9] Legros, F. (2002) “The mobility of long runout landslides”, Eng. Geology 63 301-331 [10] Leveque, R. J. (2002) Finite Volume Methods for Hyperbolic Problems Cambridge University Press [11] Louge, M.Y. and Keast, S. “On dense granular flows down flat frictional inclines” (2001)Phys. Fluids 13 1213-1233 [12] National Research Council (1991) “The eruption of Nevado del Ruiz Volcano, Colombia, South America, November 13, 1985”, The Committee on Natural Disasters [13] National Research Council (1994) “Mount Ranier: Active Cascade Volcano”, The US Geodynamics Committee [14] Patra, A. Laszloffy, A. and Long, J. “Data Structures and Load Balancing for Parallel Adaptive hp Finite Element Methods”, submitted Comp. Method. in App. Mech. and Engineering. [15] Patra, A.K., Bauer, A.C., Nichita, C.C., Pitman, E.B., Sheridan, M.F., Bursik, M., Rupp, B., Webber, A., Stinton, A., Namikawa, L., and Renschler, C. “Parallel Adaptive Numerical Simulation of Dry Avalanches over Natural Terrain” submitted J. Volcanology and Geothermal Research [16] Pitman, E.B. , Patra, A.K., Bauer, A.C., Nichita, C.C., Sheridan, M. and Bursik, M. “Computing Debris Flows” accepted, Physics of Fluids [17] Pitman, E.B. Patra, A.,. Bauer, A., Sheridan, M. and Bursik, M. “Computing Debris Flows and Landslides: Towards a Tool for Hazard Risk Evaluation”, to appear, Proc. Hyp ’02 T. Hou and E. Tadmor, eds. [18] Pitman, E. B. and Tulyakov, S. (1998) “Granular debris flows in 1 dimension” UB Math Department report [19] Rankine, W.J.M. (1857) “On the stability of loose earth”, Philos. Trans. Roy. Soc. London 147 9-27. [20] Rodriguez-Elizarraras, S., Siebe, C., Komorowski, J.-C., Espindola, J.M., and Saucedo, R. (1991) “Field observations of pristine block-and-ash flow deposits emplaced April 16-17, 1991 at Volcan de Colima, Mexico” J. Volcanology and Geothermal Research 48 399-412. [21] Savage, S.B. and Hutter, K. (1989) “The motion of a finite mass of granular material down a rough incline”, J. Fluid Mech. 199 177-215 [22] Schaeffer, D.G. (1987) “Instability in the evolution equations describing granular flow”, J. Diff. Eqn. 66 19-50. [23] Schofield, A.N. and Wroth, C.P. (1968), Critical State Soil Mechanics McGraw-Hill, London. [24] Sheridan, M.F., Stinton, A.J., Patra, A.K., Pitman, E.B., Bauer,A.C., and Nichita, C.C. “Evaluating TITAN2D mass-flow model using the 1963 Little Tahoma Peak avalanches, Mount Rainier, Washington”, submitted Journal of Volcanology and Geothermal Research. [25] Terzaghi, K. (1936) “The shearing resistance of saturated soils and the angle between planes of shear”, Proc. 1st Int. Conf. Soil Mech. 54-56 [26] Tiling, R. I., “Volcanic Hazards”, American Geophysical Union Short Course, Volume I, 1989. [27] Toro, E. F. (2001) Shock Capturing Methods for Free Surface Shallow Flows Wiley and Sons [28] Vallance, J. W. (1999) “Post-glacial lahars and potential hazards in the White Salmon River system on the Southwest flank of Mt. Adams, Washington” USGS Bulletin 2161 [29] Wang, C. and Liu, J.-G. (2003) “Positivity property of second order flux-splitting schemes for the compressible Euler equations”, Disc. Cont. Dynam. Sys. 3 201-228
Received March 2003; revised July 2003. E-mail address: [email protected], [email protected]; [email protected], [email protected], [email protected]; [email protected]
