Hazard Assessment Comparison of Tazhiping Landslide Before 2 and After Treatment

3 Dong Huang 1*, JianPing Qiao1, Meng Wang1 4 5 1. Key Laboratory of Mountain hazards and Surface process, Institute of Mountain hazards 6 and Environment, Chinese Academy of Science, Chengdu 610041, China 7 *Corresponding author (dhuang@imde.ac.cn). 8 9 Abstract: Through investigation and analysis of geological conditions and 10 mechanical parameters of the Taziping landslide, the finite volume method was 11 adopted, and, the rheological model was adopted to simulate the landslide and 12 avalanche entire mass movement process. The present paper adopted the GIS platform 13 to simulate the mass movement process before and after treatment. This paper also 14 provided the conditions and characteristic parameters of soil deposits (thickness, 15 speed, and stresses) during the landslide mass movement process and mapped the 3D 16 division of hazard zones before and after landslide treatment. Results indicated that 17 the scope of hazard zones contracted after engineering treatment of the landslide. The 18 extent of high-hazard zones was reduced by about 2/3 of the area before treatment, 19 and characteristic parameters of the mass movement process after treatment decreased 20 to 1/3 of those before treatment. Despite engineering treatment, the Taziping landslide 21 still poses significant hazard to nearby settlements. Therefore, we propose that houses 22 located in high-hazard zones be relocated or reinforced for protection. 23 24 25


Introduction
The hazards of a landslide include scope of influence (i.e., source area, possible path area, and backward and lateral expansion area) and secondary disasters (i.e., reservoir surge, blast, and landslide-induced barrier lake).A typical landslide hazard assessment aims to propose a systematic hazard assessment method with regard to a given position or a potential landslide.Current research on typical landslide hazard assessment remains immature, and there are multiple methods for interpreting landslide hazards.To be specific, the scope of influence prediction of a landslide refers to deformation and instability characteristics such as sliding distance, movement speed, and bulking thickness range.The movement behavior of a landslide mass is related to its occurrence, sliding mechanisms, mass characteristics, sliding path, and many other factors.Current landslide movement prediction methods include empirical prediction and numerical simulation.
The empirical prediction method involves analyzing landslide flow through the collection of landslide parameters in the field.It further consists of the geomorphologic method (Costa, 1984;Jackson et al., 1987;Scott and Vallance, 1993), the geometric change method (Finlay et al., 1999;Michael-Leiba et al., 2003), and the volume change method (Fannin and Wise, 2001).Empirical models are commonly simple and easy to apply, and the required data are easy to obtain as well.Numerical simulation methods are further divided into the continuous deformation analysis method (Hungr, 1995;Evans et al., 2009;Wang, et al., 2016), the discontinuous deformation analysis method (Shi, 1988), and the simplified analytical simulation method (Christen et al., 2010a;Sassa et al., 2010;Bartelt et al., 2012;Du et al., 2015).The numerical simulation method expresses continuous physical variables using the original spatial and temporal coordinates with geometric values of discrete points.Numerical simulations follow certain rules to establish an algebraic equation set in order to obtain approximate solutions for physical variables.
Empirical prediction models only provide a simple prediction of the sliding path.Due to the differences in geological environments, empirical prediction models commonly have low generality.Landslides move downslope in many different ways (Varnes, 1978).In addition, landslides can evolve into rapidly traveling flows, which exhibit characteristics of debris flows on un-channelized or only weakly channelized hillslopes.The geomorphic heterogeneity of rapid shallow landslides, such as hillslope debris flows, is larger than observed in channelized debris flows; however many of these flows can be successfully modeled using the Voellmy fluid friction model (Christen et al., 2012).The selection of model parameters remains one of the fundamental challenges for numerical calculations of natural hazards.
The continuous deformation method has the advantage of an extremely strong replication capability, but it is not recommended when analyzing flow-type landslides, lahars, or debris flows because of complicated rheological behaviors (Iverson et al., 1997;Iverson and Vallance, 2001;Hungr et al., 2001;Glade, 2005;Portilla et al., 2010;Chen et al., 2014).The fluid mechanics-based discontinuous deformation method has several shortcomings such as great computational burden, difficult parameter selection, and difficult 3-D implementation.The simplified analytical simulation method fully takes into account the flow state properties of landslides before introducing a rheological model and can easily realize 3-D implementation on the GIS platform.On that account, this paper adopted the continuous fluid mechanicsbased finite-volume method (simplified analytical simulation method).We introduce a rheological model on the basis of using mass as well as momentum and energy conservation to describe the movement of landslides.We also employed GIS analysis to simulate the entire movement process of Tazhiping landslide and map the 2-D division of hazard zones.

Kinetic analysis method
Adopting the continuous fluid mechanics-based finitevolume method, this paper took into account erosion action on the lower surface of the sliding mass and the change in frictional resistance within the landslide debris flow in order to establish a computational model.The basic idea is to divide the calculation area into a series of non-repetitive control volumes, ensuring that there is a control volume around each grid point.Each control volume is then integrated by the unresolved differential equation in order to obtain a set of discrete equations.The unknown variable is the numerical value of the dependent variable at each grid point.To solve the integral of a control volume, we make a hypothesis about the change rule of values among grid points, that is, about their piecewise distribution profile.The finite-volume method can satisfactorily overcome the finiteelement method's weakness of slow calculation, and solve the problem of complex region processing.Thus, we adopted the finite-volume method to establish the kinematic model for the landslide flow process.
The core of the finite-volume method is domain discretization.The finite-volume method uses discrete points as a substitute for continuous space.The physical meaning of the discrete equation is the conservation of the dependent variable in a finite control volume.Establishment of the conservation equation is based on the continuous movement model, that is, the continuity hypothesis about landslide substances.We divided the landslide mass into a series of units and made the hypothesis that each unit has consistent kinematic parameters (speed at a depth, density, etc.), and physical parameters (Fig. 1).We also established an Eulerian coordinate systembased conservation equation with regard to each control volume.

Control equation
The computational domain is defined as directions x and y, and the topographic elevation is given the coordinate z(x, y).H (x, y, t) is assumed as the change relationship of landslide thickness with time; U x (x, y, t) and U y (x, y, t) respectively represent the mean movement speeds along directions x and y at moment t; n x = U x / U 2 x + U 2 y and n y = U y / U 2 x + U 2 y represent the cosinoidal and sinusoidal flow vectors of the landslide on the plane x-y.The mean flow speed of substances is defined as U = U 2 x + U 2 y .Thus, the mass balance equation becomes wherein Q(x, y, t) represents the change rate (entrainment rate) of landslide volume with time.
Assuming that l(x, y, t) represents the movement distance of the landslide with time, we can obtain wherein h i represents the thickness of the ith layer of the landslide in the movement process; ρ i represents the density of the ith layer of the landslide in the movement process; ρ a represents the density of the landslide; the dimensionless parameter k i represents the entrainment rate.
The momentum balance equation is wherein S gx = g x H and S gy = g y H represent the dynamic components of the acceleration of gravity in directions x and y; g = (g x g y g z ) represents the vector of the acceleration of gravity; k a/p represents the pressure coefficient of soil; ρ a represents the density of the landslide; the dimensionless parameter k i represents the entrainment rate; S f (R) represents the frictional resistance.
The kinetic energy balance equation is wherein R(x, y, t) represents the random mean kinetic energy of the landslide; Ṗ (x, y, t) and Ḋ(x, y, t) represent the random increased kinetic energy and decreased kinetic energy of the landslide.

Constitutive relationship
The improved Voellmy rheological model is applied in the computational simulation of the landslide.See the computational formula below: wherein u i / U represents the unit vector in the movement direction of the landslide; µ represents the Coulomb friction coefficient, and is related to R(x, y, t), the random mean kinetic energy of the landslide; R t represents the gravityrelated frictional force coefficient; K represents the substrate surface curvature; ζ represents the viscous friction coefficient of the "turbulent flow".

HLLE-Heun numerical solution
Synthesizing control Eqs. ( 1), ( 3), (4), and (5), we can obtain the simplified form of the nonlinear hyperbola equation: wherein V(x, y, t) represents a vector equation consisting of four unknown vector variables; F (V) represents the flux function; G (V) represents the source term.Based on the HLLE (Einfeldt, 1988) equation of the finite-volume method and the quadrilateral grid, the node layout can adopt the grid center pattern, and the normal flux along one side of the control volume can be represented by the flux at the center of the side.The finite-volume discretization adopting the control volume as unit is depicted in Fig. 1; the Gauss theorem can be followed for the integration of Eq. ( 8), wherein C i represents the unit volume.After converting the volume integral flux function F (V) into the curved surface integral, we can obtain wherein n i represents the outward normal direction vertical to unit C i at the boundary; through adopting the HLL (Harten et al., 1983) format for the discretization of surface integral, the following simplified form can be obtained: (n) i represents the mean value of unit variables at moment t (n) ; V (n) represents the mean value of the entire grid at moment t (n) ; t := t (n−1) − t (n) wherein n ij represents the outward normal direction of the ith unit at boundary j ; the flux calculation term F (HLL) ij V (n) represents the approximate solution mode of the Riemann problem of the ith unit at boundary j ; see the computational formula below: R respectively represent the approximate values of V (n) on both sides of boundary j of the ith unit; S L and S R respectively represent the wave speeds on the left and right sides.Refer to the computational method described by Toro (1992).In addition, the gradient magnitude in the original second-order difference equation can be limited through multiplication with the flux limiter, and the second-order format of the total variation diminishing (TVD) property can be constructed to avoid the occurrence of numerical oscillation.Refer to the specific method described by LeVeque (2002).
In this paper a numerical solver within RAMMS is used, which was specifically designed to provide landslide (avalanche) engineers with a tool that can analyze problems with two-dimensional depth-averaged mass and momentum equations on three-dimensional terrain using both first-and second-order finite-volume methods (Christen et al., 2010b).Therefore, the finite-volume method is adopted to analyze the flow type (high mobility, high velocity, large scope of risks, etc.) of the landslide mass movement process.The present paper adopts the numerical approach of RAMMS and the GIS platform to simulate the mass movement process before and after treatment.The landslide depositional characteristics and the mass movement conditions can be combined to provide a scientific basis for engineering prevention , control, and forecast risk assessments for these kinds of disasters.
3 Study area and data

Tazhiping landslide
The Tazhiping landslide is located southeast of the Hongse village, Hongkou town, Dujiangyan city of Sichuan Province.The site is located at (103  urban district (Fig. 2).Its geomorphic unit is a middlemountain tectonic erosional area on the north bank of the Baisha River valley.The Tazhiping landslide is a large-scale colluvial layer landslide triggered by the Wenchuan earthquake (Fig. 3).It has a gradient of 25-40 • with an average gradient of 32 • .The landslide has an apparent round-backed armchair contour with a steep rear edge, which has a gradient of 35-50 • and an elevation of about 1370 m.The front edge is located on the south side of the mountain road, and has an elevation of about 1007 m.The landslide has an elevation difference of about 363 m, and a main sliding direction of 124 • NE.The landslide mass forms an irregular semielliptical shape, and has a length of about 530 m, an average width of 145 m and an area of approximately 7.68 × 10 4 m 2 .The landslide mass is composed of gravelly soil and is covered on by silty clay mingled with gravel.In terms of spatial distribution, the landslide is thick in the middle and thin on the lateral edges, has a thickness of 20-25 m and a volume of approximately 1.16 × 10 6 m 3 .During the earthquake, the landslide mass slid to cover the northern mountain slope of the Hongse village Miaoba settlement.The landslide has an apparent front edge boundary, and there is also a swelling deformation (Fig. 4).After the Wenchuan earthquake, the massive colluvial deposits covered the mountain slope.The colluvium is 0.5-5.0m thick at the top of the slide and is composed of rubble and gravel.The mass consists of a small amount of fine gravel, which is composed of gray or grayish-green andesite with a clast of 20-150 cm.Field surveys indicate that the rubble in the surface layer has a maximum diameter exceeding 2 m, and that fine gravel is loosely intercalated with the rubble.A small amount of yellowish-brown and gray-brown silty clay mixed with 5-40 % of non-uniformly distributed rubble composed the first 5-10 m of the slide.From 10 to 25 m deep, there is a wide distribution of gravelly soil.The soil is grayish-green or variegated in color, is slightly compact and non-uniform, and has a rock fragment content of about 50 %.The parent rock of the rock fragments is andesite, filled with silty clay or silt (Fig. 5).Table 1 shows the parameters of the surface gravelly soil of the landslide mass based on the field sampling.The landslide is an unconsolidated mass containing relatively large amounts of crushed stones and silty clay (Fig. 6).Its loose structure and strong permeability facilitate infiltration of surface water.The Wenchuan earthquake aggravated the deformation of the landslide making deposits more unconsolidated, further reducing the stability of the landslide mass.During persistent rainfall, surface water infiltrates the landslide slope resulting in increased water pressure within the landslide mass and reduced shear strength on the sliding surface.Thus, rainfall constitutes the primary inducing factor of the upper Tazhiping landslide.After infiltrating the loose layer, water saturates the slope increasing the dead weight of the sliding mass and reducing the shear strength of soil in the sliding zone.Infiltration into the landslide mass also increases the infiltration pressure of perched water, drives deformation, and poses a great threat to villages located at the front of the landslide.Slide-resistant piles and backfill were place at the toe of the slope in order to reduce the hazards of future slides.The slide-resistant piles have enhanced the overall stability of the slope; however, under heavy rainfall the upper unconsolidated landslide deposits may cut out from the top of the slide-resistant piles.
Therefore we simulate possible movement states of the Tazhiping landslide before and after treatment with slideresistant piles, comparatively analyzed the kinetic parame- ters in the movement process, and mapped the 2-D division of hazard zones.

Hazard prediction before treatment
It was assumed that the landslide was damaged before engineering treatment.According to field investigation, the sliding mass had an estimated starting volume of about 600 000 m 3 and a mean thickness of 8 m.Based on the survey report and field investigation (Hydrologic Engineering and Geological Survey Institute of Hebei Province, 2010), we adopted the survey parameters of Table 2 for the simulated calculation.These parameters were obtained from laboratory or small-scale experiments and back-analyses of relatively well-documented landslide cases.The unit weight γ = 20.8kN m −3 is from small-scale conventional triaxial test experiments in laboratory.In addition, we selected the coulomb friction coefficient µ = 0.45 and viscous friction coefficient ζ = 500 m s −2 in accordance with back-analyses of well-documented landslide cases (Cepeda et al., 2010;Du et al., 2015).The erosional entrainment rate selected was the minimum value k i = 0.0001 in the RAMMS program.See the kinematic characteristic parameters of the landslide deposits in Fig. 7.The colored bar shows the maximum values of the kinematic process for a given time step.As shown by the calculation results, deposits accumulated during the landslide movement process had a maximum flow height of 23.85 m, located around the surface gully of the middle and upper slope.The middle and lower section of the landslide deposit had a flow height of about 5-10 m; the middle and lower movement velocity of the landslide ranged from 3 and 7 m s −1 ; the landslide had a mean pressure of about 500 kPa, and the pressure of the middle and lower deposits was about 200 kPa.Thus, houses of three stories or less within the deposition range might be buried (buildings with 3 m high floors), and it was further suggested that the design strength of the gable walls of houses on the middle and upper parts of the deposit be increased above 300 kPa.

Hazard prediction after treatment
After fully accounting for the slide-resistant piles and mounds, we introduced the Morgenstern-Price method (Morgenstern and Price, 1965) to calculate the stability coefficient of Tazhiping landslide after treatment.The method was determined with an iterative approach by changing the position of the sliding surface until failure of the dump site (Fig. 8).The physico-mechanical parameters under a saturated state (Hydrologic Engineering and Geological Survey Institute of Hebei Province, 2010) were adopted to search for the sliding plane of the landslide.
Based on numerical analysis, the Tazhiping landslide stability coefficient is 0.998.Under rainfall conditions, the middle area of the Tazhiping landslide was unstable.Loose deposits in the middle part of the landslide might convert into a high-water landslide and cut out from the top of the slideresistant piles.In the damaged area, the slope had a rear edge wall elevation of about 1170 m.Its front edge was located on the south side of the mountain road, with an elevation of 1070-1072 m and a length of 182 m.Thus, the scale of the rainfall-damaged area is estimated to be about 250 000 m 3 , with a mean thickness of about 6 m.The parameters in Table 2 were again adopted for the simulated calculation.
Provided in Fig. 9 are the kinematic characteristics of the landslide deposit.The colored bar shows the maximum values of the kinematic process for a given time step.Deposits accumulated during the landslide movement process had a maximum flow height of 18.37 m, located around the surface gully of the middle and upper slope.The middle and lower portions of the landslide deposit had a flow height of approximately 3-5 m.The middle and lower movement velocity of the landslide deposits ranged between 3 and 5 m s −1 .The landslide had a mean pressure of about 330 kPa, and the pressure of the middle and lower deposits was about 100 kPa.Thus, it could be held that houses of two stories or less within the deposition range might be buried.It was further suggested that the design strength of the gable walls of houses on the middle and upper parts of the deposits be increased above 150 kPa.
After treatment, the accumulation flow height and pressure of the deposits were reduced by about one-half, and the kinematic speed is reduced by about 1/3.However, the Miaoba residential area of Red Village was still partially at risk.

Results
Landslides reflect landscape instability that evolves over meteorological and geological timescales, and they also pose threats to people, property, and the environment.The severity of these threats depends largely on landslide speed and travel distance.There may be examples where entire houses on a landslide mass are moved but not destroyed because of stable base plates.In any case, velocity plays a more important role regarding kinetic energy acting on an obstacle.However, the Miaoba residential area of Red Village is located at the frontal part of Tazhiping landslide.During landslide movement, the spatial scale indexes of a landslide mass include area, volume, and thickness.The maximum thickness of the landslide is one of the direct factors influencing the building's deformation failure status.A large landslide displacement may lead to burial, collapse, or deformation failure of the building, and thus influences its safety and stability.Thus, landslide thickness constitutes an important index for assessing the hazards of a landslide disaster, and for influencing the consequences faced by disaster-affected bodies (Fell et al., 2008;DZ/T 0286-2015DZ/T 0286- , 2015)).Provided in Table 3 is a landslide thickness-based division of the predicted hazard zones of Tazhiping landslide, in which the thickness of the landslide mass correlates with the ability of a building to withstand a landslide disaster (Hungr et al., 1984;Petrazzuoli et al., 2004;Glade et al., 2006;GB 50010-2010GB 50010- , 2010;;Hu et al., 2012;Zeng et al., 2015).After treatment with slide-resistant piles, the hazard of a future slide was reduced by about onethird overall and by two-thirds in high-hazard zones.
The hazard zones of Tazhiping landslide was given by 2-D divisions before and after engineering treatment (Fig. 10).The size of the hazard zones changed after engineering treat-   ment, particularly in the high-hazard zones.Before treatment with slide-resistant piles, the landslide posed a great hazard to eight houses on the left side of the upper Miaoba residential area, with a high-hazard zone associated with landslide mass height over 5 m and a red zone.After treatment, the number of effected houses was reduced to four.We defined outside the colored area as hazard-free.

Conclusions and discussion
The hazard assessment of landslides using numerical models is becoming more and more popular as new models are developed and become available for both scientific research and practical applications.There is some confusion about the mass movement process that is discussed by the rheological model presented in this contribution.
Landslides move downslope in many different ways (Varnes, 1978).In addition, landslides can evolve into rapidly traveling flows, which exhibit characteristics of debris flows on un-channelized or only weakly channelized hillslopes.The geomorphic heterogeneity of rapid shallow landslides, such as hillslope debris flows, is larger than observed in channelized debris flows; however many of these flows can be successfully modeled using the Voellmy fluid friction model (Christen et al., 2012).Results presented in this paper support the conclusion that Voellmy fluid rheological model can be used to simulate flow-type landslides.
The selection of model parameters remains one of the fundamental challenges for numerical calculations of natural hazards.At present, there are numerous empirical parameters obtained from 30 years of monitoring data.Such as in RAMMS, we can automatically generate the friction coefficient of an avalanche for our calculation domain based on topographic data analysis, forest information, and global parameters (WSL, 2013).The friction parameters for debris flows can found in some literature (Fannin and Wise, 2001;Iovine et al., 2003;Hürlimann et al., 2008;Scheidl and Rickenmann, 2010;Huang et al., 2015).However, there is little research regarding friction parameters of flow-type landslides.Therefore, we tested different coulomb friction coefficient µ values in the range of 0.1 ≤ µ ≤ 0.6 and viscous friction coefficient ζ values ranging between 100 ≤ µ ≤ 1000 m s −2 .Finally, we selected the coulomb friction coefficient µ = 0.45 and viscous friction coefficient ζ = 500 m s −2 in accordance with back-analyses of well-documented landslides (Cepeda et al., 2010;Du et al., 2015).Simulation results are consistent with field observations of topography and sliding path.
Based on the finite-volume method and the RAMMS program, simulation results of Tazhiping landslide were consistent with the sliding path predicted by the field investigation.This correlation indicates that numerical simulation is an effective method for studying the movement processes of flowtype landslides.The accumulation flow height and pressure of landslide deposits were reduced by about one-half, and the kinematic speed was reduced by about one-third after treatment.However, the Miaoba residential area of Red Village is still partially at risk.Considering that houses of two stories or less within the deposition range might be buried, it was further suggested that the design strength of the gable walls of houses on the middle and upper parts of the deposit be increased above 150 kPa.
By utilizing a GIS platform in combination with landslide hazard assessment indexes, we mapped the 2-D division of the Tazhiping landslide hazard zones before and after engineering treatment.The results indicated that overall hazard zones contracted after engineering treatment and, the area of high-hazard zones was reduced by about two-thirds.After engineering treatment, the number of at-risk houses on the left side of the upper Miaoba residential area was reduced from eight to four.It was thus clear that some zones are still at high hazard despite engineering treatment.Therefore, it was proposed that houses located in high-hazard zones be relocated or reinforced for protection.

Figure 6 .
Figure 6.(a) Material on the landslide surface.(b) Material in the shear zone.Photographs showing colluvial deposit cover on the mountain slope.

Figure 8 .
Figure 8. Search for the sliding plane of the Tazhiping landslide (before treatment).
represents the calculated time step; A C i represents the area of unit C i ; F represents the approximate value of the curved surface integral, as www.nat-hazards-earth-syst-sci.net/17/1611/2017/Nat.Hazards Earth Syst.Sci., 17, 1611-1621, 2017 D. Huang et al.: Hazard assessment comparison of Tazhiping landslide shown below: • 37 46 E, 31 • 6 29 N), 68 km west Chengdu city and 20 km from the Dujiangyan

Table 1 .
Parameters of surface soil of Tazhiping landslide.

Table 3 .
Division table of the predicted hazards of Tazhiping landslide (unit: m 2 ).