A MacCormack-TVD finite difference method to ... - Semantic Scholar
Computers & Geosciences 52 (2013) 1–10
Contents lists available at SciVerse ScienceDirect
Computers & Geosciences journal homepage: www.elsevier.com/locate/cageo
A MacCormack-TVD finite difference method to simulate the mass flow in mountainous terrain with variable computational domain Chaojun Ouyang a,b,c, Siming He a,b,n, Qiang Xu c, Yu Luo a,b, Wencheng Zhang a,b a
Key Laboratory of Mountain Hazards and Surface Process, Chinese Academy of Science, Chengdu 610041, China Institute of Mountain Hazards and Environment (IMHE), Chinese Academy of Sciences, Chengdu 610041, China c National Laboratory of Geohazard Prevention and Geoenvironment Protection, Chengdu University of Technology, Chengdu 610059, China b
a r t i c l e i n f o
abstract
Article history: Received 28 March 2012 Received in revised form 25 August 2012 Accepted 27 August 2012 Available online 25 September 2012
A two-dimensional mountainous mass flow dynamic procedure solver (Massflow-2D) using the MacCormack-TVD finite difference scheme is proposed. The solver is implemented in Matlab on structured meshes with variable computational domain. To verify the model, a variety of numerical test scenarios, namely, the classical one-dimensional and two-dimensional dam break, the landslide in Hong Kong in 1993 and the Nora debris flow in the Italian Alps in 2000, are executed, and the model outputs are compared with published results. It is established that the model predictions agree well with both the analytical solution as well as the field observations. & 2012 Elsevier Ltd. All rights reserved.
Keywords: Finite difference method MacCormack scheme TVD scheme Dam break Landslide Debris flow
1. Introduction Global climate change has brought along with it a fair share of disastrous weather conditions that have also increased the potential for mountainous hazards such as landslides, debris flows, dam-breaks, and floods in natural terrains. Such events can cause substantial damage to human life and engineering facilities. Therefore, hazard prediction and risk evaluation of hazardous zones are of considerable importance and have rightfully attracted the attention of many scientists. The common computational methods used to evaluate the disaster risk range for a potential hazard can be classified into three categories: empirical formula methods, discrete element methods and continuum modeling methods. The major advantage of any empirical method is its simplicity, while the disadvantage is the high demand for field or experimental data and ignorance of the guiding mechanisms (Prochaska et al., 2008; Rickenmann, 1999; Schilling and Iverson, 1997). The discrete element method allows for the quantitative reproduction of a particle trajectory, but by the very nature of such calculations, the method is timeconsuming and requires micromechanical input parameters that are often difficult to determine (Cundall, 2001). The continuum modeling method considers the moving mass to be a fluid and can
n Corresponding author at: Institute of Mountain Hazards and Environment (IMHE), Chinese Academy of Sciences, Chengdu 610041, China. Fax: þ 86 28 85222258. E-mail address: [email protected] (S. He).
0098-3004/$ - see front matter & 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.cageo.2012.08.024
therefore be depicted by the incompressible Navier–Stokes equations. These equations can then be further reduced to the two-dimensional (2D) shallow water equations by a depthaveraged procedure. Many numerical methods, including the finite volume method (Berger et al., 2011; Christen et al., 2010; Leveque, 2002; Pitman et al., 2003), the finite difference method (Beguerı´a et al., 2009; Liang et al., 2006; Wang and Sassa, 2010) and the finite element method (Chen and Lee, 2003; Quecedo et al., 2004), have been proposed to resolve the sharp discontinuities associated with the hyperbolic problems, such as the 2D shallow water equations. Nevertheless, simulation of a large-range and long-duration mountainous hazard on a natural terrain with less time and higher resolution remains a challenge. The MacCormack finite difference scheme with total variation diminishing (TVD) has been investigated by many authors and has been shown to be a feasible and robust computational solution (Tseng and Chu, 2000; Vincent et al., 2001; Wang et al., 2000). This method was further modified by Liang et al. (2006) by adding an extra term after the corrector step to save on significant computational expenses at the cost of a marginal reduction of accuracy. In this study, the Massflow-2D code based on the MacCormackTVD scheme with variable computational domain is proposed as a solution for the shallow water equations associated with the mountainous hazard dynamic procedure in natural terrains. The scheme has the following attributes: 1. Being able to simulate discontinuous flows such as those associated with shock propagation.
2
C. Ouyang et al. / Computers & Geosciences 52 (2013) 1–10
2. Suitability to natural terrain by taking into account the additional acceleration and deceleration of the topography. 3. Time efficiency along with second-order accuracy in both time and space.
in which, 8 9 h > < > = X ¼ qx ; > :q > ; y
2. Model descriptions 2.1. Governing equations The characteristics of the mass movement can be described by the Reynolds-averaged Navier–Stokes equations, which are written in the differential conservation form as @r þ rUru~ ¼ 0 @t
ð1Þ
@ru~ þ rUru~ u~ ¼ rUs~ þ rg @t
ð2Þ
Here, r ¼mass density; g¼gravitational acceleration; u¼mass velocity vector with components in the three axes of the global Cartesian coordinates x, y, and z, with z aligned with the gravity and vertical. The symbol ( ) represents the parameters varying with vertical direction. Eqs. (1) and (2) are complicated and difficult to solve in a three dimensional space. To solve these equations, they can be reduced to a two dimensional shallow water form by integrating over the vertical depth. This simplification is valid if the movement of the mass in the horizontal direction far outweighs the movement in the vertical direction. For more details of this integration, the reader is directed to the papers (Denlinger and Iverson (2004), Steffler and Yee-Chung (1993)). Consequently, the 2D depthaveraged mass and momentum balance equations, written in the vector format, reduce to @X @F @G þ þ ¼ SþT @t @x @y
ð3Þ
9 > > = kx g 0 h2 x þ ; F¼ h 2 > > > > ; : qx qy 8 qx > > < q2
G¼
h
8 9 0 > > < = @z =@x f 0 S ¼ g h 1 þ ð@zf =@xÞ2 þ ð@zf =@yÞ2 þSf x ; > > : ; 0 9 8 0 > > > > = < 0 T¼ @zf =@y > > 0 > ; : g h 1 þ ð@z =@xÞ2 þ ð@z =@yÞ2 þ Sf y > f
8 q > > > y
> > =
> > 2 0 2 > > > ; : qy þ ky g2 h > h
f
In the above equations, h, qx ¼hu and qy ¼ hv are the total mass height from the base zf to the free surface zs, the discharges per unit width in the x-direction, and the discharges per unit width in the ydirection, respectively; u and v are the depth-averaged velocity components in the x- and y-axial directions, respectively; k¼(kx,ky) is the lateral earth pressure coefficient; Sf ¼(Sfx,Sfy) represents the frictional resistance. It is worth emphasizing that we adopt the enhanced gravity g0 in place of the earth gravity, g, to take into account the acceleration variation due to the vertical velocity component. It is necessary to use non-hydrostatic accelerations and centripetal forces when flow occurs over a rugged terrain. The evolution of the vertical velocity, w, can be determined from the flow velocity at the bed, wf, and the surface, ws, as well as from their relative distribution. Furthermore, a linear variation of w through the depth of the flow is assumed as w¼
1 ðw þ ws Þ 2 f
ð4Þ
where wf and ws can be obtained from the kinematic conditions of no mass flux through the bed and no flow through the surface as follows: wf ¼ ½wz ¼ zf ¼ u
@zf @zf þv @x @y
Fig. 1. Sketch of the finite difference scheme with the variable computational domain.
ð5Þ
C. Ouyang et al. / Computers & Geosciences 52 (2013) 1–10
ws ¼ ½wz ¼ zs ¼
@zs @zs @zs þu þv @t @x @y
ð6Þ
The total vertical acceleration, g0 , can be deduced from the earth gravity, g, by accounting for the effects of dynamic acceleration as follows: g0 ¼ g þ
@w @w @w dw þu þv ¼ gþ @t @x @y dt
ð7Þ
where dw/dt is a material derivative that accounts for the effects of additional acceleration or deceleration because of the topography. A similar approach has also been used by Denlinger and Connell (2008), Denlinger and Iverson (2004) and Walters (2005).
3
The lateral earth pressure coefficients, when the Coulomb friction model is adopted, are given by 8 ð@u=@x 4 eÞ k > < active ð@u=@x r eÞ ; kx ¼ 1 ð8Þ > :k ð@u=@x o eÞ passive
8 k > < active ky ¼ 1 > :k
passive
ð@v=@y 4 eÞ ð@v=@y r eÞ
ð9Þ
ð@v=@y o eÞ
where e is small and used as a threshold to avoid the machine precision error. The lateral earth pressure coefficients in the active
Fig. 2. Comparisons of the numerical and analytical solutions of (a) the flow height and (b) the velocity, for h1 ¼ 1.0 m and h2 ¼ 0.1 m at t ¼0.1 s.
4
C. Ouyang et al. / Computers & Geosciences 52 (2013) 1–10
elongation (kactive) and passive compression (kpassive) depend on the strain rate of the moving material columns, as suggested by Savage and Hutter (1989); kactive kpassive
)
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 18 1ð1 þtan2 dÞcos2 f 1 ¼ cos2 f
friction model, (Sfx,Sfy) are expressed as
S¼ ð10Þ
8 mgh u > < Sf x ¼ ð1 þ ð@zf =@xÞ2 þ ð@zf =@yÞ2 Þ1=2 U U > : Sfy ¼ ð1 þ ð@z
mgh 2 f =@xÞ
þ ð@zf =@yÞ2 Þ1=2
U Uv
,
ð11Þ
while for the Voellmy friction model, (Sfx,Sfy) are where f and d are the internal and bed friction angles. For water flow or dilute debris flows, the lateral pressure coefficients are taken as kx ¼ky ¼1.0. Various friction models, such as the Coulomb and Voellmy friction model, have been implemented in the code. For the Coulomb
8 2 mgh > > > Sf x ¼ ð1 þ ð@z =@xÞ2 þ ð@z =@yÞ2 Þ1=2 þg Uhx Uu < f f S¼ 2 > mgh > > Sfy ¼ ð1 þ ð@z =@xÞ2 þ ð@z =@yÞ2 Þ1=2 þ g Uhx Uv : f
f
Fig. 3. Comparisons of the numerical and analytical solutions of (a) the flow height and (b) the velocity, for h1 ¼ 1.0 m and h2 ¼10 3 m at t¼ 0.1 s.
ð12Þ
C. Ouyang et al. / Computers & Geosciences 52 (2013) 1–10
where m and x are the effective friction and turbulent coefficients, respectively, and U is the moving velocity of the mass. The choice of the lateral pressure coefficient equation and the friction model depend on the physical property of the mass movement under investigation. The earth pressure coefficients are set to be one when the Voellmy friction model is used.
2.2. Numerical implementation In this study, the MacCormack-TVD finite difference scheme is adopted to solve Eq. (3) as follows. First, based on the operator-splitting technique, Eq. (3) is divided into two separate
5
one-dimensional problems namely, @X @F þ ¼ S; @t @x
ð13aÞ
@X @G þ ¼T @t @y
ð13bÞ
Next, the solution at the (n þ1)Dt is obtained by the following step: Dt Dt Dt Dt n Ly 2 Ly 1 Lx 1 Xi,j ð14Þ Xni,jþ 1 ¼ Lx 2 2 2 2 2 where Lx and Ly represent a predictor-corrector-averaged computational procedure, each executed twice to obtain the solution at
Fig. 4. Comparisons of the numerical and analytical solutions of (a) the flow height and (b) the velocity, for h1 ¼1.0 m and h2 ¼10 6 m at t ¼0.1 s.
6
C. Ouyang et al. / Computers & Geosciences 52 (2013) 1–10
Fig. 5. (a) Three-dimensional and (b) two-dimensional water level contour diagrams of the two-dimensional partial dam break at t ¼7.2 s.
C. Ouyang et al. / Computers & Geosciences 52 (2013) 1–10
3. Numerical tests
the next step. For Lx1, the solution procedure is given by Xpi,j ¼ Xni,j ðFni,j Fni1,j ÞU
Dt þ Sni,j UDt=2 2Dx
ð15Þ
Dt þ Spi,j UDt=2 2Dx
ð16Þ
Xci,j ¼ Xni,j ðFpiþ 1,j Fpi,j ÞU n þ 1=2
Xi,j
n þ þ ¼ ðXpi,j þXci,j Þ=2 þ ½Gðri,j Þ þGðr i þ 1,j ÞDXi þ 1=2,j ½Gðri1,j Þ n þ Gðr i,j ÞDXi1=2,j
ð17Þ
where the superscripts p and c indicate the predictor and corrector steps, respectively, and
DXniþ 1=2,j ¼ ðXniþ 1,j Xni,j Þ;
þ ri,j ¼
ðDXni1=2,j , DXniþ 1=2,j Þ
DXni1=2,j ¼ ðXni,j Xni1,j Þ
; ðDXniþ 1=2,j , DXniþ 1=2,j Þ
r i,j ¼
ðDXni1=2,j , DXniþ 1=2,j Þ ðDXni1=2,j , DXni1=2,j Þ
ð18Þ
ð19Þ
The function G() is expressed as GðxÞ ¼ 0:5 C ½1jðxÞ
ð20Þ
in which function j() is a minmod flux limiter function and can be written as
jðxÞ ¼ maxð0,minð2x,1ÞÞ
ð21Þ
In Eqs. (20) and (21), x is a function variable. The variable C is expressed as Cr ð1CrÞ Cr r0:5 C¼ ð22Þ 0:25 Cr 40:5 where the Cr is the local Courant number defined as pffiffiffiffiffiffi ð9qx =h9 þ ghÞDt Cr ¼ 2Dx
7
ð23Þ
3.1. One-dimensional dam break simulation In the first numerical test, the 1D dam break problem is solved by using the Massflow-2D code. Such dam break problem has been adopted to verify shock-capturing capacity and used extensively by many authors after an analytical solution was proposed by Stoker (1957). In the current framework, the Massflow-2D code can be reduced to solve the one-dimensional problem with some minor modifications. The initial upstream water depth h1 is set to 1.0 m and the downstream water depth h2 is set to 10 1, 10 3, 10 6 m. The present configuration is the same as that used by Louaked and Hanich (1998). In this case, the constant earth gravity g¼9.8 m2/s is adopted. Figs. 2–4 show the comparison of the computational and analytical solution for h2 ¼10 1,10 3,10 6 for time t¼0.1 s. It is quite evident that the computational solution shows a good agreement with the analytical solution, even for the case of h1/h2 ¼106, as shown in Fig. 3. Therefore, the current numerical scheme can be considered sufficiently robust for the case of very low initial flow depth. The results show that the solution is correct when a low initial flow depth is used—a condition that is usually adopted by most authors.
3.2. Two-dimensional dam break simulation In the two-dimensional scenario, the Massflow-2D code is applied to simulate the 2D frictionless partial dam break problem studied by Fennema and Chaudhry (1990). This problem has also been investigated by several other researchers as well
Table 1 Computational parameters for the Coulomb model. kactive
kpassive
tan1 ðmÞ ðdeg:Þ
g (kN/m3)
Dt (s)
0.8
2.5
20
22
0.005
2.3. Variable computational domain For the dam break, landslide and debris flow hazards, the area of the mobile mass is generally much smaller than the potentially affected area. However, the initial computational domain in the finite difference method must be allocated a large enough space to cover the total potentially affected area. This causes the matrices involved in computing Eqs. (15)–(17) to be very large and the corresponding calculations become quite time-consuming. The variable computational domain method, as illustrated in Fig. 1, is proposed to circumvent this shortcoming. Assuming the total number of nodes in the whole domain to be m and n, in the x- and y-axial directions respectively, a matrix X (m,n,3) can be conceived to represent the discretization of the X matrix in Eq. (1). F_m (m,n,3) is a logical matrix, which is then assigned to be true in the domain of the mobile mass and false for all other cases. Matrix F_ada (m,n,3) is another logical matrix, which is obtained by the addition of two extra nodes in all the directions, to the already existing nodes in the domain of the mobile mass. In the front of the predictor-corrector-averaged step expressed in Eqs. (15)–(17), a Matlab command X_¼X (F_ada) is executed to convert the three-dimensional matrix X into a one-dimensional (1D) array. The 1D array X can then be used in the predictor-correctoraveraged step. At the end of the predictor-corrector-averaged step, the command X (F_ada)¼X_ updates the corresponding data in the variable computational domain. In this way, the computation will be adaptively focused on the mobile mass area.
Fig. 6. Snapshots of the deposition height contour of the Hong Kong landslide at t¼ 0, 6.5 s, 12 s, and 100 s. The contour interval is 0.2 m.
8
C. Ouyang et al. / Computers & Geosciences 52 (2013) 1–10
(Alcrudo and Garcia-Navarro, 1993; Fagherazzi et al., 2004; Lin et al., 2003; Xia et al., 2010). The computational domain for this problem consists of a 200 m 200 m region and is subdivided
into 80 80 rectangular meshes. As the initial conditions, the upstream water depth is set to 10 m and the downstream water depth is set to 5 m. The breach, located at x ¼100 m and
Fig. 7. Snapshots of the deposition height contour of the Nora debris flow at t¼ 100 s, 200, 400, 600, 800, and 1000 s (a–f). (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)
C. Ouyang et al. / Computers & Geosciences 52 (2013) 1–10
95 mryr175 m, is assumed to fail instantaneously. Fig. 5(a) and (b) show the 3D view of the water surface elevation, contour map of water depths and velocity vector arrows at t ¼7.2 s. The whirlpools at both sides of the dam site are evident as well. Additionally, the current computational results show a good agreement with previous results for the same problem (Alcrudo and Garcia-Navarro, 1993; Fagherazzi et al., 2004; Fennema and Chaudhry, 1990; Lin et al., 2003; Xia et al., 2010). 3.3. Hong Kong landslide simulation The rainfall-induced landslide, which occurred in Hong Kong in 2000 and was investigated by Chen and Lee (2003), constitutes the third example in this study. The DEM data of the ground surface and the activated initial source zone are created from the terrain elevation map provided by Chen. The Coulomb friction model and the lateral earth pressure expressed in Eq. (10) are adopted for this purpose and shown in Table 1. The grid size is set to 0.5 m and the time increment is taken as 0.005 s. The successive depth contours of the mobilized debris at t ¼0, 6.5, 12 and 100 s are shown in Fig. 6. It is quite evident that the debris trajectory is well constrained in the area where field measurement is available. The estimated final front of the debris is located at 126 m, which agrees very well with the field observation of 120 m and Chen’s results of 124 m. The computed deposition depths at the front and in the body are 0.7 m and 0.95 m, respectively, which match the corresponding field observations of 0.6 and 1.2 m. Therefore, as the results tally well with the field observations and Chen’s work, it is concluded that the current computational code is able to capture the dynamic nature of real landslide hazards. 3.4. Nora debris flow simulation The dynamic procedure of the Nora debris flow in the Italian Alps in 2000, which has been carefully researched by Pirulli and Marco (2010), is computed by the Massflow-2D program to test its robustness. A high-resolution pre-event DEM with 5 m grid spacing is generated from the topographic map in Pirulli’s paper. An initial source with a volume of approximately 10,000 m3 is placed in the upstream channel. The Voellmy parameters are taken as m ¼ 0.1 and x ¼200 m/s2, which have been demonstrated to give the best-fit numerical results (Pirulli and Marco, 2010). The time step is set to 0.005 s. The sequential snapshots of the debris flow height as well as the field deposition area are shown with a red line in Fig. 7. Evidently, the debris flow initially moves down along the valley. It then bifurcates into two parts inside an alluvial fan. The final deposition zone, shown in Fig. 7(f), is consistent with the field observations and Marina’s numerical results. Moreover, the variable computational domain method reduces the computational cost significantly because the moving mass of debris occupies a very small part of the whole computational domain.
4. Conclusions and remarks In this contribution, a robust and high-resolution MacCormack-TVD finite difference scheme with variable computational domain is implemented into Massflow-2D and is verified by a series of numerical tests. The excellent agreement with previous results ensures the robustness and validity of the current calculations. To enhance and strengthen the Massflow-2D code, additional factors that include surface erosion, sediment load and porous pressure need to be incorporated and will be the focus of future studies.
9
Acknowledgments The authors would like to thank Richard M. Iverson and Roger P. Denlinger’s suggestions, which greatly improved this manuscript. We also thank two anonymous reviewers for their helpful comments. Financial support from the NSFC (Grant no. 41101008), the Youth Talent Team Program of the Institute of Mountain Hazards and Environment, CAS, the Opening Fund of the State Key Laboratory of Geohazard Prevention and Geoenvironment Protection (Chengdu University of Technology) (Grant no. SKLGP2011K010) are acknowledged.
References Alcrudo, F., Garcia-Navarro, P., 1993. A high-resolution Godunov-type scheme in finite volumes for the 2D shallow-water equations. International Journal for Numerical Methods in Fluids 16, 489–505. ¨ Beguerı´a, S., Van Asch, T.W.J., Malet, J.P., Grondahl, S., 2009. A GIS-based numerical model for simulating the kinematics of mud and debris flows over complex terrain. Natural Hazards and Earth System Sciences 9, 1897–1909. Berger, M.J., George, D.L., LeVeque, R.J., Mandli, K.T., 2011. The GeoClaw software for depth-averaged flows with adaptive refinement. Advances in Water Resources 34, 1195–1206. Chen, H., Lee, C.F., 2003. A dynamic model for rainfall-induced landslides on natural slopes. Geomorphology 51, 269–288. Christen, M., Kowalski, J., Bartelt, P., 2010. RAMMS: numerical simulation of dense snow avalanches in three-dimensional terrain. Cold Regions Science and Technology 63, 1–14. Cundall, P.A., 2001. A discontinuous future for numerical modelling in geomechanics? Geotechnical Engineering 149, 41–47. Denlinger, R.P., Connell, D.R.H.O., 2008. Computing nonhydrostatic shallow-water flow over steep terrain. Journal of Hydraulic Engineering 134, 1590–1602. Denlinger, R.P., Iverson, R.M., 2004. Granular avalanches across irregular threedimensional terrain: 1. Theory and computation. Journal of Geophysical Research 109, F01014. Fagherazzi, S., Rasetarinera, P., Hussaini, M.Y., Furbish, D.J., 2004. Numerical solution of the dam-break problem with a discontinuous Galerkin method. Journal of Hydraulic Engineering 130, 532–539. Fennema, R.J., Chaudhry, M.H., 1990. Explicit methods for 2-D transient freesurface flows. Journal of Hydraulic Engineering 116, 1013–1034. Leveque, R.J., 2002. Finite Volume Methods for Hyperbolic Problems. Cambridge University Press, Cambridge, UK. Liang, D., Falconer, R.A., Lin, B., 2006. Comparison between TVD-MacCormack and ADI-type solvers of the shallow water equations. Advances in Water Resources 29, 1833–1845. Lin, G.-F., Lai, J.-S., Guo, W.-D., 2003. Finite-volume component-wise TVD schemes for 2D shallow water equations. Advances in Water Resources 26, 861–873. Louaked, M., Hanich, L., 1998. TVD scheme for the shallow water equations. Journal of Hydraulic Research 36, 363–378. Pirulli, M., Marco, F., 2010. Description and numerical modelling of the October 2000 Nora debris flow, Northwestern Italian Alps. Canadian Geotechnical Journal 47, 135–146. Pitman, E.B., Nichita, C.C., Patra, A., Bauer, A., Sheridan, M., Bursik, M., 2003. Computing granular avalanches and landslides. Physics of Fluids 15, 3638–3646. Prochaska, A.B., Santi, P.M., Higgins, J.D., Cannon, S.H., 2008. Debris-flow runout predictions based on the average channel slope (ACS). Engineering Geology 98, 29–40. Quecedo, M., Pastor, M., Herreros, M.I., Ferna´ndez Merodo, J.A., 2004. Numerical modelling of the propagation of fast landslides using the finite element method. International Journal for Numerical Methods in Engineering 59, 755–794. Rickenmann, D., 1999. Empirical relationships for debris flows. Natural Hazards 19, 47–77. Savage, S.B., Hutter, K., 1989. The motion of a finite mass of granular material down a rough incline. Journal of Fluid Mechanics 199, 177–215. Schilling, S.P., Iverson, R.M., 1997. Automated, reproducible delineation of zones at risk from inundation by large volcanic debris flows, Proceedings of the 1997 1st International Conference on Debris-Flow Hazards Mitigation: Mechanics, Prediction, and Assessment, August 7, 1997–August 9, 1997. ASCE, San Francisco, CA, USA, pp. 176–186. Steffler, P.M., Yee-Chung, J., 1993. Depth averaged and moment equations for moderately shallow free surface flow. Journal of Hydraulic Research 31, 5–17. Stoker, J.J., 1957. Water Waves. InterScience, New York. Tseng, M.H., Chu, C.R., 2000. Two-dimensional shallow water flows simulation using TVD-MacCormack scheme. Journal of Hydraulic Research 38, 123–131. Vincent, S., Caltagirone, J.P., Bonneton, P., 2001. Numerical modelling of bore propagation and run-up on sloping beaches using a MacCormack TVD scheme. Journal of Hydraulic Research 39, 41–49.
10
C. Ouyang et al. / Computers & Geosciences 52 (2013) 1–10
Walters, R.A., 2005. A semi-implicit finite element model for non-hydrostatic (dispersive) surface waves. International Journal for Numerical Methods in Fluids 49, 721–737. Wang, F., Sassa, K., 2010. Landslide simulation by a geotechnical model combined with a model for apparent friction change. Physics and Chemistry of the Earth, Parts A/B/C 35, 149–161.
Wang, J.S., Ni, H.G., He, Y.S., 2000. Finite-difference TVD scheme for computation of dam-break problems. Journal of Hydraulic Engineering 126, 253–262. Xia, J., Lin, B., Falconer, R.A., Wang, G., 2010. Modelling dam-break flows over mobile beds using a 2D coupled approach. Advances in Water Resources 33, 171–183.
Contents lists available at SciVerse ScienceDirect
Computers & Geosciences journal homepage: www.elsevier.com/locate/cageo
A MacCormack-TVD finite difference method to simulate the mass flow in mountainous terrain with variable computational domain Chaojun Ouyang a,b,c, Siming He a,b,n, Qiang Xu c, Yu Luo a,b, Wencheng Zhang a,b a
Key Laboratory of Mountain Hazards and Surface Process, Chinese Academy of Science, Chengdu 610041, China Institute of Mountain Hazards and Environment (IMHE), Chinese Academy of Sciences, Chengdu 610041, China c National Laboratory of Geohazard Prevention and Geoenvironment Protection, Chengdu University of Technology, Chengdu 610059, China b
a r t i c l e i n f o
abstract
Article history: Received 28 March 2012 Received in revised form 25 August 2012 Accepted 27 August 2012 Available online 25 September 2012
A two-dimensional mountainous mass flow dynamic procedure solver (Massflow-2D) using the MacCormack-TVD finite difference scheme is proposed. The solver is implemented in Matlab on structured meshes with variable computational domain. To verify the model, a variety of numerical test scenarios, namely, the classical one-dimensional and two-dimensional dam break, the landslide in Hong Kong in 1993 and the Nora debris flow in the Italian Alps in 2000, are executed, and the model outputs are compared with published results. It is established that the model predictions agree well with both the analytical solution as well as the field observations. & 2012 Elsevier Ltd. All rights reserved.
Keywords: Finite difference method MacCormack scheme TVD scheme Dam break Landslide Debris flow
1. Introduction Global climate change has brought along with it a fair share of disastrous weather conditions that have also increased the potential for mountainous hazards such as landslides, debris flows, dam-breaks, and floods in natural terrains. Such events can cause substantial damage to human life and engineering facilities. Therefore, hazard prediction and risk evaluation of hazardous zones are of considerable importance and have rightfully attracted the attention of many scientists. The common computational methods used to evaluate the disaster risk range for a potential hazard can be classified into three categories: empirical formula methods, discrete element methods and continuum modeling methods. The major advantage of any empirical method is its simplicity, while the disadvantage is the high demand for field or experimental data and ignorance of the guiding mechanisms (Prochaska et al., 2008; Rickenmann, 1999; Schilling and Iverson, 1997). The discrete element method allows for the quantitative reproduction of a particle trajectory, but by the very nature of such calculations, the method is timeconsuming and requires micromechanical input parameters that are often difficult to determine (Cundall, 2001). The continuum modeling method considers the moving mass to be a fluid and can
n Corresponding author at: Institute of Mountain Hazards and Environment (IMHE), Chinese Academy of Sciences, Chengdu 610041, China. Fax: þ 86 28 85222258. E-mail address: [email protected] (S. He).
0098-3004/$ - see front matter & 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.cageo.2012.08.024
therefore be depicted by the incompressible Navier–Stokes equations. These equations can then be further reduced to the two-dimensional (2D) shallow water equations by a depthaveraged procedure. Many numerical methods, including the finite volume method (Berger et al., 2011; Christen et al., 2010; Leveque, 2002; Pitman et al., 2003), the finite difference method (Beguerı´a et al., 2009; Liang et al., 2006; Wang and Sassa, 2010) and the finite element method (Chen and Lee, 2003; Quecedo et al., 2004), have been proposed to resolve the sharp discontinuities associated with the hyperbolic problems, such as the 2D shallow water equations. Nevertheless, simulation of a large-range and long-duration mountainous hazard on a natural terrain with less time and higher resolution remains a challenge. The MacCormack finite difference scheme with total variation diminishing (TVD) has been investigated by many authors and has been shown to be a feasible and robust computational solution (Tseng and Chu, 2000; Vincent et al., 2001; Wang et al., 2000). This method was further modified by Liang et al. (2006) by adding an extra term after the corrector step to save on significant computational expenses at the cost of a marginal reduction of accuracy. In this study, the Massflow-2D code based on the MacCormackTVD scheme with variable computational domain is proposed as a solution for the shallow water equations associated with the mountainous hazard dynamic procedure in natural terrains. The scheme has the following attributes: 1. Being able to simulate discontinuous flows such as those associated with shock propagation.
2
C. Ouyang et al. / Computers & Geosciences 52 (2013) 1–10
2. Suitability to natural terrain by taking into account the additional acceleration and deceleration of the topography. 3. Time efficiency along with second-order accuracy in both time and space.
in which, 8 9 h > < > = X ¼ qx ; > :q > ; y
2. Model descriptions 2.1. Governing equations The characteristics of the mass movement can be described by the Reynolds-averaged Navier–Stokes equations, which are written in the differential conservation form as @r þ rUru~ ¼ 0 @t
ð1Þ
@ru~ þ rUru~ u~ ¼ rUs~ þ rg @t
ð2Þ
Here, r ¼mass density; g¼gravitational acceleration; u¼mass velocity vector with components in the three axes of the global Cartesian coordinates x, y, and z, with z aligned with the gravity and vertical. The symbol ( ) represents the parameters varying with vertical direction. Eqs. (1) and (2) are complicated and difficult to solve in a three dimensional space. To solve these equations, they can be reduced to a two dimensional shallow water form by integrating over the vertical depth. This simplification is valid if the movement of the mass in the horizontal direction far outweighs the movement in the vertical direction. For more details of this integration, the reader is directed to the papers (Denlinger and Iverson (2004), Steffler and Yee-Chung (1993)). Consequently, the 2D depthaveraged mass and momentum balance equations, written in the vector format, reduce to @X @F @G þ þ ¼ SþT @t @x @y
ð3Þ
9 > > = kx g 0 h2 x þ ; F¼ h 2 > > > > ; : qx qy 8 qx > > < q2
G¼
h
8 9 0 > > < = @z =@x f 0 S ¼ g h 1 þ ð@zf =@xÞ2 þ ð@zf =@yÞ2 þSf x ; > > : ; 0 9 8 0 > > > > = < 0 T¼ @zf =@y > > 0 > ; : g h 1 þ ð@z =@xÞ2 þ ð@z =@yÞ2 þ Sf y > f
8 q > > > y
> > =
> > 2 0 2 > > > ; : qy þ ky g2 h > h
f
In the above equations, h, qx ¼hu and qy ¼ hv are the total mass height from the base zf to the free surface zs, the discharges per unit width in the x-direction, and the discharges per unit width in the ydirection, respectively; u and v are the depth-averaged velocity components in the x- and y-axial directions, respectively; k¼(kx,ky) is the lateral earth pressure coefficient; Sf ¼(Sfx,Sfy) represents the frictional resistance. It is worth emphasizing that we adopt the enhanced gravity g0 in place of the earth gravity, g, to take into account the acceleration variation due to the vertical velocity component. It is necessary to use non-hydrostatic accelerations and centripetal forces when flow occurs over a rugged terrain. The evolution of the vertical velocity, w, can be determined from the flow velocity at the bed, wf, and the surface, ws, as well as from their relative distribution. Furthermore, a linear variation of w through the depth of the flow is assumed as w¼
1 ðw þ ws Þ 2 f
ð4Þ
where wf and ws can be obtained from the kinematic conditions of no mass flux through the bed and no flow through the surface as follows: wf ¼ ½wz ¼ zf ¼ u
@zf @zf þv @x @y
Fig. 1. Sketch of the finite difference scheme with the variable computational domain.
ð5Þ
C. Ouyang et al. / Computers & Geosciences 52 (2013) 1–10
ws ¼ ½wz ¼ zs ¼
@zs @zs @zs þu þv @t @x @y
ð6Þ
The total vertical acceleration, g0 , can be deduced from the earth gravity, g, by accounting for the effects of dynamic acceleration as follows: g0 ¼ g þ
@w @w @w dw þu þv ¼ gþ @t @x @y dt
ð7Þ
where dw/dt is a material derivative that accounts for the effects of additional acceleration or deceleration because of the topography. A similar approach has also been used by Denlinger and Connell (2008), Denlinger and Iverson (2004) and Walters (2005).
3
The lateral earth pressure coefficients, when the Coulomb friction model is adopted, are given by 8 ð@u=@x 4 eÞ k > < active ð@u=@x r eÞ ; kx ¼ 1 ð8Þ > :k ð@u=@x o eÞ passive
8 k > < active ky ¼ 1 > :k
passive
ð@v=@y 4 eÞ ð@v=@y r eÞ
ð9Þ
ð@v=@y o eÞ
where e is small and used as a threshold to avoid the machine precision error. The lateral earth pressure coefficients in the active
Fig. 2. Comparisons of the numerical and analytical solutions of (a) the flow height and (b) the velocity, for h1 ¼ 1.0 m and h2 ¼ 0.1 m at t ¼0.1 s.
4
C. Ouyang et al. / Computers & Geosciences 52 (2013) 1–10
elongation (kactive) and passive compression (kpassive) depend on the strain rate of the moving material columns, as suggested by Savage and Hutter (1989); kactive kpassive
)
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 18 1ð1 þtan2 dÞcos2 f 1 ¼ cos2 f
friction model, (Sfx,Sfy) are expressed as
S¼ ð10Þ
8 mgh u > < Sf x ¼ ð1 þ ð@zf =@xÞ2 þ ð@zf =@yÞ2 Þ1=2 U U > : Sfy ¼ ð1 þ ð@z
mgh 2 f =@xÞ
þ ð@zf =@yÞ2 Þ1=2
U Uv
,
ð11Þ
while for the Voellmy friction model, (Sfx,Sfy) are where f and d are the internal and bed friction angles. For water flow or dilute debris flows, the lateral pressure coefficients are taken as kx ¼ky ¼1.0. Various friction models, such as the Coulomb and Voellmy friction model, have been implemented in the code. For the Coulomb
8 2 mgh > > > Sf x ¼ ð1 þ ð@z =@xÞ2 þ ð@z =@yÞ2 Þ1=2 þg Uhx Uu < f f S¼ 2 > mgh > > Sfy ¼ ð1 þ ð@z =@xÞ2 þ ð@z =@yÞ2 Þ1=2 þ g Uhx Uv : f
f
Fig. 3. Comparisons of the numerical and analytical solutions of (a) the flow height and (b) the velocity, for h1 ¼ 1.0 m and h2 ¼10 3 m at t¼ 0.1 s.
ð12Þ
C. Ouyang et al. / Computers & Geosciences 52 (2013) 1–10
where m and x are the effective friction and turbulent coefficients, respectively, and U is the moving velocity of the mass. The choice of the lateral pressure coefficient equation and the friction model depend on the physical property of the mass movement under investigation. The earth pressure coefficients are set to be one when the Voellmy friction model is used.
2.2. Numerical implementation In this study, the MacCormack-TVD finite difference scheme is adopted to solve Eq. (3) as follows. First, based on the operator-splitting technique, Eq. (3) is divided into two separate
5
one-dimensional problems namely, @X @F þ ¼ S; @t @x
ð13aÞ
@X @G þ ¼T @t @y
ð13bÞ
Next, the solution at the (n þ1)Dt is obtained by the following step: Dt Dt Dt Dt n Ly 2 Ly 1 Lx 1 Xi,j ð14Þ Xni,jþ 1 ¼ Lx 2 2 2 2 2 where Lx and Ly represent a predictor-corrector-averaged computational procedure, each executed twice to obtain the solution at
Fig. 4. Comparisons of the numerical and analytical solutions of (a) the flow height and (b) the velocity, for h1 ¼1.0 m and h2 ¼10 6 m at t ¼0.1 s.
6
C. Ouyang et al. / Computers & Geosciences 52 (2013) 1–10
Fig. 5. (a) Three-dimensional and (b) two-dimensional water level contour diagrams of the two-dimensional partial dam break at t ¼7.2 s.
C. Ouyang et al. / Computers & Geosciences 52 (2013) 1–10
3. Numerical tests
the next step. For Lx1, the solution procedure is given by Xpi,j ¼ Xni,j ðFni,j Fni1,j ÞU
Dt þ Sni,j UDt=2 2Dx
ð15Þ
Dt þ Spi,j UDt=2 2Dx
ð16Þ
Xci,j ¼ Xni,j ðFpiþ 1,j Fpi,j ÞU n þ 1=2
Xi,j
n þ þ ¼ ðXpi,j þXci,j Þ=2 þ ½Gðri,j Þ þGðr i þ 1,j ÞDXi þ 1=2,j ½Gðri1,j Þ n þ Gðr i,j ÞDXi1=2,j
ð17Þ
where the superscripts p and c indicate the predictor and corrector steps, respectively, and
DXniþ 1=2,j ¼ ðXniþ 1,j Xni,j Þ;
þ ri,j ¼
ðDXni1=2,j , DXniþ 1=2,j Þ
DXni1=2,j ¼ ðXni,j Xni1,j Þ
; ðDXniþ 1=2,j , DXniþ 1=2,j Þ
r i,j ¼
ðDXni1=2,j , DXniþ 1=2,j Þ ðDXni1=2,j , DXni1=2,j Þ
ð18Þ
ð19Þ
The function G() is expressed as GðxÞ ¼ 0:5 C ½1jðxÞ
ð20Þ
in which function j() is a minmod flux limiter function and can be written as
jðxÞ ¼ maxð0,minð2x,1ÞÞ
ð21Þ
In Eqs. (20) and (21), x is a function variable. The variable C is expressed as Cr ð1CrÞ Cr r0:5 C¼ ð22Þ 0:25 Cr 40:5 where the Cr is the local Courant number defined as pffiffiffiffiffiffi ð9qx =h9 þ ghÞDt Cr ¼ 2Dx
7
ð23Þ
3.1. One-dimensional dam break simulation In the first numerical test, the 1D dam break problem is solved by using the Massflow-2D code. Such dam break problem has been adopted to verify shock-capturing capacity and used extensively by many authors after an analytical solution was proposed by Stoker (1957). In the current framework, the Massflow-2D code can be reduced to solve the one-dimensional problem with some minor modifications. The initial upstream water depth h1 is set to 1.0 m and the downstream water depth h2 is set to 10 1, 10 3, 10 6 m. The present configuration is the same as that used by Louaked and Hanich (1998). In this case, the constant earth gravity g¼9.8 m2/s is adopted. Figs. 2–4 show the comparison of the computational and analytical solution for h2 ¼10 1,10 3,10 6 for time t¼0.1 s. It is quite evident that the computational solution shows a good agreement with the analytical solution, even for the case of h1/h2 ¼106, as shown in Fig. 3. Therefore, the current numerical scheme can be considered sufficiently robust for the case of very low initial flow depth. The results show that the solution is correct when a low initial flow depth is used—a condition that is usually adopted by most authors.
3.2. Two-dimensional dam break simulation In the two-dimensional scenario, the Massflow-2D code is applied to simulate the 2D frictionless partial dam break problem studied by Fennema and Chaudhry (1990). This problem has also been investigated by several other researchers as well
Table 1 Computational parameters for the Coulomb model. kactive
kpassive
tan1 ðmÞ ðdeg:Þ
g (kN/m3)
Dt (s)
0.8
2.5
20
22
0.005
2.3. Variable computational domain For the dam break, landslide and debris flow hazards, the area of the mobile mass is generally much smaller than the potentially affected area. However, the initial computational domain in the finite difference method must be allocated a large enough space to cover the total potentially affected area. This causes the matrices involved in computing Eqs. (15)–(17) to be very large and the corresponding calculations become quite time-consuming. The variable computational domain method, as illustrated in Fig. 1, is proposed to circumvent this shortcoming. Assuming the total number of nodes in the whole domain to be m and n, in the x- and y-axial directions respectively, a matrix X (m,n,3) can be conceived to represent the discretization of the X matrix in Eq. (1). F_m (m,n,3) is a logical matrix, which is then assigned to be true in the domain of the mobile mass and false for all other cases. Matrix F_ada (m,n,3) is another logical matrix, which is obtained by the addition of two extra nodes in all the directions, to the already existing nodes in the domain of the mobile mass. In the front of the predictor-corrector-averaged step expressed in Eqs. (15)–(17), a Matlab command X_¼X (F_ada) is executed to convert the three-dimensional matrix X into a one-dimensional (1D) array. The 1D array X can then be used in the predictor-correctoraveraged step. At the end of the predictor-corrector-averaged step, the command X (F_ada)¼X_ updates the corresponding data in the variable computational domain. In this way, the computation will be adaptively focused on the mobile mass area.
Fig. 6. Snapshots of the deposition height contour of the Hong Kong landslide at t¼ 0, 6.5 s, 12 s, and 100 s. The contour interval is 0.2 m.
8
C. Ouyang et al. / Computers & Geosciences 52 (2013) 1–10
(Alcrudo and Garcia-Navarro, 1993; Fagherazzi et al., 2004; Lin et al., 2003; Xia et al., 2010). The computational domain for this problem consists of a 200 m 200 m region and is subdivided
into 80 80 rectangular meshes. As the initial conditions, the upstream water depth is set to 10 m and the downstream water depth is set to 5 m. The breach, located at x ¼100 m and
Fig. 7. Snapshots of the deposition height contour of the Nora debris flow at t¼ 100 s, 200, 400, 600, 800, and 1000 s (a–f). (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)
C. Ouyang et al. / Computers & Geosciences 52 (2013) 1–10
95 mryr175 m, is assumed to fail instantaneously. Fig. 5(a) and (b) show the 3D view of the water surface elevation, contour map of water depths and velocity vector arrows at t ¼7.2 s. The whirlpools at both sides of the dam site are evident as well. Additionally, the current computational results show a good agreement with previous results for the same problem (Alcrudo and Garcia-Navarro, 1993; Fagherazzi et al., 2004; Fennema and Chaudhry, 1990; Lin et al., 2003; Xia et al., 2010). 3.3. Hong Kong landslide simulation The rainfall-induced landslide, which occurred in Hong Kong in 2000 and was investigated by Chen and Lee (2003), constitutes the third example in this study. The DEM data of the ground surface and the activated initial source zone are created from the terrain elevation map provided by Chen. The Coulomb friction model and the lateral earth pressure expressed in Eq. (10) are adopted for this purpose and shown in Table 1. The grid size is set to 0.5 m and the time increment is taken as 0.005 s. The successive depth contours of the mobilized debris at t ¼0, 6.5, 12 and 100 s are shown in Fig. 6. It is quite evident that the debris trajectory is well constrained in the area where field measurement is available. The estimated final front of the debris is located at 126 m, which agrees very well with the field observation of 120 m and Chen’s results of 124 m. The computed deposition depths at the front and in the body are 0.7 m and 0.95 m, respectively, which match the corresponding field observations of 0.6 and 1.2 m. Therefore, as the results tally well with the field observations and Chen’s work, it is concluded that the current computational code is able to capture the dynamic nature of real landslide hazards. 3.4. Nora debris flow simulation The dynamic procedure of the Nora debris flow in the Italian Alps in 2000, which has been carefully researched by Pirulli and Marco (2010), is computed by the Massflow-2D program to test its robustness. A high-resolution pre-event DEM with 5 m grid spacing is generated from the topographic map in Pirulli’s paper. An initial source with a volume of approximately 10,000 m3 is placed in the upstream channel. The Voellmy parameters are taken as m ¼ 0.1 and x ¼200 m/s2, which have been demonstrated to give the best-fit numerical results (Pirulli and Marco, 2010). The time step is set to 0.005 s. The sequential snapshots of the debris flow height as well as the field deposition area are shown with a red line in Fig. 7. Evidently, the debris flow initially moves down along the valley. It then bifurcates into two parts inside an alluvial fan. The final deposition zone, shown in Fig. 7(f), is consistent with the field observations and Marina’s numerical results. Moreover, the variable computational domain method reduces the computational cost significantly because the moving mass of debris occupies a very small part of the whole computational domain.
4. Conclusions and remarks In this contribution, a robust and high-resolution MacCormack-TVD finite difference scheme with variable computational domain is implemented into Massflow-2D and is verified by a series of numerical tests. The excellent agreement with previous results ensures the robustness and validity of the current calculations. To enhance and strengthen the Massflow-2D code, additional factors that include surface erosion, sediment load and porous pressure need to be incorporated and will be the focus of future studies.
9
Acknowledgments The authors would like to thank Richard M. Iverson and Roger P. Denlinger’s suggestions, which greatly improved this manuscript. We also thank two anonymous reviewers for their helpful comments. Financial support from the NSFC (Grant no. 41101008), the Youth Talent Team Program of the Institute of Mountain Hazards and Environment, CAS, the Opening Fund of the State Key Laboratory of Geohazard Prevention and Geoenvironment Protection (Chengdu University of Technology) (Grant no. SKLGP2011K010) are acknowledged.
References Alcrudo, F., Garcia-Navarro, P., 1993. A high-resolution Godunov-type scheme in finite volumes for the 2D shallow-water equations. International Journal for Numerical Methods in Fluids 16, 489–505. ¨ Beguerı´a, S., Van Asch, T.W.J., Malet, J.P., Grondahl, S., 2009. A GIS-based numerical model for simulating the kinematics of mud and debris flows over complex terrain. Natural Hazards and Earth System Sciences 9, 1897–1909. Berger, M.J., George, D.L., LeVeque, R.J., Mandli, K.T., 2011. The GeoClaw software for depth-averaged flows with adaptive refinement. Advances in Water Resources 34, 1195–1206. Chen, H., Lee, C.F., 2003. A dynamic model for rainfall-induced landslides on natural slopes. Geomorphology 51, 269–288. Christen, M., Kowalski, J., Bartelt, P., 2010. RAMMS: numerical simulation of dense snow avalanches in three-dimensional terrain. Cold Regions Science and Technology 63, 1–14. Cundall, P.A., 2001. A discontinuous future for numerical modelling in geomechanics? Geotechnical Engineering 149, 41–47. Denlinger, R.P., Connell, D.R.H.O., 2008. Computing nonhydrostatic shallow-water flow over steep terrain. Journal of Hydraulic Engineering 134, 1590–1602. Denlinger, R.P., Iverson, R.M., 2004. Granular avalanches across irregular threedimensional terrain: 1. Theory and computation. Journal of Geophysical Research 109, F01014. Fagherazzi, S., Rasetarinera, P., Hussaini, M.Y., Furbish, D.J., 2004. Numerical solution of the dam-break problem with a discontinuous Galerkin method. Journal of Hydraulic Engineering 130, 532–539. Fennema, R.J., Chaudhry, M.H., 1990. Explicit methods for 2-D transient freesurface flows. Journal of Hydraulic Engineering 116, 1013–1034. Leveque, R.J., 2002. Finite Volume Methods for Hyperbolic Problems. Cambridge University Press, Cambridge, UK. Liang, D., Falconer, R.A., Lin, B., 2006. Comparison between TVD-MacCormack and ADI-type solvers of the shallow water equations. Advances in Water Resources 29, 1833–1845. Lin, G.-F., Lai, J.-S., Guo, W.-D., 2003. Finite-volume component-wise TVD schemes for 2D shallow water equations. Advances in Water Resources 26, 861–873. Louaked, M., Hanich, L., 1998. TVD scheme for the shallow water equations. Journal of Hydraulic Research 36, 363–378. Pirulli, M., Marco, F., 2010. Description and numerical modelling of the October 2000 Nora debris flow, Northwestern Italian Alps. Canadian Geotechnical Journal 47, 135–146. Pitman, E.B., Nichita, C.C., Patra, A., Bauer, A., Sheridan, M., Bursik, M., 2003. Computing granular avalanches and landslides. Physics of Fluids 15, 3638–3646. Prochaska, A.B., Santi, P.M., Higgins, J.D., Cannon, S.H., 2008. Debris-flow runout predictions based on the average channel slope (ACS). Engineering Geology 98, 29–40. Quecedo, M., Pastor, M., Herreros, M.I., Ferna´ndez Merodo, J.A., 2004. Numerical modelling of the propagation of fast landslides using the finite element method. International Journal for Numerical Methods in Engineering 59, 755–794. Rickenmann, D., 1999. Empirical relationships for debris flows. Natural Hazards 19, 47–77. Savage, S.B., Hutter, K., 1989. The motion of a finite mass of granular material down a rough incline. Journal of Fluid Mechanics 199, 177–215. Schilling, S.P., Iverson, R.M., 1997. Automated, reproducible delineation of zones at risk from inundation by large volcanic debris flows, Proceedings of the 1997 1st International Conference on Debris-Flow Hazards Mitigation: Mechanics, Prediction, and Assessment, August 7, 1997–August 9, 1997. ASCE, San Francisco, CA, USA, pp. 176–186. Steffler, P.M., Yee-Chung, J., 1993. Depth averaged and moment equations for moderately shallow free surface flow. Journal of Hydraulic Research 31, 5–17. Stoker, J.J., 1957. Water Waves. InterScience, New York. Tseng, M.H., Chu, C.R., 2000. Two-dimensional shallow water flows simulation using TVD-MacCormack scheme. Journal of Hydraulic Research 38, 123–131. Vincent, S., Caltagirone, J.P., Bonneton, P., 2001. Numerical modelling of bore propagation and run-up on sloping beaches using a MacCormack TVD scheme. Journal of Hydraulic Research 39, 41–49.
10
C. Ouyang et al. / Computers & Geosciences 52 (2013) 1–10
Walters, R.A., 2005. A semi-implicit finite element model for non-hydrostatic (dispersive) surface waves. International Journal for Numerical Methods in Fluids 49, 721–737. Wang, F., Sassa, K., 2010. Landslide simulation by a geotechnical model combined with a model for apparent friction change. Physics and Chemistry of the Earth, Parts A/B/C 35, 149–161.
Wang, J.S., Ni, H.G., He, Y.S., 2000. Finite-difference TVD scheme for computation of dam-break problems. Journal of Hydraulic Engineering 126, 253–262. Xia, J., Lin, B., Falconer, R.A., Wang, G., 2010. Modelling dam-break flows over mobile beds using a 2D coupled approach. Advances in Water Resources 33, 171–183.
