Journal of the Institute of Engineering, Vol. 8, No. 1, pp.

71–92
TUTA/IOE/PCU © TUTA/IOE/PCU
All rights reserved. Printed in Nepal
Fax: 977-1-5525830

Finite Element Modeling and Decoupled Seismic Stability Analysis of a

Zoned Rockfill Dam Designed By Traditional Empirical Methods

Shanker Dhakal1, N. P. Bhandary2, R. Yatabe2, P. L. Pradhan3, R. C. Tiwari4

1
Department of Civil Engineering, Nepal Engineering College, Pokhara University, Nepal
2
Department of Civil and Environmental Engineering, Ehime University, Japan
3
Civil Engineering Department, Pulchowk Campus, Institute of Engineering, Tribhuvan University, Nepal
4
Graduate School of Science and Engineering, Ehime University, Japan
Corresponding email: shankeronly@gmail.com

Abstract: Despite the recommendation made by ICOLD a few decades ago, the high and important
zoned rockfill dams, either existing or proposed, located in the region of high seismicity like Nepal,
have been found still designed simply by following the traditional empirical methods. The authors,
having identified this gap, became interested to carry out a research to know: how such dams would
perform when their seismic safety was evaluated by dynamic analysis in terms of crest settlement,
which might lead to the loss of freeboard and hence the overtopping (collapse) of the structure. Also, as
a secondary objective, the authors wished to address: how the effective geometry of the dam foundation
in the structural model may be determined when flexibility of foundation is considered by trial and error
‘vertical stress criterion’.
This paper puts forth the summary of the research. The seismic safety is evaluated by adopting a
simplified technique (‘decoupled seismic deformation analysis’) where the peak crest acceleration
response obtained from the time history analysis of the 2D plain strain finite element model (both with
and without foundation flexibility considered) is utilized to predict the potential permanent deformation
(settlement) of crest using the method based on Newmark’s Analysis Procedure.

Keywords: Zoned rockfill dam, Seismic deformation analysis, Finite element method, Displacement
decoupling, Time history analysis

1. Introduction
1.1 Zoned Rockfill Dam and Its Prospects
Dams are the structures mainly meant to store water for its controlled use (regulated flow).
Reservoir projects and the associated dams are spread all around the world. Embankment (Fill)
dams cover the majority of them; more than 80% [1]. International Commission of Large Dams
(ICOLD) differentiates the rockfill dam from the earthfill dam such that the former has its
embankment consisting of rocks more than 50 % of its volume. Rockfill dam can further be
classified as membrane (deck) type rockfill dam such as concrete faced rockfill dam (CFRD), and
zoned type rockfill dam (‘ZRD’) such as the central clay core rockfill dam (‘CCRD’) with
inverted filters in between the inner clay core and outer rockfill shells.
72 Journal of the Institute of the Engineering

Availability of local construction materials, high inherent stability and higher capacity to absorb
shocks in extreme loading [2,3] such as earthquakes mainly determine the selection of a zoned
rockfill dam. Nepal, for example, (a small developing country but very rich in its hydro-potential)
favors the construction of zoned rockfill dams. The existing dam of Kulekhani Hydroelectric
Project (the only one seasonal storage project in Nepal) and that of many of the proposed storage
projects such as Bagmati Multipurpose Project (an attractive project in the Eastern Region of the
country) are of ‘CCRD’ type.

1.2 Seismic Stability Analysis Concerns in a Zoned Rockfill Dam Design

Assuring the seepage and piping related problems are under control, the next important concern
for the stability of a zoned rockfill dam in an earthquake event appears to be the stress-
deformation analysis. Again, despite the fact that stress analysis may have its own significance,
deformations especially those leading to the loss of freeboard, should be given a prime
importance [4-6]. In this connection, the prediction of potential displacement (settlement) of the
embankment crest seems to be the most important response parameter in assessing the stability of
a fill dam. Fig.1 is drawn to illustrate the phenomena.

Fig. 1: A schematic illustration of loss of freeboard and overtopping in fill dams

Conventionally, the fill dams are designed utilizing the empirical relations and slope stability is
checked by limit equilibrium method (LEM) employing a number of trial failure surfaces to
assure that the factor of safety (FoS) meets the requirement of minimum specified value. Therein,
the earthquake effects are incorporated in terms of seismic coefficients. However, this method
has several flaws; particularly it does not give any information about the deformations in the
body of embankment to assess the dam stability against the loss of freeboard.
A few researches such as [7,8] are available which have used finite element analysis results to
calculate the factors of safety by following the basics of LEM. But, again they are not free from
the flaws of LEM.
Further, the conclusion of the 13th international conference on large dams in 1979 was that
“During earthquake resistant design of large dams or of important dams built in areas with high
seismicity (PGA > 0.2 g), an appropriate dynamic analysis should be carried out with data as
detailed as possible” [9]. This is, in fact, true because the pseudo-static method when applied to a
zoned rockfill dam, for example, concerns only a particular aspect of the behavior of the dam
with respect to an earthquake. The problem of the seismic behavior of this type of dam is,
however, more complex because of the very nature of the material. Therefore it seems essential
to design a zoned rockfill dam by the dynamic method.

1.3 NAP and Decoupled Seismic Analysis Method

Newmark (1965) first attempted to analyze seismic slope stability by quantifying the amount of
displacement due to a failure. He made the analogy that the sliding mass on the failure surface is
 Finite Element Modeling and Decoupled Seismic Stability . . . 73

similar to a block resting on an inclined plane [10] and that the contact between the potential
sliding mass and the underlying slope was assumed as rigid-plastic [11].
He found that increasing the amount of shaking decreased the factor of safety, as expected and
that for a given frictional resistance, there was one particular seismic coefficient (k) or
acceleration (a) that produced a factor of safety of one. This threshold coefficient or acceleration
at which the factor of safety is 1 was named the yield coefficient (k y) or the yield acceleration (ay)
respectively. When a block was subjected to acceleration greater than the yield acceleration, the
block moved. The permanent displacement of the block is found by double integrating the
acceleration time history where the acceleration exceeded the yield acceleration [12,13].
Later, few other scientists worked based on his procedure (commonly known as Newmark’s
Analysis Procedure or NAP, in short) and derived further relationships to facilitate the use of
NAP framework. For example, Hynes-Griffin and Franklin (1984) related the estimated
deformations of a sliding block to the ratio of yield acceleration to maximum acceleration
through upper-bound, median, and lower bound-bound correlations [14]. Being in very high
seismicity zone, the upper bound relationship proposed by Hynes-Griffin and Franklin (Fig. 2)
has been used for the case of zoned rockfill dam studied in this research.
As a matter of fact, more accurate and reliable predictions of permanent deformations can be
obtained by dynamic analysis using the elasto-plastic nonlinear codes. However, in the absence
of such rigorous codes, fairly approximate estimation can be obtained by some simplified
method. Research have revealed that results in terms of acceleration time histories of linear or
equivalent linear analyses performed by shear-beam or FEM can be used as the input for the
semi-empirical methods based on NAP such as Hynes-Griffin and Franklin (1984) as discussed
earlier to assess the permanent deformation indirectly[15,16].

Fig. 2: The upper-bound correlation between ay/amax and permanent deformation (modified from
Hynes-Griffin and Franklin, 1984)
Therefore the method may be named as ‘decoupled seismic stability analysis’ or ‘decoupled
seismic deformation analysis’. The overview of the seismic analysis and the decoupled
74 Journal of the Institute of the Engineering

displacement method may be depicted from Fig. 3. Here, the bold letters represent the methods
that the authors have exercised in this research.

1.4 Finite element modeling of zoned rockfill dam

By virtue of its irregular geometrical and material compositions (from shell to core, from
pervious to impervious, from dry to saturated, etc.), the zoned rockfill embankment cannot be
modeled appropriately with approaches other than Finite Element Method (FEM).

Fig. 3: Decoupled seismic safety evaluation model

Researchers have developed many finite element computer programs for the dynamic analysis of
embankment dams. Some of the codes are such that they can fully incorporate the complex
nonlinear plasticity constitutive models and hence more accurate and reliable estimates of
stresses and deformations can be obtained [17-20].
Researchers have also contributed much to answer whether 2D modeling is sufficient or one
should go for 3D solid modeling, in the earthquake analysis of dam structures [21-23].
Regarding the effect of foundation on structure: conventionally, the bottom nodes of the finite
elements in the mesh are fixed (completely restrained) assuming the foundation (subsurface) to
be sufficiently rigid, such as that done by Alberro (1972), and Knight and Colleagues (1985)
[24]. However, some researchers consider the flexibility of the foundation in the model by
incorporating some effective domain of the subsurface (foundation) continuously into the
structural model of embankment [25-27]. Again, many have also incorporated the foundation-
structure interactions by modeling the phenomena via spring dashpot system, such as [28].

1.5 Research problem statement and objectives

Already, an explanation on the increasing understanding over the importance of dynamic analysis
for fill dams and need to evaluate their safety in terms of loss of freeboard, has been put down
somewhere before. Despite this, also being provoked by ICOLD right from the end of seventies,
plenty of evidences are found that even in the regions of very high seismicity such as in Nepal
[29-31], the high and important (zoned rockfill) dams either built [32] or proposed [33], have
been designed simply following the traditional empirical methods, and the seismic stability is
checked only in terms of factor of safety by pseudo-static limit equilibrium method. It appeared
surprising to the authors and it urged them to address this very gap.
 Finite Element Modeling and Decoupled Seismic Stability . . . 75

They became interested to know: how the zoned rockfill dams that were designed by traditional
empirical methods would perform when their seismic safety was evaluated in terms of the
potential permanent settlement of the crest which might lead to the loss of freeboard and
complete collapse followed by overtopping.
In addition to the ultimate objective of evaluating the seismic safety in terms of permanent crest
settlement, finite element modeling is another inherent objective of this research. This research
must be the first of its kind to introduce in Nepal.
Here the response is evaluated by considering the foundation both to be rigid and flexible. Rather
than following the conventional guidelines, the authors have tried to go from the first principle. A
trial and error procedure aiming to obtain negligible vertical stresses on all the three boundaries
of the subsurface domain is adopted. This may be named as the ‘vertical stress criterion’.
As a reference problem for the study, the zoned rockfill dam of the proposed Bagmati
Multipurpose Project (BMP) of Nepal, mentioned somewhere earlier, is chosen because the first
author has some relation with it; sometimes in past, he had chosen this storage project for his
undergraduate degree dissertation work, where emphasis was put on the hydrological and the
hydraulic design [34].

1.6 Comment on the adopted method

The authors were encouraged to conduct the present research not only from the point of view of
scientific aspects of zoned rockfill dam but also to introduce a practice of rigorous research on
the area of zoned rockfill dam in the country like Nepal, and bring awareness on the need of
dynamic analysis of zoned rockfill dam in the high seismic Himalayan region like Nepal. Having
realized the situation that a designer in such area is likely not to have access to the advanced and
rigorous computer codes and it may not be possible for him to develop them himself, effort was
made to work under the existing constraints, for the time being. Therefore, the simplified method
of decoupled seismic analysis has been adopted to evaluate the seismic safety in terms of
permanent displacement (crest settlement). The available software package of SAP 2000
(V10.0.1), probably the most commonly used program in Nepal, is used as the numerical tool.
However, the authors do give due emphasis in introducing yet rigorous and nonlinear dynamic
analysis design tool (rockfill dam specific, but not general) in future which may simulate the real
behavior of a zoned rockfill dam including interface interaction modeling, so that the safety may
be obtained directly in terms of excessive deformations and stresses..

2. Previous Research
As summarized in [16,17], the numerical modeling techniques such as the finite element method
were first applied to the dynamic analysis of embankment dams by Clough and Chopra (1966).
This was followed by major improvements by Gaboussi (1967), Schnabel et al. (1972), Gaboussi
and Wilson (1973), Idriss et al. (1974), Martin et al. (1975), Finn et al. ( 1977), Lee and Finn
(1978), White et al. (1979), Zienkiewicz and Shiomi (1984), Finn et al. (1986), Medinna et al
(1990), Li et al (1992) and so on. Scientists have developed from preliminary finite element
codes such as SHAKE, QUAD-4, FLUSH to tools such as DIANA, ANSYS, FLAC, ADINA,
etc. capable of handling rigorous constitutive models.
As per the seismic displacement evaluation, researchers reveal that the analysis can be carried out
using the decoupled approach that generally provides a conservative estimate of the seismic-
76 Journal of the Institute of the Engineering

induced permanent displacements [35-38]. The displacement-based approach has been applied in
a number of analyses of ideal or real earth dams, rockfill dams and even concrete dams with base
sliding [39-42] and, in a few documented case histories, it has been successfully used to back-
calculate the measured seismic-induced permanent displacements [43-45].

3. Finite Element Modeling

3.1 Embankment section geometry
BMP high dam is a zoned type rockfill dam with central earth (clay) core. The dam is 117 m high
above the lowest foundation. The crest is 10 m wide and about 700 m long. So the ratio of length
between canyon walls and height of the maximum cross section, L/H is about 6. Fig. 4 shows the
simplified structural framework (with acceptable approximation) of the maximum cross section
of the dam, designed by the classical empirical methods.

3.2 Finite element mesh, restraints and models

The 2-D plane strain assumption was made based on the rationale that the L/H ratio is greater
than 5 [21, 22]. The finite element mesh consists of Isoperimetric quadrilateral and triangular
elements. To simulate the lift construction and incorporate the inhomogeneity [46]; with the
assumption that core is more affected, as adopted by Knight and colleagues (1985) [24], the core
was divided into ten different layers.

T W L = 2 3 0 . 0 0 , F S L = 2 2 3 .0 0 , M O L = 2 0 0 . 0 0

I m p e r v io u s C o r e
R o c k F il l
G r a v e l F i ll
R a n d o m F il l
R i v e r B e d A l lu v i u m
R o c k L in e

Fig. 4: Maximum cross section of BMP dam to be modeled (with acceptable approximations)
Two models were adopted; one with rigid foundation (all nodes at the base fully restrained) and
other considering the foundation flexibility by incorporating some effective part of foundation
into the model. The boundary conditions for this model are such that all the bottom nodes are
restrained in all directions while vertical nodes at either side of foundation are restrained against
translation in horizontal direction in addition to the common plane strain restraints.
 Finite Element Modeling and Decoupled Seismic Stability . . . 77

Regarding the effective part of foundation to be taken in the model, the trial and error ‘vertical
stress criterion’ aforementioned somewhere was used. For this, a number of trials with an
increase in the dimensions of foundation on either side and underneath the embankment were
made until the vertical stresses at the boundaries came out negligible (less than 10 to 20 percent
of the maximum). The foundation beyond this effective part is neglected. The final effective
geometry of the foundation (as the multiple of base width of embankment) was obtained to be
twice the base width on either side and four times of the base width underneath.
Fig. 5 and Fig. 6 respectively show the two dimensional finite element models assuming the
foundation to be rigid (RF model) and considering the flexibility of foundation (FF model).
For the purposes of analysis, the reservoir was assumed at the highest flood level. The
corresponding free surface of water (phreatic line) in the core is located using the conventional
procedure described in standard text books such as [47]. The required extra meshing along the
line was done manually.
Below this phreatic line, the elements are saturated and elements above this line are assumed dry
(unsaturated). The model demonstrating the division of core into saturated and dry parts is
presented in Fig. 7.

3.3 Material properties

BMP high dam consists of impervious clay core covered by rock, gravel and random fill shells
resting on a stiff foundation layer of gravel and sand alluvium deposit down to bedrock. The rock
types available nearby are reported to consist mainly of shale and sandstone [48]. The clay is
assumed to be of rolled stiff type.

Fig 5: Rigid foundation finite element model

Fig. 6: Flexible foundation finite element model

78 Journal of the Institute of the Engineering

Fig. 7: Division of core by the phreatic line.

Based on this general information available in the design report, literatures such as [49-52, 9] and
other open source links were surveyed and the material characteristics were assumed. In fact,
confidence on the reliability of these would only be achieved if the arrangements for the tests
could be made which is beyond the scope of this research.
Table 1 presents the information on the densities of homogeneous materials under normal
(unsaturated) and saturated conditions. The shells are assumed free-draining, the core to have
very low permeability and the rock foundation to be fully impervious.
Table 1: Unit weights (kN/m3)

Material Type Unsaturated state Saturated state

Rock Fill- Shell 19.90 22.10

Gravel Fill- Shell 18.50 20.00
Random Fill- Shell 18.50 20.00
Rolled Clay- Core 19.80 20.60
Alluvial- Foundation 19.50 21.50
Rock- Foundation 21.00 -
Regarding the elastic properties (modulus of elasticity and Poisson’s ratio) required for the linear
elastic analysis, the materials were assumed to be homogeneous, elastic and isotropic except for
core which, as already mentioned before, was divided into ten different layers but each layer
again assumed homogenous, elastic and isotropic. Due to the lack of experimental data, but with
acceptable approximations for structural modeling, all the shells and filter were assigned the
same elastic parameters. The elastic material data used are assumed to be from the suitable
drained triaxial test results in the case of shells and the undrained triaxial test results for the case
of core. These data for the unsaturated materials are presented in Table 2.
For saturated materials, again due to the lack of data, an assumption was made. With the
assumption that the rock skeleton are intact and shell body is free draining, the elastic parameters
are not likely to change appreciably; the only effect could be the very small contribution made by
water as lubricant.
However, due to the effect of water, the clay minerals are likely to react and the core material
may suffer degradation of strength.
 Finite Element Modeling and Decoupled Seismic Stability . . . 79

Table 2: Unsaturated materials elastic properties

Material Type E (kN/m2) υ
Rock Foundation 50,000 0.35
Alluvial Foundation 150 0.38
Shells 3,000 0.25
1/10 20 0.42
2/10 33.63 0.42
3/10 45.59 0.42
4/10 56.56 0.42
Core in the order of depth 5/10 66.87 0.42
below crest 6/10 76.67 0.42
7/10 86.06 0.42
8/10 95.13 0.42
9/10 103.91 0.42
10/10 112.46 0.42
To incorporate this effect and also the enhanced pore pressure which ultimately reduce the
strength, the elastic modulus of the clay was approximately assumed to reduce by 15 % [51, 9].

3.4 Input ground motion

Due to the lack of records of any strong earthquake ground motion in the vicinity of the BMP
dam site, the famous imperial valley earthquake of May 18, 1940 (El Centro, California, USA)
was chosen as the prescribed input earthquake. Both the N-S and vertical components are used.
Fig. 8 shows the input ground acceleration time histories of the excitations in the upstream-
downstream (N-S component) and vertical directions (vertical component) of the dam.

N-S component
80 Journal of the Institute of the Engineering

Vertical Component

Fig. 8: Input ground motion (El-Centro, 1940)

V1: time in second, V2: acceleration (in unites of g)
(Source: USGS)

3.5 Modeling the water pressures

In the case of a region with high permeability such as a rock zone, it may be safely assumed that
no steady state infiltration force is developed. When the rate of change in level of stored water is
high, so that the impervious zone is not in a steady state infiltration condition, an unsteady
seepage force develops. But since it is extremely complicated to consider the unsteady seepage
force in dynamic analysis, the situation is sometimes simplified by assuming that the static
hydraulic force acts on the upstream side of the impervious zone [9].The present research has
followed this basic idea.
The hydrostatic pressure at any depth d below the free surface of water is given by γw.d normal to
the surface; γw being the unit weight of water. Thus if the height of water surface above the base
on the core is H, then the static water pressure at a height z above the base in given by γw (H-z).
This linear pressure distribution of water pressure can be assigned to the nodes on the u/s face as
the area surface pressure loads by assigning joint pattern in SAP 2000. The joint pattern available
is in the form A x + B y + C z + D; so for z measured above the bottom of core, A = 0, B = 0,
C = -1 and D = H and the multiplier for surface pressure load = γw = 9.81 kN/m3.
The hydrodynamic added pressures were computed based on Zangar’s formulation simply
following the recommendation by Indian design standard [13]. In pseudo-static analysis, the
pressures are converted into nodal forces on the basis of tributary area by integrating the definite
integral of regression curve representing the hydrodynamic pressures. But for the case of time
history analysis, the constant values of hydrodynamic pressure forces at the nodes may not be
logical. Therefore, the concept of equivalent virtual added masses was utilized to model
hydrodynamic pressures. The hydrodynamic pressure forces were first equated to the equivalent
inertial forces added by the reservoir and hence the shape of the body of reservoir, that can be
assumed to move with the dam, was worked out. Again, the added masses to be lumped as the
nodal masses were computed on the basis of tributary areas. Table 3 shows the calculated
equivalent hydrodynamic added masses for the cores of main dam and that of the upstream coffer
dam, over which the main dam has been designed to be constructed.
 Finite Element Modeling and Decoupled Seismic Stability . . . 81

Table 3: Hydrodynamic added masses (kN-s2/m)

Relative ht. above base For Main Dam For Coffer Dam
0.1 774.620 185.139
0.2 757.354 182.239
0.3 734.997 175.738
0.4 703.558 166.579
0.5 659.046 155.125
0.6 597.468 141.161
0.7 514.834 123.896
0.8 407.152 101.961
0.9 270.432 73.410
1.0 100.681 35.716
The virtual hydrodynamic added masses thus computed were lumped at the u/s face of the core of
the main dam and both u/s and d/s faces of the core of the coffer dam..
It is believed that not considering the reservoir-embankment interaction would have negligible
effects in the present case [17].

3.6 Model verification

Because the authors were working with the general purpose finite element program, they wished
to verify a few aspects of the modeling before actually using the models to carry out the time
history analyses. Two simple verifications were done: one for foundation stresses and other for
gravity turn-on deformations of embankment.

3.6.1 Foundation stresses

In fact, the distribution of vertical stress contours in the foundation in each trial were obtained
comparing well with the pressure bulb diagram obtained by the classical theory of soil mechanics
[54]. However, to verify the magnitude too, a simple case of a long wall was modeled the plain
strain finite element results (modeled in SAP 2000) of which could be readily compared
analytically with that of the line load formulation. It was a concrete wall (1 m thick, 3 m height;
E = 2 x 104 MPa, ν = 0.20) founded on a clayey soil (E = 20 MPa, ν = 0.42). The vertical stress
contours thus obtained under the gravity load of the wall are shown in Fig. 9.
The stresses at different locations in the foundation were computed using line load formulation of
Boussinesq’s equation, and found comparable with those obtained by the FE analysis. As an
example, at a depth of 2 m vertically below the centre, FE solution gives a vertical compression
stress equal to 21.81 kN/m2 while the soil mechanics formula for line load case [37] for the same
point gives:
q 2
σz = = 22.92 kN/m2
z Π
(Since z = 2 m , q = (1 x 3) x 24 = 72 kN/m)
The average variation is below 5 %.
82 Journal of the Institute of the Engineering

Fig. 9: Pressure bulbs under a long wall by FEM

Also the parametric studies were made which revealed that changing the elastic parameters of the
foundation has negligible effects on the vertical stresses in foundation, complying with the
formula and also in accordance with a finding by Clough and Woodward in 1967 [24].
Thus, the first aspect of verification is done which gives confidence on the worked out geometry
of effective foundation of the FF model.

3.6.2 Embankment deformations

As revealed from Fig. 10, the vertical displacement (settlement) profile of the zoned rockfill dam
within the core for the linear static analysis under the self weight load at the end-of-construction
(gravity turn-on) stage shows an increasing trend from base towards crest until some depth below
the crest where maximum settlement occurs and again decreases for the portion over it. The
horizontal deformation profile is however symmetrical (Fig. 11).
 Finite Element Modeling and Decoupled Seismic Stability . . . 83

120.0

Height above the base (m)

110.0
100.0
90.0
80.0
70.0
RF
60.0
FF
50.0
40.0
30.0
20.0
10.0
0.0
0 100 200 300 400

Displacement (mm)

Fig. 10: Vertical deformation profile of core

(through the centre of the core)
400
350
Displacement (mm)

300
250
RF
200
FF
150
100
50
0
-20.00 -10.00 0.00 10.00 20.00

Distance from core centre ( m)

Fig. 11: Horizontal deformation profile of core

(through maximum vertical displacement point)
As expected the displacements in FF model are larger than that in RF model.
However, the vertical displacement profile is different from the general perception that the
maximum elastic settlement of an embankment under its self weight should occur at the crest.
The finite element investigation of a homogeneous fill dam carried out by Clough and Woodward
also reveals this. Nevertheless, the section under study here is not homogenous but it is zoned.
Therefore, it required further investigation. Parametric studies were performed taking a small
symmetrical rockfill dam model (equivalent to the upstream coffer dam of the BMP dam). First,
this zoned rockfill dam was analyzed and the settlement profile was plotted (Fig. 12). Then, the
dam is assigned with homogeneous material properties, analyzed and again the settlement profile
was plotted (Fig. 13). Figure 13 compares exactly with the finding by Clough and Woodward
(and as expected) that “for the homogeneous fill dam, under gravity loads without reservoir, the
settlement increases from base toward crest; with steep increase at the bottom portion while
difference is small toward crest”.
84 Journal of the Institute of the Engineering

1.2

Relative ht above the base

1.0

0.8

0.6

0.4

0.2

0.0
0.000 50.000 100.000 150.000 200.000 250.000

Settlement (mm)

Fig. 12: Deformation profile for the arbitrary embankment with central core zone
Investigation with varying material properties too lead to the same trends.
For the purposes of analysis, even the embankment materials of the BMP embankment were all
assigned with a single hypothetical material and also the core of the zoned dam was assigned
with some hypothetical material stiffer than the shell. But, in either case, the profile similar to
that of the homogenous dam is resulted, with maximum vertical deformation at the crest.
1.2
Relative ht above base

1.0

0.8

0.6

0.4

0.2

0.0
0 1 2 3 4 5 6 7

settlement (mm)

Fig. 13: Deformation profile for the arbitrary homogeneous embankment

Thus, it may be interpreted that it is the relative stiffness of the weak core and strong shell of a
zoned rockfill dam that results the vertical deformation profile in the core of BMP zoned rockfill
dam as already discussed, with maximum value not at the crest but at some depth below the crest.
It may be wise to remind the readers here that modeling the likely interactions at the core-shell
interfaces of a zoned dam is assumed beyond the scope of this paper, and is left as a scope for
future studies.
Having discussed so far, it seems that the rigid and flexible foundation finite element models
prepared on the general purpose platform of SAP 2000 are valid and may be utilized to work
further.

4. Seismic Analysis and Results

4.1 Time history analysis
The uncoupled modal equations of motion for MDOF system could be solved in closed form if
the excitations were a simple function but the time-stepping numerical methods are necessary for
complex excitations such as earthquake ground motion. [55]
 Finite Element Modeling and Decoupled Seismic Stability . . . 85

Time-history analysis is a step-by-step analysis of the dynamical response of a structure to a

specified loading that may vary with time. The analysis may be linear or non-linear.
The dynamic equilibrium equations to be solved are given by
(t )} + [C] {u (t )} + [K] {u (t )} = {r(t)}; {r(t)} being the response vector relative to
[M] {u
ground motion. [56]
Wilson-θ method is one of such solution techniques. This method, developed by E. L. Wilson, is
a modification of the conditionally stable linear acceleration method that makes it
unconditionally stable. This modification is based on the assumption that the acceleration varies
linearly over an extended time step δt= θ.Δt. The accuracy and stability properties of the method
depend on the value of the parameter θ, which is always greater than 1. If θ = 1, this method
reverts to the linear acceleration method, which is stable if Δt < 0.551 TN, where TN is the shortest
natural period of the system. If θ ≥ 1.37, Wilson’s method is unconditionally stable, making it
suitable for direct solution of the equation of motion; θ = 1.42 gives optimal accuracy. [55]
In dynamic analysis of structures and foundations, damping plays an important role. However,
due to the limitation in our knowledge about damping, the most effective way to treat damping
within modal analysis framework is to treat the damping value of a MDOF system as an
equivalent Rayleigh damping in the form:
[C] = α [M] + β [K]
in which [C], [M], [K] are respectively the damping-, mass- and stiffness- matrix of the physical
system; α and β are the predefined constants called Rayleigh damping coefficients.
As such, in most of the practical engineering analysis, the analyst makes the simplifying
assumptions in selecting damping ratios (constant for all significant modes). It is a fact that
modal mass participation decreases with increase in modes. Based on above, one can infer that,
as mass participation decreases with higher modes, the frequency increases and it is indeed an
observed phenomenon. With reduction in modal mass for successive modes, critical damping
will decrease with increase in mode. Overall damping of a system being a constant (since total
mass and stiffness are constant for a system), the damping ratio will increase with increasing
modes. To incorporate this reality, one option may be to compute the Rayleigh coefficients using
the following basic formulation incorporating range values for the first m significant modes: [57]
2ξ 1ω1 − 2ξ m ω m
β= ; α = ( 2ξ1ω1 – β ω12 )
ω1 − ω m
2 2

The modal damping ratio ξ for the soil system would typically be much higher than other civil
engineering structures, say 15 to 20 % .[55]
In this research, the damping variation, as commonly prevails for the geotechnical system, is
taken between 10 to 20 percent corresponding to the first and the last significant modes to
consider.
Now, in order to obtain the natural frequencies and mode shapes, the global mass matrix and the
global stiffness matrix need to be computed to solve the Eigen value problem. The Eigen values
yield natural frequencies and the Eigen vectors corresponding to the Eigen values give the mode
shapes. This job was carried out by the inbuilt modal analysis facility available in SAP2000. The
first twelve modes, by default in the program, were considered as they contributed to more than
99% of modal mass participation. Even to achieve the 90 % mass participation, as generally
86 Journal of the Institute of the Engineering

recommended by design codes, it was found that at least 10 modes are needed for BMP dam
section. Table 4 shows the range of natural periods between the first and twelfth significant
modes of the RF and FF models for the BMP dam.
Table 4: Time period of the significant modes first and twelfth significant vibration modes
Mode no. For RF model For FF model
1 0.8773 1.1957
12 0.4244 0.4981
Figure 14 shows the first and twelfth significant mode shapes of the BMP dam with the FF
model, for example. In both models, it is found that the modes are basically constituted by core
deformations.

Mode 1

Mode 12

Fig. 14: 1st and 12th vibration modes of FF model

The damping coefficients thus computed are presented in Table 5.
Table 5: Rayleigh damping coefficients
Model α β
RF 0.060938 0.026737
FF 0.212111 0.030379
Finally, the RF and FF models were each subjected to the linear time history analysis with the El
Centro Earthquake ground motion as the prescribed earthquake and using the Wilson-θ
incremental algorithm, inbuilt in SAP-2000 program, with θ = 1.42.
The computed time histories of the peak acceleration at the crest of the dam for the models are
Fig. 15 and Fig. 16 respectively.
The peak horizontal crest accelerations of the history as obtained from the time history analysis
are presented in Table 6. The acceleration response amplification is almost double for RF model,
in which case the response acceleration is nearly equal to “g”.
 Finite Element Modeling and Decoupled Seismic Stability . . . 87

Fig. 15: Horizontal crest acceleration response history for RF model

Fig. 16: Horizontal crest acceleration history of FF model

Table 6: Peak horizontal crest acceleration response
Model definition Acceleration, in units of g
RF 0.96
FF 0.51
1.2 Calculation of Permanent Crest Settlement
As per the chart shown already in Fig. 2, the permanent crest settlement of the crest of a fill dam
is obtained as a function of the ratio of yield acceleration to the peak crest acceleration during
earthquake excitation. The displacement thus obtained may be assumed equivalent to be the
seismic crest settlement of the dam [14].
Now, the value of yield acceleration was assumed to be 0.15 g based on the suggestion by Seed
(1979) that a seismic coefficient of 0.1g for magnitude 6.5 earthquakes and 0.15g for magnitude
8.25 earthquakes may be taken to obtain a factor of safety of 1.15 [15, 58].
Having known the values of yield acceleration, and already evaluated peak crest response
accelerations from time history analyses, the permanent settlement of the zoned rockfill
embankment could be determined. Table 7 shows the final results.
88 Journal of the Institute of the Engineering

Table 7: Permanent displacement results

Model definition Permanent crest settlement (m)
RF 1.20
FF 0.56

It is found that the potential permanent displacement (settlement) of the crest predicted by RF
model is higher than that by FF model. It is unlike the case of static analysis, as mentioned
somewhere before, that the FF model is more conservative.
Further the value exceeds even 1m for RF model, the model which is conventionally adopted in
design. The model appears conservative in seismic analysis and design. Nevertheless, again
conventionally speaking, some design guidelines of embankment dams such as Indian standards,
stipulate 1 m as the limit of permanent settlement. This seems to conclude that the BMP dam
designed by the traditional empirical methods and method of limit equilibrium is insufficient
when evaluated dynamically in terms of permanent crest settlement; for example, by the so
worked out decoupled seismic analysis method.

6. Concluding Remarks
In this paper, the seismic safety of a zoned rockfill dam designed by traditional empirical
methods was evaluated in terms of permanent settlement by a more detailed but yet simplified
method of ‘decoupled seismic stability analysis’. BMP dam was taken as a problem case.
Two types of finite element models were prepared and investigated: one with rigid foundation
assumption (RF model) and other incorporating some effective geometry of foundation into the
model to include foundation flexibility (FF model). The effective geometry of foundation was
worked out from the first principle by trial and error ‘vertical stress criterion’. A large domain
with twice the base width on either side and four times the base width underneath was obtained
as effective foundation to achieve the negligible stresses (below 10%) at the boundaries. But it
can be reduced little bit to allow for a higher limit such as 20% in the stress criterion, for
practical design purposes, as recommended by scientists specializing soil mechanics such as
Terzaghi. This would reduce the degrees of freedom and save some computational cost.
Few secondary studies were made on gravity turn-on foundation stress analysis and embankment
deformation analysis while validating the models prepared on the general purpose platform of
SAP 2000. Also, few important observations on vibration modes and vibration period of zoned
rockfill dam structure were made while computing the Rayleigh damping coefficients.
Having seen the result of BMP dam evaluated as a case in this research, it appears that the zoned
rockfill dams designed by the traditional empirical methods may not be safe against the
prescribed design earthquake if evaluated by rigorous methods. So application of dynamic
analysis and rigorous tools is encouraged for the design and construction of high and important
dams in the high seismicity region like Nepal, as already recommended by ICOLD few decades
ago.
Further, it is always advised to be conservative in designing the freeboard of fill dams, with due
consideration of the potential settlement during an earthquake.
In fact, the displacement-based approach provides a compromise between the rather inadequate
pseudo-static approach and the more refined numerical analyses; it has indeed the advantage of
 Finite Element Modeling and Decoupled Seismic Stability . . . 89

giving a quantitative assessment of the earthquake-induced displacement using a rather simple

analytical procedure. The authors do encourage the introduction of research and practice of more
rigorous methods which can give yet more accurate and detailed results on seismic settlements.
For important structures like dams, every possible effort should be made to ensure its safety. The
potential downstream disaster due to breaching of dam is always very fatal. Therefore, although
the authors’ present research emphasized more on seismic displacement analysis with respect to
loss of freeboard, due importance shall be given in future for evaluation of stresses and strains
leading to the cracking conditions in rockfill dam cores, such as [59].
Further, the research on advancement or innovation in material and structural composition to
improve the performance of the zoned rockfill dam are also encouraged, such as [60,61].
In addition, the authors have focused themselves on the dynamic analysis of the dam structure
assuming that the seepage related designs of the dam are under control and assuming no seismic
liquefaction of foundation would take place. However, they do wish the numerical simulation
regarding stability analysis against seepage would also be introduced amongst the design
researchers of dams in Nepal.
Finally the studies and research on comparing the seismic safety evaluation of zoned rockfill dam
structures via the rigorous effective stress nonlinear dynamic analyses, displacement-based block
analysis and through forced-based pseudo-static methods, such as [62], may be the scope for
future research in Nepal.

Acknowledgements
The first author would like to thank Dr. Sunil Kumar Lama of the Institute of Engineering,
Tribhuvan University, Nepal, his undergraduate project supervisor, for inspiring him to enter into
the world of dam and hydropower structures. He is also thankful to Professor Nelson Lam of the
University of Melbourne, Australia, for his repeated encouragement to publish the refereed
journal papers for the subsequent international research career. The first author also would like to
thank IIT-K, India for giving access to their library.

REFERENCES
[1] Water Encyclopedia (Retrieved in November 2010). http://www.waterencyclopedia.com/Da-
En/Dams.html
[2] Foster, M. A., Fell, R. and Spannagle, M. (1998). Analysis of Embankment Dam Incidents,
UNICIV
Report No. R374, The University of New South Wales, Australia
[3] ICOLD (2001). Guidelines on Design Features of Dams to Resist Seismic Ground Motion, Bulletin
120, International Commission on Large Dams
[4] Swaisgood, J. R. (2003). Embankment Dam Deformations Caused by Earthquakes, Proceedings of
the 2003 Pacific Conference on Earthquake Engineering, Christchurch NZ, Paper no. 14
[5] Cetin, H., Laman, M. and Ertunc, A. (2000). Settlement and Slaking Problems in the World’s
Fourth Largest Rock-fill Dam, the Ataturk Dam in Turkey, Engineering Geology, 56, 225-242
[6] Gikas, V. and Sakellariou, M. (2008). Settlement Analysis of the Mornos Earth Dam (Greece):
Evidence from Numerical Modeling and Geodetic Monitoring, Engineering Structures, 30, 3074-
3081
90 Journal of the Institute of the Engineering

[7] Meskouris, K., Konke, C. and Chudoba, R. (1999). Seismic Safety of Rockfill Dams, European
Conference on Computational Mechanics, München, Germany
[8] Hammouri, N. A., Malkawi, A. I. H. and Yamin, M. M. A. (2008). Stability Analysis of Slopes
Using the Finite Element Method and Limiting Equilibrium Approach, Bull Eng Geol Environ, 67,
471–478
[9] Japanese Society of Civil Engineers (2000). Dynamic Analysis and Earthquake Resistant Design,
Vol. 3: Dynamic Analysis of Special Structures, A. A. Balkema, Rotterdam, Netherlands, 2001
[10] Kramer, S.L. (1995). Geotechnical Earthquake Engineering, Prentice Hall, Upper Saddle River,
New Jersey
[11] Newmark, N. M. (1965). Effects of Earthquakes on Dams and Embankments, Geotechnique, 15,
139-160
[12] Wilson, R. C. and Keefer, D. K. (1983). Dynamic Analysis of a Slope Failure from the 6 August,
1979 Coyote Lake, California Earthquake, Bulletin of the Seismological Society of America, 73 (3),
863-877
[13] Jibson, R. W., Harp, E. L. and Michael, J. A. (1998). A Method for Producing Digital Probabilistic
Seismic Landslide Hazard Maps: An Example from the Los Angeles, California Area, Open-File
Report 98-113, USGS. http://pubs.usgs.gov/of/1998/ofr-98-113/ofr98-113.html
[14] Hynes-Griffin, M. E. and Franklin, A. G. (1984). Rationalizing the Seismic Coefficient Method,
Miscellaneous Paper GL-84-13, US Army Corps of Engineers, Waterways Experiment Station,
Vicksburg, Mississippi
[15] IIT-K, GSDMA (2005). Guidelines for Seismic Design of Earth Dams and Embankments, IIT-K
GSDMA
[16] Baldovin, E. and Paoliani, P. (1994). Dynamic Analysis of Embankment Dams, The Numerical
Analysis and Modeling of Soil Structure Interaction (Editor: J.W. Bull), E & FN Spon, London
[17] ANCOLD (1998). Guidelines for Design of Dams for Earthquakes, Australian National Committee
on Large Dams, Melbourne, Australia.
[18] Lacy, S. J. and Prevost, J. H. (1987). Nonlinear Seismic Response Analysis of Earth Dams, Soil
Dynamics and Earthquake Engineering, 6(1), 48-63
[19] Mircevska, V. and Bickovski, V. (1998). Two Dimensional Nonlinear Dynamic Analysis of
Rockfill Dam, Proceeding of the International Symposium on New Trends and Guidelines on Dam
Safety, Barcelona, Spain
[20] Uddin, N. (1997). A Single Step Procedure for Estimating Seismically-induced Displacements in
Earth Structures, Computers and Structures, 64(5-6), 1175-1182
[21] Cascone, E. and Rampello, S. (2003). Decoupled Seismic Analysis of an Earth Dam, Soil Dynamics
and Earthquake Engineering, 23, 349-365
[22] Jafarzadeh, F. and Yaghubi, A. (1998). Three Dimensional Dynamic Behaviour of Zoned Rockfill
Dams with Emphasis on a Real Case Study, Proceeding of the International Symposium on New
Trends and Guidelines on Dam Safety, Barcelona, Spain.
[23] Yu, Y., Xie, L. and Zhang, B. (2005). Stability of Earth-rockfill Dams: Influence of Geometry on
the Three-dimensional Effect, Computers and Geotechnics, 32, 326-339
[24] Singh, B. and Varshney, R. S. (1994). Engineering for Embankment Dams, Oxford and IBH
Publishing Co. Pvt. Ltd., New Delhi, India
[25] Jafarzadeh, F. and Soleimanbeigi A. (2007). Effect of Geometry and Foundation Conditions on the
Accuracy of the Steady State Seepage Analysis Results for Rockfill Dams, Scientia Iranica, 14 (3),
212-220
 Finite Element Modeling and Decoupled Seismic Stability . . . 91

[26] Clough, R. W. and Woodward, R. J. (1967). Analysis of Embankment Stresses and Deformations,
Journal of Soil Mechanics and Foundation Division, ASCE, 93 (SM4), 529-49
[27] Gui, M.-W. and Chiu, H.-T. (2009). Seismic Response of Renyitan Earth-fill Dam, Journal of
GeoEngineering, 4 (2), 41-50
[28] Swoboda, G. and Lei, X. Y. (1994). Simulation of Arch Dam–Foundation Interaction with a New
Friction Interface Element, International Journal for Numerical and Analytical Methods in
Geomechanics, 18(9), 601–617
[29] Nepal Electricity Authority (xxxx). Upper Karnali Hydroelectric Project, Chapter 7: Seismology
[30] Nepal Electricity Authority (xxxx) West-Seti Hydroelectric Project, Local Seismicity Study Report
[31] Maskey, P. N. (2005). Generation of Site Dependent Earthquake Ground Motion Parameters, The
Nepalese Journal of Engineering, 1 (1), 84-91
[32] Consulting Engineers, Nippon Koei Co. Ltd., Tokyo (1983). Kulekhani Hydroelectric Project,
Project Completion Report, HMG/Nepal, Kulekhani Hydroelectric Development Board
[33] Hydro-Engineering Services (P) Ltd. (1999). Final Report on Feasibility Study Level Analysis
(Updating Parameters) of Bagmati Multipurpose Project, A report submitted to WECS, HMG/Nepal
[34] Dhakal, S. and colleagues (2005). Analysis and Design of Storage Type Hydro-Project (BMP), B.E.
Final Year Project Report, Department of Civil Engineering, Pulchowk Campus, Institute of
Engineering (IOE), Tribhuvan University (TU), Nepal
[35] Lin, J. S. and Whitman R. V. (1983). Decoupling Approximation to the Evaluation of Earthquake-
induced Plastic Slip in Earth Dams, Earthquake Engineering and Structural Dynamics, 11, 667–78
[36] Gazetas, G. and Uddin, N. (1994). Permanent Deformation on Preexisting Sliding Surfaces in
Dams, Journal of Geotechnical Engineering, ASCE 120(11), 2041–2061
[37] Rathje, E. M. and Bray, J. D. (1999). An Examination of Simplified Earthquake-induce
Displacement Procedures for Earth Structures, Canadian Geotechnical Journal , 36, 72–87
[38] Makdisi, F. I. and Seed H. B. (1978). Simplified Procedure for Estimating Dam and Embankment
Earthquake-induced Deformations, Journal of Geotechnical Engineering, ASCE 1978, 104(7), 849–
867
[39] Nejad, B. G., Soden, P., Taiebat, H. and Murphy. S. (2010). Seismic Deformation Analysis of a
Rockfill Dam with a Bituminous Concrete Core, IOP Conference Series, Materials Science and
Engineering, 10, 012106
[40] Singh, R. Roy, D. and Jain, S. K. (2005). Analysis of Earth Dams Affected by the 2001 Bhuj
Earthquake, Engineering Geology, xx, xxx–xxx (Article in press)
[41] Gazetas, G. and Dakoulas, P. (1992). Seismic Analysis and Design of Rockfill Dams: state-of-the-
art. Soil Dynamics and Earthquake Engineering, 11, 27–61
[42] Chopra, A. K. and Zhang, L. (1991). Earthquake-induced Base Sliding of Concrete Gravity Dams,
Journal of Structural Engineering, ASCE, 117(12), 3698–3719
[43] Elgamal, A. W. M., Scott, R. F., Succarieh, M. F. and Liping, Y. (1990). La Villita Dam Response
during Five Earthquakes Including Permanent Deformations, Journal of Geotechnical Engineering,
ASCE, 116(10), 1443–1462
[44] Bardet, J.P. and Davis, C. A. (1994). Performance of San Fernanado Dams during 1994 Northridge
Earthquake, Journal of Geotechnical Engineering, ASCE, 120(7), 554–564
[45] Ozkan, M. Y., Erdik, M., Tuncer, M. A. and Yilmaz, C. (1986). An Evaluation of Surgu Dam
Response during 5 May 1986 Earthquake, Soil Dynamics and Earthquake Engineering, 15, 1–10
[46] Dakoulas, P. and Gazetas, G. (1986). Seismic Shear Strains and Seismic Coefficients in Dams and
Embankment, Soil Dynamics and Earthquake Engineering, 5(2), 75-83.
92 Journal of the Institute of the Engineering

[47] Murthy, V. N. S. (1993). A Text Book of Soil Mechanics and Foundation Engineering, UBSPD
[48] Nepal Electricity Authority (1980). Bagmati Multipurpose Project, Phase-I Feasibility Report,
Annex 18: Bagmati High Dam and Hydroelectric Development
[49] Central Material Testing Laboratory (2002). Final Report of Soil Investigation Work for Proposed
Lalitpur Bishalbazar Commercial Complex, Department of Civil Engineering, Pulchowk Campus,
IOE, TU, Nepal
[50] Kutzner, C. (1997). Earth and Rockfill Dams - Principles of Design and Construction. Oxford and
IBH Publishing Co. Pvt. Ltd., New Delhi, India
[51] Goodman, R.E. (1989). Introduction to Rock Mechanics. John Wiley and Sons
[52] UNDP, Central Soil and Materials Research Station (1986). Workshop on Seismic Analysis and
Design of Earth and Rockfill Dams, Volume 1 and Volume 2, Central Board of Irrigation and
Power, New Delhi
[53] Bureau of Indian Standards (2002). IS 1893-2002: Indian Standard, Criteria for earthquake
Resistant Design of Structures, PART 1: General Provision and Buildings, Manak Bhavan, 9
Bahadur Shah Zafar Marg, New Delhi 110002
[54] Das, B. M. (2004). Principles of Geotechnical Engineering, Thomson Books/ Cole
[55] Chopra, A. K. (2005). Dynamics of Structures, Pearson Prentice Hall, India
[56] Computers and Structures Inc. (2005). SAP 2000 Analysis Reference Manual, Berkeley, California,
USA
[57] Chowdhary, I. and Dasgupta, S. P. (2003). Computation of Rayleigh Damping Coefficients for
Large Systems, Department of Civil Engineering, Indian Institute of Technology, Kharagpur, India
http://www.ejge.com/2003/Ppr0318/Ppr0318.pdf
[58] Seed, H. B. (1979). Considerations in the Earthquake Resistant Design of Earth and Rockfill Dams,
Geotechnique, 29, 215-26
[59] Osuji, S. O. and Anyata, B. U. (2007). Susceptibility of Clay Core to Cracks in Rockfill Dams by
Finite Element Modeling, Advanced Material Research (Volume 18-19): Advances in Materials and
System Technologies
[60] Noorzad, R. and Omidvar, M. (2010). Seismic Displacement Analysis of Embankment Dams with
Reinforced Cohesive Shell, Soil Dynamics and Earthquake Engineering, 30, 1149–1157
[61] Bagherzadeh Kh. A., Mirghasemi, A. A. and Mohammadi, S. (2011). Numerical Simulation of
Particle Breakage of Angular Particles Using Combined DEM and FEM, Powder Technology, 205,
15–29
[62] Rampello, S. and Silvestri, F. (2009). Forced Based Pseudo-Static Methods versus Displacement
Based Methods for Slope Stability Analysis (Editor: E. Cosenza), Eurocode 8 Perspectives from the
Italian Standpoint Workshop, Doppiavoce, Napoli, Italy, 249-262