Propagation of seismic waves through liquefied soils

ARTICLE IN PRESS Soil Dynamics and Earthquake Engineering 30 (2010) 236–257 Contents lists available at ScienceDirect Soil Dynamics and Earthquake Engineering journal homepage: www.elsevier.com/locate/soildyn Propagation of seismic waves through liquefied soils Mahdi Taiebat a,, Boris Jeremić b, Yannis F. Dafalias b,c, Amir M. Kaynia d, Zhao Cheng e a Department of Civil Engineering, University of British Columbia, Vancouver, BC, Canada V6T 1Z4 Department of Civil and Environmental Engineering, University of California, Davis, CA 95616, USA Department of Mechanics, National Technical University of Athens, Zographou 15780, Hellas d Norwegian Geotechnical Institute, P.O. Box 3930 Ullevaal Stadion, N-0806 Oslo, Norway e Earth Mechanics Inc., Oakland, CA 94621, USA b c a r t i c l e in f o a b s t r a c t Article history: Received 23 March 2009 Received in revised form 21 November 2009 Accepted 23 November 2009 To predict the earthquake response of saturated porous media it is essential to correctly simulate the generation, redistribution, and dissipation of excess pore water pressure during and after earthquake shaking. To this end, a reliable numerical tool requires a dynamic, fully coupled formulation for solid– fluid interaction and a versatile constitutive model. Presented in this paper is a 3D finite element framework that has been developed and utilized for this purpose. The framework employs fully coupled dynamic field equations with a u–p–U formulation for simulation of pore fluid and solid skeleton interaction and a SANISAND constitutive model for response of solid skeleton. After a detailed verification and validation of the formulation and implementation of the developed numerical tool, it is employed in the seismic response of saturated porous media. The study includes examination of the mechanism of propagation of the earthquake-induced shear waves and liquefaction phenomenon in uniform and layered profiles of saturated sand deposits. & 2009 Elsevier Ltd. All rights reserved. Keywords: Numerical analysis Wave propagation Earthquake Liquefaction Constitutive modeling 1. Introduction Performance-based design (PBD) of geotechnical structures is gaining popularity in professional practice and is fostering research in the academic community. PBD relies heavily on modeling and simulation tools. One of the challenges in PBD of geotechnical and geophysical problems is analysis of dynamic transient phenomena in fluid-saturated porous media and, in particular, modeling and simulations of seismic wave propagation. This subject, which can be related to liquefaction, continues to challenge engineering research and practice. In addition, lateral movement of sloping ground as a result of pore pressure generation poses other challenges in PBD. Pore pressure generation phenomenon commonly occurs in loose to medium dense sands that are fully saturated and may induce two related phenomena: the flow liquefaction which leads to flow slides and the cyclic mobility which leads to lateral spreads. Pore pressure generation phenomenon commonly occurs in loose to medium dense sands that are fully saturated and may induce two related phenomena: the flow liquefaction that leads to flow slides and the cyclic mobility that leads to lateral spreads. Flow liquefaction occurs when shear stresses required for static equilibrium exceed residual shear  Corresponding author. Tel.: + 1 604 822 3279; fax: + 1 604 822 6901. E-mail address: mtaiebat@civil.ubc.ca (M. Taiebat). 0267-7261/$ - see front matter & 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.soildyn.2009.11.003 strength, as opposed to liquefaction, which implies a zero effective stress. In this phenomenon an earthquake brings soil to the point of instability, and deformations are driven by static stresses. The associated failure occurs rapidly with little warning and produces large deformations. On the other hand, cyclic mobility occurs when shear stresses required for static equilibrium are less than residual shear strength. These deformations are driven by dynamic stresses, occur incrementally, and can be large or small [1]. Accumulation of permanent deformations, degradation of soil moduli, increase of hysteretic damping, and change of soil fabric as a function of imposed cyclic shear strains require advanced models and implementations that are tricky but doable. In order to systematically investigate this subject it is necessary to predict the generation, redistribution, and dissipation of excess pore pressures during and after earthquake shaking and their impact on the transmitted waves. A fundamental approach requires a dynamic coupled stress-field analysis. Fully coupled transient response of solid–pore fluid interaction and constitutive behavior of the soil skeletonsolid skeleton play equally important roles in successful numerical simulation of response in saturated granular soil medium. The mechanical model of this interaction, when combined with a suitable constitutive description of the solid phase and with efficient, discrete, computation procedures, allows most transient and static problems involving deformations to be properly modeled and accurately simulated. ARTICLE IN PRESS M. Taiebat et al. / Soil Dynamics and Earthquake Engineering 30 (2010) 236–257 An advanced numerical simulation tool has been developed and utilized for this purpose. The issues of interests here are the capabilities of the advanced constitutive model and the rigorous framework that is used for modeling the fully coupled solid skeleton–pore fluid interaction. Many constitutive models with different levels of complexity have been formulated to describe the response of soils in cyclic loading, e.g., [2–7]. The bounding surface models have generally proved efficient and successful in simulations of cyclic loading of sands [8]. Recently, greater consideration has been given to the role of the progressive decay in soil stiffness with increasing pore pressure, accumulation of deformation, stress dilatancy, and hysteretic loops in formulating constitutive models for liquefiable soils [9,10], among others. The bounding surface constitutive model in this study predicts with accuracy the soil response of noncohesive soils during loading under different paths of monotonic and cyclic loading, drained and undrained conditions, and for various soil densities, stress levels, and loading conditions. From a practical point of view, it is probably most important that for all aforementioned conditions a single soil-specific set of constants is needed. Besides the superiority of this constitutive model in modeling of displacement proportional damping, the other attractive feature of the present framework is proper modeling of velocity proportional damping. This is done by taking into account the interaction of pore fluid and solid skeleton. In this way the present formulation models in a realistic way the physical damping, which dissipates energy during any wave propagation. This paper summarizes the key modeling features of our modeling and simulation tool, including the soil stress–strain behavior and the finite element formulation. A detailed calibration of the employed constitutive model for three different sands has been presented; this can be used as a guideline for future applications of the model. Using these model parameters, performance of the constitutive model in reproducing the material response for different types and conditions of loading, as well as different densities and stress states, is illustrated. This is followed by verification of the implemented fully coupled element in modeling of a complex step loading on a saturated elastic medium. The rest of the paper deals with modeling and simulation of earthquake-induced soil liquefaction and its effects on performance of layered soil systems. Specifically, the seismic response of a uniform soil column with relatively dense sand is compared with that of a layered soil column with a deep liquefiable loose layer. This liquefiable layer could operate as an isolating layer, causing reduction of accelerations in shallow layers. In addition, an example for illustration of the possible effects of a loose interlayer and the magnitude of the base acceleration on the resulting lateral displacements as well as details of the wave propagation in a relatively dense sloping soil column on the resulting lateral displacements as well as details of the wave propagation have been numerically studied and discussed. The simulations presented rely on verified formulation and implementation of behavior of fully coupled porous media and on validated constitutive material modeling. The numerical simulations presented in this paper illustrate the ability of advanced computational geomechanics tools to provide valuable detailed information about the underlying mechanisms. 2. Mathematical formulation The material behavior of soil skeletonsolid skeleton is interrelated to the pore fluid pressures. The behavior of pore fluid is 237 assumed to be elastic, and thus all the material nonlinearity is concentrated in the soil skeletonsolid skeleton. The soil behavior (mix of soil skeletonsolid skeleton and pore fluid) can thus be described using a single-phase constitutive analysis approach for the skeleton combined with the full coupling with pore fluid. These are the two major parts of the formulation in the present study. 2.1. Constitutive model The SANISAND constitutive model is used here for modeling of soil response. SANISAND is the name used for a family of Simple ANIsotropic SAND constitutive models within the frameworks of critical state soil mechanics and bounding surface plasticity. Manzari and Dafalias [7] constructed a simple stress-ratiocontrolled constitutive model for sand in a logical sequence of simple steps. The model is fully compatible with critical state soil mechanics principles; it renders the slope of the dilatancy stress ratio (also known as the phase transformation line), a function of the state parameter c, such that at the critical state where c ¼ 0 the dilatancy stress ratio coincides with the critical state failure stress ratio. In addition, softening of dense samples is modeled within a collapsing peak stress ratio bounding surface formulation. The peak stress ratio is again made a function of c such that at the critical state where c ¼ 0 it becomes the critical state stress ratio, following an original suggestion by Wood et al. [11]. The bounding surface feature enables reverse and cyclic loading response simulation. Dafalias and Manzari [10] extended the model to account for the fabric changes during the dilatant phase of deformation on the subsequent contractant response upon load increment reversals. An additional mechanism was also proposed by Dafalias et al. [12] for the inherent anisotropy effects. Using the concept of limiting compression curve [13] and a proper closed yield surface, Taiebat and Dafalias [14] eliminated the inability of the previous versions of the model to induce plastic deformation under constant stress-ratio loading, which is especially important at high confining pressures causing grain crushing. In the present paper the focus is on wave propagation in granular media. To involve fewer model parameters and for simplicity, the version of the SANISAND model with fabric change effects [10] has been considered as the constitutive model for the soil. The inherent anisotropy [12] and the plastic strains under constant-stress ratios [14] have not been accounted for in the present work. For quick reference an outline of the constitutive model in its generalized form for multiaxial stress space is given in Appendix A. 2.2. Fully coupled solid skeleton–pore fluid interaction The mechanical model of the interaction between solid skeleton and pore fluid, when combined with a suitable constitutive description of the solid phase and with efficient, discrete, computation procedures, allows one to solve most transient and static problems involving deformations. The modeling framework described here, based on the concepts originally outlined by Biot [15], is appropriate for saturated porous media. For modeling of the fully coupled solid–pore fluid interaction, Zienkiewicz and Shiomi [16] proposed three general continuum approaches, namely (a) u–p, (b) u–U, and (c) u–p–U formulations (u, solid displacement; p, pore pressure; and U, fluid displacement). In this study the coupled dynamic field equations with u–p–U formulation have been used to determine pore fluid and soil ARTICLE IN PRESS 238 M. Taiebat et al. / Soil Dynamics and Earthquake Engineering 30 (2010) 236–257 skeletonsolid skeleton responses. This formulation uses pore fluid pressure as an additional (dependent) unknown field to stabilize the solution of the coupled system. The pore fluid pressures have been connected to (dependent on) the displacements of pore fluid so that with known volumetric compressibility of the pore fluid pressure can be calculated. Despite its power, this formulation has rarely been implemented into finite element codes. The formulation takes into account velocity proportional damping (usually called viscous damping) by proper modeling of coupling of pore fluid and solid skeleton, while the displacement proportional damping is appropriately modeled using elastoplasticity with the powerful material model chosen. No additional (and artificial) Rayleigh damping has been used in the finite element model. The main components of the u–p–U formulation for porous media are outlined in Appendix B for quick reference. Detailed descriptions of the u–p–U formulation, finite element discretization, and time integration are presented in [17]. Fig. 1. Schematic description of Verification and Validation procedures and the computational solution (after [28,37]). 3.2. Verification and validation 3. Numerical tool 3.1. Program implementation In order to study the dynamic response of saturated soil systems as a boundary value problem, an open-source 3D numerical simulation tool has been developed. Parts of OpenSees framework [18] have been used to connect the finite element domain. In particular, FEM classes from OpenSees (namely, class abstractions for Node, Element, Constraint, Load, and Domain) have been used to describe the finite element model and to store the results of the analysis performed on the model. In addition, Analysis classes were used to drive the global level finite element analysis, i.e., to form and solve the global system of equations. As for the Geomechanics modules, a number of elements, algorithms, and material models from UC Davis Computational Geomechanics toolset have been used. In particular, sets of NewTemplate3Dep [19] numerical libraries have been used for constitutive level integrations, nDarray numerical libraries [20] were used to handle vector, matrix, and tensor manipulations, while FEMtools libraries (u–p–U element) were used in implementation of the coupled solid–fluid interaction at the finite element level. Finally, solvers from the uMfPACK set of libraries [21] were used to solve the nonsymmetric global (finite element level) system of equations. The SANISAND constitutive model and the fully coupled u–p–U element are implemented in the above-mentioned numerical simulation tool [17,19,22]. The constitutive model is integrated using a refined explicit integration method with automatic error control for yield surface drift based on the work by Sloan et al. [23]. In order to develop integration of the dynamic finite element equation in the time domain, the weak form of the finite element formulation is rewritten in a residual form [24], and the resulting set of residual (nonlinear) dynamic equations is solved using the Hilber–Hughes–Taylor (HHT) amethod [25–27]. The developed finite element model for this formulation uses eight-node brick elements. Because the pore fluid is compressible, there are no problems with locking, particularly if good equation solvers are used, therefore there is no need for lower order of interpolation for pore fluid pressures. We use u–p–U formulation where an additional unknown for fluid is used (pore pressure in addition to fluid displacements); this helps in stabilizing the system. Pore fluid is compressible but it is many orders of magnitude more stiff than solid skeleton. All of the above libraries and implementations are available either through their original developers or through the second author’s website (http://geomechanics.ucdavis.edu). Prediction of mechanical behavior comprises use of computational model to foretell the state of a physical system under conditions for which the computational model has not been validated [28]. Confidence in predictions relies heavily on proper verification and validation processes. Verification is the process of determining that a model implementation accurately represents the developer’s conceptual description and specification. It is a Mathematics issue. Verification provides evidence that the model is solved correctly. Verification is also meant to identify and remove errors in computer coding and verify numerical algorithms. It is desirable in quantifying numerical errors in computed solution. Validation is the process of determining the degree to which a model is an accurate representation of the real world from the perspective of the intended uses of the model. It is a Physics issue. Validation provides evidence that an appropriate model is used. Validation serves two goals: (a) a tactical goal in identification and minimization of uncertainties and errors in the computational model and (b) a strategic goal in increasing confidence in the quantitative predictive capability of the computational model. Verification and validation procedures are hence the primary means of assessing accuracy in modeling and computational simulations. Fig. 1 depicts relationships between verification, validation, and the computational solution. In order to verify the u–p–U formulation, a number of closed form or very accurate solutions should be used. To this end and also to illustrate the performance of the formulation and versatility of the numerical implementation, several elastic 1D problems have been studied; these include drilling of a well, the case of a spherical cavity, consolidation of a soil layer, line injection of fluid in a reservoir, and shock wave propagation in saturated porous medium. The verification process for shock wave propagation in porous medium, which is the most important and difficult example among these five, has been presented here. Validation has been carried out on an extensive set of physical tests on Toyoura, Nevada, and Sacramento river sands and has been presented along with detailed explanation on calibration of the constitutive model parameters. 3.2.1. Verification: u–p–U element The analytic solution developed by Gajo [29] and Gajo and Mongiovi [30] for 1D shock wave propagation in ARTICLE IN PRESS 239 M. Taiebat et al. / Soil Dynamics and Earthquake Engineering 30 (2010) 236–257 permeability k ¼ 106 cm=s, which creates very high coupling between porous solid and pore fluid. The second model had permeability k ¼ 102 cm=s, which creates a low coupling between porous solid and pore fluid. Good agreement was obtained between the numerical simulations and the analytical solution as presented in Fig. 3. Fig. 2. Viscous coupling finite-element analysis of infinite half-space subjected at the surface to a step displacement boundary condition. Table 1 Simulation parameters used for the shock wave propagation verification problem. Parameter Symbol Value Poisson ratio Young’s modulus n 0.3 E Solid particle bulk modulus Ks 1:2  106 kN=m2 Fluid bulk modulus Kf Solid density Fluid density rs rf Porosity HHT parameter n 3:6  107 kN=m2 2:17  106 kN=m2 2700 kg=m3 1000 kg=m3 0.4  0.2 a 3.2.2. Validation: SANISAND model Material model validation can be performed by comparing experimental (physical) results and numerical (constitutive) simulations for different type of sands in different densities, confining pressures, loading paths, and drainage conditions. The SANISAND constitutive model has been specifically tested for its ability to reproduce a series of monotonic and cyclic tests on Toyoura sand [31,32], Nevada sand [33], and Sacramento river sand [34] in a wide range of relative densities and confining pressures and in different drainage conditions. The 14 material parameters required by the model for these sands are listed in Table 2; they are divided into different groups according to the particular role they play. The material parameters can be selected mainly from standard types of laboratory tests. Here the calibration has been carried out using drained and undrained triaxial compression and extension tests at different values of initial void ratio and confining pressure in order to calibrate different features of hardening/softening and Table 2 Material parameters of the SANISAND constitutive model for three types of sands. Parameter elastic porous medium was used in the verification study. A model was developed consisting of 1000 eight-node brick elements, with boundary conditions that mimic 1D behavior (Fig. 2). In particular, no displacement of solid (ux ¼ 0, uy ¼ 0) and fluid (Ux ¼ 0, Uy ¼ 0) in x and y directions, i.e., in horizontal directions, was allowed along the height of the model. Bottom nodes have full fixity for solid ðui ¼ 0Þ and fluid ðUi ¼ 0Þ displacements, while all the nodes above base are free to move in the z direction for both solid and fluid. Pore fluid pressures are free to develop along the model. Loads to the model consist of a unit step function (Heaviside) applied as (compressive) displacements to both solid and fluid phases of the model, with an amplitude of 0.001 cm. The u–p–U model dynamic system of equations was integrated using the HHT algorithm. Table 1 presents relevant parameters for this verification. Two set of permeability of material were used. The first model had Elasticity Toyoura Nevada Sacramento G0 125 0.05 150 0.05 200 0.2 1.25 0.712 0.934 0.019 0.7 1.14 0.78 0.83 0.027 0.45 1.35 0.65 0.96 0.028 0.7 2.1 1.05 2.0 0.704 0.81 0.8 1.25 2.56 1.2 n CSL M c e0 l x Dilatancy 1.5 nd A0 Kinematic n Hardening h0 ch 7.05 0.968 9.7 1.02 5.0 1.03 Fabric dilatancy zmax cz 2.0 600 5.0 800 – – b 1.5 1.25 Fluid Displ. (x10−3cm) Solid Displ. (x10−3cm) Symbol K=10−6cm/s 1 0.75 K=10−2cm/s 0.5 0.25 K=10−2cm/s 1.25 1 K=10−6cm/s 0.75 0.5 0.25 Lines: Closed−form Symbols: FEM 0 0 4 6 8 10 12 time (µsec) 14 4 6 8 10 time (µsec) 12 Fig. 3. Compressional wave in both solid and fluid, comparison with closed form solution. 14 ARTICLE IN PRESS 240 M. Taiebat et al. / Soil Dynamics and Earthquake Engineering 30 (2010) 236–257 dilatancy/contractancy of the model. A step-by-step calibration process for the model constants for Toyoura sand is shown in Figs. 4–7. The G0 parameter that defines the elastic shear modulus of sand in Eq. (2)1 can be calibrated using the stress–strain curves or from the elastic wave propagation tests in the field or laboratory. In this paper, a good estimation of G0 for monotonic shearing in Toyoura sand is obtained by fitting the triaxial versions of Eqs. (1)1 and (2)1 to the initial stages of the deviatoric stress– strain ðeq qÞ data of triaxial drained tests in Fig. 4. The critical state parameters consist of Mc;e , the critical stress ratio in triaxial compression and extension, and the parameters ðe0 Þ, l, and x of the equation ec ¼ e0 lðpc =pat Þx , which defines the position of the critical state line (CSL) in the space of void ratio vs. mean effective stress [35]. Calibration of these parameters requires monotonic tests that approach critical state. Results of the drained and undrained triaxial compression tests on Toyoura sand are presented in Fig. 5 for fitting these parameters. These bounding and dilatancy surfaces in the SANISAND model are generalizations of the bounding (or peak), dilatancy (or phase transformation) lines in triaxial space, expressed analytically by Z ¼ Mexpðnb cÞ and Z ¼ Mexpðnd cÞ, respectively. The parameters nb and nd can be determined by evaluating these equations at the peak and phase transformation states. The calibration process is b shown in Fig. 6(a and b), where c and Zb are the values of c and Z at the peak stress-ratio state, and cd and Zd are the values of c 20 Go=125 15 10 5 0 0 0.2 0.4 0.6 Deviatoric strain, (%) 0.8 Fig. 4. Calibration of the G0 constant using monotonic triaxial tests on Toyoura sand [31]. 4000 1 CIUC−after shearing CIDC−after shearing CIUC−initial CIUC−after shearing CIDC−initial CIDC−after shearing 0.95 q (kPa) 3000 2000 M=(q/p)critical=1.25 1000 Void ratio, e 100q(1+e)/[3(patp)1/2(2.97−e)2] 25 and Z at the phase transformation state. Note that the phase transformation state corresponds to the peak of the volumetric strain ev in drained tests and the peak of the excess pore pressure u2u0 in undrained tests. Direct estimation of the dilatancy parameter A0 requires good quality stress-dilatancy data, where A0 can then be calibrated based on ev 2eq curves. Here the parameters A0 (related to dilatancy) and h0 and ch (related to the effect of distance from bounding surface) are estimated by trial and error, as shown in Fig. 6(c–f). Finally, fabric-dilatancy constants zmax and cz require trial and error fitting of loading–unloading reverse loading, or cyclic data, preferably undrained, where Z must exceed M d in the process so that evolution of z is activated. Fig. 7 shows calibration of the zmax parameter using the loading–unloading undrained triaxial tests on Toyoura sand. Figs. 8 and 9 compare the experimental data [31] and model simulations for drained and undrained triaxial compression tests followed by unloading on isotropically consolidated samples of Toyoura sand. These experiments/simulations cover a wide range of confining pressures and densities and show variations of response from highly dilatant to highly contractant. This variety of response depends on the combination of density and confining pressure that is reflected in the value of state parameter c. Figs. 10 and 11 compare the experimental data [32] and model simulations for drained constant-p cyclic triaxial tests on isotropically consolidated loose and dense samples of Toyoura sand. In these tests the amplitude of shear strain (g ¼ ea er , as defined in [32]) has been increased with cyclic loading. Accumulation of compressive volumetric strain with cyclic loading can be observed in Fig. 10, and clearly when the stress ratio exceeds a certain value, which varies by the state, the specimen begins to dilate. As expected, the dilatancy feature is more pronounced for the dense sample in Fig. 11. Experimental results [33] and numerical simulation for constant-p triaxial compression tests on Nevada sand in different densities and confining pressures are presented in Fig. 12. In addition, the stress path and shear response of an undrained cyclic triaxial test on anisotropically consolidated Nevada sand, at relative density of about 40%, are presented in Fig. 13. The last set of simulations for the purpose of model validation are presented in Fig. 14, where experimental results of drained triaxial compression tests on Sacramento river sand [34] have been compared with the model simulations. The experiments/ simulations are on both loose and dense samples of Sacramento river sand in different initial confining pressures. It should be noted that the fabric-dilatancy parameters, needed in correct 0.9 0.85 0.8 0.7 0.75 ec=0.934−0.019(p/pat) 0.7 0 0 1000 2000 p (kPa) 3000 0 2 4 6 8 10 (p/pat)0.7 Fig. 5. Calibration of the CSL constants using monotonic triaxial tests on Toyoura sand [31]. 12 ARTICLE IN PRESS 241 0 0 −0.02 −0.1 ln(ηd/M) ln(M/ ηb) M. Taiebat et al. / Soil Dynamics and Earthquake Engineering 30 (2010) 236–257 −0.04 nb=1.25 −0.06 −0.2 nd=2.1 −0.3 −0.08 −0.06 −0.04 −0.02 −0.4 −0.2 0 −0.15 ψb 300 −0.1 −0.05 0 ψd 3000 pin=100 kPa e=0.833, Dr=37.9% A0=0.7 200 A0=0.2 q (kPa) q (kPa) 250 150 100 A0=2.0 2000 A0=0.2 A0=0.7 1000 A0=2.0 50 0 0 0.9 0.92 0.94 0.96 Void ratio, e 0.98 0 1 1000 2000 p (kPa) 3000 1.6 1600 ho(1−che)=3.0 1.2 1.4 e=0.833 0.6 800 400 1.6 0.9 0.4 0 0 5 7.05(1−0.968e) ho=7.05 ch=0.968 0.8 0.4 e=0.907 10 15 20 Axial strain (%) ho(1−che) q (kPa) 1200 0 25 0.8 0.84 0.88 0.92 Void ratio, e 0.96 Fig. 6. Calibration of the bounding and dilatancy surface constants using monotonic triaxial tests on Toyoura sand [31]. 3000 Close match of physical testing data with constitutive predictions represents a satisfactory validation of the material models. This validation with previous verification gives confidence in simulations of real, prototype behavior. 2000 4. Numerical simulation examples q (kPa) 4000 1000 zmax=2 zmax=0 0 0 1000 2000 p (kPa) 3000 Fig. 7. Calibration of the zmax parameter using the loading–unloading triaxial tests on Toyoura sand [31]. simulation loading–unloading cycles, have not been calibrated for Sacramento river sand, as the authors could not find enough related experimental data in the literature. Liquefaction of level and sloping ground represents a common situation during earthquakes. It is therefore of interest to use the developed numerical tool for simulating the earthquake wave propagation in such cases. This section presents results of numerical simulation for different cases of a 1D site response, using the available 3D elements in the developed numerical simulation tool. Two general examples have been studied. The first example presents results of a study on the isolating effect of a liquefied sand layer in propagation of seismic waves. In the second example the seismicinduced shear deformation of a gentle slope in the presence or absence of a liquefiable layer has been studied. It should be noted that the simulations are based on the small deformation assumption, which might introduce a large ARTICLE IN PRESS 242 M. Taiebat et al. / Soil Dynamics and Earthquake Engineering 30 (2010) 236–257 Experiment Simulation 1200 1200 900 pin=500 kPa pin=100 kPa 600 q (kPa) q (kPa) 900 300 600 300 0 0 0 5 10 15 20 Axial strain (%) 1200 25 30 0 5 10 15 20 Axial strain (%) 1200 pin=500 kPa pin=100 kPa 25 30 pin=500 kPa pin=100 kPa 900 q (kPa) 900 q (kPa) pin=500 kPa pin=100 kPa 600 300 600 300 0 0 0.8 0.85 0.9 0.95 Void ratio, e 1 0.8 0.85 0.9 0.95 Void ratio, e 1 Fig. 8. Simulations vs. experiments in drained triaxial compression tests on isotropically consolidated samples of Toyoura sand [31]. Simulation 4000 3000 3000 q (kPa) q (kPa) Experiment 4000 e=0.735, Dr=63.7% e=0.833, Dr=37.9% e=0.907, Dr=18.5% 2000 1000 1000 0 0 0 4000 5 10 15 20 Axial strain (%) 25 30 0 4000 e=0.735, Dr=63.7% e=0.833, Dr=37.9% e=0.907, Dr=18.5% 2000 5 10 15 20 Axial strain (%) 25 30 e=0.735, Dr=63.7% e=0.833, Dr=37.9% e=0.907, Dr=18.5% 3000 q (kPa) 3000 q (kPa) e=0.735, Dr=63.7% e=0.833, Dr=37.9% e=0.907, Dr=18.5% 2000 2000 1000 1000 0 0 0 1.000 2.000 p (kPa) 3.000 0 1000 2000 p (kPa) 3000 Fig. 9. Simulations vs. experiments in undrained triaxial compression tests on isotropically consolidated samples of Toyoura sand [31]. ARTICLE IN PRESS 243 M. Taiebat et al. / Soil Dynamics and Earthquake Engineering 30 (2010) 236–257 Experiment Simulation 2 2 p=100 kPa (const.) ein=0.845 Stress ratio, q/p Stress ratio, q/p p=100 kPa (const.) 1 0 −1 −1 0 1 Shear strain, γ (%) 2 −2 2.5 −1 0 1 Shear strain, γ (%) 2 2.5 Volumetric strain, εv (%) Volumetric strain, εv (%) 0 −1 −2 2 1.5 1 0.5 0 −3 −2 −1 0 1 Shear strain, γ (%) 2 2 1.5 1 0.5 0 −3 3 2.5 −2 −1 0 1 Shear strain, γ (%) 2 3 2.5 Volumetric strain, εv (%) Volumetric strain, εv (%) ein=0.845 1 2 1.5 1 0.5 0 2 1.5 1 0.5 0 −1 −0.5 0 0.5 Stress ratio, q/p 1 1.5 −1 −0.5 0 0.5 Stress ratio, q/p 1 1.5 Fig. 10. Simulations vs. experiments in constant-p cyclic triaxial tests on a relatively loose sample of Toyoura sand [32]. error in the strain tensor (by neglecting the quadratic portion of displacement derivatives) when soil undergoes large deformations. 4.1. Example 1: isolation of ground motions in level ground A 10 m vertical column of saturated sand consisting of 20 u–p– U eight-node brick elements was considered, subjected to an earthquake shaking at the base. Two particular cases have been studied in this example. In the first case, all 20 elements of the soil column are assumed to be at a medium dense state at the beginning of analysis ðein ¼ 0:80Þ. In the second case, the deepest 2 m of the soil column are assumed to be at looser initial states compared with the first case. In particular, two of the deepest elements are assumed to be at ein ¼ 0:95 (loose) and the next two elements at ein ¼ 0:875 (medium loose), to make a smoother transition between the loose and dense elements. The rest of the soil column, i.e., the upper 16 elements (8 m) are at the same density as the first case ðe ¼ 0:80Þ. A schematic illustration of these two cases is presented in Fig. 15. The elements are labeled from E01 (bottom) to E20 (surface) in Fig. 15. The boundary conditions are such that the vertical displacement degrees of freedom (DOF) of the soil and water at the base ðz ¼ 0Þ are fixed, while the pore pressure DOFs are free. The soil and water vertical displacement DOFs at the surface ðz ¼ 10 mÞ are free to simulate the upward drainage. The pore pressure DOFs are fixed at the surface. On the sides, solid skeleton and water are prevented from moving in the y direction, whereas movements in the x and z directions are free. It is emphasized that the displacements of solid skeleton and pore fluid are different. In order to simulate the 1D behavior, all of the related DOFs of the nodes at each depth are connected in a master–slave fashion, i.e., the nodes at the same level have the same ui ði ¼ x; y; zÞ. The same thing is true for p and Ui . The values of ui and Ui , however, may be different at a node, allowing for relative movement of soil and fluid. The soil is assumed to be Toyoura sand and is characterized by the validated SANISAND model explained in previous sections ARTICLE IN PRESS 244 M. Taiebat et al. / Soil Dynamics and Earthquake Engineering 30 (2010) 236–257 Experiment Simulation 2 2 p=100 kPa (const.) ein=0.653 1 Stress ratio, q/p Stress ratio, q/p p=100 kPa (const.) 0 −1 −1 0 1 Shear strain, γ (%) 2 −2 −1 0 1 Shear strain, γ (%) −2 −1 0 1 Shear strain, γ (%) 2 0.6 Volumetric strain, εv (%) 0.6 Volumetric strain, εv (%) 0 −1 −2 0.3 0 −0.3 −0.6 0.3 0 −0.3 −0.6 −2 −1 0 1 Shear strain, γ (%) 2 2 0.6 Volumetric strain, εv (%) 0.6 Volumetric strain, εv (%) ein=0.653 1 0.3 0 −0.3 −0.6 −1.5 −1 −0.5 0 0.5 1 Stress ratio, q/p 1.5 0.3 0 −0.3 −0.6 −1.5 −1 −0.5 0 0.5 1 Stress ratio, q/p 2 1.5 2 Fig. 11. Simulations vs. experiments in constant-p cyclic triaxial tests on a relatively dense sample of Toyoura sand [32]. 1 Dr=60% 40% 100 60% 40% 60% q (kPa) 200 p=160 kPa p=80 kPa p=40 kPa 40% Volumetric strain (%) 300 0 Dr=40% −1 p=160 kPa −2 80 40 p=160 kPa −3 Dr=60% −4 80 40 −5 0 0 2 4 6 Axial strain (%) 8 10 0 2 4 6 Axial strain (%) 8 10 Fig. 12. Simulations (solid lines) vs. experiments (symbols) in constant-p triaxial compression tests on isotropically consolidated samples of Nevada sand [33]. ARTICLE IN PRESS 245 M. Taiebat et al. / Soil Dynamics and Earthquake Engineering 30 (2010) 236–257 Simulation 80 60 60 40 40 q (kPa) q (kPa) Experiment 80 20 0 0 −20 −20 −40 −40 0 40 80 p (kPa) 120 160 0 80 80 60 60 40 40 q (kPa) q (kPa) 20 20 0 −20 −20 0 0.5 1 Axial strain (%) 1.5 2 80 p (kPa) 120 160 20 0 −40 −0.5 40 −40 −0.5 0 0.5 1 Axial strain (%) 1.5 2 Fig. 13. Simulations vs. experiments in undrained cyclic triaxial tests on Nevada sand with Dr  40% [33]. Loose Dense 6 pin=1240 kPa q (kPa) 3000 1030 590 2000 290 440 200 100 100 1000 0 0 2 4 6 Axial strain (%) 8 Volumetric strain (%) 4000 Loose Dense 3 440 200 100 100 290 0 −3 590 −6 pin=1240 kPa −9 10 1030 0 2 4 6 Axial strain (%) 8 10 Fig. 14. Simulations (solid lines) vs. experiments (symbols) in drained triaxial compression tests on isotropically consolidated loose and dense samples of Sacramento river sand [34]. with the model parameters given in Table 2. The other parameters related to the boundary value problem are given in Table 3. The permeability is assumed to be isotropic and equal to k ¼ 5:0 104 m=s. A self-weight analysis was performed before the base excitation. The resulting fluid hydrostatic pressures and soil stress states along the soil column serve as initial conditions for the subsequent dynamic analysis. An input acceleration time history (Fig. 15), taken from the recorded horizontal acceleration of Model No. 1 of VELACS project at RPI [36], was considered in which amax ¼ 0:2g and main shaking lasts for about 12 s. In addition to the (physical) velocity and displacement proportional damping from the element and the material model, respectively, a small amount of numerical damping has been introduced through constants of the HHT algorithm in order to stabilize the dynamic time stepping. No Rayleigh damping has been used. Results of the analysis during and after the shaking phase for both cases (i.e., the uniform and the layered soil profiles) are presented in Figs. 16–19. In these figures the plots on the left and right correspond to the cases of uniform and layered soil profiles, respectively. In particular, Figs. 16 and 17 show variations of the shear stress sxz vs. vertical effective stress sz and shear strain g, respectively. Figs. 18 and 19 show time histories of void ratio e and horizontal ðxÞ component of acceleration in the solid part of the mixture au;x , respectively. Note that different rows of plots in Figs. 16–18 correspond to different elements. Similarly, different rows in Fig. 19 correspond to different elevations of the nodes. Fig. 16 a shows the typical mechanism of cyclic decrease in vertical effective stress for the case of the uniform soil column. ARTICLE IN PRESS 246 M. Taiebat et al. / Soil Dynamics and Earthquake Engineering 30 (2010) 236–257 Fig. 15. Illustration of the problem in Example 1 in terms of the soil layering, the finite element mesh, and the input base acceleration. Table 3 Parameters used in the boundary value problem simulations (in addition to the material model parameters). Parameter Symbol Value 3:6  107 kN=m2 Solid particle bulk modulus Ks Fluid bulk modulus Kf Solid density Fluid density rs rf a HHT parameter 2:17  106 kN=m2 2700 kg=m3 1000 kg=m3  0.2 From early stages of shaking, signs of the so-called butterfly loops in the effective stress path can be partially observed at the upper levels, as these layers are in lower confining pressures compared with the deeper layers and thus are more dense than critical in the critical state soil mechanics terminology, i.e., with less contractive tendency. In later stages of shaking, when the confining pressures are reduced to smaller values, the butterfly shapes of the stress paths are more pronounced, and almost all depths of the soil column experience the cyclic mobility mechanism of liquefaction. In this stage the soil layers momentarily experience small values of effective stress followed by a dilative response in which the elements recover their strength in cycles of loading owing to their denser than critical state. As a result, the seismic-induced shear waves propagate all the way to the surface of the soil column. The response can be observed as an average degradation of stiffness and accumulation of shear strains in all levels of the soil column, as shown in Fig. 17 a. In the case of the layered soil profile, however, Fig. 16 b shows that the first cycle of shaking degrades the vertical effective stress in the loose layer of element 1, E01, to very small values in a flow liquefaction mechanism. The shear wave in this first cycle of shaking propagates to the surface and causes some reduction in vertical effective stress in upper soil layers. Because element 1, E01, is in loose states and has a strong contractive tendency, it does not recover its strength in cyclic loading, and therefore it remains in this liquefied state after the first cycle with negligible shear resistance. Therefore this layer acts as an isolating layer, and no shear stress can be transmitted to the upper layers after the first cycle. This mechanism can also be observed in the shear stress–shear strain plot of Fig. 17 b, as the first cycle completely liquefies the E01; after this cycle no shear stress is propagated to the upper layers. Of course the very loose element E01 shows considerable shear strain due to liquefaction. As shown in Fig. 16 b, in the upper layers after the first cycle of shaking no reduction in vertical effective stress occurs due to shearing (induced by shaking); however, one can observe that the vertical effective stress keeps decreasing in upper layers during (and even for a while after) shaking. This is because the dissipated pore pressure in deeper layers travels upward and reduces the effective stress in upper layers. The effective stress will then start to recover as the excess pore pressure dissipates from these layers. Fig. 18 shows the redistribution of void ratio during and long after the shaking. The horizontal dashed line in each plot shows ein;shaking , that is, the value of the void ratio for the corresponding layer at the beginning of shaking (at the end of self-weight analysis). For the uniform soil column, Fig. 18 a shows a general trend of consolidation in all of the soil layers. This consolidation starts later in upper layers because of the water pumped from the lower layers, which is to be dissipated by time. In the layered soil column (Fig. 18 b) the loose layer in element 1, E01, shows considerable consolidation due to shaking. The dissipated volume of water from this layer is pumped to the upper layers, and the result is some dilation (increase in void ratio) in these layers for some time after the end of shaking. By the time this pumped volume of water in the upper layers dissipates, all the layers show a consolidation response. In other words, in the layered soil column pumping of water causes initial loosening of soil in the upper layers, which is then followed by densification. Fig. 19 shows propagation of the horizontal acceleration through the soil column during shaking. In the uniform soil column Fig. 19 a shows that the base acceleration is transmitted to the surface of the soil column. It shows some spikes of acceleration with large amplitude at the surface of the soil resulting from the dilative phases of soil response during the shaking phase. In the layered soil column, however, Fig. 19 b shows that only the first cycle of shaking is transmitted to the surface (although with reduced amplitude) and subsequently the upper layers are isolated because of the fully liquefied state of the loose layer in element 1, E01. The transmitted accelerations to the surface of the soil column after the first cycle of shaking have very small amplitudes. Finally, Figs. 20 and 21 show contours of excess pore pressure (ue) and excess pore pressure ratio (Re) in the uniform and layered soil columns at different times. Note that ue =p pin =(sv)in  sv and Ru =ue/(sv)in, where pin and p are the initial and current levels of the pore pressure DOF in the u–p–U element, and the (sv)in and sv are the initial and current levels of the vertical effective stress. The initial state here refers to the beginning of the shaking stage when the excess pore pressure due to self-weight loading is completely dissipated and ue =Ru = 0. The excess pore pressure develops during the shaking stage (the first 15 s), propagates from deeper layers upward and dissipates long after the end of shaking. While excess pore pressure develops in the soil column as result of shaking, the value of sz reduces from its initial value (sv)in. The extreme case of sz = 0 results in Ru =1, ARTICLE IN PRESS 247 M. Taiebat et al. / Soil Dynamics and Earthquake Engineering 30 (2010) 236–257 20 20 E19 0 −20 20 E17 0 E17 0 −20 20 −20 20 E15 0 E15 0 −20 20 −20 20 E13 0 E13 0 −20 20 E11 0 −20 20 E09 0 −20 20 σxz (kPa) −20 20 σxz (kPa) E19 0 −20 20 E11 0 −20 20 E09 0 −20 20 E07 0 −20 20 E07 0 −20 20 E05 0 −20 20 E05 0 −20 20 E03 0 −20 20 E03 0 −20 20 E01 0 −20 E01 0 −20 0 25 50 σz (kPa) 75 100 0 25 50 σz (kPa) 75 100 Fig. 16. Variation of shear stress sxz vs. vertical effective stress sz for Elements E01–E19 during and after shaking. (a) Uniform. (b) Layered. which corresponds to the full loss of effective stress in the soil element and is known as flow liquefaction. For the case of uniform soil column Fig. 16a shows that sv drops to zero almost in all depths of the soil column. As a result the corresponding Ru during shaking approaches to 1 as shown in Fig. 21a. The dissipated pore pressure keeps the upper layers in a liquefied state even for long after the end of shaking before it completely dissipates from the surface. The situation in slightly different in the layered soil column. Fig. 21b shows that the value of excess pore pressure ratio Ru remains in the safe region (low values) during the shaking. Liquefaction is observed only in the deepest layers during shaking. The dissipated pore pressure from these layers propagates toward the surface of the soil column and causes some decrease in the magnitude of Ru in upper layers but does not liquefy those layers in this example. 4.2. Example 2: earthquake-induced shear deformations in gentle slopes A 10 m gently inclined ð1:53 Þ column of saturated sand consisting of 20 u–p–U eight-node brick elements subjected to a base shaking was considered. In terms of initial densities two particular cases have been studied in this example. In the first, the soil column is assumed to be uniform with e ¼ 0:7 (relatively dense, Dr C 75%), while in the other case a similar soil column is considered that includes a relatively loose layer with e ¼ 0:9 ðDr C 25%Þ at depth 2.5–3.5 m (Fig. 22). The boundary conditions are exactly the same as those presented in Example 1. The soil parameters and the other parameters related to the boundary value problem are also the same as those in Example 1; these are given in Tables 2 (Toyoura sand) and 3. The permeability is assumed to be isotropic and equal to k ¼ 103 m=s. A static application of gravity analysis was performed before the base excitation. The resulting fluid hydrostatic pressures and soil stress states along the soil column serve as initial conditions for the subsequent dynamic analysis. The input acceleration time history (Fig. 22) is considered with amax ¼ 0:2g, and shaking lasts for about 8 s. The following three cases have been studied in this paper: (I) Uniform soil profile with amax ¼ 0:2g. (II) Layered soil profile with amax ¼ 0:2g. (III) Layered soil profile with amax ¼ 0:4g. ARTICLE IN PRESS 248 M. Taiebat et al. / Soil Dynamics and Earthquake Engineering 30 (2010) 236–257 20 20 E19 0 −20 20 −20 20 E17 0 −20 20 E15 0 E15 0 −20 20 −20 20 E13 0 E13 0 −20 20 −20 20 E11 −20 20 E09 0 −20 20 E11 0 σxz (kPa) 0 σxz (kPa) E17 0 −20 20 −20 20 E09 0 −20 20 E07 0 −20 20 E07 0 −20 20 E05 0 −20 20 E05 0 −20 20 E03 0 −20 20 E03 0 −20 20 0 −20 −4.5 E19 0 −3 −1.5 γ (%) 0 E01 0 1.5 −20 −4.5 E01 −3 −1.5 γ (%) 0 1.5 Fig. 17. Variation of shear stress sxz vs. shear strain g for Elements E01–E19 during and after shaking. (a) Uniform. (b) Layered. The predicted behavior for cases I–III in terms of the vertical effective stress, shear stress, shear strain, excess pore pressure ratio ðRu Þ, and lateral displacement are presented and discussed in the rest of this section. The response is shown for selected elements along the soil profile (for positions of the elements refer to Fig. 22). In particular, Figs. 23a, 24a, and 25a show changes in the shear stress and vertical effective stress. Figs. 23b, 24b, and 25b show changes in the shear stress and shear strain. Finally, Figs. 23c, 24c, and 25c show time histories of the excess pore pressure ratio, Ru . Contours of excess pore pressure ðRu Þ, shear strain, and lateral displacements in all depths vs. time are presented in parts d, e and f of Figs. 23–25 for the above three cases. Note that Ru ¼ ue =ðsv Þin , where ue is the excess pore pressure and ðsv Þin is the initial effective vertical stress. Clearly Ru cannot be larger than 1 since ue r ðsv Þin (in cohesionless soils the stress state cannot go to the tensile side, in which ue 4 ðsv Þin ). The typical mechanism of response observed in different elements during the shaking is explained in the following. Generally, at the early stages of loading in a cycle the element shows a contractive response in shearing that is associated with an increase in excess pore pressure and a decrease in vertical effective stress and stiffness (in particular the shear stiffness). As loading continues to larger stress ratios the mechanism of response in the element might change to a dilative regime. This means that as shearing continues the excess pore pressure decreases and the vertical effective stress and the soil stiffness increase. In other words soil regains its stiffness during the dilative phase. Upon and after reversal of shear stress increment, again the sample shows a contractive regime with an increase in excess pore pressure and a decrease in vertical effective stress and soil stiffness. As the absolute value of the stress ratio exceeds a certain level (phase transformation) in the reverse loading (in negative stress ratios), the element again shows the dilative response. The final stress-ratio increment reversal in the loading ARTICLE IN PRESS 249 M. Taiebat et al. / Soil Dynamics and Earthquake Engineering 30 (2010) 236–257 0.7981 0.7981 E19 0.7969 E19 0.7969 0.7957 0.797 0.7957 0.7967 E17 E17 0.7955 0.7959 0.7941 0.7961 0.7952 0.7961 E15 0.7942 0.7952 0.7924 0.7952 E13 0.7929 E13 0.7947 0.7906 0.7942 0.7936 0.7955 E11 0.789 0.7942 E09 0.7914 0.7886 0.7939 E11 0.7942 void ratio 0.7916 void ratio E15 0.7943 0.7958 0.7929 0.7954 E09 0.7939 0.7923 0.7953 E07 0.7906 0.7874 0.793 E07 0.7934 0.7916 0.7955 E05 0.7894 0.7859 0.792 E05 0.7932 0.791 0.8688 E03 0.7878 0.7835 0.7914 E03 0.8649 0.861 0.9349 E01 0.7854 0.7793 E01 0.9166 0.8983 0 50 100 150 time (sec) 200 0 50 100 150 time (sec) 200 Fig. 18. Variation with time of void ratio e for Elements E01–E19 during and after shaking. (a) Uniform. (b) Layered. cycle is again associated with a contractive response. The above mechanism typically happens in all cycles. The overall contractive response could be stronger or weaker than the dilative response depending on the state of soil, i.e., its void ratio and effective stress. For looser samples or at higher confining pressures the contractive response is more pronounced, and therefore shearing at such states shows increase in excess pore pressure and decrease in vertical effective stress and stiffness. In denser samples or at lower confining pressures the response is the opposite. Fig. 23 shows detailed response of the soil profile in case I. Based on the previous general explanation the response starts with a contractive regime which decreases the vertical effective stress and the stiffness and increases the excess pore pressure ratio. Since a relatively dense state is assumed for this analysis ðe ¼ 0:7; Dr C 75%Þ, the soil has a general tendency of dilation as the stress ratio exceeds the dilation line (phase transformation). On the other hand, the inclination of the soil column poses asymmetric horizontal shear stresses (toward the direction of inclination) during cycles of shaking. This essentially makes the sample more prone to dilation at parts of the shaking cycles that are toward the direction of inclination (positive directions of shear stress and strain in this case), and the soil regains its stiffness and strength ðsv Þ in the dilative parts of the loading cycle. The dilative tendency of response and the offset shear stress in cyclic loading lead to asymmetric butterfly loops in stress space during the dilative and the subsequent intense contractive phases of response. The asymmetric horizontal shear stresses also result in asymmetric shear strains that accumulate toward the inclination of the soil profile and create permanent horizontal displacements in this direction. The dilative phases lead to instantaneous drops in the excess pore pressure ratio Ru . These drops are more pronounced in upper elements that have lower confining pressures and therefore have a stronger dilative tendency. Fig. 24 shows the response of the soil profile in case II. In this case the presence of the relatively loose layer at depth 2.5–3.5 m ðe ¼ 0:9; Dr C 25%Þ affects the response of the system. This layer has a stronger contractive response in cycles of loading. Such a ARTICLE IN PRESS 250 M. Taiebat et al. / Soil Dynamics and Earthquake Engineering 30 (2010) 236–257 0.5 0.5 z=10m 0 −0.5 0.5 z=9m 0 z=9m 0 −0.5 0.5 −0.5 0.5 z=8m 0 z=8m 0 −0.5 0.5 −0.5 0.5 z=7m 0 z=7m 0 −0.5 0.5 −0.5 0.5 z=6m −0.5 0.5 z=5m 0 −0.5 0.5 z=4m 0 −0.5 0.5 z=6m 0 au,x (g) 0 au,x (g) z=10m 0 −0.5 0.5 −0.5 0.5 z=5m 0 −0.5 0.5 z=4m 0 −0.5 0.5 z=3m 0 −0.5 0.5 z=3m 0 −0.5 0.5 z=2m 0 −0.5 0.5 z=2m 0 −0.5 0.5 z=1m 0 −0.5 0.5 z=1m 0 −0.5 0.5 z=0m 0 −0.5 z=0m 0 −0.5 0 3.75 7.5 11.25 time (sec) 15 0 3.75 7.5 11.25 time (sec) 15 Fig. 19. Time history of horizontal component of acceleration in solid part of the mixture au;x for nodes at different elevations during shaking. (a) Uniform. (b) Layered. Fig. 20. Contours of excess pore pressure Ru, kPa in the soil column at different times. (a) Uniform. (b) Layered. tendency causes faster increase in excess pore pressure and decrease in vertical effective stress and stiffness. Degradation of stiffness leads to faster accumulation of the permanent shear strain in the symmetric loops of stress–strain. The time histories of surface lateral displacement for these two cases are shown in Fig. 26 (results of case III presented in this figure will be discussed ARTICLE IN PRESS M. Taiebat et al. / Soil Dynamics and Earthquake Engineering 30 (2010) 236–257 251 Fig. 21. Contours of excess pore pressure ratio (Ru ) in the soil column at different times. (a) Uniform. (b) Layered. Fig. 22. Illustration of the problem in Example 2 in terms of the soil layering, the finite element mesh, and the input base acceleration. later). It can be observed that the permanent surface lateral displacement in case II ð C 0:40 mÞ is almost twice as much as case I ð C 0:19 mÞ because of the presence of the relatively loose layer. The excess pore pressure ratio in the loose layers shows instantaneous spikes to Ru ¼ 1 and drops mainly to two values of C0:35 or C 0:65, which are related to the end of dilative phases in positive and negative shear stresses, respectively. The pore water expelled from this layer during shaking is pumped to the upper layers and shows a considerable increase in the excess pore pressure ratio of element 19 (E19) after the end of shaking. Fig. 25 shows response of the soil profile in case III. The geometry and parameters in case III are similar to those in case II. The only difference is the intensity of the input base acceleration, which is 0:4g in case III. The interesting observation in the simulation results is that although the magnitude of the base acceleration in case III is twice that of case II, the calculated residual shear strains are not considerably different in the two cases. This could be attributed to the fact that in case III the larger amplitude of acceleration (and shear stress) soon brings the soil state to the dilative phase during which the vertical effective stress increases and soil regains its strength. Comparing the stress path of element 03 (E03) in cases II and III makes this phenomenon more clear. Owing to increase of the vertical effective stress during the intense dilation in case III the effective stress state remains larger in this case compared with case II, and therefore less degradation is observed in the subsequent cycles of loading. The strongly pronounced dilative response in all layers of soil profile in case III clearly affects the resulting excess pore pressure ratios (Fig. 25 c) comparing with what was observed in case II (Fig. 24 c). More specifically, case III experiences more intense drops in plots of Ru even to negative values, especially at upper elements, which are even more dilative because of lower surcharge (confining pressure). The time histories of the surface lateral displacement and variations with depth of the permanent lateral displacement for different cases are presented and compared in Fig. 26. Although case III shows larger amplitudes of cyclic motions, its permanent surface lateral displacement is about 20% smaller than that of case II. That is quite interesting and emphasizes the importance of interaction of dynamic characteristics of the earthquake and the response of the soil (and its components) as described by the SANISAND constitutive model. Fig. 26 b provides more detailed information about the lateral displacements of soil layers in different cases. Owing to the aforementioned stronger dilative response in case III, the dense layers adjacent to the loose layer appear to get stronger in case III and show lesser lateral displacements. This could also be partially explained by smaller surface lateral displacement in case III compared with case II. Note that the permanent displacement is equivalent to integration of shear strains in the whole soil column. 5. Summary and conclusions An efficient finite element formulation and a numerical tool for analysis of wave propagation phenomena in fluid-saturated porous media has been developed. The saturated porous medium is modeled as a two-phase, fully coupled system consisting of a porous solid and a fluid phase. The numerical tool includes the u– p–U formulation for fully coupled behavior of soils and the SANISAND model as an advanced elastic–plastic constitutive model for modeling of stress–strain response of the solid phase. The formulation takes into account velocity proportional ARTICLE IN PRESS 252 M. Taiebat et al. / Soil Dynamics and Earthquake Engineering 30 (2010) 236–257 50 E19 0 0.5 E19 0 −50 50 E15 0.5 E15 0 −50 50 −50 50 E11 0 E07 −50 50 E07 0 E03 −50 0 20 40 60 σz (kPa) 80 100 E03 0 −50 −5 E11 0.5 E07 0 1 −50 50 0 0.5 0 1 −50 50 0 Ru σxz (kPa) E11 E15 0 1 −50 50 0 E19 0 1 −50 50 0 σxz (kPa) 1 50 0.5 0 0 5 10 γ (%) 15 20 0 3 6 9 12 E03 15 time (sec) Fig. 23. Results of response of the uniform soil column with amax ¼ 0:2g (case I) in selected elements along soil profile (E03, E07, E11, E13, E19): (a) shear stress vs. vertical effective stress, (b) shear stress vs. shear strain, (c) excess pore pressure ratio histories, (d) excess pore pressure ratio vs. time and depth, (e) shear strain vs. time and depth, (f) lateral displacement vs. time and depth. damping (usually called viscous damping) by proper modeling of coupling of pore fluid and solid skeleton, while the displacement proportional damping is appropriately modeled by the dissipation mechanism of the constitutive model of elasto-plasticity used. Numerical results that demonstrate the accuracy and versatility of the formulations and capabilities of the model are presented. The verified and validated models and computational simulations tool are used to predict behavior of a layered soil system during seismic loading. Two general problems are studied. Firstly, a numerical study is conducted on isolating effects of a liquefied sand layer in propagation of seismic waves. Secondly, seismically induced shear deformation of a gentle slope in the presence or absence of liquefiable layer is studied. Attention is given to propagation of shear waves through soil layers and to variation of stresses and strains during such propagation. Spatial evolution of liquefaction in loose and medium dense soil layers is also investigated. In the first example liquefaction of loose layers in depth is shown to prevent transmission of earthquake-induced motions and shear stresses to the upper layers. While prevention of transmission of motions may be considered as a positive feature, the water flux coming from the liquefied bottom layers might force the upper layers towards lower effective stress and possibly instability, even in the absence of earthquake-induced shear stresses. More importantly, attention needs to be paid to the fact that liquefaction of the deep loose layers could result in large strains in these layers, which might result in (rigid block) translation of the upper layers. The second example is used to present complexities arising from variation of soil density in a layered system as well as the input base acceleration. For example, it is shown that larger shaking produces smaller lateral spread, because of close interaction of earthquake loading and soil constitutive response (with its solid and fluid components). The capabilities of the u–p–U element, the SANISAND constitutive model, and the nonlinear dynamic finite element numerical tool are demonstrated in handling complex phenomena. These studies provide new insight into the mechanisms of wave propagation and ARTICLE IN PRESS 253 M. Taiebat et al. / Soil Dynamics and Earthquake Engineering 30 (2010) 236–257 50 E19 0 0.5 E19 0 −50 50 E15 0.5 E15 0 −50 50 −50 50 E11 0 E07 −50 50 E07 0 E03 −50 0 20 40 60 σz (kPa) 80 100 E03 0 −50 −5 E11 0.5 E07 0 1 −50 50 0 0.5 0 1 −50 50 0 Ru σxz (kPa) E11 E15 0 1 −50 50 0 E19 0 1 −50 50 0 σxz (kPa) 1 50 0.5 0 0 5 10 γ (%) 15 20 0 3 6 9 time (sec) 12 E03 15 Fig. 24. Results of response of the layered soil column with amax ¼ 0:2g (case II) in selected elements along soil profile (E03, E07, E11, E13, E19): (a) shear stress vs. vertical effective stress, (b) shear stress vs. shear strain, (c) excess pore pressure ratio histories, (d) excess pore pressure ratio vs. time and depth, (e) shear strain vs. time and depth, (f) lateral displacement vs. time and depth. seismic behavior of saturated soil systems. Fully coupled nonlinear dynamic numerical simulations provide accurate means for detailed investigation of complex loading conditions and geometries and different material states. Such high fidelity simulations provide detailed input to PBD methods that most existing empirical approaches used in liquefaction modeling fail to provide. In particular one must emphasize the importance of using a robust and realistic constitutive model such as SANISAND that can account for the various important features of the soil response in relation to its density and loading conditions. For example, it is important that depending on the level of density and confining pressure the constitutive model with a single set of parameters can show the dilative or contractive response. This is accomplished by using the state parameter c. This parameter in the model controls the peak stress ratio and by consequence the plastic modulus, as well as the dilatancy stress ratio (phase transformation stress ratio) and by consequence the dilatancy. Similarly, by using a fabric-dilatancy tensor the model has the capability to account for intense contractive response upon reverse loading increments after a dilative phase of loading. This feature has an important effect on the predicted displacement and pore water pressure built up during cyclic loading. Had these constitutive properties not been realistically described, the overall description of the phenomenon would have been much poorer. Appendix A. Constitutive platform The main components of the SANISAND model are summarized here. Both stress and strain quantities are assumed positive in compression (as is common in mechanics of materials), and the effect of this sign convention has been considered on the model equations. All stress components in this paper should be considered as effective stress. Finally, in terms of notation, tensor quantities are denoted by bold-faced symbols and operations explained accordingly. ARTICLE IN PRESS 254 M. Taiebat et al. / Soil Dynamics and Earthquake Engineering 30 (2010) 236–257 0.5 E19 0 E19 0 E11 0 E07 0 E07 0 E03 0 −50 0 20 40 60 80 σz (kPa) 100 E11 0.5 E07 0 1 −50 50 −50 50 0.5 0 1 −50 50 −50 50 E15 0 1 −50 50 Ru E11 0 σxz (kPa) −50 50 0.5 E15 0 E15 0 E19 0 1 −50 50 −50 50 σxz (kPa) 1 50 50 E03 0 −50 −5 0.5 0 0 5 10 γ (%) 15 20 0 3 6 9 12 time (sec) E03 15 Fig. 25. Results of response of the layered soil column with amax ¼ 0:4g (case III) in selected elements along soil profile (E03, E07, E11, E13, E19): (a) shear stress vs. vertical effective stress, (b) shear stress vs. shear strain, (c) excess pore pressure ratio histories, (d) excess pore pressure ratio vs. time and depth, (e) shear strain vs. time and depth, (f) lateral displacement vs. time and depth. 0 0.4 II. Layered, amax=0.2g III. Layered, amax=0.4g 0.2 I. Uniform, amax=0.2g 0.1 III II Loose Layer Depth (m) ux (m) 0.3 Case I −2 −4 −6 −8 0 −10 0 3 6 9 time (sec) 12 15 0 0.1 0.2 0.3 permanent ux (m) [@15 sec] 0.4 Fig. 26. (a) Time histories of surface lateral displacement, and (b) variations with depth of the permanent lateral displacement for cases I–III. ARTICLE IN PRESS M. Taiebat et al. / Soil Dynamics and Earthquake Engineering 30 (2010) 236–257 Ad ¼ A0 ð1 þ/z : nSÞ Elasticity: Isotropic hypoelasticity is defined by s_ e e_ ¼ ; 2G p_ e_ ev ¼  K ð1Þ with s is the stress deviator; p the mean effective pressure; e the strain deviator; ev the volumetric strain; superscript e denoting elastic; superposed dot denoting the rate; and G, K the hypoelastic shear and bulk moduli, respectively, given by   ð2:97eÞ2 p 1=2 2ð1 þ nÞ G ¼ G0 pat G ð2Þ ; K¼ ð1þ eÞ pat 3ð12nÞ where G0 is a dimensionless material constant, n is a constant Poisson’s ratio, e is the void ratio, and pat is the atmospheric pressure used for normalization. Yield surface: The yield surface is defined by pffiffiffiffiffiffiffiffiffi ð3Þ f ¼ ½ðspaÞ : ðspaÞ1=2  2=3mp ¼ 0 in terms of a deviatoric back-stress ratio a, and can be visualized as a cone. The gradient of f ¼ 0 in stress space is obtained based on Eq. (3) as @f 1 ¼ n ðn : rÞI; @r 2 ra n ¼ pffiffiffiffiffiffiffiffiffi 2=3m ð4Þ in terms of deviator unit tensor n ðn : n ¼ 1Þ shown along the radius ra of the yield surface circular trace, where I is the identity second-order tensor. The n of Eq. (4) is used to define an effective Lode angle y as pffiffiffi cos3y ¼  6 tr n3 ð5Þ The critical stress ratio for a given Lode angle y is denoted by M and has been obtained by interpolation between values of Mc at y ¼ 0 (compression) and Me at y ¼ p=3 (extension) according to M ¼ Mc gðy; cÞ; gðy; cÞ ¼ 2c ; ð1 þcÞð1cÞcos3y c¼ Me Mc ð6Þ Bounding, dilatancy and critical surfaces: Three concentric and homologous surfaces, the bounding, dilatancy, and critical are considered for the model in the pplane. Their tensor points aby , ady , and acy along the direction n emanating from the origin at a Lode angle y are defined analytically by pffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffi aby ¼ 2=3½Mexpðnb cÞmn ¼ 2=3aby n ð7aÞ pffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffi ady ¼ 2=3½Mexpðnd cÞmn ¼ 2=3ady n acy ¼ pffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffi 2=3½Mmn ¼ 2=3acy n ð7bÞ ð7cÞ x where c ¼ eec is the state parameter; ec ¼ e0 lðpc =pat Þ as the analytical expression of the critical state line with constants e0 , l and x [35]; and nb and nd are material constants. Plastic flow: In general the model employs a non-associative flow rule. The plastic strain rate direction R has deviatoric and volumetric part according to R ¼ R0 ð1=3ÞDI, where D is the dilatancy coefficient to be specified in the sequel. The deviatoric part R is defined as the normal to the critical surface at point acy , Eq. (7c). The plastic strain is given by e_ p ¼ /LSR ð8aÞ R ¼ Bnþ Cðn2 13 IÞ13DI ð8bÞ B ¼ 1þ 3 1c gcos3y; 2 c rffiffiffi 3 1c g 2 c C¼3 ð8cÞ where L is the loading index; and / S the Macauley brackets (/LS ¼ L if L 40 and /LS ¼ 0 if L r0). The dilatancy coefficient D in Eq. (8a) is defined by D ¼ Ad ðady aÞ : n ð9aÞ 255 ð9bÞ where A0 is a material constant. The so-called fabric-dilatancy internal variable z is introduced in order to account for the effect of fabric changes during the dilatant phase of deformation on the subsequent contractant response upon load increment reversals. The evolution law for the fabric-dilatancy tensor z will be presented later. Evolution laws: This model has two tensorial internal variables, namely, the back-stress ratio tensor a and the fabric-dilatancy tensor z. The evolution law for the back-stress ratio a is a function of the distance between bounding and current back-stress ratio in the form of a_ ¼ /LS23hðaby aÞ ð10Þ While constant values of h yield reasonable simulations, better results are obtained if h is made function of the current state variables as h¼ b0 ðaa~ Þ : n ð11aÞ b0 ¼ G0 h0 ð1ch eÞðp=pat Þ1=2 ð11bÞ where a~ is the initial value of a at initiation of a new loading process and is updated to the new value when the denominator of Eq. (11a) becomes negative. The h0 and ch are material constants. Finally the evolution law for the fabric-dilatancy z is introduced as z_ ¼ cz /LDSðzmax nþ zÞ ð12Þ with cz and zmax as the material constants that control the maximum value of z and its pace of evolution. Eqs. (9b) and (12) are introduced to enhance the subsequent contraction in unloading after the dilation in loading in order to generate the wellknown butterfly shape in the undrained cyclic stress path. Appendix B. Finite element formulation The overall equilibrium or momentum balance equation for the soil–fluid mixture can be written as sij;j ru€ i rf w€ i þ rbi ¼ 0 ð13Þ where u€ i is the acceleration of the solid part, bi is the body force € i is the fluid acceleration relative to the solid part. per unit mass, w For fully saturated porous media (no air trapped inside), density is equal to r ¼ nrf þ ð1nÞrs , where n is the porosity, rs and rf are the soil particle and water density, respectively. For the pore fluid, the equation of momentum balance can be written as € i p;i Ri rf w rf w€ i n þ rf bi ¼ 0 ð14Þ where R is the viscous drag force. According to Darcy’s seepage law, the viscous drag forces R between soil matrix and pore fluid _ (water) can be written as Ri ¼ k1 ij w j , where kij is the tensor of anisotropic Darcy permeability coefficients. For simple case of isotropic permeability, scalar value of permeability k is used. The permeability k used here with dimension of ½L3 TM1  is different from the permeability used in the usual soil mechanics (K) which has the same dimension of velocity, i.e., ½LT 1 . Their values are related by k ¼ K=g rf , where g is the gravitational acceleration and the permeability K is measured in an experiment. The final equation is the mass conservation of the fluid flow expressed by _ i;i þ ae_ ii þ w p_ ¼0 Q ð15Þ ARTICLE IN PRESS 256 M. Taiebat et al. / Soil Dynamics and Earthquake Engineering 30 (2010) 236–257 where bulk stiffness of the mixture Q is expressed as 1=Q ¼ n=Kf þ ðanÞ=Ks and Ks and Kf are the bulk moduli of the solid and fluid phases, respectively. In the above governing equations, convective and terms of lower order are omitted [38]. A change of variable is performed by introducing an alternative variable Ui , defined as Ui ¼ ui þ UiR ¼ ui þwi =n, that represents absolute displacement of the pore fluid. The basic set of unknowns is then comprised of the soil skeleton displacements ui , the water pore pressure p, and the water displacements Ui . These unknown vector (ui and Ui ) and scalar (p) fields can be approximated using shape functions and the corresponding (unknown) nodal values (solid displacements u Ki , pore fluid pressures p K , and pore fluid displacements U Ki ), following standard finite element discretization [39,40] p ¼ NKp p K ; ui ¼ NKu u Ki ; Ui ¼ NKU U Ki ð16Þ where NKu , NKp and NKU are (same) shape functions for solid displacement, pore pressure and fluid displacement, respectively. Each node of the ðu2p2UÞ element has thus seven degrees of freedoms in three dimensions (three for solid displacements, one for pore fluid pressures and three or pore fluid displacements). After some tensor algebra and manipulations, the final u2p2U form can be written as 2 32 € 3 0 ðMs ÞKijL 0 u Lj 7 6 0 76 € 7 0 0 4 56 4 pN 5 0 0 ðMf ÞKijL U€ Lj 32 _ 3 2 ðC1 ÞKijL 0 ðC2 ÞKijL u Lj 7 76 6 _ 7 0 0 0 þ4 56 4 pN 5 ðC2 ÞLjiK 0 ðC3 ÞKijL U_ Lj 2 32 3 u Lj 0 ðG1 ÞKiM 0 6 6 7 PMN ðG2 ÞLjM 7 þ 4 ðG1 ÞLjM 54 p M 5 U Lj 0 ðG2 ÞKiL 0 2R u 0 3 2 u 3 f Ki O NK;j sij dO 7 6 7 6 6 ð17Þ þ4 ¼ 0 5 4 0 7 5 U 0 f Ki The left hand side components of the above matrix equation are given as Z Z ðMs ÞKijL ¼ NKu ð1nÞrs dij NLu dO; ðMf ÞKijL ¼ NKU nrf dij NLU dO O ðC1 ÞKijL ¼ Z ðC3 ÞKijL ¼ Z ðG2 ÞKiM ¼ Z O O O Z u NKu n2 k1 ij NL dO; ðC2 ÞKijL ¼ U NKU n2 k1 ij NL dO; ðG1 ÞKiM ¼ O p nN U K;i NM dO; PNM ¼ Z O NNp O Z U NKu n2 k1 ij NL dO O p u NK;i ðanÞNM dO 1 p N dO Q M ð18Þ while the right hand side components are given as Z Z Z ðf s ÞKi ¼ NKu sij0 nj dG NKu ðanÞpni dG þ NKu ð1nÞrs bi dO Gt Gp O Z Z ðf f ÞKi ¼  NKU npni dG þ NKU nrf bi dO Gp O ð19Þ References [1] Kramer SL. Geotechnical earthquake engineering. Upper Saddle River, NJ: Prentice-Hall; 1996. [2] Martin GR, Seed HB, Finn WDL. Fundamentals of liquefaction under cyclic loading. Journal of the Geotechnical Engineering Division 1975;101(5): 423–38. [3] Finn WDL, Martin GR, Lee KW. An effective stress model for liquefaction. Journal of the Geotechnical Engineering Division 1977;103(6):517–33. [4] Zienkiewicz OC, Leung KH, Pastor M. Simple model for transient soil loading in earthquake analysis. I: basic model and its application. International Journal for Numerical and Analytical Methods in Geomechanics 1985;9(5): 453–76. [5] Prevost JH. A simple plasticity theory for frictional cohesionless soils. Soil Dynamics and Earthquake Engineering 1985;4(1):9–17. [6] Iai S, Matsunaga Y, Kameoka T. Strain space plasticity model for cyclic mobility. Soils and Foundations 1992;32(2):1–15. [7] Manzari MT, Dafalias YF. A critical state two-surface plasticity model for sands. Géotechnique 1997;47(2):255–72. [8] Dafalias Y. Overview of constitutive models used in VELACS. In: Proceedings of the international conference on the verification of numerical procedures for the analysis of soil liquefaction problems. Rotterdam, Netherlands: Balkema; 1994. p. 1293–303. [9] Elgamal A, Yang Z, Parra E, Ragheb A. Modeling of cyclic mobility in saturated cohesionless soils. International Journal of Plasticity 2003;19(6):883–905. [10] Dafalias YF, Manzari MT. Simple plasticity sand model accounting for fabric change effects. Journal of Engineering Mechanics 2004;130(6):622–34. [11] Wood DM, Belkheir K, Liu DF. Strain softening and state parameter for sand modeling. Géotechnique 1994;42(2):335–9. [12] Dafalias YF, Papadimitriou AG, Li XS. Sand plasticity model accounting for inherent fabric anisotropy. Journal of Engineering Mechanics 2004;130(11): 1319–33. [13] Pestana JM, Whittle AJ. Compression model for cohesionless soils. Géotechnique 1995;45(4):611–31. [14] Taiebat M, Dafalias YF. SANISAND: simple anisotropic sand plasticity model. International Journal for Numerical and Analytical Methods in Geomechanics 2008;32(8):915–48. [15] Biot MA. Theory of propagation of elastic waves in a fluid–saturated porous solid. The Journal of the Acoustical Society of America 1956;28:168–78. [16] Zienkiewicz OC, Shiomi T. Dynamic behaviour of saturated porous media; the generalized Biot formulation and its numerical solution. International Journal for Numerical Methods in Engineering 1984;8:71–96. [17] Jeremić B, Cheng Z, Taiebat M, Dafalias YF. Numerical simulation of fully saturated porous materials. International Journal for Numerical and Analytical Methods in Geomechanics 2008;32(13):1635–60. [18] McKenna FT. Object oriented finite element programming: Framework for analysis, algorithms and parallel computing. PhD thesis, Department of Civil and Environmental Engineering, University of California, Berkeley, CA, USA; 1997. [19] Cheng Z. Computational solid and coupled nonlinear geomechanics in small and large deformation. PhD thesis, Department of Civil and Environmental Engineering, University of California, Davis, CA, USA; 2006. [20] Jeremić B, Sture S. Tensor data objects in finite element programming. International Journal for Numerical Methods in Engineering 1998;41: 113–26. [21] Davis TA, Duff IS. A combined unifrontal/multifrontal method for unsymmetric sparse matrices. Technical Report TR97-016 (UofF) and TR-97-046 (RAL), University of Florida and Rutherford Appleton Laboratory; 1997. [22] Taiebat M. Advanced elastic–plastic constitutive and numerical modeling in geomechanics. PhD thesis, Department of Civil and Environmental Engineering, University of California, Davis, CA, USA; 2008. [23] Sloan SW, Abbo AJ, Sheng D. Refined explicit integration of elastoplastic models with automatic error control. Engineering Computations 2001;18(1/2):121–54. [24] Argyris J, Mlejnek H-P. Dynamics of structures. Amsterdam, USA: Elsevier, North-Holland; 1991. [25] Hilber HM, Hughes TJR, Taylor RL. Improved numerical dissipation for time integration algorithms in structural dynamics. Earthquake Engineering and Structure Dynamics 1977;5(3):283–92. [26] Hughes TJR, Liu WK. Implicit-explicit finite elements in transient analysis: implementation and numerical examples. Journal of Applied Mechanics 1978;45:375–8. [27] Hughes TJR, Liu WK. Implicit–explicit finite elements in transient analysis: stability theory. Journal of Applied Mechanics 1978;45:371–4. [28] Oberkampf WL, Trucano TG. Hirsch C. Verification, validation and predictive capability in computational engineering and physics. In: Proceedings of the foundations for verification and validation on the 21st century workshop. Laurel, MD: Johns Hopkins University, Applied Physics Laboratory; 2002. p. 1–74. [29] Gajo A. Influence of viscous coupling in propagation of elastic waves in saturated soil. ASCE Journal of Geotechnical Engineering 1995;121(9): 636–44. [30] Gajo A, Mongiovi L. An analytical solution for the transient response of saturated linear elastic porous media. International Journal for Numerical and Analytical Methods in Geomechanics 1995;19(6):399–413. [31] Verdugo R, Ishihara K. The steady state of sandy soils. Soils and Foundations 1996;36(2):81–91. [32] Pradhan TB, Tatsuoka F, Sato Y. Experimental stress-dilatancy relations of sand subjected to cyclic loading. Soils and Foundations 1989;29(1): 45–64. ARTICLE IN PRESS M. Taiebat et al. / Soil Dynamics and Earthquake Engineering 30 (2010) 236–257 [33] Arulmoli K, Muraleetharan KK, Hosain MM, Fruth LS. VELACS laboratory testing program, soil data report. Technical Report 90-0562, The Earth Technology Corporation, Irvine, CA, report to the National Science Foundation, Washington, DC; March 1992. [34] Lee KL, Seed HB. Drained strength characteristics of sands. Journal of the Soil Mechanics and Foundations Division, Proceedings of the American Society of Civil Engineers 1967;93(SM6):117–41. [35] Li X-S, Wang Y. Linear representation of steady-state line for sand. Journal of Geotechnical and Geoenvironmental Engineering 1998;124(12):1215–7. [36] Arulanandan K, Scott RF. (Eds.), Proceedings of the international conference on the verification of numerical procedures for the analysis [37] [38] [39] [40] 257 of soil liquefaction problems. Rotterdam, Netherlands: A. A. Balkema; 1993. Schlesinger S. Terminology for model credibility. Simulation 1979;32(9): 103–4. Zienkiewicz OC, Chan AHC, Pastor M, Schrefler BA, Shiomi T. Computational geomechanics with special reference to earthquake engineering. New York: Wiley; 1999. Zienkiewicz OC, Taylor RL. The finite element method, vol. 1, 4th ed. New York: McGraw-Hill Book Company; 1991. Zienkiewicz OC, Taylor RL. The finite element method, vol. 2, 4th ed. New York: McGraw-Hill Book Company; 1991.