# Inverse Parameter Identification Technique Using PSO Algorithm Applied to Geotechnical

Authors: Joerg Meier; Winfried Schaedler; Lisa Borgatti; Alessandro Corsini; Tom Schanz
Journal: Journal of Artificial Evolution and Applications
 (2008)

---

## Abstract

This paper presents a concept for the application of particle swarm optimization in geotechnical engineering.
For the calculation of deformations in soil or rock,
numerical simulations based on continuum methods are widely used.
The material behavior is modeled using constitutive relations that require sets of material parameters to be specified.
We present an inverse parameter identification technique, based on statistical analyses and a particle swarm optimization algorithm,
to be used in the calibration process of geomechanical models.
Its application is demonstrated with typical examples from the fields of soil mechanics and engineering geology.
The results for two different laboratory tests and a natural slope clearly show that particle swarms are an efficient and fast tool for finding improved parameter sets to represent the measured reference data.


---

## Body

## 1. Introduction

When compared to most other engineering tasks, geotechnical problems are often characterized by the following peculiarities. The materials involved are geological materials, that is, soils and rocks, which are inhomogeneous and consist of various phases in
different states of aggregation. Initial and boundary conditions tend to be
complex and heterogeneous. Furthermore, in real geotechnical field problems,
the exact geometry is usually not known, with the available geometry
information being limited to topographical surface data and punctual outcrops
or soundings. For this reason, in geotechnics, there is always a need for a
high level of simplification and abstraction. Frequently, continuum methods are
used to calculate deformations in soil or rock, and the material behavior is
simulated by means of constitutive models, which require a certain set of
material parameters.Normally, in geotechnical engineering, the values of these parameters are set based on the results of
laboratory experiments, literature data, or even just experience values are
used. The results of the calculation, a forward calculation, are then compared to measurement data obtained in the laboratory or in the field. Provided that the simplifi- cations made and the constitutive model chosen are appropriate and provided that the performed
calculation gives plausible results, parameter values are then varied by trial
and error in order to reach an improved fit of the calculation results to the
measured data, the reference data.Though this is done based on the
experience of the geoscientist, the procedure remains to a certain extent
arbitrary or at least subjective.In recent years, due to the availability of sufficiently fast computer hardware, there has been a growing interest in the application of inverse parameter identification
strategies and optimization algorithms to geotechnical modeling in order to
make this procedure automated [1–5] and thus more traceable and
objective. Furthermore, this approach provides statistical information, which
can be used to quantify the calibration quality of the developed geotechnical
model.Applications of optimization procedures in geotechnics were described by many authors, for
example, in the calibration process of geotechnical models [1, 2], or to identify hydraulic parameters from field drainage tests [6]. Already in 1996,
Ledesma et al. [7] and Gens et al. [8] applied gradient methods to a synthetic and a real example of a tunnel drift simulation. Also during the
excavation of a cavern in the Spanish Pyrenees the above mentioned group of
authors was applying gradient methods for the identification of geotechnical
parameters [9]. Malecot et al. [10] used inverse parameter identification
techniques for analyzing pressure-meter tests and finite-element simulations of
excavation problems. For the identification of soil parameters, also genetic algorithms
were studied [11, 12]. Feng et al. [13] used an inverse technique for the
determination of the parameters of viscoelastic constitutive models for rocks,
based on genetic programming and a particle swarm optimization algorithm. In
the field of geoenvironmental engineering, Finsterle [14] examined the potential use of standard optimization algorithms for the solution of aquifer remediation
problems in three-phase and three-component flow and transport simulations of
contamination plumes. As a different aspect of parameter identification, Cui
and Sheng [15] determined the minimum parametric distance to the limit state of
a strip foundation by optimizing a reliability index. In 2006, J. Meier and T.
Schanz in [5] applied particle swarm optimization
techniques to geotechnical field projects and laboratory tests, namely, a multistage
excavation and the desaturation of a sand column.All cited references agree on the fact that back-calculation
of model parameters by means of optimization routines is possible in the field
of geotechnics, if an appropriate forward calculation depending on adequately
realistic model assumptions is provided, for example, Calvello and Finno [1, 2].
In this context, particle swarms represent a powerful tool for finding
parameter sets that best represent the reference data, with acceptable
calculation effort and time consumption.


## 2. Working Scheme of the Adopted Parameter Identification Strategy

The starting point of the parameter
identification strategy presented in this study is given by an ordinary
geotechnical modeling-task, the so-called forward calculation. This forward
problem consists of a specified geometry with given initial and boundary
conditions and a material model, which requires a set of material parameters to
be determined. It is generally also possible to identify geometrical parameters
[4], but this issue will not be discussed in this article. For the first run of the forward calculation, the user presets the unknown parameters, for example,
entering estimated values, or the preset of the parameter vector is done by a
random generator within values margins specified by the user. The relevant
results of the forward calculation are then read out and their deviation from a
set of reference data, usually measured data, is determined by means of an
objective function. This procedure is repeated many times, while at any one
time, an optimization algorithm, based on the parameter combinations and the
values of the objective function during the previous forward calculations,
identifies an improved parameter set to be used in the next forward
calculation.This sequence of cycles, illustrated in Figure1, is interrupted when one of the following stop criteria
is fulfilled:Figure 1
Flowchart of the adopted iterative
procedure.(i)
a maximum
number of runs or maximum calculation time is reached;(ii)
the
deviation from the reference dataset, described by the value of the objective
function, falls below a specified limit;(iii)
the
deviation could not be lowered during a certain number of cycles.Hence we use a direct approach as
described by Cividini et al. [16] to solve a back-analysis problem. In an iterative procedure, the trial values of the unknown parameters are corrected
by minimizing an error function. It is therefore not necessary to formulate the
inverse problem itself, the desired solution is obtained by combining the
results of numerous forward calculations with an optimization routine.To quantify the deviation between the
reference data and the modeling results, we chose the frequently used and
relatively simple method of least squares. In this method, the objective
functionf(x) for more than one reference dataset is defined as
(1)f(x)=∑g[wgfg′(x)]
with
(2)fg′(x)=1m∑h=1mwh(yh,gcalc(x)−yh,gmeas)2.
In (1) and (2), x denotes the parameter
vector to be estimated and wg are positive weighting factors
associated correspondingly with the error measure fg′(x).
Via the weights   wg,
the different series  g can be
scaled to the same value range and different precisions can be merged, for
example, a series of measuring data possessing a higher precision is included
with a higher weighting factor compared to more uncertain data. The particular
numbers for the weights have to be given manually respecting the engineer’s
experience and they have to be specified depending on the optimization problem.
The weighting factors wh are used to provide a possibility for
considering different precisions and measurement errors within one and the same
data series. The dimensions of the weighting factors can be taken in the way to
obtain a dimensionless objective function quantity. For minimizing the
objective function, we use a particle swarm optimization algorithm described by
Eberhart and Kennedy [17, 18].A computer program developed by the first author of this paper implements this algorithm
and disposes of interfaces to several commercial finite-element packages used
in geotechnical engineering. A short pseudocode of the implemented PSO used in
this study is presented in Algorithm1.Algorithm 1:Pseudocode of the
used particle swarm optimizer.
“initialization
of swarm”create
particle listfor
each particleset
initial particle position and velocitynext
particle“processing
loop”do“get
global best particle position”determine
best position in past of all particles“determine
current particle positions”for
each particlecalculate
and set new velocity based on a corresponding equationcalculate
and set new position based on a corresponding equationnext
particle“parallelized
calculation of objective function values”for
each particlestart
forward calculation and calculate objective function valuenext
particle“join
all calculations”wait
for all particle threads“post-processing”for
each particleif
current objective function is less than own best in past:save
position as own bestend
ifnext
particleloop until stop criterion is met


## 3. Concept for the Application to Geotechnical Problems

If the parameter identification
strategy described above is to be applied to geotechnical tasks, the
geotechnical model of the forward problem usually has to be adapted for its use
in the optimization routine. The runtime of a single forward calculation has to
be minimized in order to allow for a high number of calls. The number of calls
needed depends furthermore on the number of parameters to be identified and, of
course, on the used optimization algorithm.While the reduction of calculation time demands simplification and abstraction, the model still should be sufficiently complex to reproduce the reference data with the required
accuracy. Furthermore, the number of required forward calculations can be
reduced by applying hypersurface approximation methods [4].As a next step, it is essential to select the parameters to be identified and to decide on the upper and lower
limits of their plausible values margins. The values of some of the parameters
might be fixed with the aid of previous knowledge. These specifications must be
done with care and require the experience of the geoscientist, since they
influence the obtained results. However, due to the application of a particle swarm
optimizer instead of, for example, a gradient method, no initial guess for the
parameter set is necessary because initial positions are generated randomly
within the parameter value margins.Due to the inhomogeneities of the
geological materials involved and the uncertainties related to the initial
conditions and the geometrical boundary conditions, geotechnical problems tend
to be underdetermined. In order to improve and ensure the efficiency of the
back-analysis, it is of significant importance to check if the set of
parameters to be identified may be reduced and if for the prescribed trusted
zone the optimization problem is well posed. For this purpose, a statistical
analysis is done based on the results of a Monte Carlo procedure including a sufficient parameter set. Figure2 shows the principle
scheme of the matrix plot used here to visualize the results of the Monte Carlo simulation. A standard mathematical tool for
examining multidimensional datasets is the scatter plot matrix (see [19]), whic is included in the
matrix plot presented in Figure 2, where each nondiagonal element shows the
scatter plots of the respective parameters. The matrix is symmetric. Matrix
element D-B, for example, may suggest that the involved
parameters B and D are not independent, but strongly correlated. The diagonal matrix
plot elements (A-A,…,D-D) show plots where the
value of the objective function is given over the parameter which is associated
with the corresponding column. These plots are called hereafter objective
function projections. If the problem is well posed, each of these plots of the
objective function projections has to present one firm extreme value as it is
the case in the diagram D-D. Otherwise, the respective parameter
could not be identified reliably. By filtering out data points that have
objective function values larger than a certain threshold level, the distribution
of the remaining points gives a rough idea of the size and shape of the extreme
value (solution) space. For further statistical analyses, the well-known linear
2D correlation coefficient can be calculated from the individual scatter plots.
Referring to [20] in the analysis of our case studies, we consider variables
with a correlation coefficient of less than 0.5 as “noncorrelated”.Figure 2
Statistical analyses via
matrix plot, principle scheme.


## 4. Applications

### 4.1. Description of the Studied Geological Material and the Adopted
Constitutive Model

In the following examples of applications of the presented parameter identification technique using PSO, we are modeling the mechanical behavior of a natural soil. It is of geotechnical interest because it favors the development of numerous landslides, namely,
rotational soil slips, earth slides, and earthflows. The studied material is clay that results from the weathering of structurally complex geological formations, the San Cassiano formation (Kassianer Schichten), and the La Valle formation (Wengener Schichten) of the Alpine
Trias. These rock formations are made up by interbedded strata of marls, tuffites,
claystones, limestones, dolomites, and sandstones. Like its source rocks, the
soil is characterized by a high clay and silt content. In the field, it
consists of a clayey matrix with coarser components, of diameters from
centimeters up to meters, floating in it without mutual support. Therefore, it
is not possible to sample the material as a whole representatively. As it has
been completely remolded by earthflow phenomena, no preferred orientation of
the components can be observed. For this study, only the fraction smaller than
2 mm was taken, as it is considered to determine the relevant
mechanical properties of the entire soil and it allows for the manufacturing of
reproducible samples and test specimens. Some properties of the studied
material are listed in Table1. Unless otherwise expressly stated, all tests and classifications were carried out according to the German standard DIN.Table 1
Properties of the studied soil.
Percentage
of clay particles (diameter<2 μm)26–34%Bulk
density (g/cm³)1.82Dry
density (g/cm³)1.29Density
of grains (g/cm³)2.76Porosityn53%Lime
content36%Loss
of ignition6.5%Liquid
LimitWL0.50Plastic
LimitWP0.27Water
content0.41From this description of the material, it becomes clear that its mechanical behavior is expected to be very complex and therefore it is only possible to model some important aspects of this behavior. Like many soils with a high clay and silt content, the studied
material is highly compressible and exhibits a significant amount of creep
deformations, thus its behavior is strongly time-dependent. As a constitutive
model, we chose the soft soil creep model, which was developed by Vermeer and
Neher [21] to account particularly for these phenomena. The soft soil creep model requires the following material parameters to be specified (see Table 2).Table 2
Parameters of the soft soil
creep model.
ParameterDescription(Unit)cEffective
cohesion(kPa)φEffective
friction angle(°)ψDilatancy
angle(°)λ*Modified
compression index(dimensionless)κ*Modified
swelling index(dimensionless)μ*Modified
creep index(dimensionless)A set of three parameters(c,φand,ψ) is
needed to model failure according to the Mohr-Coulomb criterion. Two further
parameters are used to model the amount of elastic and plastic strains and
their stress dependency. The modified compression index (λ*) represents the
slope of the normal consolidation line during one-dimensional or isotropic
logarithmic compression. In the same manner, the modified swelling index (κ*)
is related to the unloading or swelling line. The modified creep index (μ*)
serves as a measure to simulate the development of volumetric creep
deformations with the logarithm of time.In the modeling examples of this article, we will not take into account the development and the influence of water pressures, which would be also possible but imply a considerable increase in calculation effort; and in the case of the slope example in Section4.4, more reference data would be needed. All forward calculations were carried out applying the finite-element method using the commercial code PLAXIS (Version 8.2, professional, update-pack 8, build 1499) and considering the effect of
large deformations by means of an updated Lagrangian formulation (updated mesh
analysis) [22].


### 4.2. Oedometer Test

A one-dimensional compression test
was conducted by MFPA Weimar (Germany) [23] in a fixed oedometer ring with an inner diameter of around 7 cm (71.45 mm) and a height of around 2 cm (20.21 mm). Drainage was allowed on the top and at the bottom of the soil sample. All load steps were applied vertically, while the sample was held radially,
impeding horizontal displacements. First, the sample was preloaded with 9 kPa
during two days and with 13 kPa during one day. Then the load was doubled
successively, with each load step lasting 24 hours, loading the sample with 25,
50, 100, 200, 400, and 800 kPa. After that, it was unloaded at 400, 200, 100,
and 50 kPa, and finally it was reloaded again with 100, 200, 400, and 800 kPa
(last step took 43 hours). The displacements of the sample top were recorded
continuously.For the numerical model of the test setting, we used an axisymmetric geometrical configuration with the exact dimensions of the test specimen. In order to minimize calculation runtime, the discretization was done with two six-node triangular elements only, which is
the minimum possible number, as the software offers only triangular elements.
Thereby, the duration of a forward calculation could be reduced to less than
one minute on an ordinary personal computer. The accuracy of the deformation
results was checked by carrying out comparative analyses with finer meshes.
Horizontal fixities were assigned to the lateral boundary and to the rotation
axis, simulating the stiff oedometer ring and vertical fixities were attributed
to the basal boundary, representing the fix filter plate at the bottom. After
generating the initial stress state by applying the soil self-weight (gravity
loading procedure), distributed loads were applied perpendicular to the top
boundary analog to the laboratory conditions.Three parameters of the material model(λ*,κ*,andμ*) can be determined directly from the oedometer test. As
the material is known to show no dilatancy, ψ can be set to 0° in all
calculations. Laboratory data from shear tests on similar soil samples reported
by Panizza et al. [24] is shown in Table 3.Table 3
Shear-test data reported
by Panizza et al. 2006 [24].
Effective friction angleφ′(°)Effective cohesionc′(kPa)Direct shear tests18201810204918391649146918252020Triaxial tests1972814We averaged these values, giving
double weight to the triaxial test data, which we assumed to be more precise,
coming out with an average friction angle of20° and an average cohesion of 27
kPa. The set of experimental parameter values is shown by Table 4 and the results of a forward calculation using these parameters are presented in Figure 3 comparing them with the reference data of the oedometer test.Table 4
Parameters obtained from
experiments and identified parameters using PSO.
Experimental parameter valuesIdentified values PSOλ*0.0640.082κ*0.0350.051μ*1.46E–03—φ (°)20—c (kPa)2725ψ (°)0—Figure 3
Oedometer test calculation
results versus reference data.The graph shows that the deformations are underestimated by the simulation. In order to test the ability
of the PSO algorithm to find good parameter combinations, wide search areas
were chosen for the five parameters. A statistical analysis (see Section3) comprising 2000 calls of the forward calculation was then performed varying these parameters. Their value margins are displayed in Table 5.Table 5
Search intervals for the
parameters of the oedometer test.
MaximumMinimumln (max.)ln (min.)φ308——c1000.001——λ*10.0020.00−6.21κ*0.50.001−0.69−6.91μ*0.750.00001−0.29−11.51The scatter plot matrix of all data
points with objective function values lower than 10−7 is given in
Figure 4. For the parameters λ*,κ*,   and   μ*, logarithmized margins of the
search intervals were used in order to avoid overrepresentation of high
parameter values. Furthermore, a parameter constraint was prescribed, demanding
for λ*>κ*, which has to apply for all materials. The objective function
projections of c,κ*, and λ* indicate that the data points showing good model
fits seem to concentrate in quantifiable value ranges of these parameters,
whereas μ* and φ cannot be identified. The modified compression
index (λ*) and the cohesion (c) appear to be correlated (correlation
coefficient of 0.88), to a lesser extent, this holds true also for λ* and κ*
(correlation coefficient of 0.64). According to these findings, λ*,κ*, and c
were selected for the optimization procedure. The friction angle (φ) and the
modified creep index (μ*) were fixed on their experimental values.Figure 4
Scatter plot matrix for
the oedometer test.After 159 cycles (1590 calls) the particle
swarm optimizer had reduced the deviation to5*10−9, which was found
to be a sufficiently low value to stop the optimization routine. The identified
best parameter set and the corresponding calculation results can be seen also
in Table 4 and Figure 3. It becomes clear that the identified parameter set represents the measurement data much better than the available laboratory parameters.


### 4.3. Isotropic Compression Test

An isotropic compression test was
performed in the triaxial apparatus on a cylindrical soil sample with a
diameter of 5 cm (49.88 mm) and a height of 10 cm (98.55 mm). Drainage was allowed on top and at the bottom of the specimen. All load steps were performed
isotropically, applying a hydrostatic cell pressure. The sample was preloaded
at 30 kPa for three hours and, after that, at 50 kPa for 17 hours. Then it was
gradually loaded to 800 kPa in one hour, increasing the load by steps of 100 kPa. After reaching this target load, the stress level was left constant for 3
weeks. Top displacements and volume change of the sample were recorded during
the whole test. A reference dataset for the horizontal displacements was
calculated from the vertical displacements and the volume change, assuming the
shape of the specimen to remain exactly cylindrical until the end of the test.For the numerical model of the test
setting, an axisymmetric geometrical configuration with the exact dimensions of
the test specimen was used. Again, the model was discretized only with two
six-node triangular elements, to save calculation time. The accuracy of the deformation
results was checked by carrying out several comparative analyses with finer
meshes. Horizontal fixities were assigned to the rotation axis. Vertical
fixities were attributed to the basal boundary, representing the fix filter
plate at the bottom. After generating the initial stress state by applying the
soil self-weight (gravity loading procedure), two independent and identical
distributed loads were applied, one perpendicular to the upper boundary
(vertically) and the other one perpendicular to the lateral boundary
(radially). Loading was carried out the same way as in the laboratory, but
instead of the stepwise application of the 800 kPa target load, this load was
applied directly after the 50 kPa load step. For this reason, the 50 kPa load
step in the model was prolonged in such a manner that the integral of the load
as a function of time equals the test conditions.In the example of the isotropic
compression test, except for the relatively short phases before reaching the
target load, only one load step (800 kPa) is applied. Therefore, of the model
parameters onlyμ* can be determined directly from the test. This value (1.3*10−3) is very similar to the one obtained from the oedometer test. Figure 5 shows the results of a forward calculation using this value together with the laboratory values of Section 4.2.Isotropic compression
test calculation results versus reference data.
(a)(b)Also in this example, the deformations are underestimated by the simulation. A statistical analysis with 1820 calls was carried out varying the parametersc (cohesion), λ* (modified
compression index), and μ* (modified creep index) within the boundaries given
in Table 7, logarithmic values were used for the search intervals of the latter two parameters.As the modeled test contains no
unloading phases, it makes no sense to identify the modified swelling indexκ*.
Its value was therefore linked to the value of λ* by multiplying this parameter
by 0.5, which is
the typical κ*/λ* ratio, we observed in our laboratory tests performed on this material and similar materials.The results of the statistical
tests presented in Figure6 suggest that good fits can be obtained for cohesion
values between 20 kPa and 90 kPa, but apart from this, the cohesion value seems
to have no influence on the quality of the model calibration. Whereas for the
modified compression index (λ*) and the modified creep index (μ*), the objective function projections suggest a preferred value range for the data points with low deviation values.Figure 6
Scatter plot matrix for
the isotropic compression test.Furthermore, it can be concluded thatλ* and c might be quite closely related to each other (correlation coefficient 0.96), this means, for a given λ*, an appropriated cohesion value
could be computed by the equation given also in Figure 6. This linear relationship seems to be valid for cohesion values between 20 and 90 kPa and λ*
values between 0.064 and 0.165.Therefore, only the parametersλ* and μ* were selected for the optimization procedure via PSO. We stopped this procedure after 500 calls (50 cycles). The identified values are shown in Table 6. In Figure 5, the calculation results are compared to the reference dataset
and the results obtained by using only laboratory data.Table 6
Parameters from laboratory
tests and identified parameters using PSO.
Experimental valuesIdentified values (PSO)λ*0.0640.089κ*0.035—μ*1.30E–032.33E–03φ (°)20—c
(kPa)27—ψ (°)0—Table 7
Isotropic compression test
search intervals for the varied parameters.
MaximumMinimumln (max.)ln (min.)c
(kPa)1000.001——λ*1.006.74E–030.00−5.00μ*0.750.00001−0.29−11.51Like for the oedometer test, an improved parameter set could be found also for the isotropic compression test. It can be observed that the results are much better for the vertical displacements, although a weighting factor of one had been assigned to both
datasets.This may be due to the fact that the precision of the reference data is lower for the horizontal displacements, since the volume change of the sample could not be measured with the same
accuracy as the top displacements. In addition to this, small inhomogeneities
of the material or a small frictional resistance at the sample top could have
caused a slight distortion of the cylindrical shape of the specimen which was
assumed to calculate the horizontal displacements.


### 4.4. Deformations Along a Shear Zone in a Natural Slope

#### 4.4.1. Analyzed Section and Reference Data

The presented parameter identification technique was also applied to the2D-model of a natural slope which is located in the municipality of Corvara in the Dolomites (Italy). It
shows continuous creep deformations at the basis of a 20–40 m thick soil
cover consisting mainly of the above studied material and very similar
materials. As slopes of this type are known to show potential acceleration
phases that can endanger human settlements, a detailed study of the geology and
geomorphology as well as comprehensive monitoring were carried out by the
University of Modena and Reggio Emilia and by the National Research Council’s
Group for Hydrogeological Catastrophes Defense (CNR-GNDCI) [25] by order of the Autonomous Province of Bolzano/South-Tyrol [24].In this context, the displacements of several surface points were observed regularly by means of global positioning
system (GPS) measurements. In our study, a subarea of the slope was modeled
along a representative2D section. A site plan with the location of the section
and the GPS points is given in Figure 7, whereas Figure 8 depicts a field
survey of this section and shows the abstracted geometry model.Figure 7
Site map with the
location of the section and the GPS points.Field survey along the
section and derived geometry model.
(a)(b)


#### 4.4.2. Geotechnical Model

The geometry was determined using
all the information available, that is, a core drilling near section-pointB, the local geomorphology, refraction seismics, and direct current resistivity (DC-resistivity) [24]. As exposed in
Figure 9, the vertical profile of the slope was divided into three layers
interpreting various inclinometer profiles reported by Corsini et al. [25]; the illustrated one is located near section-point B. The uppermost layer, the soft soil cover, which is showing little internal deformations, was only considered
in the form of its weight acting on the intermediate layer, the shear zone.
Therefore, the displacement vectors of section-points C and F are presumed to
be equal to those of section-points B and E, respectively. The shear zone is
assumed to be a thin, soft, and highly plastic layer, exhibiting a pronounced
time dependency in its mechanical behavior. The third layer is given by the
underlying weathered bedrock, which is supposed to be stable. The earth
pressure at the foot of the slope was assumed to be slightly lower than the
earth pressure at rest; this means that there is no support obtained from the
soil layer further downslope. A description of the assumptions made for the
foot load is shown in Figure 10.Figure 9
Interpretation of the
inclinometer data from Corsini et al. 2005 [25].Figure 10
Earth pressure
assumptions made for the foot of the slope.The actual main detachment zone is modeled as an open crack. No tensile forces are acting across it onto the
downslope section of the sliding body, which is moving as a whole. Around section-pointA, the soil body below the secondary shear zone is assumed to move at the same
velocity as point B. The material properties of the secondary shear zone, which
is not subject of this study, were fixed to comply with this criterion. At
present, the displacement rates observed along the slope are more or less
constant, being superposed only by seasonal variations attributed to
fluctuations of the groundwater conditions which are not modeled in this study.
Therefore, our geotechnical model features displacement rates remaining constant
with time.


#### 4.4.3. Numerical Model

Figure11 shows the characteristics
of the numerical model of the studied slope. A plane strain geometrical
configuration with the real dimensions of the slope was used. The model was
discretized with 1070 triangular six-node elements. To save calculation time,
the number of elements was reduced by modeling only the uppermost 20 m of the
bedrock layer. The upper and the lower layers were meshed with the automatic
meshing procedure of the software and using a very coarse setting. A linear
elastic material model was assigned to them. Finally, one forward calculation
took approximately three minutes on an ordinary personal computer. The
intermediate layer was meshed manually by predefining geometry points in order
to assure a sufficiently fine mesh and a suitable orientation of the triangles
in order to go against an excessive distortion of the element shapes by the
calculated deformations. For the same reason, the updated mesh option was only
used in the last three years of the simulation (which are compared to the
reference dataset). The detachment zones were modeled by means of interfaces on
both sides of their geometry lines. The interfaces were modeled with a
Mohr-Coulomb material model, with negligible cohesion, the same friction angle
as the basal shear zone and with a constant reference stiffness in the order of
magnitude of the shear zone stiffness. The accuracy of the calculated
deformation results was checked by carrying out several comparative analyses
with finer meshes and also with a horizontal basal boundary. Horizontal
fixities were assigned to the lateral edges of the model, which extends over a
total length of 1080 m. Vertical fixities were attributed to the basal boundary,
representing the stable bedrock.Discretization of the slope
model: foot zone, vertical profile, and detachment zone.
(a)(b)(c)


#### 4.4.4. Calculation Phases

In the first calculation phase, all
three layers are made up by the bedrock material. An initial stress state is
generated by applying the self-weight of this material (gravity loading
procedure). In the second calculation phase, the two upper layers are replaced
by the weaker material of the soil cover. The third calculation phase marks the
starting point  of the slope instability in the model. The shear zone material
with time-dependent mechanical behavior is inserted and the horizontal load at
the foot is set. After that, the model is left creeping with unchanged boundary
conditions during a period of 33 years. As the loading history of the shear
zone material is unknown, this time period had to be chosen arbitrarily to
reach constant displacement rates as they are presently observed along the
natural slope.


#### 4.4.5. Results of Initial Model Using Parameters Derived from Experiments

In a first trial forward calculation, forλ*,κ*, and μ*, the parameters calculated from the laboratory experiments were used as input
values. As the deformations along the shear zone are known to persist since hundreds
or thousands of years [26], the shear strength of this zone has decreased to a
residual value that is characteristic for the soils originated by the
weathering of the San Cassiano and La Valle beds outcropping in the whole slope
area. Therefore, cohesion was assumed to be negligible (0.01 kPa) and a
friction angle of 10° was adopted, according to the average slope inclination
observed in nearby areas which were formed since the Late Glacial by the
studied processes (earth slides and earthflows) and covered by comparable soil
covers [27]. The stiffness of the uppermost layer was set equal to the stiffness
modulus observed in the oedometer test during unloading and reloading between
the load steps 400 kPa and 800 kPa. For Poisson’s ratio of this layer, we used
0.35, a value that is considered to be characteristic of clayey soils. The
experimental parameter values are shown in Table 8 and the deformations calculated on their basis for the last three years of the creep phase are
presented by Figure 12. The latter are only in the range of millimeters, and
thus not representing the actual situation in the field, where between
September 2001 and September 2004, displacements from several centimeters to
several decimeters were measured.Table 8
Laboratory values of the
parameters used for the slope example.
Experimental parameter valuesλ*0.064κ*0.035μ*0.00146c
(kPa)0.01φ (°)10G
(kPa)5560v0.35Figure 12
Results of a forward
calculation using laboratory values of the parameters.


#### 4.4.6. Results of Statistical Analysis and Optimization Procedure

A statistical analysis was carried
out (which will not be reported in detail here). One interesting finding of
this analysis was that the friction angle and the modified creep index appeared
to be closely correlated (coefficient of 0.92). The parametersλ*,κ*, and μ*;
the friction angle; and the stiffness of the uppermost layer (represented by
its shear modulus G) were chosen for the optimization procedure during which
they were varied within the intervals specified in Table 9.Table 9
Search intervals and
identified parameters for the slope example.
ParameterFixed parametersVaried parametersIdentified values (PSO)MaximumMinimumv0.35———G (kPa)—200002005160c (kPa)0.01———ψ (°)0———κ*—10.0050.60λ*—20.011.42μ*—1.50.0010.145φ (°)—16810.7After 82 cycles, each of them consisting of 10 forward calculations, the procedure was stopped because, from then on, the deviation could no longer be reduced significantly. The resulting parameter set is also given in Table9. Figure 13 depicts the calculated deformations
using the identified parameter combination, together with the displacement
vectors of the GPS measurement points.Figure 13
Calibrated slope model—comparison of modeled displacement
vectors (blue)
with reference data derived from the GPS measurements
(brown).It can be observed that the identified parameter set is able to reproduce the field measurements
qualitatively. Because of the simplifications made in the model, no exact fit
of the displacement vectors is possible. The presented back analysis procedure
gives one of a number of possible approximate solutions to the geotechnical
problem and the result returned by the particle swarm optimizer can be seen as
a parameter set that best represents the reference data.


### 5. Conclusions

A back analysis procedure for the identification of material parameters of constitutive models applied to geotechnical problems was presented. This procedure represents a direct
approach based on the method of least squares, correlation analyses, and a particle
swarm optimization algorithm. The applicability and suitability of the
technique was demonstrated by means of three examples from the fields of soil
mechanics and engineering geology. The studied material was a natural soil.
Besides being way more objective and less arbitrary than the conventional
trial and error procedure, the outlined method provides valuable information on
the quality of the model calibration, the uniqueness of an obtained solution,
or the determinateness of the problem. In all three examples, the particle swarm
optimizer was able to identify an improved parameter set after a justifiable
amount of forward calculations. Further research should also concentrate on the
identification of the geometrical parameters of geotechnical problems.


## 4.1. Description of the Studied Geological Material and the Adopted
Constitutive Model

In the following examples of applications of the presented parameter identification technique using PSO, we are modeling the mechanical behavior of a natural soil. It is of geotechnical interest because it favors the development of numerous landslides, namely,
rotational soil slips, earth slides, and earthflows. The studied material is clay that results from the weathering of structurally complex geological formations, the San Cassiano formation (Kassianer Schichten), and the La Valle formation (Wengener Schichten) of the Alpine
Trias. These rock formations are made up by interbedded strata of marls, tuffites,
claystones, limestones, dolomites, and sandstones. Like its source rocks, the
soil is characterized by a high clay and silt content. In the field, it
consists of a clayey matrix with coarser components, of diameters from
centimeters up to meters, floating in it without mutual support. Therefore, it
is not possible to sample the material as a whole representatively. As it has
been completely remolded by earthflow phenomena, no preferred orientation of
the components can be observed. For this study, only the fraction smaller than
2 mm was taken, as it is considered to determine the relevant
mechanical properties of the entire soil and it allows for the manufacturing of
reproducible samples and test specimens. Some properties of the studied
material are listed in Table1. Unless otherwise expressly stated, all tests and classifications were carried out according to the German standard DIN.Table 1
Properties of the studied soil.
Percentage
of clay particles (diameter<2 μm)26–34%Bulk
density (g/cm³)1.82Dry
density (g/cm³)1.29Density
of grains (g/cm³)2.76Porosityn53%Lime
content36%Loss
of ignition6.5%Liquid
LimitWL0.50Plastic
LimitWP0.27Water
content0.41From this description of the material, it becomes clear that its mechanical behavior is expected to be very complex and therefore it is only possible to model some important aspects of this behavior. Like many soils with a high clay and silt content, the studied
material is highly compressible and exhibits a significant amount of creep
deformations, thus its behavior is strongly time-dependent. As a constitutive
model, we chose the soft soil creep model, which was developed by Vermeer and
Neher [21] to account particularly for these phenomena. The soft soil creep model requires the following material parameters to be specified (see Table 2).Table 2
Parameters of the soft soil
creep model.
ParameterDescription(Unit)cEffective
cohesion(kPa)φEffective
friction angle(°)ψDilatancy
angle(°)λ*Modified
compression index(dimensionless)κ*Modified
swelling index(dimensionless)μ*Modified
creep index(dimensionless)A set of three parameters(c,φand,ψ) is
needed to model failure according to the Mohr-Coulomb criterion. Two further
parameters are used to model the amount of elastic and plastic strains and
their stress dependency. The modified compression index (λ*) represents the
slope of the normal consolidation line during one-dimensional or isotropic
logarithmic compression. In the same manner, the modified swelling index (κ*)
is related to the unloading or swelling line. The modified creep index (μ*)
serves as a measure to simulate the development of volumetric creep
deformations with the logarithm of time.In the modeling examples of this article, we will not take into account the development and the influence of water pressures, which would be also possible but imply a considerable increase in calculation effort; and in the case of the slope example in Section4.4, more reference data would be needed. All forward calculations were carried out applying the finite-element method using the commercial code PLAXIS (Version 8.2, professional, update-pack 8, build 1499) and considering the effect of
large deformations by means of an updated Lagrangian formulation (updated mesh
analysis) [22].


## 4.2. Oedometer Test

A one-dimensional compression test
was conducted by MFPA Weimar (Germany) [23] in a fixed oedometer ring with an inner diameter of around 7 cm (71.45 mm) and a height of around 2 cm (20.21 mm). Drainage was allowed on the top and at the bottom of the soil sample. All load steps were applied vertically, while the sample was held radially,
impeding horizontal displacements. First, the sample was preloaded with 9 kPa
during two days and with 13 kPa during one day. Then the load was doubled
successively, with each load step lasting 24 hours, loading the sample with 25,
50, 100, 200, 400, and 800 kPa. After that, it was unloaded at 400, 200, 100,
and 50 kPa, and finally it was reloaded again with 100, 200, 400, and 800 kPa
(last step took 43 hours). The displacements of the sample top were recorded
continuously.For the numerical model of the test setting, we used an axisymmetric geometrical configuration with the exact dimensions of the test specimen. In order to minimize calculation runtime, the discretization was done with two six-node triangular elements only, which is
the minimum possible number, as the software offers only triangular elements.
Thereby, the duration of a forward calculation could be reduced to less than
one minute on an ordinary personal computer. The accuracy of the deformation
results was checked by carrying out comparative analyses with finer meshes.
Horizontal fixities were assigned to the lateral boundary and to the rotation
axis, simulating the stiff oedometer ring and vertical fixities were attributed
to the basal boundary, representing the fix filter plate at the bottom. After
generating the initial stress state by applying the soil self-weight (gravity
loading procedure), distributed loads were applied perpendicular to the top
boundary analog to the laboratory conditions.Three parameters of the material model(λ*,κ*,andμ*) can be determined directly from the oedometer test. As
the material is known to show no dilatancy, ψ can be set to 0° in all
calculations. Laboratory data from shear tests on similar soil samples reported
by Panizza et al. [24] is shown in Table 3.Table 3
Shear-test data reported
by Panizza et al. 2006 [24].
Effective friction angleφ′(°)Effective cohesionc′(kPa)Direct shear tests18201810204918391649146918252020Triaxial tests1972814We averaged these values, giving
double weight to the triaxial test data, which we assumed to be more precise,
coming out with an average friction angle of20° and an average cohesion of 27
kPa. The set of experimental parameter values is shown by Table 4 and the results of a forward calculation using these parameters are presented in Figure 3 comparing them with the reference data of the oedometer test.Table 4
Parameters obtained from
experiments and identified parameters using PSO.
Experimental parameter valuesIdentified values PSOλ*0.0640.082κ*0.0350.051μ*1.46E–03—φ (°)20—c (kPa)2725ψ (°)0—Figure 3
Oedometer test calculation
results versus reference data.The graph shows that the deformations are underestimated by the simulation. In order to test the ability
of the PSO algorithm to find good parameter combinations, wide search areas
were chosen for the five parameters. A statistical analysis (see Section3) comprising 2000 calls of the forward calculation was then performed varying these parameters. Their value margins are displayed in Table 5.Table 5
Search intervals for the
parameters of the oedometer test.
MaximumMinimumln (max.)ln (min.)φ308——c1000.001——λ*10.0020.00−6.21κ*0.50.001−0.69−6.91μ*0.750.00001−0.29−11.51The scatter plot matrix of all data
points with objective function values lower than 10−7 is given in
Figure 4. For the parameters λ*,κ*,   and   μ*, logarithmized margins of the
search intervals were used in order to avoid overrepresentation of high
parameter values. Furthermore, a parameter constraint was prescribed, demanding
for λ*>κ*, which has to apply for all materials. The objective function
projections of c,κ*, and λ* indicate that the data points showing good model
fits seem to concentrate in quantifiable value ranges of these parameters,
whereas μ* and φ cannot be identified. The modified compression
index (λ*) and the cohesion (c) appear to be correlated (correlation
coefficient of 0.88), to a lesser extent, this holds true also for λ* and κ*
(correlation coefficient of 0.64). According to these findings, λ*,κ*, and c
were selected for the optimization procedure. The friction angle (φ) and the
modified creep index (μ*) were fixed on their experimental values.Figure 4
Scatter plot matrix for
the oedometer test.After 159 cycles (1590 calls) the particle
swarm optimizer had reduced the deviation to5*10−9, which was found
to be a sufficiently low value to stop the optimization routine. The identified
best parameter set and the corresponding calculation results can be seen also
in Table 4 and Figure 3. It becomes clear that the identified parameter set represents the measurement data much better than the available laboratory parameters.


## 4.3. Isotropic Compression Test

An isotropic compression test was
performed in the triaxial apparatus on a cylindrical soil sample with a
diameter of 5 cm (49.88 mm) and a height of 10 cm (98.55 mm). Drainage was allowed on top and at the bottom of the specimen. All load steps were performed
isotropically, applying a hydrostatic cell pressure. The sample was preloaded
at 30 kPa for three hours and, after that, at 50 kPa for 17 hours. Then it was
gradually loaded to 800 kPa in one hour, increasing the load by steps of 100 kPa. After reaching this target load, the stress level was left constant for 3
weeks. Top displacements and volume change of the sample were recorded during
the whole test. A reference dataset for the horizontal displacements was
calculated from the vertical displacements and the volume change, assuming the
shape of the specimen to remain exactly cylindrical until the end of the test.For the numerical model of the test
setting, an axisymmetric geometrical configuration with the exact dimensions of
the test specimen was used. Again, the model was discretized only with two
six-node triangular elements, to save calculation time. The accuracy of the deformation
results was checked by carrying out several comparative analyses with finer
meshes. Horizontal fixities were assigned to the rotation axis. Vertical
fixities were attributed to the basal boundary, representing the fix filter
plate at the bottom. After generating the initial stress state by applying the
soil self-weight (gravity loading procedure), two independent and identical
distributed loads were applied, one perpendicular to the upper boundary
(vertically) and the other one perpendicular to the lateral boundary
(radially). Loading was carried out the same way as in the laboratory, but
instead of the stepwise application of the 800 kPa target load, this load was
applied directly after the 50 kPa load step. For this reason, the 50 kPa load
step in the model was prolonged in such a manner that the integral of the load
as a function of time equals the test conditions.In the example of the isotropic
compression test, except for the relatively short phases before reaching the
target load, only one load step (800 kPa) is applied. Therefore, of the model
parameters onlyμ* can be determined directly from the test. This value (1.3*10−3) is very similar to the one obtained from the oedometer test. Figure 5 shows the results of a forward calculation using this value together with the laboratory values of Section 4.2.Isotropic compression
test calculation results versus reference data.
(a)(b)Also in this example, the deformations are underestimated by the simulation. A statistical analysis with 1820 calls was carried out varying the parametersc (cohesion), λ* (modified
compression index), and μ* (modified creep index) within the boundaries given
in Table 7, logarithmic values were used for the search intervals of the latter two parameters.As the modeled test contains no
unloading phases, it makes no sense to identify the modified swelling indexκ*.
Its value was therefore linked to the value of λ* by multiplying this parameter
by 0.5, which is
the typical κ*/λ* ratio, we observed in our laboratory tests performed on this material and similar materials.The results of the statistical
tests presented in Figure6 suggest that good fits can be obtained for cohesion
values between 20 kPa and 90 kPa, but apart from this, the cohesion value seems
to have no influence on the quality of the model calibration. Whereas for the
modified compression index (λ*) and the modified creep index (μ*), the objective function projections suggest a preferred value range for the data points with low deviation values.Figure 6
Scatter plot matrix for
the isotropic compression test.Furthermore, it can be concluded thatλ* and c might be quite closely related to each other (correlation coefficient 0.96), this means, for a given λ*, an appropriated cohesion value
could be computed by the equation given also in Figure 6. This linear relationship seems to be valid for cohesion values between 20 and 90 kPa and λ*
values between 0.064 and 0.165.Therefore, only the parametersλ* and μ* were selected for the optimization procedure via PSO. We stopped this procedure after 500 calls (50 cycles). The identified values are shown in Table 6. In Figure 5, the calculation results are compared to the reference dataset
and the results obtained by using only laboratory data.Table 6
Parameters from laboratory
tests and identified parameters using PSO.
Experimental valuesIdentified values (PSO)λ*0.0640.089κ*0.035—μ*1.30E–032.33E–03φ (°)20—c
(kPa)27—ψ (°)0—Table 7
Isotropic compression test
search intervals for the varied parameters.
MaximumMinimumln (max.)ln (min.)c
(kPa)1000.001——λ*1.006.74E–030.00−5.00μ*0.750.00001−0.29−11.51Like for the oedometer test, an improved parameter set could be found also for the isotropic compression test. It can be observed that the results are much better for the vertical displacements, although a weighting factor of one had been assigned to both
datasets.This may be due to the fact that the precision of the reference data is lower for the horizontal displacements, since the volume change of the sample could not be measured with the same
accuracy as the top displacements. In addition to this, small inhomogeneities
of the material or a small frictional resistance at the sample top could have
caused a slight distortion of the cylindrical shape of the specimen which was
assumed to calculate the horizontal displacements.


## 4.4. Deformations Along a Shear Zone in a Natural Slope

### 4.4.1. Analyzed Section and Reference Data

The presented parameter identification technique was also applied to the2D-model of a natural slope which is located in the municipality of Corvara in the Dolomites (Italy). It
shows continuous creep deformations at the basis of a 20–40 m thick soil
cover consisting mainly of the above studied material and very similar
materials. As slopes of this type are known to show potential acceleration
phases that can endanger human settlements, a detailed study of the geology and
geomorphology as well as comprehensive monitoring were carried out by the
University of Modena and Reggio Emilia and by the National Research Council’s
Group for Hydrogeological Catastrophes Defense (CNR-GNDCI) [25] by order of the Autonomous Province of Bolzano/South-Tyrol [24].In this context, the displacements of several surface points were observed regularly by means of global positioning
system (GPS) measurements. In our study, a subarea of the slope was modeled
along a representative2D section. A site plan with the location of the section
and the GPS points is given in Figure 7, whereas Figure 8 depicts a field
survey of this section and shows the abstracted geometry model.Figure 7
Site map with the
location of the section and the GPS points.Field survey along the
section and derived geometry model.
(a)(b)


### 4.4.2. Geotechnical Model

The geometry was determined using
all the information available, that is, a core drilling near section-pointB, the local geomorphology, refraction seismics, and direct current resistivity (DC-resistivity) [24]. As exposed in
Figure 9, the vertical profile of the slope was divided into three layers
interpreting various inclinometer profiles reported by Corsini et al. [25]; the illustrated one is located near section-point B. The uppermost layer, the soft soil cover, which is showing little internal deformations, was only considered
in the form of its weight acting on the intermediate layer, the shear zone.
Therefore, the displacement vectors of section-points C and F are presumed to
be equal to those of section-points B and E, respectively. The shear zone is
assumed to be a thin, soft, and highly plastic layer, exhibiting a pronounced
time dependency in its mechanical behavior. The third layer is given by the
underlying weathered bedrock, which is supposed to be stable. The earth
pressure at the foot of the slope was assumed to be slightly lower than the
earth pressure at rest; this means that there is no support obtained from the
soil layer further downslope. A description of the assumptions made for the
foot load is shown in Figure 10.Figure 9
Interpretation of the
inclinometer data from Corsini et al. 2005 [25].Figure 10
Earth pressure
assumptions made for the foot of the slope.The actual main detachment zone is modeled as an open crack. No tensile forces are acting across it onto the
downslope section of the sliding body, which is moving as a whole. Around section-pointA, the soil body below the secondary shear zone is assumed to move at the same
velocity as point B. The material properties of the secondary shear zone, which
is not subject of this study, were fixed to comply with this criterion. At
present, the displacement rates observed along the slope are more or less
constant, being superposed only by seasonal variations attributed to
fluctuations of the groundwater conditions which are not modeled in this study.
Therefore, our geotechnical model features displacement rates remaining constant
with time.


### 4.4.3. Numerical Model

Figure11 shows the characteristics
of the numerical model of the studied slope. A plane strain geometrical
configuration with the real dimensions of the slope was used. The model was
discretized with 1070 triangular six-node elements. To save calculation time,
the number of elements was reduced by modeling only the uppermost 20 m of the
bedrock layer. The upper and the lower layers were meshed with the automatic
meshing procedure of the software and using a very coarse setting. A linear
elastic material model was assigned to them. Finally, one forward calculation
took approximately three minutes on an ordinary personal computer. The
intermediate layer was meshed manually by predefining geometry points in order
to assure a sufficiently fine mesh and a suitable orientation of the triangles
in order to go against an excessive distortion of the element shapes by the
calculated deformations. For the same reason, the updated mesh option was only
used in the last three years of the simulation (which are compared to the
reference dataset). The detachment zones were modeled by means of interfaces on
both sides of their geometry lines. The interfaces were modeled with a
Mohr-Coulomb material model, with negligible cohesion, the same friction angle
as the basal shear zone and with a constant reference stiffness in the order of
magnitude of the shear zone stiffness. The accuracy of the calculated
deformation results was checked by carrying out several comparative analyses
with finer meshes and also with a horizontal basal boundary. Horizontal
fixities were assigned to the lateral edges of the model, which extends over a
total length of 1080 m. Vertical fixities were attributed to the basal boundary,
representing the stable bedrock.Discretization of the slope
model: foot zone, vertical profile, and detachment zone.
(a)(b)(c)


### 4.4.4. Calculation Phases

In the first calculation phase, all
three layers are made up by the bedrock material. An initial stress state is
generated by applying the self-weight of this material (gravity loading
procedure). In the second calculation phase, the two upper layers are replaced
by the weaker material of the soil cover. The third calculation phase marks the
starting point  of the slope instability in the model. The shear zone material
with time-dependent mechanical behavior is inserted and the horizontal load at
the foot is set. After that, the model is left creeping with unchanged boundary
conditions during a period of 33 years. As the loading history of the shear
zone material is unknown, this time period had to be chosen arbitrarily to
reach constant displacement rates as they are presently observed along the
natural slope.


### 4.4.5. Results of Initial Model Using Parameters Derived from Experiments

In a first trial forward calculation, forλ*,κ*, and μ*, the parameters calculated from the laboratory experiments were used as input
values. As the deformations along the shear zone are known to persist since hundreds
or thousands of years [26], the shear strength of this zone has decreased to a
residual value that is characteristic for the soils originated by the
weathering of the San Cassiano and La Valle beds outcropping in the whole slope
area. Therefore, cohesion was assumed to be negligible (0.01 kPa) and a
friction angle of 10° was adopted, according to the average slope inclination
observed in nearby areas which were formed since the Late Glacial by the
studied processes (earth slides and earthflows) and covered by comparable soil
covers [27]. The stiffness of the uppermost layer was set equal to the stiffness
modulus observed in the oedometer test during unloading and reloading between
the load steps 400 kPa and 800 kPa. For Poisson’s ratio of this layer, we used
0.35, a value that is considered to be characteristic of clayey soils. The
experimental parameter values are shown in Table 8 and the deformations calculated on their basis for the last three years of the creep phase are
presented by Figure 12. The latter are only in the range of millimeters, and
thus not representing the actual situation in the field, where between
September 2001 and September 2004, displacements from several centimeters to
several decimeters were measured.Table 8
Laboratory values of the
parameters used for the slope example.
Experimental parameter valuesλ*0.064κ*0.035μ*0.00146c
(kPa)0.01φ (°)10G
(kPa)5560v0.35Figure 12
Results of a forward
calculation using laboratory values of the parameters.


### 4.4.6. Results of Statistical Analysis and Optimization Procedure

A statistical analysis was carried
out (which will not be reported in detail here). One interesting finding of
this analysis was that the friction angle and the modified creep index appeared
to be closely correlated (coefficient of 0.92). The parametersλ*,κ*, and μ*;
the friction angle; and the stiffness of the uppermost layer (represented by
its shear modulus G) were chosen for the optimization procedure during which
they were varied within the intervals specified in Table 9.Table 9
Search intervals and
identified parameters for the slope example.
ParameterFixed parametersVaried parametersIdentified values (PSO)MaximumMinimumv0.35———G (kPa)—200002005160c (kPa)0.01———ψ (°)0———κ*—10.0050.60λ*—20.011.42μ*—1.50.0010.145φ (°)—16810.7After 82 cycles, each of them consisting of 10 forward calculations, the procedure was stopped because, from then on, the deviation could no longer be reduced significantly. The resulting parameter set is also given in Table9. Figure 13 depicts the calculated deformations
using the identified parameter combination, together with the displacement
vectors of the GPS measurement points.Figure 13
Calibrated slope model—comparison of modeled displacement
vectors (blue)
with reference data derived from the GPS measurements
(brown).It can be observed that the identified parameter set is able to reproduce the field measurements
qualitatively. Because of the simplifications made in the model, no exact fit
of the displacement vectors is possible. The presented back analysis procedure
gives one of a number of possible approximate solutions to the geotechnical
problem and the result returned by the particle swarm optimizer can be seen as
a parameter set that best represents the reference data.


## 4.4.1. Analyzed Section and Reference Data

The presented parameter identification technique was also applied to the2D-model of a natural slope which is located in the municipality of Corvara in the Dolomites (Italy). It
shows continuous creep deformations at the basis of a 20–40 m thick soil
cover consisting mainly of the above studied material and very similar
materials. As slopes of this type are known to show potential acceleration
phases that can endanger human settlements, a detailed study of the geology and
geomorphology as well as comprehensive monitoring were carried out by the
University of Modena and Reggio Emilia and by the National Research Council’s
Group for Hydrogeological Catastrophes Defense (CNR-GNDCI) [25] by order of the Autonomous Province of Bolzano/South-Tyrol [24].In this context, the displacements of several surface points were observed regularly by means of global positioning
system (GPS) measurements. In our study, a subarea of the slope was modeled
along a representative2D section. A site plan with the location of the section
and the GPS points is given in Figure 7, whereas Figure 8 depicts a field
survey of this section and shows the abstracted geometry model.Figure 7
Site map with the
location of the section and the GPS points.Field survey along the
section and derived geometry model.
(a)(b)


## 4.4.2. Geotechnical Model

The geometry was determined using
all the information available, that is, a core drilling near section-pointB, the local geomorphology, refraction seismics, and direct current resistivity (DC-resistivity) [24]. As exposed in
Figure 9, the vertical profile of the slope was divided into three layers
interpreting various inclinometer profiles reported by Corsini et al. [25]; the illustrated one is located near section-point B. The uppermost layer, the soft soil cover, which is showing little internal deformations, was only considered
in the form of its weight acting on the intermediate layer, the shear zone.
Therefore, the displacement vectors of section-points C and F are presumed to
be equal to those of section-points B and E, respectively. The shear zone is
assumed to be a thin, soft, and highly plastic layer, exhibiting a pronounced
time dependency in its mechanical behavior. The third layer is given by the
underlying weathered bedrock, which is supposed to be stable. The earth
pressure at the foot of the slope was assumed to be slightly lower than the
earth pressure at rest; this means that there is no support obtained from the
soil layer further downslope. A description of the assumptions made for the
foot load is shown in Figure 10.Figure 9
Interpretation of the
inclinometer data from Corsini et al. 2005 [25].Figure 10
Earth pressure
assumptions made for the foot of the slope.The actual main detachment zone is modeled as an open crack. No tensile forces are acting across it onto the
downslope section of the sliding body, which is moving as a whole. Around section-pointA, the soil body below the secondary shear zone is assumed to move at the same
velocity as point B. The material properties of the secondary shear zone, which
is not subject of this study, were fixed to comply with this criterion. At
present, the displacement rates observed along the slope are more or less
constant, being superposed only by seasonal variations attributed to
fluctuations of the groundwater conditions which are not modeled in this study.
Therefore, our geotechnical model features displacement rates remaining constant
with time.


## 4.4.3. Numerical Model

Figure11 shows the characteristics
of the numerical model of the studied slope. A plane strain geometrical
configuration with the real dimensions of the slope was used. The model was
discretized with 1070 triangular six-node elements. To save calculation time,
the number of elements was reduced by modeling only the uppermost 20 m of the
bedrock layer. The upper and the lower layers were meshed with the automatic
meshing procedure of the software and using a very coarse setting. A linear
elastic material model was assigned to them. Finally, one forward calculation
took approximately three minutes on an ordinary personal computer. The
intermediate layer was meshed manually by predefining geometry points in order
to assure a sufficiently fine mesh and a suitable orientation of the triangles
in order to go against an excessive distortion of the element shapes by the
calculated deformations. For the same reason, the updated mesh option was only
used in the last three years of the simulation (which are compared to the
reference dataset). The detachment zones were modeled by means of interfaces on
both sides of their geometry lines. The interfaces were modeled with a
Mohr-Coulomb material model, with negligible cohesion, the same friction angle
as the basal shear zone and with a constant reference stiffness in the order of
magnitude of the shear zone stiffness. The accuracy of the calculated
deformation results was checked by carrying out several comparative analyses
with finer meshes and also with a horizontal basal boundary. Horizontal
fixities were assigned to the lateral edges of the model, which extends over a
total length of 1080 m. Vertical fixities were attributed to the basal boundary,
representing the stable bedrock.Discretization of the slope
model: foot zone, vertical profile, and detachment zone.
(a)(b)(c)


## 4.4.4. Calculation Phases

In the first calculation phase, all
three layers are made up by the bedrock material. An initial stress state is
generated by applying the self-weight of this material (gravity loading
procedure). In the second calculation phase, the two upper layers are replaced
by the weaker material of the soil cover. The third calculation phase marks the
starting point  of the slope instability in the model. The shear zone material
with time-dependent mechanical behavior is inserted and the horizontal load at
the foot is set. After that, the model is left creeping with unchanged boundary
conditions during a period of 33 years. As the loading history of the shear
zone material is unknown, this time period had to be chosen arbitrarily to
reach constant displacement rates as they are presently observed along the
natural slope.


## 4.4.5. Results of Initial Model Using Parameters Derived from Experiments

In a first trial forward calculation, forλ*,κ*, and μ*, the parameters calculated from the laboratory experiments were used as input
values. As the deformations along the shear zone are known to persist since hundreds
or thousands of years [26], the shear strength of this zone has decreased to a
residual value that is characteristic for the soils originated by the
weathering of the San Cassiano and La Valle beds outcropping in the whole slope
area. Therefore, cohesion was assumed to be negligible (0.01 kPa) and a
friction angle of 10° was adopted, according to the average slope inclination
observed in nearby areas which were formed since the Late Glacial by the
studied processes (earth slides and earthflows) and covered by comparable soil
covers [27]. The stiffness of the uppermost layer was set equal to the stiffness
modulus observed in the oedometer test during unloading and reloading between
the load steps 400 kPa and 800 kPa. For Poisson’s ratio of this layer, we used
0.35, a value that is considered to be characteristic of clayey soils. The
experimental parameter values are shown in Table 8 and the deformations calculated on their basis for the last three years of the creep phase are
presented by Figure 12. The latter are only in the range of millimeters, and
thus not representing the actual situation in the field, where between
September 2001 and September 2004, displacements from several centimeters to
several decimeters were measured.Table 8
Laboratory values of the
parameters used for the slope example.
Experimental parameter valuesλ*0.064κ*0.035μ*0.00146c
(kPa)0.01φ (°)10G
(kPa)5560v0.35Figure 12
Results of a forward
calculation using laboratory values of the parameters.


## 4.4.6. Results of Statistical Analysis and Optimization Procedure

A statistical analysis was carried
out (which will not be reported in detail here). One interesting finding of
this analysis was that the friction angle and the modified creep index appeared
to be closely correlated (coefficient of 0.92). The parametersλ*,κ*, and μ*;
the friction angle; and the stiffness of the uppermost layer (represented by
its shear modulus G) were chosen for the optimization procedure during which
they were varied within the intervals specified in Table 9.Table 9
Search intervals and
identified parameters for the slope example.
ParameterFixed parametersVaried parametersIdentified values (PSO)MaximumMinimumv0.35———G (kPa)—200002005160c (kPa)0.01———ψ (°)0———κ*—10.0050.60λ*—20.011.42μ*—1.50.0010.145φ (°)—16810.7After 82 cycles, each of them consisting of 10 forward calculations, the procedure was stopped because, from then on, the deviation could no longer be reduced significantly. The resulting parameter set is also given in Table9. Figure 13 depicts the calculated deformations
using the identified parameter combination, together with the displacement
vectors of the GPS measurement points.Figure 13
Calibrated slope model—comparison of modeled displacement
vectors (blue)
with reference data derived from the GPS measurements
(brown).It can be observed that the identified parameter set is able to reproduce the field measurements
qualitatively. Because of the simplifications made in the model, no exact fit
of the displacement vectors is possible. The presented back analysis procedure
gives one of a number of possible approximate solutions to the geotechnical
problem and the result returned by the particle swarm optimizer can be seen as
a parameter set that best represents the reference data.


## 5. Conclusions

A back analysis procedure for the identification of material parameters of constitutive models applied to geotechnical problems was presented. This procedure represents a direct
approach based on the method of least squares, correlation analyses, and a particle
swarm optimization algorithm. The applicability and suitability of the
technique was demonstrated by means of three examples from the fields of soil
mechanics and engineering geology. The studied material was a natural soil.
Besides being way more objective and less arbitrary than the conventional
trial and error procedure, the outlined method provides valuable information on
the quality of the model calibration, the uniqueness of an obtained solution,
or the determinateness of the problem. In all three examples, the particle swarm
optimizer was able to identify an improved parameter set after a justifiable
amount of forward calculations. Further research should also concentrate on the
identification of the geometrical parameters of geotechnical problems.



---

*Source: 574613-2008-05-04.xml*