Seismic response of earth dams considering dynamic properties of unsaturated zone

It is conventionally assumed in the analysis and design of earth dams that the soil located above the phreatic line, i.e. the uppermost seepage flow line, is completely dry. However, there is often an unsaturated flow of water through an unsaturated zone above this borderline and variation in moisture content in this zone results in variation of matric suction throughout this region. Variation of matric suction, in turn, results in variation of effective stresses in this zone. In this research, the seismic response of earth dams in terms of the displacement and acceleration at the crown of the dam as well as the stress distribution in the dam body is investigated. Taking into account the effect of unsaturated zone, a comparison is made to investigate the effect of conventional simplification in ignoring the dynamic characteristics of the unsaturated zone above the phreatic line and the more complicated analysis which includes the unsaturated zone. A function for the soil-water retention curve (SWRC) was assigned to the soil in the unsaturated zone to determine the variation of matric suction in this zone and analyses were made using finite difference software (FLAC). Results are then compared to the conventional method for homogeneous dams. In these analyzes the soil shear modulus was assumed to vary with the mean effective stress both for saturated and unsaturated zones. Among various results, it was notable that the history of crest x-displacement, and acceleration show higher values in models accounting for the unsaturated region. It was attributed to the considerably lower values of damping ratio in the crest region in the unsaturated models.


Introduction
For designing earth dams, analyzing the dynamic behavior of the dam against earthquakes is a major step.Common practice for analyzing the earth dams is to use the available numerical codes and commercial softwares.But these softwares mostly do not consider the unsaturated sate of soils.
Severe failures for dams subjected to earthquakes have been reported due to the various mechanisms such as liquefaction (e.g. the lower San Fernando Dam during 1971 San Fernando earthquake and damage of Mochikoshi Dam in 1978 Izu-Ohshim-Kinkar earthquake), slope failure and cracks leading to slope instability (i.e.Austrian Dam during 1989 Loma Prieta earthquake) and so on [1].Although the multiphase physics of unsaturated region in dams can be very important for such failure mechanisms, there are very few studies on the seismic behavior of dams taking into account the unsaturated region above the phreatic line.
The existence of more than one fluid phase in unsaturated soils results in the capillarity phenomenon.The behavior of unsaturated soils is, therefore, a function of capillary action and the state of matric suction.Zienkiewicz et al. (1990) were first to study the dynamic behavior of unsaturated soils [2].Their pioneering work has been followed by a couple of other studies in which different aspects of the unsaturated soils' response to the seismic excitations have been investigated [3][4][5][6][7].Different constitutive equations and assumptions in these studies have been considered.For instance, Muraleetharan and Wei and Ravichandran and Muraleetharan developed the governing equations for the seismic behavior of unsaturated soils and used a two stress state approach [3,4].Whereas, Uzuoka and Borja (2012), Shahbodagh-  2016) have used the effective stress approach and single stress framework [5][6][7].It has been stated that single stress approach can eliminate the need for multiple constitutive parameters for nonlinear response of unsaturated soils.On the other hand, in the effective stress approach, the definition of the effective stress parameter has faced a debate in the recent years for which multiple formulations have been introduced.Thus, in dynamic analysis of the response of unsaturated soils, various studies have employed different assumptions for the effective stress parameter.For instance, Uzuoka and Borja (2012), and Ghorbani et al. (2016) have used degree of saturation as the effective stress parameter whereas Shahbodagh-Khan et al. (2015) have considered it to be a function of matric suction which takes the hydraulic hysteresis into account.
One importance question when applying such constitutive equations to large scale engineering problems is that whether considering the unsaturated state of the soil would affect the results of the analysis or not and how influential such considerations would be.One may ask how important is to take the unsaturated region into account.The importance of considering the difference in saturated and unsaturated soil behavior for designing earth dams has been stated before in the previous studies by Alonso and Pinyol (2008) and Alonso and Cardoso (2010) [8,9].But the dynamic response of the earth dams considering the unsaturated state of the soil has not yet been comprehensively looked into.In particular the dynamic response of unsaturated earth dams with varying design parameters has not yet been examined in details.
Despite the slow but considerable and important progress in theoretical and numerical formulations for dynamics of unsaturated soils, the application of such models in engineering practice has faced less attention.There are only a few studies focusing on engineering projects such as slope and embankment.Mori et al. (2011) have modeled the dynamic behavior of a fill slope and have compared it to its monitored behavior and Matsumaru and Uzuoka (2015) have simulated the behavior of the shaking table tests on embankments [10][11].The focus of these studies has, however, been mainly on the effects of different assumptions for pore air pressure field on the simulation results.
In order to apply unsaturated models for analysis of earth dams, more studies are required and different characteristics of dams such as different material and geometric parameters have to be considered.
To investigate the influence of considering unsaturated sate in dynamic analysis of earth dams, in this study the results of a numerical study have been presented and discussed.Dams with different heights have been considered and the differences between the dams' response to seismic excitation with and without assumption of unsaturated region above the phreatic line have been examined.In the subsequent sections, the details of the model, the theoretical background of the model and the results have been presented.

Theoretical background
Numerical analysis were carried out using the finite difference program FLAC (Fast Lagrangian Analysis of Continua) which employs a continuum finite difference discretization based on the Langrangian approach.
In the analysis and design of earth dams, it is conventionally assumed that the soil located above the phreatic line, i.e. the uppermost seepage flow line, is completely dry.However, there is often an unsaturated flow of water through an unsaturated zone above this borderline and variation in moisture content in this zone results in the variation of matric suction throughout this region.Variation of matric suction, in turn, results in variation of effective stresses in this zone.
In order to consider dynamic properties of unsaturated zones in earth dams, the required constitutive equations were implemented in the numerical finite difference code FLAC2D using its built-in FISH language.The introduced functions were then assigned to the dam materials and seismic behavior of the dam was analysed in the plane strain condition.
Two different scenarios were used for modeling.At first similar to most of conventional studies, it was assumed that the uppermost seepage flow line, is completely dry (namely SAT-DRY condition) and in the second modeling scenario, unsaturated zone above this borderline was taken into account (SAT-UNSAT).The dynamic analysis was done in two steps.At first equivalent linear method was employed to specify shear modulus and damping ratio in each zone, then non-linear analysis was utilized to study the seismic behavior of the considered earth dams.
The Equivalent Linear method (E.L.) starts a dynamic analysis with the specified soil stiffness.In the first stage, FLAC processes the entire earthquake record and identifies the peak shear strains at each Gauss numerical integration point in each element.The shear modulus is then modified according to a specified G reduction function and the process is repeated.This iterative procedure continues until the required G modifications are within a specified tolerance.Soils subjected to dynamic stresses tend to 'soften' in response to cyclic shear strain.In the Equivalent Linear soil model, this softening is described as a ratio relative to maximum shear modulus.This is called a G-reduction function.The cyclic shear strain is obtained from the finite element analysis.The computed shear strain together with the specified maximum shear modulus are used to compute new shear modulus values in each iteration.
Ishibashi and Zhang (1993) developed an expression for estimating the G-reduction function [12].Where, the two main variables are PI (plasticity index) and confining pressure.Equation 1 describes the equation for G max which is used in G-reduction procedure, the involved parameter in this equation are described in equations 2 to 5.   is mean effective stress,  specifies the shear strain (0.65  max ).
The damping ratio,  , is a function of the cyclic shear strain, the same as the G-reduction function.Previous studies of Ishibashi and Zhang (1993) have led to the development of an expression [12], which can be used to estimate the Damping Ratio function.The variables in the expression are the plasticity index (PI), the G modulus reduction ratio G/G max and indirectly the confining pressure.The expression is as follows: In order to account for the presence of unsaturated parts in dam, the following relationship for the effective stress has been used: There are various formulations for the effective stress parameter.As discussed by Nikooee et al. ( 2013) and Pereira et al. (2010), degree of saturation is not a good candidate for the effective stress parameter, especially in fine grained soils [13,14].Recently Lu et al. (2010) suggested the use of effective degree of saturation for the effective stress parameter to account for the fact that the immobile water does not contribute to the mechanical behavior of unsaturated soils [15].Furthermore, Oh and Lu (2014) have shown that the use of effective degree of saturation results in a reasonable approximation of suction stress under different hydraulic and mechanical states [16].Therefore, in this study the effective stress parameter,  , was assumed to be effective degree of saturation.For the Gmax function, the required value of the mean stress was computed based on the equations ( 7) to (10).
Where earth pressure coefficient, 0 K , is defined based on Poisson's ratio,  .
Shear strength expressed as a function of effective normal stress, and the shear strength parameters, namely, the cohesion intercept ( c ) and the angle of shearing resistance ( ).

Initial condition:
In the numerical analyses, the initial condition of stress and pore pressure were considered to be a steady state condition at which earthquakes are more possible.The steady state condition was reached after a layer by layer construction analysis followed by an impounding analysis.
In fact, to have numerical stability, the element size (spatial element size), ΔL , has been suggested to be smaller than one-tenth to one-eighth of the wavelength associated with the highest frequency component of the input wave [17].In first step of dynamic analysis, the initial value of Rayleigh damping was considered to be 3%.Then by means of fish language and programming in FLAC in each step, the damping ratio was updated based on the shear modulus ratio and plasticity index (Equation 6).Furthermore in each step, the effective stress value for the unsaturated case was calculated and used to find the value of maximum shear modulus, G max .Then shear strain and G max were used to update the value of the shear modulus.For this purpose, programming was performed using fish language and the required subroutine were implemented in the FLAC.

Dam characteristics:
The heights of the simulated dams were considered as 35m, 50m and 70m, respectively.Their upstream and downstream slopes were invariably 1:3.5, and foundation was considered to be 12m deep.
As mentioned earlier, to investigate the seismic response of earth dams considering dynamic properties of unsaturated zone, dam modes were analyzed in two different cases.First the dams were modeled without considering dynamic properties of unsaturated zone, therefore above water surface, soil elements were dry.In the second case the dam models were investigated while the more complicated analysis which includes the unsaturated zone was implemented.This could help to understand the effect of conventional simplification in ignoring the dynamic characteristics of the unsaturated zone above the phreatic line.For this purpose, the above mentioned function for the soil-water retention curve (SWRC) was assigned to the soil elements in the unsaturated zone to determine the variation of matric suction in this zone.The modified G-reduction for unsaturated case which considers the values of skeleton stress in unsaturated cases was also implemented.

Results
Figure 2 presents the variation of the shear stress in two aforementioned scenarios and for dams of different heights.As it can be seen the shear stress in unsaturated case is higher due to the increase in the density of soil elements above water table which results in an increase in total stress.The increase in the total stress as well as the presence of the negative pore water pressure lead to the increase of the effective stress and accordingly an increase in shear stress values.
Water table (phreatic line) in dams with 35m, 50m and 70m is located in the elevation of 24m, 36m and 54m, respectively.Fig       The influence of suction and pore pressure on the soil shear modulus is noticeable.The shear modulus which is directly proportional to the effective stress and has a reverse relationship to the shear strain increases above water table in the unsat-sat models  The damping ratio for unsat-sat models above phreatic line is reduced by height.
 The history of crest x-displacement, and acceleration show higher values in unsat-sat models following the considerably lower values of damping ratio in the crest region.
Following two recommendations can be considered for future studies:  In order to be able to have a good insight on seismic behavior of dams considering unsaturated region therein, it is recommended to explore various design variables and parameters.The effect of foundation stiffness, cross sectional shape and different dam body material parameters can be considered for the future studies.
 Each study carries specific assumptions which can be released or investigated in future researches.The effect of different formulations for soil water retention curve as well as the effective stress parameter can be looked into in the future studies.
during the 1978 Tabas earthquake is used in all examples.Its duration is 30 sec and its maximum acceleration is about 0.8g.To minimize errors, this accelerogram has been adjusted according to the base line corrections.In addition this record was normalized to 0.3g and only the downstream-upstream component has been considered.Fig.1shows the corrected accelerogram.It has been applied uniformly along the base of the foundation.

Fig 4 .
Figure2presents the variation of the shear stress in two aforementioned scenarios and for dams of different heights.As it can be seen the shear stress in unsaturated case is higher due to the increase in the density of soil elements above water table which results in an increase in total stress.The increase in the total stress as well as the presence of the negative pore water pressure lead to the increase of the effective stress and accordingly an increase in shear stress values.Water table (phreatic line) in dams with 35m, 50m and 70m is located in the elevation of 24m, 36m and 54m, respectively.Fig 3.depicts the shear modulus profiles of the earth dams with height.Under the ground water table in both cases, shear modulus decrease with increasing elevation but above water table (phreatic line), shear modulus increases in saturated-unsaturated scenario.The reason is that shear modulus has a direct relation with the effective stress.Fig 4.Presents the damping ratio profiles.As it can be seen damping ratio has a reverse concordance with soil stiffness represented by G/Gmax ratio.Above phreatic line effective stress increases in unsaturated condition due to suction effect leading to an increase in G/Gmax, and hence, reduction in damping ration of the soil.For the saturated dry condition, there is no suction effect and the decrease in effective stress with elevation due to surchage reduction controls the dam behavior (see also equation 6 for further insight).

Figure 2 .
Figure 2. The variation of shear stress for 3 different heights for two modeling scenarios (dry-sat and unsat-sat models)

Figure 7 .Figure 8 .Figure 9 .Figure 10 .
Figure 7. Crest acceleration history for the dam with 35 m height in two modeling scenarios conventional modelling which neglects the unsaturated zone above phreatic line.
Where G is shear modulus, G max denotes max shear modulus, K 2,max is a material coefficient which depends on porosity of soil (in our study, it was set as 60), P a is

Table 1
presents the dam soil properties.The hydraulic properties (Van Genuchten parameters) were obtained from SoilVision database, They are presented in Table2.Where a, nvg, ad mvg are material coefficients (constant parameters), ws and wr are UNSAT saturated and residual volumetric water content, respectively and  denotes matric suction.
3. depicts the shear modulus profiles of the earth dams with height.Under the ground water table in both cases, shear modulus decrease with increasing elevation but above water table (phreatic line), shear modulus increases in saturated-unsaturated scenario.The reason is that shear modulus has a direct relation with the effective stress.