Contribution of seismic tomography for taking soil variability into

account in landslide stability numerical modelling

J. Gance (1), S. Bernardie (1), R. Vandromme (1) & G. Grandjean (1)
(1) BRGM – Natural Hazard and C02 Security Storage Department. E-mail: g.grandjean@brgm.fr, Tel. +33-
(0)238643524

ABSTRACT: The Ballandaz landslide (Le Planay, France) has been studied for years. Geotechnical charac-
terization, geophysical surveys and geomechanical modelling have been carried out, but contributions these
different approaches can provide each others were not deeply analysed. The idea of this work is to use geo-
physical and borehole data obtained during previous campaigns, for constraining a finite element modelling
taking into account soil variability. Instead of constructing the slope model by geological interpretations of
geophysical sections and boreholes, as it is generally done, we use a velocity distribution derived from seismic
tomography computations to fill the mechanical mesh with different mechanical properties. Simulations are
carried out with the GEFDYN code. A parametric study is also performed for determining the optimal number
of materials, the optimal element size, the velocity classification algorithm, etc. Finally, our results are com-
pared to a classical 3-layers model for discussion.

The work is divided into two parts. The first one

presents how the numerical mesh modeling the slope
1 INTRODUCTION is generated automatically from the P-wave velocity
(Vp) tomogram (Grandjean & Sage, 2004). For this
Nowadays, geomechanical finite-element codes aspect, the numerical mesh is composed of a given
running on modern computers allow complex calcu- number of materials, each material being character-
lations and produce coherent results, particularly ized by its mechanical parameters. Each parameter is
when homogeneous media are considered. Taking estimated from its petrophysic relationship with Vp.
into account the spatial heterogeneity of geological The second part presents a parametric study where
materials is becoming a new challenge in numerical different algorithm options and parameters entering
modeling, since simulating natural processes con- in the computations are tested with the GEFDYN
form as much as possible to reality is not straight- code (Aubry et al., 1986). The main goal of this
forward. Some scientific publications already deal study being dedicated to measure the reliability of
with this issue, but generally focalized on stochastic our approach in the context of heterogeneous slopes,
analyses of the problem. These methods give rela- we finally compare the final section of the failure
tively good results for slope stability models (Cho, criterion to some field observations.
2007; Griffiths & Fenton, 2004; Vanmarcke 1980;
Meo et al., 2008) but remains in complex cases
strongly time consuming (Rohmer 2009). Montgom- 2 AUTOMATIC CREATION OF THE MODEL
ery (1994) showed that the simple introduction of
topography in the model allows increasing the poten-
2.1 Creation of the mesh
tiality for shallow landslides occurrence.
A transformation program is used for creating the
We propose here a model that takes into account mechanical mesh. It uses directly the Vp data grid -
the spatial variability of the different constitutive the seismic tomogram - to build a numerical mesh
materials of a slope, and the changes of their me- taking into account the topography, each element be-
chanical properties, from seismic tomographic imag- ing characterized by a Vp value. Then, a k-means
es. To measure the pertinence of such approach, we classification algorithm is used to attribute a material
compare it with a simple 3-layers model that is gen- number to elements, the total number of material be-
erally used for standard numerical studies. ing arbitrary fixed.
 Each material is attributed to a group by the same 2 1 1 2
algorithm so that the materials belonging to a group E vp (2)
1
have some mechanical properties in common (Figure
1). The number of group is also arbitrary fixed. - The permeability k is calculated with the relation
proposed by Berg (1970):
6 5.1
2.2 Determination of the hydro-mechanical k 5,1 10 d 2 exp( 1.385 ) (3)
parameters Where k is the permeability in darcy; is porosity in
percent, d is median diameter of grains in mm, φ is
After the mesh has been created, the program at- the standard deviation and equals to P90 – P10.
tributes to each material the values of 9 hydro-
mechanical parameters necessary to define the We choose, after calibration, φ = 0 and d = 0.12
Drucker-Prager constitutive law in GEFDYN. mm. The coefficient of earth pressure at rest is cal-
culated from the friction angle:
K0 1 sin (4)
The friction angle will be taken under 35.3° for the
good accuracy of the results (Desrues 2002).

Figure 2 (top) shows the mesh constituted by the

three different groups, each group containing a cer-
tain number of materials (non-visible here), and each
material being defined by 9 hydro-mechanical prop-
erties. Figure 2 (bottom) shows the classical 3-layers
model to compare.

Figure 1. Creation of a mesh from the Vp data grid.

Four of these parameters are identical to all mate-

rials composing a group. These properties are the
solid density, the cohesion, the friction angle and the
Poisson ratio. If the model is constituted of 3 groups,
three different sets of parameters need to be speci-
fied.
The other parameters are supposed to vary spa-
tially within a specific material; they are calculated
for each of them from the Vp values and the constant
parameters of the group. Although the petrophysical
relationships used for this estimation are rather valid
for rocks, they have been selected because of their
simplicity and their coherent results in most of geo-
logical materials. By this mean, continuous varia-
tions of realistic parameters can be modelled, even if
the reality is anyway more complex to describe. The
parameters derived from Vp are the following: Figure 2. Model created automatically from the Vp data
grid (top). 3-layers model used for the comparison.
- The porosity n: it is calculated with the relation
proposed by Castagna (1985):
2.3 Boundary conditions
Vp 5810 9420 n 2210 Vcl (1)
The calculations were carried out with the Druck-
where n refers to the porosity in percent, Vcl, the clay er-Prager constitutive law defined in the GEFDYN
volume in percent. Vcl is taken equal to zero to keep finite-elements code, using the plane strain approxi-
a good range of porosity in the model. mation. Vertical displacements and horizontal ones
were respectively fixed at the lower edge of the
- The Young modulus E: it is calculated with a rela- mesh and at the left one. The right edge was as-
tion restricted to isotropic homogeneous medium: sumed to be free of displacements for the strength
not to be constrained (Figure 3).
 3.2 Failure criterion
In the case of hydro-mechanical calculation, the In order to visualize the potential failure of the
pore pressure is imposed at each node. We model the slope, we need to define a failure criterion. The most
presence of water using a linear water table. All the
commonly used is the Mohr-Coulomb criterion, but
pressures imposed are hydrostatic ones (Figure 4). it can be out of range in the case of low normal
The simulation then refers to a non-coupled hydro- stresses. We thus decided to use a parabolic criterion
mechanical calculation. (Lade, 2010) that corrects this problem:
b
'
a Pa (5)
Pa
With τ is the tangential stress, σ’ the effective nor-
mal stress, Pa the atmospheric pressure, a and b two
dimensionless coefficients.

This criterion is simply based on the assumption

Figure 3. Conditions on displacements. that it fits the Mohr-Coulomb criterion for maximum
stresses (Figure 3). The a et b coefficients are chosen
so that the two curves are superimposed for normal
stresses higher than 50 kPa. The failure is studied
only on horizontal facets. For avoiding to mix two
criteria (parabolic and Drucker-Prager criteria), the
computations are realized in the elastic domain. The
following equation represents the distance to the cri-
terion:

criterion if criterion
(6)
0 if criterion
Figure 4. Imposed pore pressure.

3 RESULTS OF CALCULATIONS

3.1 Preliminary tests on the size of the elements and

the number of materials
A first simple study shows that the number of ma-
terials doesn’t change very much calculation times, Figure 3. The parabolic and Mohr-Coulomb criteria repre-
contrary to the size of the elements, as clearly shown sented in the Mohr plane.
in Figure 2.
3.3 Comparing the 3-layers model and the
automatic one.
We compare here the impact of the two following
models on the simulation results:
- The model 1: a 3-layers model, created by hand,
that contains 3 materials;
- The model 2: generated automatically from Vp
values that contains 3 materials too.
For each material, the properties are the same; they
were calculated from Vp values and equations 1 to 4
(Table 1).
Figure 2. CPU time in function of the size of the elements
for different numbers of materials. Table 1. Parameter fixed for each group.
______________________________________________
Property Material 1 Material 2 Material 3
______________________________________________
On the other hand, the size of the elements plays a Poisson ratio 0.3 0.2 0.17
role in the quality of the results. This value was ad- Solid density 1800 2200 , 2400 (Kg/m3)
justed to 1.3 m, for avoiding the numerical artifacts Friction angle 35° 29° 29°
due to the mesh roughness without increasing the Cohesion 10 300 high (kPa)
computation times too much.
 Let’s analyze stress variations along a vertical 4 PARAMETRIC STUDY
section for the two models (Figure 4). We can notice 4.1 The number of materials
that results are quite similar in the homogeneous In this part, we study the number of materials that
zones for both of them, but some differences are vis-
can be used in the model. By setting a large number
ible at the interfaces. Indeed, the model 2 presents of materials, we expect to take into account the
some high shear stress picks at the interface of het- smooth variations of hydro-mechanical properties, as
erogeneities that are not present in the model 1. shown on the Vp tomogram, avoiding consequently
Moreover, the higher stress values observed at the sharp variations between two consecutive materials.
interfaces of materials are deeper for the model 2
than for the model 1. Six different automatically generated models have
been tested with a number of materials comprised
between 3 and 100. The results show that the places
of the potential failure are affected by this number.
We can see for example that the failure zones are
deeper (1 to 5 meters) with a large number of mate-
rials as shown on Figure 6. We will see later that this
representation of the failure plane fits quite well with
the field data. The amplitude of the strength is also
higher in the 3-layers model than in the others. This
could be due to strength artifacts located near the in-
terfaces separating too contrasted materials as it is
Figure 4. Shear stress on the cross section X=110m. the case for the 3-layers model. Moreover, the com-
parison between the Figure 8 and Figure 6 shows
As a first conclusion, we can highlight that the that a high number of material smoothes the result
simple consideration of the materials spatial variabil- and then give more coherent stresses. Indeed, in the
ity - that represents the only difference between the classical 3-layers model, the failure appears as guid-
two models - can strongly modify the slope behav- ed by the layered structure, leading to artifacts in
ior. Instead of giving a large zone of failure, the cal- some places, whereas in the 100 materials model the
culation with our automatic model shows three lo- failure zones are larger and more regular.
calized zones situated in the second material at the
places where some heterogeneities are observed
(Figure 8).

Figure 6. Failure criterion on a 100 material model.

Evidently, a high number of materials seems more

appropriate to model this kind of slope. For the rest
of the studies, we will use 100 materials.

4.2 Variation of the cohesion inside a group.

Because we have observed that the continuity of
parameters within the model was a good assumption
in our case, we can try to smooth the cohesion line-
arly from the first material of a group to the first ma-
terial of the following one. Because we can’t set our
Figure 8: Localized failures in the case of model 1 (top) and cohesion values to huge values in order to stay in the
model 2 (bottom) elastic hypothesis, we selected one of the two crite-
ria, respectively when adapted to the first (weak co-
hesion) or other (strong cohesion) groups. This op-
tion is supposed to be valid since the Drucker-Prager
criterion is reached only in the first group and the
failure is located in the second one.
 minimize those artifacts, we built a model with a
As shown on Figure 10, the strengths are now shorter size of elements and with a higher number of
smoothed at the interfaces and the contrast of stress- materials. These new conditions made the resulting
es between two materials is lower. In addition, the images smoother, with a good integration of the spa-
convergence of computations is better. tial variability of properties and a good convergence
of the computations.

5 COMPARISON WITH THE FIELD DATA

Two boreholes have been drilled close to the stud-

ied profile. The first one (SC1) give a surface of
failure at 22 meters of depth. The second one (I4) is
not located exactly on the profile, but we expected a
Figure 7. Shear strength on the cross section X=191m for surface of failure between 10 and 15 meters of depth
fixed cohesion (solid line) and linear variation of cohesion (dot just under the road.
line)
We performed a calculation on a 100 materials
4.3 Test on the classification algorithms model, created with all the properties described in
the previous paragraph. Then, we calculated the par-
In order to understand the impact of the classifica-
abolic failure criterion and from this result, we iden-
tion of materials on the results, we have also tested 3
tified a likely surface of failure on the model (Figure
different classification algorithms. The k-means al-
8). This result was compared to the one obtained
gorithm was compared to two algorithms that (i)
from the 3-layers model.
separates the Vp in classes of the same size and (ii)
creates classes that contain the same number of ele-
ments.

These algorithms have been chosen to be very dif-

ferent in order to highlight their impact on the re-
sults. After the tests, we can note that this impact is
quite negligible. The shear strength is slightly differ-
ent only in the downhill part of the slope of about
6%. This difference seems not to be due to the
change of the geometry since we uses a 100 materi-
als model. It is more related to the change of the Vp
values attributed to each material, and then to the Figure 8. Comparison between two likely surfaces of fail-
change of the mechanical parameters that are calcu- ures.
lated from this Vp.
We can see that the two identified surfaces are
FInally, the effects of the classification algorithm quite different. The surface found with our model is
on the simulation are minor in the case of a 100 ma- deeper than with the 3-layer one. It seems that the
terials model. This test also shows that a variation of surface of failure of our model fits better with the
mechanical properties in the materials has not al- field data than the other one. Moreover, the simple
ways an effect on the resulting strength. We will thus model doesn’t give information under I4, whereas
keep the k-means algorithm for the rest of the study. the other one permits to draw a surface around 13m.

4.4 Conclusions on the automatically created model

6 HETEROGENEITIES IN THE MODEL
We saw that the automatic constructed model al-
lows considering all information present in the seis-
To identify the surface of failure as we did in the
mic tomogram since it takes into account the spatial
previous paragraph, we used the parabolic failure
variability of the medium heterogeneities after they
criterion, post processed from the strength calculated
have been translated into mechanical properties.
by the finite elements code GEFDYN. We stated that
the automatic model take into account the heteroge-
The principal drawback of the method concerns
neities, since it is able to consider inclusions of a
the artifacts due to crenellated boundaries separating
strong material in a softer one. The parameters of
the materials: these artifacts impact the simulations,
this heterogeneity are then stronger that those of the
particularly when the size of the elements is large
soils all around. Nevertheless, the finite elements
and when the number of material is low. In order to
code handles the inclusions and the bulk material as
one. The solution in displacements has to be contin- Some improvements on the method are still possible.
uous, and then, doesn’t admit the slide of a material Indeed, we have seen that the role of heterogeneities
on another. That case is however probable in reality. in the failure were not fully taken into account, and
that it was not recommended to model two very con-
This problem can be solved in GEFDYN by using trasted materials, because of the smoothing tech-
a mechanical interface element that allows the slide nique proposed to stabilize the convergence. These
of a volume on another one, but the difficulty is to aspects should be addressed in the future steps of our
automate this process, and in particular to detect work.
such a situation all over the model: not all heteroge-
neities are supposed to slide on the bulk material sit- Finally it would have been interesting to test some
uated below. dynamic conditions. For example, the sudden raising
of the water table level could have modified tempo-
To finalize our interpretation, we drew a new sur- rarily the stresses on the model because of the spatial
face of failure considering the slide of the heteroge- variability of the permeability.
neity was possible, and make it go just under the
downhill block (Figure 9).
8 BIBLIOGRAPHIE

Aubry, D., D. Chouvet, A. Modaressi, and H.

Modaressi. GEFDYN - Notice Scientifique. Ecole
Centrale PARIS, 1986.
Berg, 1970
Castagna, 1985
Cho, Sung Eun. "Effects of spatial variability of soil
propertiies on slope stability." Engineering
Figure 9. Surface of failure drawn from the result of a calcu- Geology, 2007.
lation on our model. The depths of the surface is written at SC1 Desrues, Jacques. "Limitations du choix de l'angle
and I4 de frottement pour le critère de plasticité de
We then compared our results to the field data and Drucker-Prager." Revue Française de
concluded that the depths of the sliding plane in SC1 Géotechnique, 2002: 853-862.
and I4 was well correlated to the simulated maxi- Grandjean, G., and S. Sage. "JaTS: a fully portable
mum of failure. seismic tomography software based on Fresnel
wavepaths and a probabilistic reconstruction
approach." Computers and Geosciences 30, 2004:
7 CONCLUSIONS 925–935.
Griffiths, D.V., and Gordon A. Fenton.
We tested and developed a method for building a "Probabilistic slope stability analysis by finite
finite elements model from the P-wave velocities elements." Journal of the Geotechnical
value obtained by seismic tomography. This ap- Engineering Division, 2004, ASCE 130 (5), 507-
proach preserve in the modeling the sliding domains 518.
of the slope compared to the common used method Meo, M., U. Tammaro, and P. Capuano. "Influence
that orient them along contrasts of mechanical prop- of topography on ground deformation at Mt.
erties. We have demonstrated that the simple contri- Vesuvius (Italy) by finite element modelling."
bution of the geometry of the soil provide rather ac- Non-linear Mechanics, 2008.
curate results. Montgomery, David R. "A physically based model
for the topographic control on shallow
The first drawback of this method is that the cren- landsliding." Water Ressources Research, 30, 4,
ellated interface between materials creating artifacts 1994 : 1153-1171.
on strength can be limited if a high number of mate- Rohmer, Jeremy. Développement d'une
rials is used, and if the contrast of properties be- méthodologie de propagation de champ aléatoires
tween two consecutive materials is low. The study dans les modèles éléments finis. Journées
also highlights the difficulty to model properly a Vulnerisc 2009, Orléans, France, 2009.
medium constituted of several materials with quite Lade, V. P. “The mechanics of surficial failure in
different properties. The proposed method solves soil slopes.” Engineering Geology, 2010.
this issue by smoothing linearly properties belonging Vanmarcke, E.H. "Probabilistic stability analysis of
to two contiguous materials. This approach allows a earth slopes." Engineering Geology, 1980.
good convergence for high contrasted models, with
good accuracy of resulting strength.