# IJEGE-11_BS-Martinez-et-alii

*Italian Journal of Engineering Geology and Environment - Book www.ijege.uniroma1.it © 2011 Casa Editrice Università La Sapienza*

*DOI: 10.4408/IJEGE.2011-03.B-051*

**QUASI-THREE DIMENSIONAL TWO-PHASE DEBRIS FLOW MODEL**

**ACOUNTING FOR BOULDER TRANSPORT**

**K**

**ey**

**words***: debris flow, boulders accumulation, finite*

*element method, discrete element method, lagrangian*

formulation

formulation

**INTRODUCTION**

of poorly sorted sediments, rocks and fine material,

agitated and mixed with water, surge down slopes in

response to water flow and gravitational attraction. A

typical surge of debris flow has a steep front or “head”

with the densest slurry, the highest concentration of

boulders and the greatest depth. A progressively more

dilute and shallower tail follows this head.

tion and clearly divide previous debris flow research

into two distinct categories. The first, based upon the

pioneering work of J

model describes a single-phase material that remains

rigid unless stresses exceed a threshold value: the

plastic yield stress. Various rheological models have

been proposed, derived from experimental results or

from theoretical considerations, such as the Bingham

model (b

*et alii*, 1989) and the quadratic model pro-

**ABSTRACT**

a continuum non-Newtonian fluid phase composed

by water and fine sediments, and a non-continuum

phase for large particles such as boulders. Particles

are treated in a Lagrangian frame of reference us-

ing the 3D Discrete Element Method. The fluid

phase flow equations are solved by the RiverFLO-

2D computational model which is based on the 2D

depth-averaged shallow water approximation and

uses the Finite Element Method on a triangular

non-structured mesh. The model considers particle-

particle and wall-particle collisions, and considers

that particles are immersed in a fluid and subject to

gravity, friction and drag forces. Bingham and Cross

rheological models are used for the continuum phase

providing very stable results, even in the range of

very low shear rates. Results show that the Bing-

ham formulation proves better able to simulate the

stopping of the fluid when applied shear stresses are

low. Results from numerical simulations comparison

with analytical solutions and data from flume-exper-

iments, show that the model is capable of reasonable

approximating the motion of large particles moving

in the fluid flow. An application to simulate a debris

flow event that occurred in Venezuela in 1999 shows

that the model can model the main boulder accumu-

lation reported for the alluvial fan in that event.

*C.E. MARTINEZ , F. MIRALLES-wILHELM & R. GARCIA-MARTINEZ*

*5th International Conference on Debris-Flow Hazards Mitigation: Mechanics, Prediction and Assessment Padua, Italy - 14-17 June 2011*

mentum equations in Cartesian coordinates can be

written as follows:

evation, and and are the depth averaged velocities in

x and y directions respectively; g is the gravitational

acceleration and is fluid density. FD represents the

fluid-solid interaction force exerted on the fluid by

particles through the fluid drag force, this force is

evaluated as:

number of particles in the cell.

*fx*

*fy*

model the slurry. Assuming a Bingham rheologi-

cal model and Manning’s formula, as proposed by

o’b

such as t

fractions and mass changes respectively.

esses of debris flows are still challenging in many

aspects. In particular, there are many factors related

to the movement and interaction of individual boul-

ders and coarse sediments that have not been fully

addressed in previous works. a

*et alii*(2003)

simulate the motion of solid particles in debris flows.

DEM is a numerical method to simulate dry granular

flows whereeach particle is traced individually in a

Lagrangian frame ofreference by solving Newton’s

equation of motion.

debris flows, considering a continuum fluid phase,

and large sediment particles, such as boulders, as a

non-continuum phase. Large particles are treated in

a Lagrangian frame of reference using DEM, and

the fluid phase composed by water and fine sedi-

ments is modeled with an Eulerian approach using

the depth-averaged Navier–Stokes equations in two

dimensions. Bingham and Cross rheological models

are used for the continuum phase. Particle’s equa-

tions of motion are fully threedimensional. Numeri-

cal simulations have been compared with analytical

solutions, with data from laboratory experiments and

with a real debris flow event.

**GOVERNING EQUATIONS**

*Fig. 1 - Schematic representation of debris flow with*

*large solid particles*

**QUASI-THREE DIMENSIONAL TWO-PHASE DEBRIS FLOW MODEL ACOUNTING FOR BOULDER TRANSPORT**

*Italian Journal of Engineering Geology and Environment - Book www.ijege.uniroma1.it © 2011 Casa Editrice Università La Sapienza*

**F**

*R*is the particle radius and

*ρ*

*p*

*is the particle*

*C*

*d*

*is the drag coefficient,*

**u**is the fluid veloc-

**v**is the particle

velocity vector.

Based on the simplified model that uses a spring-

dashpot-slider system to represent particle interac-

tions (a

*et alii*, 2003), the normal contact force

*N*is the Manning roughness coefficient.

*μ*and yield stress

*y*,

concentration

*Cv*, using the relationships proposed by

o’b

*α*

*1*

*, α*

*2*

*β*

*1*

*β*

*2*

with various sediment mixtures.

are expressed as

*m*

*eff*

tracked using Newton’s second law; considering grav-

ity, buoyancy, fluid drag and collision forces.

*C.E. MARTINEZ , F. MIRALLES-wILHELM & R. GARCIA-MARTINEZ*

*5th International Conference on Debris-Flow Hazards Mitigation: Mechanics, Prediction and Assessment Padua, Italy - 14-17 June 2011*

this direction,

*w*. Then, torque in

*x*and

*y*directions

are not detectable by the numerical model. Last com-

ponent is the torque in direction perpendicular to the

sloping surface, which is essentially negligible.

triangular elements. To solve the resulting system of

ordinary differential equation, the model applies a

four-step time stepping scheme and a selective lump-

ing method, as described by G

*et alii*.

from the particle governing equation, which is then

integrated to find velocity and displacement of each

particle.

**RESULTS**

*ANALYTICAL SOLUTIONS*

of mud flows. This modeling would account for the

representation of the fluid phase of the debris flow.

The numerical model was run using RiverFLO-2D

software, a finite modeling system for detailed analy-

sis of river hydrodynamics, sediment transport and

bed evolution (G

*et alii*, 2006). In

were implemented, the first, including Bingham the-

ory and Manning’s formula, as proposed by o’b

in m

*et alii*(2007).

This implementation provided enough data for verifi-

cation and testing of the new rheological formulations

proposed (m

*et alii*, 2007).

*LABORATORY EXPERIMENTS*

mixtures for the continuum phase and spherical mar-

bles for the discrete phase. The experiments were

performed in a 1.9 m long, 0.19 m wide, Plexiglas

walled flume, with adjustable slope. The downstream

*k*

*N*

*δ*

*N*

*i*and

*j*. The

maximum overlap is dependent on the stiffness

*k*

*N*

*.*

requiring stiffness of the order 105-107 N/m.

**F**

*ν*

*N*

*C*

*N*

*C*

*N*

*b*, defined as the ra-

tio of the normal component of the relative velocities

before and after collision.

*F**TC*

frictional limit, at which point the particles begin to

slide over each other. Prior to sliding, the tangential

contact force is calculated using a linear spring rela-

tionship,

*k*

*T*

*is the tangential stiffness coefficient, and*

*δ*

*T*

*i*and

*j*since their initial contact.

When

*k*

*t*

*δ*

*T*

*μ*

*f*

*F*

*NC*

fined as

*μ*

*f*

**F**

*TD*

curs, damping is accounted for from friction. Particle

rotation is not considered, since it its impact on boul-

der transport is assumed to be of much less impor-

tance than the friction and drag forces.

**QUASI-THREE DIMENSIONAL TWO-PHASE DEBRIS FLOW MODEL ACOUNTING FOR BOULDER TRANSPORT**

*Italian Journal of Engineering Geology and Environment - Book www.ijege.uniroma1.it © 2011 Casa Editrice Università La Sapienza*

*EXPERIMENT 1*

slope was set to 4° and the initial volume released was

6.3 L. For

*t*= 3 s the wave practically stopped flowing

as shown in Fig. 2. The propagation of the wave was

recorded for different times t to construct the spread-

ing diagram shown in Fig. 3.

lation. Numerically, the condition of stopping the fluid

is not easy to achieve; however, it is possible to appre-

ciate how the maximum velocity in the fluid decreases

with time and it becomes very close to zero about the

time the fluid must stop. This fact shows that velocity

criteria could be used numerically to stop the fluid.

mesh sizes. Figures 4 and 5 show how mesh refine-

ment contributes to improve substantially the

*at alii*(2009) was im-

tion, and then there is a well defined interface between

dry and wet elements.

platform, 0.75 m long and 0.95 m wide. A dam-break

type of flow was initiated by an abrupt removal of a

gate, releasing mixtures from a 0.40 m long reservoir

situated on the upstream part of the flume. Water-clay

mixtures were used in all the experiments, with vol-

ume sediment concentration 23.5 % and 26.5 %. For

preparation of the mixtures, kaolinite clay with spe-

cific unit weight of 2.77 was used. Fluid density was

measured in the laboratory and rheological parameters

*μ*and

*τ*

*y*

*were determined using equations (7) and (8)*

*α*

*1*

*β*

*2*

*α*

*β*

*2*

*Tab. 1 - Rheological properties of experimental fluids*

*Fig. 2 - Experiment 1, fluid stops flowing over the slop-*

*ing channel*

*Fig. 3 - Experiment 1, spreading relation*

*Fig. 4 - Final free-surface longitudinal profile and*

*U*

*max*

*- mesh size 0.03 m*

*Fig. 5 - Final free-surface longitudinal profile and*

*Umax - mesh size 0.01 m*

*C.E. MARTINEZ , F. MIRALLES-wILHELM & R. GARCIA-MARTINEZ*

*5th International Conference on Debris-Flow Hazards Mitigation: Mechanics, Prediction and Assessment Padua, Italy - 14-17 June 2011*

minimum depth parameter that determines the distinc-

tion between dry and wet elements. Best results were

found with a minimum depth parameter equal to 0.01

times the average fluid depth as shown in Fig. 5.

*EXPERIMENT 2*

was set to 9.54° and the initial volume released was

6.4 L. The objective of this test was to study the

spreading of the fluid in the fan and to study particle

movement into the fluid.

the gate. By the time the fluid was released, the piece

of wood was quickly removed, so the particles could

start their movement along the channel with the fluid

at

*t*= 0.5 s. Blue particles represent those particles

placed initially in the first row, orange particles are

those placed in the second row. In Figure 6(b) can

be noticed how particles in the center tend to move

forward to reach the front of the wave, particles in

the second row displace particles in the first row to

the sides and these are then left behind because of the

fluid velocity gradient. By the time the flow reaches

the fan, particles move to the sides of the flow as it is

shown in Figure 7.

measured at the lab) with the numerical results. In the

numerical solution can be seen that there exist some

delay on the particles positioned close to the walls;

this is due to the velocity boundary condition at the

walls. In practical applications of the model, it is

necessary to allow slip at the walls, since the no slip

condition formulated in finite elements becomes very

restrictive. However; it is not possible a total slip con-

dition, since for this case no velocity profile would be

created across the channel. In this simulation 90% of

slip at the wall was considered.

*EXPERIMENT 3*

flume bottom slope was increased to 10.7° and the

initial volume released was 11.1 L. The objective of

this test was to study the spreading of the fluid and

study particle movement into a mixture with higher

clay concentration.

*Fig. 6. - Experiment 2, numerical simulation (a) t = 0.5*

*s, (b) t = 1.6 s*

*Fig. 7 - Experiment 2: Final position of particles (a)*

*experimental data(b) numerical solution*

*Fig. 8 - Experiment 3, spreading relation*

**QUASI-THREE DIMENSIONAL TWO-PHASE DEBRIS FLOW MODEL ACOUNTING FOR BOULDER TRANSPORT**

used to map the distribution and thickness of depos-

its and to draw contours of maximum boulder size, as

shown in Figure 11.

The element characteristic size was 12 m on average.

The topography data used to define the finite element

mesh was interpolated from the original cartographic

information prior to the event (G

12, with an average volume sediment concentra-

tion of

*Cν*=0.3. The Manning coefficient considered

compared with numerical results obtained using Bing-

ham rheological model and using Cross rheological

model. Both rheological formulations produce very

similar results, they are not totally capable of resem-

bling the spreading of the flow; however, they show

a final fluid extend, when velocities in the fluid be-

come very close to zero, very similar to the real one.

Bingham formulation shows to be more effective in

decreasing the velocities oalong the fluid to zero.

ous experiment. Figure 9 compares the final particle

positions obtained numerically against final observed

particle location. Note that some particles lag behind

close to the flume wall and that the general location

of the particles on the alluvial fan is very close to the

observed locations.

AL FAN DEBRIS FLOODING EVENT

steep slopes of Cerro El Avila, north of Caracas, Ven-

ezuela, and caused flooding and massive debris flows

in the channels of major drainages that severely dam-

aged coastal communities along the Caribbean Sea.

The largest fan on this area is that of San Julián River

at Caraballeda, shown in Figure 10. This fan was one

of the most heavily damaged areas in the event. The

thickness of sediment deposition, maximum size of

transported boulders, and size of inundated area were

all notably larger in this drainage in comparison to the

other close watersheds.

*et alii*, 2001), measuring slope, de-

*Fig. 9 - Experiment 3: Final position of particles, (a)*

*experimental data(b) numerical solution*

*Fig. 10 - Caraballeda alluvial fan, Venezuela*

*Fig. 11 - Contours of maximum transported boulder size*

*on the Caraballeda alluvial fan, Venezuela.*

*From USGS, 2002*

*C.E. MARTINEZ , F. MIRALLES-wILHELM & R. GARCIA-MARTINEZ*

to take into account terrain irregularities. The same

value was used by G

ships (7) and (8) were selected for the calculation of

fluid rheological properties, using the parameters for

water-clay mixtures. As a result of the volume sedi-

ment concentration,

*Cν*= 0.3,

*ρ*= 1531 kg/m

*m*=

*t*

*y*

in the event. The boulders were placed into the fan

during the first three hours of simulation, at a rate of

50 particles every 6 min. This amount of particles was

selected to ensure a manageable computational time.

Density for the boulders is

*r*

*p*

mostly found in the area by the USGS.

Figure 12. Comparing this region with the post-event

aerial view shown in the background, it can be noted

that the model acceptably reproduces the extent of the

area affected by the debris flow.

apex, where the discharge of the river is simulated.

Velocities decrease at the urban areas, ranging from

1.0 to 6 m/s at the time of the hydrograph maximum

value. Higher velocities are developed in the avulsion

zone, at the center of the fan,, reaching 10 m/s. The ve-

locities calculated by the model are in good agreement

with those estimated by USGS, which ranged from 1.3

to 13.6 m/s.

the right side of the fan, while smaller boulders take

the central path. According to the USGS report, the

slope at the center was 4 degrees, while the concrete

channel direction, was steeper, with a slope gradient

of 5.5 to 6 degrees, then larger boulders were trans-

*Fig. 12 - Inflow hydrograph for a 500 year-return pe-*

*riod, including solid concentration by volume*

*of 0.3*

*Fig. 13 - Flooded area at time t =2.2 h. Legend indicates*

*flow depth in m*

*Fig. 14 - Velocity field at t =2.2 h. Legend indicates ve-*

*locity in m/s*

**QUASI-THREE DIMENSIONAL TWO-PHASE DEBRIS FLOW MODEL ACOUNTING FOR BOULDER TRANSPORT**

**CONCLUSIONS**

to simulate debris flows, considering large particles,

such as boulders. The continuum non-Newtonian

phase is solved by RiverFLO-2D finite element model

in two horizontal directions and the particle transport

with the Discrete Element Method in 3D.

in a laboratory flume-fan, including the effect of parti-

cle-particle and wall-particle collisions.

shear rates. In the simulation of mud dam-break prob-

lems, Bingham formulation was better able to simu-

late the stopping stage of the fluid; however, the Cross

formulation proved more accurate for early stages of

the solution.

the capability of the model to simulate large scale real

events. Results show that the model reasonably ap-

proximates the flood extent affected by the debris flow

and the observed boulder accumulation areas, including

distribution boulders sizes. Future work includes com-

parison with field events using larger number of boul-

ders and enhancement of the detecting particle contact

algorithm in order to improve computational time.

ameter and slope steepness reflect USGS observations

that for the larger transported and deposited boulders

there was a proportional relationship between boulder

size and slope steepness.

central direction alignment, some of them reached the

shoreline or entered into the sea. Larger boulders were

deposited in the avulsion zone or took right direction

to the concrete channel. None of these large boulders

reached the shoreline. In Figure 16, it can be seen that

the model predicts reasonably boulder locations as

compared with the field data given in Figure 11.

dence that strongly supports transport by debris flow.

At other sites, the largest boulders were observed iso-

lated along the concrete channel, fact that suggests

that these boulders moved sliding along the bottom of

the channel in a dilute fluid until deposition occurred

(USGS Report 01-0144). There is no indication of big

boulders close to the shoreline at this site of the fan.

in order from bigger to smaller as it proceeds down-

stream on alluvial fans. This process was better ob-

served along the central direction and it was also rep-

licated in the numerical simulation.

*Fig. 15 - Boulder positions at time t =1.8 h, Legend in-*

*dicates particle diameter in m*

*Fig. 16 - Boulder positions at time t = 6.0 h, Legend in-*

*dicates particle diameter in m*

*C.E. MARTINEZ , F. MIRALLES-wILHELM & R. GARCIA-MARTINEZ*

**REFERENCES**

*The potential of the Discrete Element Method to simulate debris flow.*

**1**: 435-445.

*Experiments on a gravity-free dispersion of large solid spheres in a Newtonian fluid under shear.*Pro-

**225**: 49-63.

*An introduction to rheology.*Elsevier. Amsterdam.

*aint a plastic material and not a viscous liquid; the measurement of its mobility and yield*

*value.*Proceedings of American Society of Testing Materials,

**19**: 640-664.

*n explicit two-dimensional finite element model to*

*simulate short and long term bed evolution in alluvial rivers.*J. of Hyd. Res.,

**44**(6): 755-766.

*ud Flow Hazard Maps for Vargas State. Final Report for the Avila Project.*Fluid Mechanics

*Dam-break flood routing.*Chapter 5 of Dam-Break Problems, Solu-

*The physics of debris flows*. Rev. of Geophysics,

**35**: 245-296.

*A 2D finite element debris flow model based on the cross*

*rheology formulation.*Fourth International Conference on Debris-Flow Hazards Mitigation: Mechanics, Prediction and

Assessment. Chendu, China.

*Physical properties and mechanics of hyperconcentrated sediment flows.*ASCE Specialty

*Laboratory analysis of mudflows properties.*J. of Hyd. Eng.,

**114**(8): 877-887.

*Debris Flows.*Balkema, Rotterdam.

*Debris-flow and flooding hazards associated*

*with the December 1999 storm in coastal Venezuela and strategies for mitigation*. U.S. Geological Survey, Open File Report

01-0144.