Seismic performance of arch dams on non-homogeneous and discontinuous foundations (a case study: Karun 4 Dam)

The present study set out to investigate the nonlinear seismic response of the dam–reservoir–rock foundation system, taking into consideration the effects of change in the material properties of discontinuous foundation. To this end, it is important to provide the proper modeling of truncated boundary conditions at the far-end of rock foundation and reservoir fluid domain and to correctly apply the in situ stresses for rock foundation. The nonlinear seismic response of an arch dam mainly depends on the opening and sliding of the dam body’s contraction joints and foundation discontinuities, failure of the jointed rock and concrete materials, etc. In this paper, a time domain dynamic analysis of the 3D dam–reservoir–foundation interaction problem was performed by developing a nonlinear Finite Element program. The results of the analysis of Karun-4 Dam revealed the essential role of modeling discontinuities and boundary conditions of rock foundation under seismic excitation.


Introduction
A large number of high arch dams in the world are built or to be built in the seismically active areas. Therefore, it is essential to carry out the seismic analysis of coupled dam-foundation system to ensure the safety and reliability of high arch dam structures. The structural strength of an arch dam under ground motion is mostly dependent on the stability and strength of its abutments. Actually, even high safety margins for unexpected ground motions do not guarantee the stability of dam if it is constructed on an uncertain foundation. In addition, due to the complex nature of rock foundation including non-homogeneous materials, the existence of joint sets and faults, etc., final judgment about the performance of dam will be difficult. Collapse of Malpasset Dam in France in 1959 is an apparent example of the lack of foundation strength.
During the past few years, extensive research has been done in the analysis and design of concrete dams, but the need for more accurate modeling of abutments in a coupling system with regard to the effects of mass, flexibility and non-homogeneity of discontinuous rock foundation are still felt. In the current study, a numerical program for 3D nonlinear dynamic analysis of concrete arch dams is developed in FORTRAN language and applies it to Karun 4 Dam. The study attempts to control two main features which directly affect the accuracy and precision of analytical results, i.e., correct in situ stresses for the bedrock and the choice of a proper boundary condition for the farend of discontinuous rock foundation.

A review of research areas and solving methods
The main factors that influence significantly the three-dimensional nonlinear dynamic analysis of arch dams are identified: (1) dam-reservoir interaction and resulted distribution of hydrodynamic pressure, (2) reservoir-foundation interaction and the related effects of reservoir bottom sediments, (3) dam-foundation interaction and the role of non-homogeneity and discontinuities in bedrock, (4) nonuniform input of the free-field motions, (5) nonlinear behavior of the quasi-brittle material of concrete and jointed rock, and contact in the contraction and peripheral joints of dam body, and (6) boundary conditions. Great effort has been made to develop the fundamentals and analytical methods of the above-mentioned subjects and a brief review of the main issues related to research are presented here.
Fluid-solid interaction is a very complicated problem which involves both structural and fluid dynamics. Several researchers have developed advanced numerical methods based on the finite element, boundary element or combination of them to model the dynamic dam-reservoir interactions. Two common approaches in fluid domain are Eulerian-and Lagrangian-based formulations (Bouaanani and Lu 2009). The former approach is known as pressureor potential-based formulations where fluid pressure or velocity potential is selected as state variable. The Lagrangian approach in fluid domain is an extension of the solid finite element formulation with nodal displacements as degrees of freedom. As a result, the fluid domain is formulated using the same shape functions as structural elements and, in this way, compatibility at the fluidstructure interface is automatically fulfilled. In such cases, fluid elements are characterized by a volumetric elastic modulus equal to the fluid bulk modulus (or fluid compressibility); with a negligible shear resistance and Poisson constant being nearly 0.5 to simulate the fluid flow more reliable. The major disadvantage of Lagrangian approach is the generation of spurious circulation modes, due to the zero shear modulus. There are different methods for eliminating these zero energy modes that have been used by a large number of researchers (Chopra et al. 1969;Shugar and Katona 1975;Hamdi et al. 1978;Zienkiewicz and Bettess 1978;Akkaş et al. 1979;Wilson and Khalvati 1983;Olson and Bathe 1983;Dogangün et al. 1996;Dogangün and Livaoglu 2006, etc.). One approach is based on assuming that the shear modulus of the fluid is numerically very small. This approach has been adopted by this study. The partial absorption of pressure waves at sediment layers of reservoir bottom and its lateral sides may significantly affect the magnitude of hydrodynamic forces while the response of dam due to the ground motion is being investigated (Mirzabozorg et al. 2003). In the present study, the Lagrangian approach is used for modeling the fluid and sediment domain. Also, the interface elements with low shear stiffness are used to model the common surfaces of the fluid and solid elements.
Dam-foundation interaction is mainly related to the bedrock's flexibility, changes of physical properties of rock foundation, existence of joints and faults in it, geometry of dam body, water leakage, uplift pressure, etc. It has been proven that in a nonlinear dynamic analysis including dam-foundation interaction, the foundation's mass, flexibility and radiation damping are important (Tan and Chopra 1995). In addition, in order to model the behavior of jointed rock masses, their strength and deformability should be expressed as a function of joint orientation, joint size, and joint frequency. Moreover, it is not possible to represent every joint individually in a constitutive model. Therefore, it is necessary to use simple techniques such as equivalent continuum method which can capture reasonably the behavior of jointed rock mass. The finite element method developed in the present study applies for the foundation rock model both: nonlinear solid element for modeling the jointed rock as an equivalent continuum whose properties represent material properties of the jointed rock, and nonlinear interface element to account for the surface roughness of discontinuities.
The best way to determine the size of foundation model is based on the ratio of deformation modulus of foundation (E f ) to the elastic modulus of concrete dam (E c ). For a flexible rock foundation with E f /E c less than , the foundation model should extend at least twice the dam height in all directions (Federal Energy Regulatory Commission Division of Dam Safety and Inspections, Washington DC 1999).
The definition of suitable boundary conditions related to its surrounding domain is another important part of modeling procedure. In the present study, the governing equations and related boundary conditions consist of water free surface and far-end truncated boundaries of the reservoir and rock foundation. By neglecting the effects of gravity waves, a zero-pressure boundary condition is prescribed at the horizontal free surface water (assuming negligible surface tension). This simplification is a common assumption in the analysis of concrete dams, particularly for deep reservoirs. Also, several studies have been carried out into the improvement of the boundary condition at farend of the reservoir and rock foundation in the dynamic analysis of coupling system. A radiation condition (such as, Sommerfeld's and Sharan's boundary conditions and Küçükarslan 2004) can be applied at the truncated boundary of the reservoir. A similar boundary for eliminating waves propagating outward from the foundation domain is Lysmer and Kuhlemeyer (1969) boundary. Boundary element method has also been widely used (Beer et al. 2008;Brebbia and Dominguez 1992). In this method, only the boundaries of the unbounded medium are discretized so that the spatial dimension is reduced by one and the radiation condition is satisfied automatically as a part of the fundamental solution (Estorff and Kausel 1989). Many researchers have successfully used infinite boundary elements to model wave propagation problems on far-end boundary (Valliappan and Zhao 1992, etc.). In this study, for the truncated boundary conditions of reservoir and rock foundation, interface elements have been used. Using these boundary conditions does not prevent the sliding at foundation discontinuities due to seismic loading, thus providing non-uniform excitation.
It has been recognized that for infrastructures with extended foundations such as concrete arch dams, the ground motion is non-uniform along the canyon due to the wave passage effects, coherency effects and site response effects (Der kiureghian 1996). Several researchers have studied on the response of concrete dams when the system was excited non-uniformly (Bayraktar et al. 1996;Mirzabozorg et al. 2007, etc.). Some studies indicate that the non-uniform seismic input can have an important role in the dam response (Alves and Hall 2006).
In addition, the nonlinear material properties of concrete dam and bedrock as well as the nonlinear effects of contraction joints of dam body and discontinuities in the dam foundation are modeled by appropriate methods which will be described later (Ahmadi and Razavi 1992;Ahmadi et al. 2001;Mojtahedi and Fenves 2000).

Description of the analysis procedure
The schematic view of dam-reservoir-foundation system is shown in Fig. 1. The full system is modeled by an assemblage of solid and interface elements. The isoparametric 8-node solid elements with 2 9 2 9 2 Gauss integration are used for all domains including dam body, reservoir, sediment and rock abutments. Also, eight-node interface elements are used in common surfaces of domains, such as the contraction and perimetral joints in dam body and discontinuities of rock abutments. The finite element formulations support both geometric and material nonlinearity. Assuming that non-linearities are limited to the concrete dam, rock blocks, contraction joints of dam body and rock discontinuities, the stiffness of these elements needs to be updated at each iteration. Interface elements are placed between continuum (solid) elements, as shown in Fig. 2. A summary flow chart of the finite element program is shown in Fig. 3. In the programing, the capabilities of several open source programs that are developed for the analysis of concrete dams such as ADAP-88, EAGD-SLIDE, EACD-3D-96 (Mojtahedi et al. 1992;Chavez and Fenves 1994;Tan and Chopra 1996) were investigated and useful subroutines and subprograms of them were used (Smith and Griffiths 2004;Bathe 1996).

Nonlinear dynamic analysis
The governing equations of motion for 3D nonlinear dynamic analysis of coupling system (subjected to earthquake loads) were discretized using Newmark's method. By adopting very large time increments, static nonlinear analysis can be accomplished as a special case of dynamic analysis. The discretized nonlinear dynamic equation of   (2017) 9: 191-203 193 motion is given by (Bathe 1996;Zienkiewicz and Taylor 2005): where [M] is the mass matrix, [C] is the damping matrix, [K] is the stiffness matrix and {R} is the nodal external forces.
and U f g are the acceleration, velocity and displacement vectors, respectively, at the (n ? 1)th time step. Also, [K T ] is the tangent stiffness matrix and [C T ] is the updated damping matrix which changes at the same time as the stiffness reduces. The dynamic equilibrium of the system at time ''n ? 1'' can be written in terms of the unknown node displacements U n?1 by substitution of (Wilson 2002): where the constants b 1 to b 6 are defined: The damping matrix is determined based on the wellknown Rayleigh damping: The parameters a and b are pre-defined constants and can be evaluated by the solution of a pair of simultaneous equations if the two damping ratios (n i ) associated with two specific frequencies ðx i Þ are known: The full Newton-Raphson iteration scheme can be used to resolve the nonlinearity. The parameters b and c determine the stability and accuracy characteristics of the algorithm. The constant acceleration method is obtained when b = 1/4 and c = 1/2.

Nonlinear models for foundation rock and Dam body
Combinations of Mohr-Coulomb (1882-1900) yield function with a tension cut-off (i.e., The Modified Mohr-Coulomb model suggested by Paul in 1961) are used for both concrete and jointed rock materials. Crook et al. (2003) presented that the Modified Mohr-Coulomb model is able to effectively model both brittle-tensile or axial splitting fractures and shear features model. The Mohr-Coulomb criterion is based on Coulomb's equation (1773). If r 11 ! r 22 ! r 33 are the principal stresses, we can write this criterion as (Lubliner 1990): where / and c are the internal friction angle and cohesion, respectively. The 3D failure surface of the Mohr-Coulomb criterion can be expressed in terms of stress invariants where the lode angle is: For the tension cutoff yield function (Rankine crack model), we have where f 0 t is tension strength. This criterion can be fully described by the following equation: Reservoir modeling and simulation The finite element formulation for fluid is based on Lagrangian approach in which the fluid strains are calculated from the linear strain-displacement equations. The pressure volume relationship for a linear fluid is expressed by where p, k and ev are pressure that is equal to the mean stress, the bulk modulus and the volumetric strains of the fluid, respectively. The volumetric strain can be expressed in terms of Cartesian displacement components as follows: where ux, uy and uz are displacement components related to axes x, y and z, respectively.

Nonlinear interface element
Different models have been developed to represent the contact between two surfaces. One is based on the cohesive law which can be expressed in a way that the local traction ðtÞ across the interface is taken as a function of displacement jump ðdÞ across the cohesive surfaces. A formula worked out by Ortiz and Pandolfi (1999) accounts for the free energy density per unit undeformed area (H) so that the traction acting on the interface: where d c denotes the value of d (Needleman 1990) at peak traction ðt max ¼ r c Þ: ; ðe ¼ expð1ÞÞ ð15Þ Figure 4 illustrates the curve of loading and unloading responses of interface elements. When the contact surfaces undergo a combination of shear deformation and normal compression, the effective separation ðdÞ is computed only from the shear components. Also, under normal compression the cohesive material behaves like a linear spring. The weighting coefficient ðbÞ defines the ratio between critical shear and normal tractions (Ruiz et al. 2001).
The value of interface stiffness depends on the roughness of contact surfaces, as well as the relevant properties of filling material and its moisture. For an initially closed interface, the normal stiffness K n and the tangential stiffness K S are assumed to have large values. These values can be estimated from the lowest Young's modulus and shear modulus of the adhesive domain around it, according to the following equations: where m i ði ¼ 1; 2Þ is a factor that controls the contact properties (only in penetration), E and G are the smaller elastic and shear modulus when considering the contact between two different materials, L e is the characteristic thickness of the adjacent solid element perpendicular to the interface, and A e is the surface area of the interface element.
For a proper numerical modeling of the far-end boundaries of massed foundation and reservoir, the viscous boundary condition must be applied to prevent wave reflection from the artificial boundary of infinite media in finite element analysis (Ghaemian et al. 2005).
Combinations of viscous and spring boundaries were used as the basis for interface element to find a simple and efficient model. Therefore, interface elements are placed around the problem boundaries with equivalent stiffness of infinite media and related damping to absorb the energy of outgoing waves in the normal and tangential directions. To calculate the damping coefficient C tb for the elements of far-end boundaries, the following equations are used: In which V P and V S are the primary and secondary wave propagation velocity within the foundation medium and are given as: where subscript f indicates the parameters pertinent to the rock foundation. The radiation damping derived from Eq. (17) is applied on the far-end boundaries of the foundation which are added to the global damping matrix of the structure, [C]. Similarly, linear viscous elements can be inserted at the upstream boundary of the reservoir that will allow the wave to pass and the strain energy in the water to ''radiate'' away from the dam.

Verification example
To verify the accuracy and validity of the finite element modeling and developed computer code, the tallest monolith with unit width of well-known Pine Flat Dam, a concrete gravity dam in California, which is 122 m high, is selected. A water depth of 116 m is considered as the full reservoir condition. The geometry and FE model of the Pine Flat dam monolith with unit width is shown in   (2017) 9: 191-203 195 concrete is taken to be 3.36 MPa which is about 12% of its compressive strength. The dynamic tensile strength must be equivalent to the direct tensile strength multiplied by a factor of 1.50 (Raphael 1984). The analysis results consist of the weight of the dam, the static pressure of impounded water and the seismic response for earthquake excitation of horizontal-x component of Taft-1952 Lincoln California ground motion with scaling to a PGA = 0.4 g. The applied proportional damping of the dam provides a critical damping ratio of 5% in the fundamental vibration mode of it. Figure 6 shows the comparison between the obtained results for crest displacement of the first 15 s excitation from the program developed for the present study and commercially available ANSYS program for the case of linear material behavior. In addition, Fig. 7 shows the analysis results obtained from developed program by considering the effects of material nonlinearity and base sliding. For transient structural analysis in ANSYS, a corresponding 8-node 3-D concrete element SOLID65 is selected.

Results of case study analysis
Karun-4 Dam is a double curvature arch dam on Karun River in the province of Chaharmahal-e Bakhtiari, Iran. Its whole crest length is divided by 20 contraction joints. The geometric characteristics of the dam shape are listed in Table 1 Table 2. These plates cause six large blocks, as shown in Fig. 10. The water in the reservoir was assumed to have a constant depth of 155 m, mass density = 10 kN/m 3 , bulk modulus = 2131 MPa and Poisson's ratio = 0.495 for nearly incompressible fluid. For the upstream/downstream sediments with the assumption of about 60/30 m depth: mass density q s = 13.6 kN/m 3 , and bulk modulus B s = q s c s 2 = 3312 MPa, where the sound speed profile is estimated through the physical sediment properties using Biot theory and assuming c s = 1560 m/s. The developed finite element mesh for the reservoir and sediments is shown in Fig. 11.
Interface elements are used for the modeling of rock discontinuities, vertical joints between the cantilevers and the intersection of the dam with canyon rock, as well as the boundary between reservoir and surrounding domain with negligible shear stiffness. At the truncated boundaries of the reservoir and rock foundation (later called as ''Moving B.C.''), the interface elements are available in the developed numerical program. The properties of several interface elements are presented in Table 3. The developed coupled model includes 11,764 nodes and 9348 elements. The position of interface elements is shown in Fig. 12.
The loads applied to the model consist of static and dynamic loading. Static loads are dead weight and the hydrostatic pressure of reservoir water at its normal level in addition to the sediment weight.
Information about the in situ stresses of the rock field is a fundamental parameter for the dam-foundation analysis and has a direct effect on dynamic design of such a coupled system. The in situ stress in a rock mass is simply equal to the weight of the overlying material; therefore, the discontinuities will control the magnitude and direction of this stress field. In this study, firstly the static load of discontinuous rock weight was applied to investigate the in situ stress. For this loading case, the dam body should remain free of stress because of canyon deformation. To overcome this problem, the numerical program has the ability to change the material properties in loading steps. Therefore, in a pre-loading step, Young's modulus of dam body and a region of the rock abutments near the dam gradually decrease and Poisson's ratio increases to 0.49. In the next dummy load step, material properties gradually change to real values.
The first 40 s of the three components of Taft Lincoln School records obtained from the 1952 earthquake happening in Kern County, California are used as input ground motion. The peak ground acceleration of x, y (horizontal components), and z (vertical component) directions are 0.156, 0.178 and 0.108g, respectively. For seismic hazard study of Karun 4 Dam site, the earthquake time histories are scaled to the maximum credible level at the middle height of canyon (PGA hor = 0.49g, PGA ver = 0.26g). A time step of 0.01 s is chosen for the analysis. The    Adv Struct Eng (2017) 9:191-203 197 displacement time histories for the three components of Taft earthquake are shown in Fig. 13.
To present the effect of foundation interaction on the seismic response, several cases of massive foundation are chosen: • C0 Continuous rock foundation-Moving B.C (without interface elements between the rock blocks and with interface elements on the truncated boundary); • C1 Continuous rock foundation-Fixed B.C (without interface elements between the rock blocks and on the truncated boundary); • C2 Discontinuous rock foundation-Moving B.C (called ''Real Case''); • C3 Condition ''C1'' with rigid and massless foundation; • C4 Condition ''C2'' with applying a reduction factor of 10% for deformation modulus and allowable bearing capacity of the rock blocks RB1, RB2, RB3 and RB4 (demonstrated in Fig. 10); • C5 Condition ''C2'' with applying a reduction factor of 10% for deformation modulus and allowable bearing capacity of rock blocks RB5 and RB6 (demonstrated in Fig. 10);   Fig. 9) in upstream-downstream direction for several cases. As can be seen, detailed modeling of foundation with high accuracy has a very important role in the coupled system analysis. Also, using the interface elements with appropriate characteristics on the far-end boundaries and major fault zones of the foundation changes the seismic response of dam significantly. It should also be noted that such boundary conditions and modeling of discontinuities in bedrock are critical for an actual response of dam as compared in Fig. 14. The foundation flexibility effects on dam response is also compared in Fig. 15.
The results of the analysis, as shown in Fig. 16, illustrate the importance of the effects of inhomogeneity and change in the material properties of discontinuous rock foundation on the seismic response of concrete dams. As can be seen, the foundation flexibility with a reduction factor of 10 and 20% for deformation modulus and allowable bearing   Contraction joints in the dam 3:0 Â 10 9 1:5 Â 10 9 Peripheral joints at the dam-foundation 4:0 Â 10 9 2:0 Â 10 9 Discontinuities in rock masses 1:0 Â 10 9 0:8 Â 10 9 Far-end boundaries of rock foundation 6:0 Â 10 9 3:0 Â 10 9 Fig. 15 191-203 199 capacity of domain parts has significantly affected the dam response. Figure 17 shows the difference between time history response of three adjacent nodes of the truncated boundary of foundation (belongs the three blocks RB1, RB2 and RB3). This is a result of non-uniform excitation of discontinuities foundation that has been composed of blocks with different mass, geometry and boundary conditions. To study the effect of contraction joint nonlinearity on the seismic response of dam, the time history comparison of joint sliding at the crown cantilever in crest and mid-height levels is considered as shown in Fig. 18. The maximum joint sliding is 10.1 cm at time 29.5 s for case C2. It should be noted that for the model with joint nonlinearity, there are two nodes at each point which are located on the interface elements. For example, nodes 64 and 53 at crest level and nodes 775 and 814 at the mid-height level are both on the crown cantilever at the same position (demonstrated in Fig. 9). Also, the contours of maximum and minimum principal stress obtained from nonlinear analysis for real case (Case C2) with a deformed scale of 20 at time 40 s are shown in Fig. 19a, b. A summary of the upstream-downstream and vertical displacements for all the analyses is provided for comparison in Table 4.

Conclusions
In this study, after conducting a survey of the literature on numerical methods for modeling of coupled concrete damreservoir water-foundation rock systems, a finite element computer program developed for the analysis and design of For the purpose of the study, the highest arch dam in Iran, i.e., ''Karun 4''is selected. To evaluate the seismic response of this arch dam more accurately, various effects of dam-foundation interaction in time domain including the effects of inertia and flexibility of non-homogeneous foundation rock as well as the proper exertion of in situ stresses have been taken into account.
As mentioned in the previous section, ten cases are considered: the linear model, and nine nonlinear models with flexibility, non-homogeneous and discontinuities  effects of foundation rock and hydrodynamic effects of reservoir on dam response. In each case, the time history response of crest displacement of the crown cantilever in upstream-downstream direction is obtained. The results demonstrate that the response of concrete arch dam-reservoir-foundation system is significantly affected by ''concrete and rock material nonlinearity'', ''the role of various discontinuities present in real system'', ''flexibility of non-homogeneous foundation'', and ''farend boundary condition''. Also, the results show that the application of material non-homogeneity in foundation for dynamic analysis is a key factor for seismic response, but neglecting the foundation discontinuities and moving and viscous properties of far-end boundary condition leads to inappropriate results.