Debris flow simulation 2D(DFS 2D):Numerical modelling of debris flows and calibration of friction parameters

Minu Trs Abrhm,Nlim Stym,Biswjt Prdhn,Hongling Tin

a Department of Civil Engineering,Indian Institute of Technology Indore,Indore,India

b Centre for Advanced Modelling and Geospatial Information Systems (CAMGIS),School of Civil and Environmental Engineering,Faculty of Engineering and Information Technology,University of Technology Sydney,Sydney,Australia

c Center of Excellence for Climate Change Research,King Abdulaziz University,Jeddah,Saudi Arabia

d Earth Observation Centre,Institute of Climate Change,University Kebangsaan Malaysia,Bangi,Malaysia

e Key Laboratory of Mountain Hazards and Surface Process,Institute of Mountain Hazards and Environment,Chinese Academy of Sciences,Chengdu,China

Keywords:Debris flows Numerical model Rheology Debris flow simulation 2D (DFS 2D)

ABSTRACT Debris flows are rapid mass movements with a mixture of rock,soil and water.High-intensity rainfall events have triggered multiple debris flows around the globe,making it an important concern from the disaster management perspective.This study presents a numerical model called debris flow simulation 2D (DFS 2D) and applicability of the proposed model is investigated through the values of the model parameters used for the reproduction of an occurred debris flow at Yindongzi gully in China on 13 August 2010.The model can be used to simulate debris flows using three different rheologies and has a userfriendly interface for providing the inputs.Using DFS 2D,flow parameters can be estimated with respect to space and time.The values of the flow resistance parameters of model,dry-Coulomb and turbulent friction,were calibrated through the back analysis and the values obtained are 0.1 and 1000 m/s2,respectively.Two new methods of calibration are proposed in this study,considering the crosssectional area of flow and topographical changes induced by the debris flow.The proposed methods of calibration provide an effective solution to the cumulative errors induced by coarse-resolution digital elevation models (DEMs) in numerical modelling of debris flows.The statistical indices such as Willmott’s index of agreement,mean-absolute-error,and normalized-root-mean-square-error of the calibrated model are 0.5,1.02 and 1.44,respectively.The comparison between simulated and observed values of topographic changes indicates that DFS 2D provides satisfactory results and can be used for dynamic modelling of debris flows.

1.Introduction

Debris flows are a moving mass of loose soil particles ranging from clays to rocks with high water content,driven by gravity.Varnes(1978)classification uses the term debris to denote material that contains a high percentage of coarse particles.The updated classification by Hungr et al.(2014)defines debris flows as surging flow of saturated debris with very rapid to extremely rapid speed in a steep channel,which entrains material from the flow path.Such flows are usually triggered by extremely heavy rainfall or snowmelt.Debris flows are one among the highly disastrous hazard in the mountainous terrains due to their high-speed and long runout distances.These catastrophic landslides are induced due to rainfall and the critical rainfall conditions that may result in such events can be estimated (Teja et al.,2019; Dikshit et al.,2020; Bulzinetti et al.,2021; Segoni et al.,2021).The debris becomes fluidised due to high water content and flows downslope with very high momentum (Hungr et al.,2014).The flows exhibit viscous and turbulent properties due to high saturation and particle interactions.As the movement progresses,sediments from the bed are often entrained by the flowing mass,increasing the volume of flow significantly (Reid et al.,2016; Simoni et al.,2020).If there are no drainage paths,the flow cuts a V-shaped channel on the hillside as it moves downwards (Varnes,1978).The coarser particles become deposited and the more fluid part moves further (Varnes,1978).Such flows are chief agents of sediment transport in hilly areas(Beguería et al.,2009).When the flow occurs in human inhabited slopes,it is a threat to life and property and has an increased impact than other landslide types for their high velocity.

Calculating runout distance of debris flows has been attempted using different approaches(Hsu,1978;Hungr,1995;Denlinger and Iverson,2004;Christen et al.,2010).The initial attempts were based on an empirical relationship between the features of landslides and the runout distance (Hsu,1978; Corominas,1996; Rickenmann et al.,2006).Later on,analytical and numerical models were developed.The calculation of the runout distance can be approached by means of empirical methods (Aulitzky,1973) or by mathematical models.Models can be empirically-statistic based(Sassa,1988; Hungr,1995; Hurlimann et al.,2007),topographic gradient-based or numerical-based (Pastor et al.,2009,2014;Crosta et al.,2009,2017; Ceccatelli et al.,2017; Braun et al.,2018).Most of the existing numerical model uses a single rheological model and can be used for calculating either mudflows or debris flows.The dynamic models provide precise information about the runout process,including the flow shape,height,velocity and travel time,which can be used to evaluate the risk due to debris flows(Melo et al.,2020).

The debris-flow hazard can be evaluated by means of numerical routing models that provide maps of velocity and flow,deposition and erosion depths,as well as an estimate of the entrained and transported volume of sediments (Christen et al.,2010).This can help in assessing the impact of debris flow on obstacles in the path,possible damage and expected runout distance(Hussin et al.,2014).The simulations help in visualising the flow along with the magnitude of flow parameters (height and velocity) for different boundary conditions.

When describing the flow behaviour,rheology is the key parameter which decides the deformation and flow.Most debris flows follow non-Newtonian behaviour and the commonly adopted rheologies are: Coulomb model (Hungr and McDougall,2009),Bingham model (Coussot,1997; Malet et al.,2005) and Voellmy-Salm model (Voellmy,1955).Studies on grain inertial debris flows follow Bagnold-like rheology(Bagnold,1954),considering the flow of cohesionless particles in inertial and viscous regimes (Brufau et al.,2000; Deangeli,2008; Gregoretti et al.,2019),where the dispersive stresses induced by the contact between particles are also considered (Takahashi,2007).The simplest case is the assumption of laminar flow regime of viscoplastic materials with a constant values of yield strength and viscosity,as in Bingham model(Coussot,1994).This model is valid for flows with higher fine fraction.When the constant yield strength assumption of Bingham model is revised by considering the effect of cohesion and friction,Coulomb-viscous model can be defined (Johnson and Rodine,1984).Voellmy-Salm rheology and its modifications are being used to consider the effect of both Coulomb friction and turbulence in the flow behaviour and has been found satisfactory while performing a back analysis of past debris flows (Pirulli and Sorbino,2008; Hussin et al.,2012; Simoni et al.,2012; Frank et al.,2015;Frey et al.,2016; Bezak et al.,2020).

A number of numerical landslide models are currently available,using different numerical schemes and rheologies (Rickenmann,2005; Beguería et al.,2009; Christen et al.,2010).However,most of them are not user-friendly and provide no guidelines about the parameters to be used (McDougall,2017).McDougall (2017) has summarised the practices and runout modelling of landslides and concluded that the limited guidance in the selection of input parameter is a major challenge in runout modelling.The research also calls for the need for improved model efficiency and user friendliness of the existing models (McDougall,2017).

This study introduces a new tool,called debris flow simulation 2D(DFS 2D)which can be used to simulate the dynamics of debris flows by considering multiple rheologies,assuming single phase flow.The model also studies the effect of bed entrainment,which is a major factor in deciding the final flow volume and its impact on obstacles along the path (Cuomo et al.,2014,2016).The entrained sediments can make sudden changes in the flow behaviour (Kang and Chan,2018) and hence it is an important aspect to be considered in the modelling of debris flows(Iverson,1997;McDougall and Hungr,2005; Takahashi,2009; Cuomo et al.,2014; Iverson and Ouyang,2015).

DFS 2D is an attempt to develop an easy-to-use tool for debris flow modelling,considering multiple rheologies and bed entrainment.The model provides a brief information about each input parameter in the user interface.The model uses the soil type,moisture content and surface cover map,which are not found in the existing models,to provide the approximate range of dry-Coulomb friction,viscosity and critical shear stress values.The process of debris flow is highly complex and the parameters vary in all the three directions (McDougall,2017).Solving the variations in three dimensions is highly complicated and hence DFS 2D assumes that the velocity of the flow is constant along the depth of flow,following the shallow water conditions (Beguería et al.,2009).Shallow water equations are used when the thickness of flow is much lesser than the length and width of the flow.For smaller thickness values,the flow properties in the direction perpendicular to the plane of flow can be considered as uniform,i.e.the gradient is a constant.Thus,the complex three-dimensional (3D) problem is simplified to a two-dimensional (2D) problem which is experimentally justified in the literature(Savage and Hutter,1989;Greve and Hutter,1993).The mass and momentum conservation equations are used to dynamically model the flow and calculate the flow parameters with respect to space and time.The availability of different rheological models in the same tool provides the option to directly compare outputs and internal variables in all the models(Patra et al.,2020).Some numerical models have already provided this option to use multiple rheologies in the same model,but are too complex due to the requirement of complex input parameters.The model is an attempt to simulate the complex phenomenon of debris flows using single phase shallow water equations.A complex problem with a simple form is considered.The main objective is to provide a user-friendly tool for modelling debris flows,which can be easily used by the authorities to simulate the debris flow runout process.The input parameters are kept to a minimum and the interface provides a range of value for each parameter,which can be correlated to the type of soil at the site,making DFS 2D an easy-touse tool.

2.Numerical modelling using shallow water equations

2.1.Governing equations

DFS 2D assumes the flow material to be a single phase and homogeneous,based on Savage-Hutter theory(Savage and Hutter,1989).Ideal behaviour cannot be expected in the case of a complex phenomenon like debris flows and,hence,the choice of rheology is critical in the definition of governing equations.(Rickenmann,2005; Bertolo and Bottino,2008).Following the classical approach of debris flow modelling,the flow is modelled as a 2D continuum,using depth integrated shallow water equations(Savage and Hutter,1989;Iverson and Denlinger,2001;McDougall and Hungr,2004; Mangeney-Castelnau et al.,2005).The variation of kinematic characteristics across the space is calculated by the model,while it is assumed to be uniform along the depth of flow.DFS 2D uses a 2D Cartesian coordinates with global x and y in a Euclidian space to define the governing equations(Mangeney et al.,2007;Beguería et al.,2009).

The governing equations are the mass and momentum conservation equations,given by

where Uxand Uyare the velocity components in x and y directions,respectively; and H is the flow thickness in the direction perpendicular to the flow.The coefficients cxand cyare used to transfer the geometry from local to the global coordinate system.They are the cosine values of the slope angles between bed and horizontal plane in the x (αx) and y (αy) directions,respectively.The term Q is the mass source term which is used to incorporate the effect of bed entrainment.The bulking flow is controlled by the bed erosion.The volume of the flow increases by entraining the sediments along the flow path.The modelling of bed erosion in DFS 2D is controlled by the parameters: erosion rate (dz/dt),potential erosion depth (dz/ds),maximum erosion depth (emax) and critical shear stress (τz).DFS 2D uses a threshold value,known as critical shear stress(τcr),which decides the onset of erosion in any cell.Erosion starts only when the shear stress at any cell(τz)exceeds τcr(Eq.(3)).Then the erosion at every time step t is controlled by the potential erosion depth as mentioned in Eq.(4).

Also,the value of e at each time step is limited by

The maximum value of erosion in each cell is limited to emax,which can be provided as an input.The value can be found using subsurface explorations,assuming that erosion is possible up to the depth of bedrock.The soil layers above the bedrock may become eroded,if the shear stress due to flow exceeds the critical shear stress of the bed material.The value of τzis calculated depending on the flow rheology and is a function of velocity of flow.Each term in Eq.(2) represents different accelerations considered.The acceleration due to gravity is denoted by g.The convective acceleration is the time rate of change occurring due to the difference in position in space.It is represented in the left-hand side of Eq.(2),using the space derivatives.The right hand side considers the time or local acceleration; the first term represents the acceleration due to gravity,and the second term indicates the rate of change due to pressure differences (pressure acceleration) within the flow.Sx=tanαxand Sy=tanαyare the gradients of bed slope in x and y directions,respectively.The term Sfis the gradient of flow resistance,which is characterised by the rheology of the flow.The term k indicates the coefficient of earth pressure,which is defined as a function of the angle of internal friction φ.For perfect fluids,hydrostatic pressure is equal in all directions and the value of k remains 1 but when the material is plastic,the value can vary(Savage and Hutter,1989; Hungr,1995),depending upon the compressions and dilations within the flow.The magnitude of k in active(ka) and passive(kp) states are given as

The value of kashould be considered regions of expansion,where lateral stress is less than the normal stress defined by

The value of kpshould be considered in regions of compression:

In the passive case,the flow heights are smaller,and the flow is in a slow regime.This is due to the higher value of earth pressure coefficients.This formulation is used to control the deposition and spreading of mass in the flow path.

The terms Uxand Uyare the coefficients derived from the velocity values given by

where the negative sign is used to make sure the direction of Sfis opposite to the direction of velocity.

2.2.Rheology model

DFS 2D considers three different rheologies for modelling debris flows.The effect of rheology and the behaviour of flow are denoted by the term Sfin the momentum conservation equations.This variable is responsible for the dissipation of energy,and the governing equations can be used for different rheological behaviours by modifying the value of Sf.The depth integrated single phase model assumes constant flow properties,and the behaviour is characterised using the rheology.Such models can easily be used for the calibration of the flow parameters and thereby the calculation of runout paths of future debris flow events.Debris flows and mud flows in a single phase were modelled in a laminar flow regime and in recent models (Coussot,1994),grain inertial and friction-turbulent regimes are also considered(Brufau et al.,2000;Takahashi,2007; Medina et al.,2008; Christen et al.,2010; Rosatti and Begnudelli,2013; Frank et al.,2015).The laminar flow assumes viscoplastic material,with a linear stress-strain relationship beyond the yield point,as in the Bingham model.The Bingham model is a special case of the Hershel-Bulkley model,which considers the shear thinning behaviour and is in better agreement with the experiments.The following equation can be used to define both Bingham and Herschel-Bulkley models:

where the shear stress at any depth z is given by τz,?v/?z is the shear rate,μ is the viscosity parameter,and τcis the constant yield strength due to the cohesion of soil in the static state.The parameter β is empirical and takes the value of 1 in Bingham model.In a series of experiments conducted by Parsons et al.(2001),it was proved that addition of sand or clay to the slurry results in Bingham-like fashion.Also,Bingham rheology satisfactorily predicts the flow boundaries of 1997 Wartschenbach debris flow event(Austria) (Beguería et al.,2009).

When the constant yield strength is replaced with the Coulomb shear stress which is the combined effect of cohesion and friction in a material,the Coulomb viscous model is derived.The shear stress is defined as

where φ is the angle of internal friction,σ is the bed normal stress,and u is the pore pressure.The term (σ-u) indicates the effective normal stress.The model uses both frictional and viscous behaviour of the flow to describe the shear stress in granular particles (Naef et al.,2006).While comparing Coulomb viscous and Bingham model for modelling the 1997 Wartschenbach event,it was reported that Bignham model satisfactorily predicts the shape of flow,while Coulomb viscous model is better in estimating the deposition thickness (Beguería et al.,2009).For another event occurring in 2003 in Faucon (France),both the rheologies provided simulated the flow shape satisfactorily,yet Bingham model underpredicted the extent of flow and Coloumb viscous model overpredicted the same (Beguería et al.,2009).The findings from one-dimensional(1D) modelling of 2003 Kamikamihori event (Japan) reported that Coulomb viscous rheology resulted in flatter flow profile with more material remaining in the upsteam side.Overall bulk volume bevaiour of the predicte flow was found to be similar to the actual observations (Naef et al.,2006).

When the particles are of larger size,the assumption of laminar flow regime does not hold true (Liu et al.,2017).Hence,the turbulence of flow,induced by the interparticle interaction should also be considered for modelling.This is evaluated using a Voellmy-Salm rheology in DFS 2D.The Voellmy-Salm rheology has been developed from the Voellmy model(Voellmy,1955),to understand the flow behaviour of snow avalanches (Salm,1993).The model is based on the assumption that the concentration of shear deformations is near the basal surface (Salm,1993).It considers particles to be spherical.The model is found to be satisfactory in estimating the flow parameters of turbulent events at multiple locations (Christen et al.,2010; Hussin et al.,2012; Frank et al.,2015).The back analysis of the debris flow occurring in Faucon catchment in 2003 concludes that the rheology accurately estimates the flow volume and velocity (Hussin et al.,2012).The modelled runout distance had a variation of 1.4%from the observed value and the flow height had a variation of 44.8% from the observed value.The force of contact between the particles is defined by the angle of internal friction.The shear resistance due to impact of collision is given by

where ρ is the density of material,vgis the velocity,and ξ is the coefficient of viscous turbulent friction or simply turbulent friction with dimension LT-2.The term turbulent coefficient was first introduced by Voellmy (1955),in analogy to the turbulent open channel flow equation by Chezy.The parameter is a function of the terrain geometry and velocity.During modelling,the value is kept as constant and it needs to be calibrated for every event.The value of turbulent friction will be larger for steep terrains and fluid like flow.It decreases with the increase in resistance to flow (Bartelt et al.,2017).The velocity vgcan be expressed as (Salm,1993):

where α is the slope angle and μfis the dry-Coulomb friction.The value of Sfhas been calculated from the shear stress equation for each rheology model.The simplified Bingham expression of shear stress has been used to derive the value of Sffollowing Coussot(1997) as

For Voellmy-Salm rheology,the equation as mentioned in Christen et al.(2010) has been used as

2.3.Numerical solution

The complicated hyperbolic partial differential equations are simplified and expressed in vector form given by

where

The difficulties in implementing the hyperbolic partial differential equations are solved by adopting a simpler solution in the vector form.This form allows a satisfactory level of stability,controlled diffusivity,and required accuracy(Beguería et al.,2009).The model uses a finite difference mesh with fixed coordinates(i,j)for each timestep ‘n’ which can be solved using a Lax-Wendroff scheme.For the definition of mesh,the resolution of digital elevation model(DEM)is used.Hence,each grid will be of size s=Δx=Δy.Using a central difference forward scheme,Eq.(13) can be solved as

Fig.1.Interface for data input: (a) ASCII files and rheology,and (b) Parametric inputs.

where Δt is the time duration for each time step.The simplification results in possible dispersion which may lead to unphysical oscillations.The chances of such oscillations are much higher in the case of larger gradients.Hence,it is required to add a certain amount of numerical regularisation to stabilise the term in the right-hand side of Eq.(14) as

The function NR represents the numerical regularisation applied to the whole area in each time step.It does a weighted space averaging over the variable A.CFL is the Courant-Friedrichs-Lewy condition,which is a convergence condition used in the solution of hyperbolic partial differential equations.The quantity of regularisation is set based on the value of CFL.When there are sudden changes in flow values,the value of CFL may go greater than 1.Hence,multiple iterations are performed to keep the value of CFL below 1,providing only the minimum amount of regularisation(Beguería et al.,2009)and stability of the solution.When the mass is experiencing acceleration or deceleration,there can be nonuniformity issues for the flow resisting gradient.This issue is overcome using inner loops in each time step.The time Δt is divided into smaller steps to reduce the non-uniformity of the resistance term.

2.4.Inputs,outputs and user interface of DFS 2D

The model is implemented in Python and has used Tkinter(Moore,2018) for developing a user-friendly interface (Fig.1) for inputting the required data.The outputs are saved in different file formats in a user-specified folder.

The interface is designed in two different tabs.The first tab deals with the input files and rheology while the second tab deals with the input parameters.Both the tabs contain an‘INFO’button,which briefly describes the required inputs.The inputs required for the simulation and the outputs obtained are listed in Table 1.

Table 1 List of inputs required,and outputs obtained from DFS 2D.

The inputs can be broadly classified into ASCII file inputs and parametric inputs.The parametric inputs are again classified into soil parameters,entrainment parameters,and rheological parameters.The soil parameters can be obtained from laboratory testing,and the entrainment and rheological parameters should be calibrated by conducting a back analysis of past debris flow events.The calibrated values can be then used for calculating the flow behaviour of future debris flows.The working of DFS 2D can be summarised and shown in Fig.2.

In the ASCII files tab,the user must provide the four input maps in.txt format.The files should be in ASCII grid form and the maximum spatial extent will be for DEM.The spatial domain of all other input files should be within that of the DEM.All input files should be in a projected coordinate system with cell size in m.The release area file provides the possible area of release with the depth of release value in each cell.This file is used to calculate the instantaneous release volume at the beginning of the flow.Each cell in this file has the value of depth of failure at the particular location.Thus,the release volume can be calculated as the sum of depth values in each cell multiplied by the unit cell size of the release area file.The domain file defines the maximum possible extent of the flow.If the flow does not stop within this boundary,the user must enlarge the domain.The vegetation file is used to provide the surface cover of the DEM,which plays a crucial role in the bed entrainment.The user can also select the rheology of flow from this tab.

Soil parameters can be entered in the ‘Parameters’ tab.For making the model more user-friendly,a range of values for each parameter based on the soil type is provided in the button next to the input column.For example,the default value of dry-Coulomb friction is set based on the angle of internal friction and moisture content of the soil.Similarly,default value for viscosity is provided,based on the density and moisture content values,and the default value of critical shear stress is provided based on the soil type and surface cover file.All these values can be modified by the user,however,the default values provide an approximate idea about the range of values to be used in each case.These values are currently set based on the limited number of literature and back calculations and will be updated in future with more case studies.The user can define the maximum number of time steps also in this tab.All flow parameters ranging from density and rheological parameters to the entrainment parameters are provided in a single tab.After entering the values for all parameters,the user has to return to ASCII tab and click on the ‘Run Simulation’ button.The user will be asked to define a directory to save the output files.

Another critical part of a numerical solution is the boundary conditions provided.The DEM,which is the primary input,defines the basal boundary of the simulation.This also defines the computational domain and the size of the grid or mesh used for analysis.The release area defines the spatial boundary for the inlet.Flow starts from here and the maximum spatial limits of the flow are defined by the domain file.No flow is allowed outside this boundary,as defined by the user.The initial discharge is given by the depth of release values provided in the release area file.With respect to time,the simulation is continued until the user-defined number of time steps or until the time when the maximum flow momentum reaches 5%of the initial momentum.If this condition is not reached within the defined time steps,the user must increase the time steps.In case if the boundary cells of domain file have reached before the maximum time steps or momentum criterion,the user must enlarge the domain for next simulation.

The outputs are provided in three different formats:

(1).txt files for each time step,

(2) Heatmap in.jpg format for each time step,and

(3).mp4 video which demonstrates the changes in each time step.

The output contains the values of entrainment,flow height and velocity in each cell in each time step in the three different formats mentioned above.While the.txt file provides exact values for deriving any necessary flow parameters and for detailed analysis,the.mp4 video and heatmaps help in visualising the flow which helps in better understanding.The.mp4 video for entrainment is not provided,as the video outputs are solely for visualisation purpose.All files will be saved in the user-defined directory.The applicability of the numerical model has been evaluated by calibrating the friction parameters for a debris flow occurring in Yindongzi gully in China on 13 August 2010.

3.Back analysis of historical debris flows

3.1.Yindongzi gully

Yindongzi gully is located in Lianhe Village,in the Sichuan province of China (Fig.3).Catastrophic landslides have become an important concern in China(Yin,2011;Chen et al.,2015).It is on the right bank of the Baisha river,the general flow direction is from northeast to southwest,the catchment area is 1.1 km2,and the valley is deep and V-shaped.The blocking of rivers due to debris flows is increasing and such debris flows lead to cascadingevents which lead to huge losses(Zou et al.,2020).The average annual rainfall in the whole basin is about 1700 mm.The main gullyis 2.5 km long,with an elevation ranging from 1979 m to 1033 m.The higher elevation region,from 1979 m to 1560 m,is the clear water area,which is a funnel-shaped terrain with three sides surrounded by mountains and one exit,with a rain collecting area of 0.45 km2.It is a mountainous landform with shallow terrain and relatively high terrain difference,which is the source area of the ditch water system.

Fig.2.Schematic diagram depicting the workflow of DFS 2D.

Fig.3.Location details of the study area: (a) China,(b) Sichuan province,and (c) Yindongzi gully.

The gully is deep and has a steep terrain.This topographical condition allowed the rapid flow of debris flows(Wu and He,2015).A large amount of collapsed deposits formed by the Wenchuan earthquake (12 May 2008) on both sides of the valley provided many sediment sources for the formation of debris flows.

The lower elevation region,from 1150 m to 1033 m,is the deposition area,with relatively flat terrain,located at the outlet of Yindongzi gully (Fig.4).The debris flow deposits throughout the ditch piled up near the outlet of the gully.From the field observations,the debris flow deposits also contain gravel material,and around 15%of the total mass is contributed by particles of diameter greater than 0.5 m.The debris flow rushed out of the mountain pass,it first crossed the road and cascaded towards the Baisha river.Due to the narrowness of this section of the river,and the limited drainage capacity,the debris flow overflowed the channel and gets deposited on both sides of the river.

The gully has strong crustal activity,with faults developing from north to south.The rocks are highly fragmented with loose slope surface.The long-term weathering and erosion due to tectonic activities have led to the formation of loose granular deposits along the slope,which easily gets entrained along with the debris flow.The Wenchuan earthquake caused five major landslides and one shallow landslide in the debris flow formation area and circulation area of Yindongzi,which provided a rich source of debris flow.After the earthquake,multiple debris flows have occurred in the gully on 17 July 2009 and 13 August 2010.In this study,we considerd the debris flow event on 13 August 2010 using the pre-and post-event DEMs from Alos Palsar (advanced land observing satellite phased array type L-band synthetic aperture radar).According to meteorological data,on 12 and 13 August 2010,the maximum daily rainfall was 183.2 mm (13 August 2010),The total rainfall reached 136.6 mm in 6 h,and the maximum rainfall in 1 h is 90.6 mm(from 1:00 a.m.to 2:00 a.m.on 13 August 2010),with an average rainfall intensity of 22.8 mm/h.The heavy rainfall event has triggered multiple debris flows in the region and had killed a person.The debris flow solids almost filled up the valley,causing the original channel to change greatly and the deposits at the outlet of the gully.The high resolution terrain corrected data from Alos Palsar with 12.5 m resolution was used for the simulation.The pre-and postevent DEMs date back to June and December 2010,respectively,and volume of release used for the simulation has been decided from the field investigation data,as 4790 m3at the crown of debris flow.The erosion and deposition values in all cells range from 0 m to 7 m,calculated from the difference in elevation between pre-and post-event DEMs,and the values were used for the calibration of the friction parameters.

3.2.Numerical modelling

The pre-event DEM from Alos Palsar was used as the input for numerical modelling.The land was barren after the earthquake and previous debris flows,with loose material.Hence,uniform surface properties were provided across the domain.The soil was cohesionless,mostly granular,with a bulk density of 19 kN/m3,cohesion of 9 kPa and angle of internal friction of 24°.The basal friction angle was considered being 2/3 of the angle of internal friction,i.e.16°,and the moisture content was considered as 100%.The flow was turbulent in nature with more gravel particles and hence Voellmy-Salm rheology was chosen and multiple trials were carried out by varying the dry-Coulomb friction(μf)and viscous turbulent friction(ξ).The entrainment properties were fixed after initial trials.The erosion rate was 0.05 m/s with a critical shear stress of 1.5 kPa and maximum erosion depth of 10 m.The erosion rate was fixed to a certain value when shear stress overpasses a critical threshold.The potential erosion depth was kept constant as 0.2 m/kPa for barren land.The friction parameters for Voellmy-Salm rheology were calibrated using the shape (Abraham et al.,2021) and volume information of debris flow,as described in the next section.

3.3.Process of calibration

The conventional process of calibration of input parameters involves comparison of total volume of flow to the modelled volume,in order to find out the best suited parameters(Hussin et al.,2012).But this method is inappropriate when a coarse resolution DEM is used(Abraham et al.,2021).Due tothe coarse-resolution of DEM,the volume of flowateach cell will be higher,the error gets added up and the final volume will be much larger than the observed one.If the total volume of flow is used for calibration,the simulated shape will largely differ from the actual shape.From the perspective of risk reduction,the shape of debris flowand estimation of possible extent is important,in order to plan necessary actions of mitigation.Hence,the method of calibration should consider the shape of flow and the erosion or deposition occurring in each cell(Abraham et al.,2021).Twodifferent methods for calibration are proposed in this study.The first one uses different cross-sections along the flow channel,in order to understand the effect of each parameter on the modelled shape(width,runout path)and topographic changes.In the second method,the modelled topographic changes in each cell are compared with the observed values,to find the simulation with a minimum variation from observed values.

3.3.1.Comparison of cross-sectional area

For the process of calibration,20 different sections along the flow path were considered (Fig.4).The cross-sectional area of the debris flow path at each section was calculated using the pre-and post-event DEMs and the results were used for calibrating the friction parameters using DFS 2D.Fig.5 describes the process of calculation of cross-sectional area.The change in elevation profile is calculated using the difference between pre-and post-event DEMs,to calculate the actual change in cross-sectional area of erosion or deposit along a section.For example,the bottom boundary in Fig.5 is the ground profile after the flow and the cross-section at ‘aa’represents the eroded profile.Similarly,deposition can occur at some sections and in some sections,both erosion and deposition may occur.An eroded section is shown in Fig.5 for representation.The area of cross-section can be calculated using the pre-and postevent DEMs and can be compared with the modelled topographic changes,for calibrating the friction parameters.

3.3.2.Cell-by-cell analysis

While comparing thecross-sectional areasin differentsections,it was observedthat the results werebetterwhenμf=0.1,yet there are chances that some data are missed while doing the section wise comparison.To overcome this limitation,the calculated change in elevationwas compared with the actual difference between the preand post-event DEMs.Using this approach,the modelled erosion/deposition values at each cell is compared with the observed values and the simulation with a minimum variation from the observed values can be identified.The comparison is done cell by cell,using different statistical indices,such as Willmott’s indexof agreement,d(Eq.(20)),mean-absolute-error,MAE (Eq.(21)),and normalizedroot-mean-square-error,RMSE (Eq.(22)).Willmott (1981) proposed an index of agreement (d) as a standardised measure of the degree of model prediction error which varies between 0 and 1.The value represents the ratio of the mean-square-error and the potential error and d value of 1 indicates a perfect match,and 0 indicates no agreement at all.MAE calculates the average magnitude of error in the modelled values and RMSE follows a quadratic scoring rule to calculate the average error magnitude.

Fig.4.Details of debris flow path on 13 August 2010: (a) Elevation profile,and (b) Shape of flow and different cross-sections considered for calibration process.

Fig.5.A schematic diagram depicting the modified channel profile and cross-sectional area considered for calibration.

where n is the total number of cells through which flow has occurred,Dobs,iis the observed topographical variation obtained from pre-and post-event DEMs,and Dmod,iis the modelled elevation difference obtained after simulation.

4.Results

The Voellmy-Salm model was chosen to simulate the turbulent debris flow event.The calibration of friction parameters has been done for constant entrainment parameters.For the process of calibration,both μfand ξ were varied and multiple trials have been conducted.The values of μfwere varied from 0.01 to 0.4 and those of ξ were varied from 100 m/s2to 2000 m/s2.The size of particles were evaluated during field investigations and it was found that more than 50% of the debris material was composed of rock fragments and granular material.The upper limit of μfwas considered as 0.4,considering the larger particle size of debris flow materials at the site.

The initial trials were conducted by varying the values of friction parameters,and those trials which largely vary from the actual shape of debris flow were discarded.The process of calibration is explained using the best fitted parameters,as the effect of one parameter is best explained when the other one is kept constant.The effect of dry-Coulomb friction was evaluated for constant value of ξ.The maximum flow height simulated at each cell for different values of μfis plotted in Fig.6.

From Fig.6,it can be observed that the flow height is increased at locations where a sudden change is observed in flow direction or elevation.The maximum height in all cases is recorded in the zone of deposition.The width of flow increases as the value of μfincreases.The maximum height of flow increased when the value of μfwas increased from 0.01 to 0.1 and further increase in μfresulted in the decrease of maximum flow height.It should also be noted that,as the value of μfis increased,the locations with larger flow height are also increased.Since both width and height are equally important in modelling debris flows,the cross-sectional areas at each section are used for quantitative comparison,to calibrate the friction parameters.

During the flow,a cell can be subjected to both erosion and deposition.The unique available observation is the final cell elevation and therefore,erosion or deposition can be associated with a cell,based on the values at last time step.The cross-sectional areas of 20 different sections were then compared with the observed values for a better understanding of the effect of varying the friction parameters.Fig.7 shows the areas of cross-section at different sections considered (in logarithmic scale),for different values of μf.

It can be observed from Fig.7 that up to section 7-7,all the simulations are in good agreement with the observed values.This is also evident from Fig.7 that the initial sections of the flow with the steep slope have similar values in all simulations.Even though the maximum flow height is reduced when the value of μfis increased beyond 0.1,the enlarged widths result in increased cross-sectional areas.The simulations with μf=0.4 have the maximum crosssectional area in all cases and are greatly varying from the observed values.The maximum variation is observed at zone of deposition,where higher values of μfcalculate much larger crosssectional areas.At section 17-17 at the beginning of deposition zone,none of the simulations are in good agreement with the observed value.

Along with the calculation of cross-sectional areas,a cell-by-cell analysis is also used for calibration.The statistical indices for different values of μfare mentioned in Table 2.It can be inferred from Table 2 that the performance is better when μf=0.1.Any deviation from this value(increase or decrease)results in increased error and reduced d value.

Once the value of μfhas been calibrated,several trials were done by varying the value of ξ to find the optimum value.The flowheights for different trials are plotted in Fig.8.The maximum flow height initially increased from 10.26 m to 11.26 m when the value of ξ was increased from 100 m/s2to 400 m/s2.When the value of ξ was further increased beyond 400 m/s2,the maximum height of flow was found to be decreasing.The variation was negligible for an increase in ξ from 1000 m/s2to 2000 m/s2.Also,the distribution of flow height along the flow path was similar in the case of all simulations.The width of flow was found to be decreasing with an increase in the value of ξ.

To understand the effect of ξ in detail,the cross-sectional areas at different sections(in logarithmic scale)were compared with the observed values and are plotted in Fig.9.Similar to the case of μ,the maximum deviation from observed values was found in the zone of deposition,from section 17-17.The cross-sectional areas calculated by simulations with ξ less than 1000 m/s2were larger than the observed values,and those with ξ=2000 m/s2were found to be lesser than the actual areas in most of the sections.

Apart from the section wise analysis,a cell-by-cell analysis was also carried out for all the simulations (Table 3) and it was found that the results are better when ξ=1000 m/s2.The maximum value of d (0.5) and the minimum values of MAE and RMSE are obtained when ξ=1000 m/s2.

Table 2 Results of cell-by-cell evaluation:Statistical indices derived for different values of μf(ξ=1000 m/s2).

Table 3 Results of cell-by-cell evaluation:Statistical indices derived for different values of ξ(μ=0.1).

Both sections’ wise and cell-by-cell analyses imply that the model with μ=0.1 and ξ=1000 m/s2performs better than the other models.The calculated flow parameters can be evaluated in detail to understand the flow behaviour.

5.Discussion

The friction parameters for the debris flow event were calibrated as μf=0.1 and ξ=1000 m/s2,based on the analysis.The d value for the best suited model (0.5) indicates the mutual agreement between the observed and modelled elevation changes.The MAE and RMSE values are also the least for the same model,and the comparison of cross-sectional area also agrees with the same calibrated values.Considering multiple approaches ensures fine calibration and better estimation of flow parameters.The major variations in the observed and modelled values were noted in the zones of deposition,where the values are overestimated by the model,and the reasons should be explored in detail.Since there is no record of the original flow velocity or the change of flow parameters with respect to time,only the final values are compared with the available data from pre-and post-event DEMs.The agreement between observed and modelled values of depth of erosion along a longitudinal section through the debris flow channel is depicted in Fig.10.The comparison of cross-sectional area along the flow path shows that the observed and calculated values are in good agreement with each other (Fig.10),except in the case of zone of deposition.The deposition mostly occurred in the river,which cannot be evaluated using the difference between pre-and postevent DEMs.Hence,the values of deposition is not considered for the performance evaluation of DFS 2D.

Fig.6.Effect of dry-Coulomb friction coefficient on the debris flow height (ξ=1000 m/s2).

Fig.7.Scatter plot of observed and simulated cross-sectional areas at different sections(ξ=1000 m/s2).

The flow height,velocity and depth of entrainment at each cell were calculated by DFS 2D and saved as a text file.DFS 2D allows for bed entrainment only if the shear stress induced by the flow is greater that the user-defined critical shear stress.In this analysis,the critical shear stress was considered as 1.5 kPa.The bed entrainment was less in sections with steep slope,even though the velocity values were higher (sections 1-1 to 8-8).The maximum entrainment value calculated by the simulation is 7.5 m,which is in good agreement with the observed value of 7 m.The release volume of 4790 m3was increased to 265,342 m3at the end of simulation,due to bed entrainment.The calculated value is much higher than the value estimated from field observations.This is due to the poor resolution of DEM used for simulation.The coarse-resolution DEM overestimates the total volume of flow,which limits the use of total volume or total erosion for the use of calibration.This has been overcome in this study using two different approaches for calibration,using the comparison of cross-sectional area and cell-by-cell analysis,instead of comparing the total volume of flow.The effect of the cumulative error in volume and erosion due to coarse-resolution DEM can be avoided using the calibration methods proposed in this study.The limitation of this method is that the proposed calibration method over total volume approach may lead to the selection of parameters that would produce a more dispersed flow with respect to the observed deposit,since the shallow distal debris flow modelled deposit would not influence much the overall grand total.However,our major concern is the area affected due to debris flows,from the perspective of risk management.Hence,this approach can be considered over the total volume approach,when the DEM resolution is poor.

Fig.8.Effect of viscous turbulent friction coefficient on the debris flow height (μf=0.1).

Fig.9.Scatter plot of observed and simulated cross-sectional areas at different sections (μ=0.1).

From Fig.11,it can be understood that the flow velocity is maximum at the middle of cross-sections.Towards the edges,the velocity decreases and eventually the lateral spread of flow is limited.Along a longitudinal section,the velocity first increases and then drops (Fig.12) when the flow changes its direction.At locations closer to crown,bed entrainment is more than the maximum flow height,while the flow height is larger than bed entrainment at locations beyond 400 m from crown.

Fig.10.Comparison between observed and modelled flows along the debris flow path: (a) Area of cross-sections,and (b) Depth of erosion.

Fig.11.Calculated flow parameters in each cell.

Also,after failure,the flow started with very high velocity,which got decreased with the decrease in slope.At flat terrains,the velocity of flow and the entrainment depth are decreased.The velocity almost reduces to zero at locations where flow changes its direction.Thus,the modelled flow almost stops at these locations,which results in larger values of deposition at these zones.The large reduction in velocity values calculated by the model results in the overestimation of deposition values.But the calculated values in all other locations are found to be satisfactory.

DFS 2D uses a simple linear model for calculating the erosion values from the case study,and the modelled values were found to be in satisfactory agreement with the observed values.The simple model helps in reducing the computational time and unphysical oscillations during the calculation.However,the agreement with observed values can be evaluated using more case studies and limitations if any can be rectified in future updates.

Fig.12.Maximum values of height,velocity and entrainment recorded along the flow path.

Hence,using DFS 2D,we can calculate the flow height,velocity and bed entrainment at each cell and from the obtained results,the impacts of flow,energy,etc.can be calculated.Such information is crucial in designing any mitigation works along debris flow path(Bernard et al.,2019).The information can also be made useful in quantitative evaluation of risk due to debris flows on any structure along the path.

6.Conclusions

A new debris flow simulation model,DFS 2D,is proposed in this study,which can be used to simulate debris flows of varying moisture contents and particle sizes,considering multiple rheologies.The model considers shallow water conditions and uses the mass and momentum conservation equations as the governing equations and the increase in volume of flow due to bed entrainment.The interactive graphical user interface of the model allows the user to input the DEM,domain,release area and surface cover maps in ASCII grid format along with the soil,entrainment and rheological parameters.

The applicability of DFS 2D was evaluated using the data from a major debris flow event occurring in China on 13 August 2010.The pre-and post-event DEMs of the debris flow were used to compare the simulated and observed results.Two different methods were adopted for the process of calibration,one by considering 20 different cross-sections along the flow path,and one by comparing the absolute error of elevation change in each cell considered for the analysis.The analysis proved that the flow parameters were satisfactorily calculated using Voellmy-Salm rheology.The rheological parameters were calibrated as μf=0.1 and ξ=1000 m/s2.The statistical indices such as Willmott’s index of agreement,MAE and RMSE of the calibrated model are 0.5,1.02 and 1.44,respectively.

From the analysis,it is understood that DFS 2D can help in understanding the runout process of debris flows and the proposed calibration techniques provided a better solution for back analysis of debris flows using coarse-resolution DEMs.The model can be used for calibrating the rheological and entrainment parameters for historical debris flows,which can help in calculating the extent of possible failures in the future.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

This research was financially supported by Department of Space,India(Grant No.ISRO/RES/4/663/18-19).The authors would also like to express their sincere gratitude to Institute of Mountain Hazards and Environment,Chinese Academy of Sciences for their support in field investigations.