Three-dimensional seismic analysis of submarine slopes

ARTICLE IN PRESS

Soil Dynamics and Earthquake Engineering 26 (2006) 870–887 www.elsevier.com/locate/soildyn

Three-dimensional seismic analysis of submarine slopes A. Aziziana,, R. Popescub a

EBA Engineering Consultants, 9th Floor, 1066 West Hastings Street, Vancouver, BC, Canada V6E 3X2 Faculty of Engineering & Applied Science, Memorial University of Newfoundland, St. John’s, NL, Canada A1B 3X5

b

Accepted 15 October 2005

Abstract Three-dimensional effects in seismic analysis of submarine slopes are assessed by comparing results of two- and three-dimensional (2D & 3D) analyses, in terms of predicted displacements, shear strains, and excess pore water pressure ratios. Limits of applicability of the 2D, plane strain analysis assumptions are quantitatively assessed. Some regression equations are also presented that express ratios of 3D vs. 2D predictions as a function of slope width/height ratio and earthquake peak acceleration. It is found that the results of 2D and 3D dynamic slope stability analysis are within a tolerance of about 15% for width/height ratios larger than 3–5, and 3D effects induced by lateral boundaries become insignificant for width/height ratios larger than 6–7. The results of the present dynamic, fully coupled, nonlinear analyses are also compared with those of static slope stability analyses. r 2006 Elsevier Ltd. All rights reserved.

1. Introduction In the framework of slope stability analysis by factor-ofsafety approach, a considerable number of studies have compared the factor of safety obtained from a full 3D analysis (F3D) with the one obtained from a simplified 2D analysis (F2D). The majority of those studies have concluded that, for a given slope, F2DoF3D [1]. These studies typically account for the difference between the 3D and 2D shapes of the failure surface (Fig. 1), and they show that with increasing width of the failing soil wedge assumed in a 3D analysis, F3D converges to F2D. Due to the limitations of the limit equilibrium method when applied to dynamic analysis, results of such studies are not applicable to the seismic analysis of submerged slopes, where soil liquefaction and other effects of dynamic earthquake loading are significant. In this study, a dynamic, fully coupled, finite element analysis implementing a multi-yield surface plasticity model (Prevost [2,3]) is used to compare the behavior of slopes with various width to height ratios and subjected to a range of acceleration levels, as predicted by 3D analyses, with the predicted behaviour of their equivalent 2D Corresponding author. Tel.: +1 604 685 0275; fax: +1 604 684 6241.

E-mail address: [email protected] (A. Azizian). 0267-7261/$ - see front matter r 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.soildyn.2005.10.008

configurations. The slope behavior is indicated by three representative responses related to displacement, shear strain, and excess pore pressure ratio. The most important objective of this study is to set reasonable limits for the applicability of the two-dimensional (2D), plane strain simplifying assumption in seismic analysis of slopes. By quantifying the 3D/2D ratios of some predicted quantities related to seismic response of slopes, the study provides means of extrapolating 2D analysis results to real 3D situations. 2. Review of 3D vs. 2D analyses Justo and Saura [4], Byrne et al. [5], Stark and Eid [6], Jeremic [7], among many others, have studied the importance of considering the 3D effects in the analysis of slopes. Duncan [1] cites over 20 studies (since 1960s) concluding that in the factor-of-safety approach, with a few exceptions, two-dimensional analysis yields conservative results compared to three-dimensional analysis, i.e. F2DoF3D provided that F2D is calculated for the most critical 2D section. Since F2DoF3D, Leshchinsky and Huang [8] and Duncan and Stark [9] emphasize that in order to obtain post-failure in situ shear strength of soils by back-analysis of slope failures, 3D analysis should be avoided so that

ARTICLE IN PRESS A. Azizian, R. Popescu / Soil Dynamics and Earthquake Engineering 26 (2006) 870–887

3D (finite curved surface)

2D (infinite cylindrical surface) Fig. 1. 3D (finite curved) versus 2D (infinite cylindrical) assumptions of the shape of failure surface.

shear strength (to be used in 2D analysis) is not overestimated. Arellano and Stark [10] present several curves showing the ratio of 3D to 2D factor of safety for different width/ height ratios and slope angles, with translational failure mechanism. These curves show that for a 3H:1 V slope, for example, the 3D factor of safety is about 40% larger than the 2D factor of safety if the width/height ratio is about 2. In this type of failure mechanism, which usually occurs in flatter slopes (see also [6]), the mobilized shear strength along back scarp and sides of the slide mass are significantly different from those along the base. In deformation analysis approach, Lefebvre et al. [11] found that 2D analysis can significantly overestimate movements of a dam in a V-shaped steep-wall valley because the effects of cross-valley arching are ignored in 2D analysis. Comparing the results of 2D and 3D analyses of the transverse sections of the dam showed that if the valley walls slopes were 1:1 or steeper, plane strain results would be significantly inaccurate. 2D/3D ratios of principal stresses, maximum shear stress, and displacements were also presented. For example, average 2D/3D ratio of horizontal displacements for 1:1 wall was 2.68, whereas that for 6:1 wall was 1.05. Some other studies, such as [1,4,12] also compare the results of 3D and 2D finite element analyses for dams (simulating reservoir filling, water rise, etc.) and generally indicate that deformations in 3D models are significantly smaller than in 2D models, for width/height ratios of about 2–3. From these studies, it is also possible to note that beyond a certain distance from the abutments, there is no significant variation in the predicted crest displacement. This distance is normally about 2.5–3 times the height. Prevost et al. [13] performed 2D and 3D total stress, dynamic FE analysis of Santa Felicia earth dam using a non-linear multi-surface plasticity model and compared the measured and computed earthquake responses. Their analysis showed that the first 10 frequencies of the 3D model all fall within the first two frequencies of the 2D

871

model, indicating that more intermediate modal configurations are generated, despite the fact that the dam is a relatively long dam and 3D effects should not be significant. The study demonstrates the importance of 3D effects being more pronounced for strong shaking, in terms of crest acceleration and permanent deformations. In case of strong shaking, the 3D horizontal crest acceleration response is significantly lower than the 2D one due to significant contributions from higher modes of vibration. The two-dimensional idealization is normally based on the following considerations: 1. Site material idealization: in most cases, it is assumed that soil layering is perpendicular to the plane of interest and a cross section represents all sections. If nonhomogeneity or anisotropy of the slope material is important, then performing a 2D analysis is not appropriate. 2. Site geometry idealization: in numerical or analytical approaches toward many geotechnical problems, three main geometric idealizations are used to simplify and speed-up the analysis significantly (see, e.g., [14]): plane strain, plane stress, and axi-symmetry. Almost all twodimensional slope stability analyses use plane strain assumption, in which the value of the strain component perpendicular to the plane of interest is zero. Analysis time as well as necessary computational resources will decrease significantly, especially in a seismic step-by-step time-domain analysis. However, the plane strain assumption is valid if one dimension is very large in comparison with the other two. It also requires all of the following: no curvature or corners exist in geometry of slope; the slope deformation is not constrained significantly by a near lateral boundary (such as a dam in a narrow rock-walled valley); no curvature exists in the shape of (postulated) failure surface in the direction perpendicular to the plane of interest (that is, failure surface is the same in any cross section). 3. Loading idealization: design recommendations usually suggest that from three components of earthquake acceleration, only the larger horizontal acceleration would suffice for analysis purposes (see, e.g., [15]). Vertical and smaller horizontal components of acceleration are then ignored. In the analysis of slopes, the horizontal acceleration is usually applied to the crosssection of the slope in its plane, that is, no instability due to transverse (out-of-plane) excitation is taken into account.

3. Numerical model 3.1. Finite element code A multi-yield surface plasticity constitutive model implemented in the finite element program, Dynaflow [3], is used in this study. Dynaflow is a state-of-the-art general-

ARTICLE IN PRESS A. Azizian, R. Popescu / Soil Dynamics and Earthquake Engineering 26 (2006) 870–887

Fig. 2. The finite element meshes used in this study: (a) 3D; (b) 2D coarse; and (c) 2D fine meshes.

000

0

z=0

y z

3

f ac e : 311

x

10 100 Sym

Corners: 311011 001

3110 11

metr yf

a ce:

Base : 310 0010

3110 11

L= 1

40 m

001

0 k: 3 Bac

Legend 0: Active degree-of-freedom 1: Prescribed displacement 3. Prescribed acceleration

010

01

F

t: 3 ron

000

00

310

311

le ri ab Va

0 01 301

: B/2

0 01 3 11

purpose finite element analysis program for linear and nonlinear, two- and three-dimensional systems. Coupled analysis of porous media is performed for simulating the behavior of saturated soil materials under partially drained conditions by means of an extension of Biot’s formulation [16] in the non-linear regime [17]. The multi-yield plasticity model is a kinematic hardening model based on a relatively simple plasticity theory [2] and is applicable to both cohesive and cohesionless soils. The concept of a ‘‘field of work-hardening moduli’’ is used by defining a collection of nested yield surfaces in the stress space. Rounded Mohr–Coulomb yield surfaces are employed in this study. The yield surfaces define regions of constant shear moduli in the stress space, and in this manner, the model discretizes the smooth elastic–plastic stress–strain curve into a number of linear segments. The outermost surface, i.e. failure surface, corresponds to zero shear modulus. The plastic flow rule is associative in its deviator component. A non-associative flow rule is used for the dilatational component to reflect the dependency of plastic dilatancy on effective stress ratio. The material hysteretic behavior and shear stress-induced anisotropic effects are simulated by a purely kinematic hardening rule. Altogether, accurate simulation of shear-induced plastic dilation, hysteretic effects, and full solid–fluid coupling leads to a more precise calculation of excess pore water pressure (EPWP) build-up/dissipation as well as modelling gradual soil softening/hardening during and after the earthquake shaking. The required constitutive model parameters can be derived from the results of conventional laboratory and/ or in situ soil tests and liquefaction strength analysis (see e.g. [18–20]). The soil constitutive model, its implementation algorithm, and the methodology for estimating the constitutive model parameters have been repeatedly validated in the past for soil liquefaction computations, based on both full scale measurements and centrifuge experimental results (e.g [21,18,19,22–24]).

0 01 30 1

872

010

Fig. 3. Boundary conditions of the 3D model. The boundary condition codes are according to the definitions given in Table 1, with the following order: x, y, z components of solid phase displacements, and x, y, z components of fluid phase displacements.

3.2. Finite element mesh The general scheme of the finite element meshes used in this study is shown in Fig. 2. In addition to the variation of the model width (B), the effect of the finite element discretization (mesh refinement) on the results is also investigated. The cross section of the 3D mesh is either coarse or fine (Figs. 2b and c). The width of elements in the z-direction is either 5 or 2 m depending on the model width. The slope is 1V:3H, however, through a limited number of analyses, the effect of slope angle changes is also addressed. This slope angle is selected according to the outline of the centrifuge experiments to be carried out in the framework of the COSTA-Canada project [25,26]. For the full 3D analyses, due to symmetry of the structure and loading with respect to a plane normal to the z-axis (see Fig. 3), only half of the structure is analyzed.

The soil is discretized into eight-node, bilinear, brick finite elements with six degrees of freedom per node, three for solid and three for fluid displacements. Depending on the model width (B) and width of elements in the zdirection, the total number of elements in 3D analysis varies from 2 to 20 times the number of elements in 2D analysis. For the 2D plane–strain analyses, the soil is discretized into four-node, bilinear, quadrilateral finite elements with four degrees of freedom per node, corresponding to vertical and horizontal components of solid and fluid displacements. The 2D coarse and fine meshes have 156 and 372 elements, respectively. The numbers of equations for the 2D coarse and fine meshes are 659 and 1541, respectively. For the 3D mesh, the ranges of the number of equations are: (1) from 2045 to 20819, in case of

ARTICLE IN PRESS A. Azizian, R. Popescu / Soil Dynamics and Earthquake Engineering 26 (2006) 870–887

coarse mesh ðB=H ¼ 0:8220Þ; and (2) from 7118 to 28655 in case of fine mesh ðB=H ¼ 10240Þ.

The analysis domain has an idealized shape, especially with regard to the lateral boundaries perpendicular to the face of the slope as they are selected to be vertical and rigid to model presence of a much stiffer soil. A realistic model would have inclined lateral boundaries (V-shaped valleys) and interface elements with Coulomb friction. In this study, it was first selected to include the most restrictive boundary condition, which is vertical walls and stick contacts. As it will be discussed later, the effect of the vertical rigid boundary is quite significant in the cases analyzed here, which necessitated the use of transmitting boundary conditions. The boundary conditions used in 3D analyses are shown in Fig. 3. The boundary condition codes are according to the definitions given in Table 1 [3], with the following order: x, y, z components of solid phase motion, and x, y, z components of fluid phase motion. The input acceleration is not applied to the fluid phase on the front and back faces of the model (Fig. 3) to avoid numerical noise in the response. To simulate impervious boundaries on the front and back faces, the horizontal component of fluid displacement in the x-direction (i.e. the fourth degree-offreedom) of each node located on those faces is slaved to that of the corresponding node located on the very next vertical plane of nodes. The nodes located at the corners of the model (shown by solid circles Fig. 3) have also this type of slaved degrees-of-freedom. As rigid boundaries are assumed for the analysis domain, the earthquake acceleration is applied at the base of the model (y ¼ 0 face) as well as all lateral boundaries, namely x ¼ 0; L, and z ¼ 0 faces. The earthquake motion is an acceleration time history generated according to the response spectrum for soil type 3 (soft to medium clays and sands) recommended by the Uniform Building Code [27]. Fig. 4 shows the time history scaled to amax ¼ 0:3 g. In this study, amax varies from 0.1 to 0.5 g. The analysis time is 20 s with time step set to 0.01 s. The dynamic analysis of the present coupled porous continuum, according to Biot’s theory of porous media [16], is treated as a Hyperbolic Initial Boundary Value problem. For time integration, the Newmark algorithm is selected, which for hyperbolic analysis requires a
0:5. Also, as the step-by-step implicit analysis is chosen, the

Table 1 Description of boundary condition codes [3] Boundary condition code

Description

0 1 3

Unspecified (i.e., active degree of freedom) Prescribed displacement Prescribed acceleration

0.4 0.3 Acceleration (g)

3.3. Boundary conditions and earthquake loading

873

0.2 0.1 0 -0.1 -0.2 -0.3 -0.4

0

5

10

15

Time (s) Fig. 4. Synthesized acceleration time history used in this study ðamax ¼ 0:3 gÞ compatible to the response spectrum recommended for soft to medium soils (Type 3—uniform building code [27]).

time step size is dictated only by the sampling time step used for the synthesized input motion. In addition, for this type of analysis the Newmark parameter b is set internally to ða þ 0:5Þ2 =4. This type of implicit time integration is unconditionally stable and does not impose any restrictions on the time step size. For non-linear iterations the quasiNewton BFGS method (see [3] for details) is selected because it provides faster solutions compared to other methods and it is sufficiently accurate. 3.4. Material properties The material properties used in this study (Table 2) are estimated as described in detail by Azizian [41] based on values reported in the literature for the Fraser River sand (BC, Canada) at Dr ¼ 40%. The properties are updated and modified slightly after obtaining the results of a series of cyclic simple shear tests performed at the University of British Columbia (UBC), within the framework of ‘‘Earthquake Induced Damage Mitigation from Soil Liquefaction’’ project, (http://www.civil.ubc.ca/liquefaction/). The procedure for updating the material properties, as explained by Azizian [41] is based on liquefaction strength analysis and the calibrated parameters are calculated using a statistical tool (Response Surface Methodology, see e.g. [42]–[44]). The updated properties are also used in a limited number of 3D analyses. Hereafter, the properties inferred from the literature are referred to as set #1, and the modified properties are referred to as set #2 (see Table 2). In Table 2, the values of low-strain shear modulus correspond to a reference mean effective normal stress p0 ¼ 100 kPa. The fluid is assumed compressible with a bulk modulus of 2000 MPa (e.g. [45]). All model parameters can be estimated using the results of traditional laboratory soil tests, excepting for the dilation parameter (Xpp), which is calibrated by trial-anderror to match the predicted liquefaction response of the sand with the observed laboratory response. This is done by fitting the predicted and observed liquefaction strength curves in terms of cyclic stress ratio vs. the number of

ARTICLE IN PRESS A. Azizian, R. Popescu / Soil Dynamics and Earthquake Engineering 26 (2006) 870–887

874

cycles to liquefaction (Fig. 5). The predicted response is calculated by simulating the laboratory tests using element tests, which consists of solving the constitutive relationships of the material model with known values of parameters, subject to a particular mechanism of cyclic loading with known magnitude and frequency, starting at a certain state of stresses (such as isotropically consolidated), and subject to certain conditions of drainage. The element tests do not involve any finite element mesh or boundary conditions and the undrained condition is simply imposed by assuming a very large (penalty) value for fluid bulk modulus.

Table 2 Material properties of Fraser River sand Soil property

Symbol

Set #1 Ref. (Literature)

Set #2 (UBC tests)

Mass density—solid (kg/m3) Porosity Low-strain shear modulus (MPa) Poisson’s ratio Power exponent Friction angle at failure Cohesion (kPa) Maximum deviator strain in compression/ extension % (in pconstant triaxial test) Dilation angle (phase transformation) Dilation parameter

rs

2720

[28]

2710

nw G0

0.45 40

[28–30] [28,31–33]

0.45 47

n n f c emax dev

0.3 0.5 36 0 5/3

e.g. [34] [35,36] [29,30,36–38] -

0.3 0.5 36 0 8/7

c

34

[29,30,39]

34

Xpp

0.2

Calibration using data given in [29,39] [40]

0.3

2.3104

Hydraulic conductivity k (m/s)

2.3104

Cyclic triaxial or simple shear tests are normally used for this purpose; however, the cyclic stress ratio has different expressions for the two tests. Vaid and Sivathayalan [39] recommend a correction factor (Cr) as the ratio of cyclic stress under simple shear (tcyc =s0vc ) and triaxial ðsdcyc =2s03c Þ stress conditions: ðtcyc =s0vc Þsimpleshear ¼ C r ðsdcyc =2s03c Þtriaxial ,

where, tcyc ¼ cyclic shear stress, s0vc ¼ initial effective vertical stress (simple shear), sdcyc ¼ cyclic deviator stress, and s03c ¼initial effective all-round pressure (triaxial). For Fraser Delta sand with Dr ¼ 40%, Vaid and Sivathayalan [39] recommend the value C r ¼ 0:78. Similar correction factors equal to 0.6 and 0.7 are recommended by Seed and Peacock [46] and Castro [47], respectively; however, since the Vaid-Sivathayalan Cr is based on the tests carried out on the Fraser sand, it is deemed more accurate for the present study. The liquefaction strength data used here for calibration of the dilation parameter (Xpp) for set #1 is obtained from the results given by [29] for Dr ¼ 40%. These are shown in Fig. 5 (solid triangles and squares). The values corresponding to solid squares resulted from those corresponding to solid triangles corrected by the Cr factor to represent simple shear tests results. The corresponding Xpp values, obtained from liquefaction strength analysis, are given in Table 3, where, s0nc is the initial consolidation pressure,
CSR is cyclic stress ratio sd;cyc =2s0nc with sd,cyc ¼ deviator stress, and NLiq is the number of cycles to liquefaction. The initial liquefaction was considered to occur when reading 5% double amplitude axial strain, or unit excess pore pressure ratio, whichever occurred first. With the above soil parameters (Table 2, set #1), the theoretical liquefaction strength curve of the sand, obtained from element tests simulating isotropically consolidated,

0.25

New UBC Data

0.2 Cyclic Stress Ratio

Vaid et al. [29] Triaxialto Simple Shear (Cr) 0.15

Element Test (Set #1) Triaxial to Simple Shear

Element Test (Set #2, UBC data)

0.1

0.05

0 1

10

(1)

100

Number of Cycles to Liquefaction Fig. 5. Liquefaction strength analysis for Fraser River sand: comparison between laboratory and numerical simulation results.

ARTICLE IN PRESS A. Azizian, R. Popescu / Soil Dynamics and Earthquake Engineering 26 (2006) 870–887

875

Table 3 Liquefaction strength data [29] and the estimated dilation parameter, Xpp Test No.

s0nc (kPa)

CSR

NLiq.

Xpp

1 2 3

100 200 400

0.145 0.145 0.136

10 10 10

0.18 0.18 0.20

undrained, cyclic triaxial test is shown in Fig. 5 by solid circles. The material properties inferred from the literature (set #1) were then updated using an experimental liquefaction strength curve, based on the results of three cyclic simple shear tests on Fraser River Sand with Dr ¼ 40%, performed at UBC (http://www.civil.ubc.ca/liquefaction/) and represented by solid diamonds in Fig. 5. The points obtained using the new set of numerically calibrated parameters (set #2) are shown in Fig. 5 by open triangles. 4. Factors and Responses Results of 2D and 3D analyses are compared to evaluate the applicability of 2D, plane strain simplifying assumption. A number of predicted responses, including displacement, shear strain, and excess pore water pressure ratio, are considered. The comparison is performed in terms of the ratios of those results obtained from 3D analyses, on the plane of symmetry, to those obtained from 2D analyses. Several possible factors are considered in the study. The most important factor is the width to height ratio of the slope ðB=HÞ, since with increase of this ratio the behaviour of a real 3D structure becomes closer to the idealized plane strain situation. The effects of B=H are studied for a range of peak ground accelerations. In addition, the influence of mesh refinement on numerical results is investigated. The factors are (Fig. 2): width to height ratio ðB=HÞ, maximum acceleration (amax) and mesh size (fine and coarse). These are the selected factors for the main part of the study; however, effects of some other factors such as slope angle, soil properties, and transverse loading are partially investigated, as presented later. The selected responses for comparison are the ratios of 3D to 2D results, predicted at the end of strong shaking ðt ¼ 15 sÞ: 1. Displacement: predicted displacement at crest on the plane of symmetry (node ‘a’ shown in Fig. 2a). 2. Shear strain index: average (weighted by element area) of predicted maximum shear strain in a zone located on the symmetry plane shown in Fig. 6a. 3. Excess pore pressure ratio (ru) index: average (weighted by element area) of the ratio of predicted excess pore pressure to initial vertical effective stress, i.e. ru ¼ ue =s0v0 , in a zone located on the symmetry plane shown in Fig. 6b.

Fig. 6. Selected areas on the symmetry plane of the 3D mesh (z ¼ B=2 in Fig. 2) for calculating: (a) shear strain index; and (b) ru index.

The ratios of above responses (3D/2D) are hereafter referred to as: 1. Rd: 3D/2D crest displacement ratio, 2. Rs: 3D/2D shear strain index ratio, 3. Rr: 3D/2D ru index ratio (ru is the ratio of excess pore water pressure with respect to the initial effective vertical stress). The selection of elements is based on the typical distributions of maximum shear strain and ru (see Figs. 12–17 discussed later) in a way that the corresponding indices represent the slope performance (deformation, instability, etc.) and the influence of the lateral boundaries is not significant. 5. Numerical analysis aspects 5.1. Screening by RSM In the first phase of the analysis, response surface methodology (RSM, see e.g. [42]–[44]) is used to identify the significant factors. Based on the results of this screening phase, further insight is gained by performing more analyses focusing on the effects of resulting significant factors, namely width/height ratio ðB=HÞ and seismic acceleration intensity (characterized in this study by the maximum value, amax). The ranges of these factors are selected as: model width/ height ratio ðB=H ¼ 228Þ, maximum acceleration ðamax ¼ 0:120:5 gÞ. The mesh size effect has been also investigated by including two different mesh sizes (Fig. 2) in the screening phase. The RSM analysis showed that: (1) B=H and amax were the factors with significant effects; (2) there was no significant interaction between factors; and (3) the mesh size effect was insignificant. As mesh size effect resulted insignificant for this study, the coarse mesh (Fig. 2) has been used for further analyses. It is worth noting that the insignificance of mesh size is in terms of both the 3D/2D ratios and the individual 3D and 2D analysis sets; that is, both coarse and fine meshes yield very close results (e.g. crest displacement) for the cases analyzed.

ARTICLE IN PRESS A. Azizian, R. Popescu / Soil Dynamics and Earthquake Engineering 26 (2006) 870–887

876

5.2. Lateral boundary effects When assuming rigid lateral boundaries, 3D analyses predict maximum response in terms of displacements and strains at a certain distance from the lateral boundaries, irrespective of the B=H ratio (see Fig. 7 for predicted crest displacements). For the slope geometry and soil properties considered here, the peak crest deformation occurs at a distance of about 25 m from the lateral boundary for amax ¼ 0:1 g. This distance slightly increases as the peak acceleration increases. Occurrence of such a peak in crest deformation is believed to be produced by the wave interference phenomenon caused by the reflection of waves from the rigid lateral boundaries (z ¼ 0 and B), such as occurring in an elastic rod subjected to a sinusoidal wave with subsequent generation of standing waves. The present situation, i.e. the 3D model of a non-linear material subjected to irregular earthquake loading, is much more complex owing to all types of wave propagation, more importantly, shear waves (SH-waves) propagating vertically from the base and the longitudinal waves (P-waves) originating from the x ¼ 0 and L faces, in addition to those propagating and being reflected horizontally from the z ¼ 0 and B faces. Nonetheless, a qualitative analogy, as explained by Azizian [41], indicates that by changing the artificial rigid boundary conditions, which may simulate real conditions in some cases, to some transmitting boundary conditions, the peak in crest deformation can be eliminated. While the results such as presented in Fig. 7 seem to reflect the true behavior of a soil deposit with rigid vertical boundaries, large predicted displacements at a distance of about 25 m from the z ¼ 0 boundary prevents meaningful result analysis in terms of the effect of B=H ratio on 3D vs. 2D results. Moreover, real-life soil deposits are unlikely to include rigid, parallel, vertical boundaries; and therefore,

wave reflection/interference phenomena, such as those simulated here, may not be of interest for practice. Hence, to obtain meaningful results in terms of 3D effects for slopes with different B=H ratios, the wave interference effects have been eliminated by using transmitting boundaries on the z ¼ 0 face. In this case, the boundary conditions on that face (see Fig. 3) are changed from prescribed acceleration to free movement (code 3 to code 0) and the earthquake acceleration is applied to the transmitting-boundary nodal elements. This is equivalent to imposing a shear stress to the nodal points at the z ¼ 0

Fig. 8. Contours of predicted displacements (m): (a) without, and (b) with transmitting boundary conditions for B=2 ¼ 70 m and amax ¼ 0:1 g at t ¼ 15 s. Deformation magnification scale ¼ 10.

0.8 0.7

B/H = 0.8 B/H = 2 B/H = 3.2 B/H = 5

Crest Disp. (m)

0.6 0.5 0.4

B/H = 8 B/H = 14 B/H = 2 2D

0.3 0.2 0.1 0 0

20

40 60 Distance from side boundary (m)

80

100

Fig. 7. Predicted deformed shape of the crest at t ¼ 15 s, for various B=H ratios and amax ¼ 0:1 g ðH ¼ 10 mÞ.

ARTICLE IN PRESS A. Azizian, R. Popescu / Soil Dynamics and Earthquake Engineering 26 (2006) 870–887

877

boundary, which is a function of the incident (input) motion and soil properties (shear wave velocity and density [3]). This type of transmitting boundary condition is similar to the one presented by Lysmer and Kuhlemeyer [48] in the sense that a shear stress is applied at the boundary, and it generally belongs to the Engquist and Majda [49] family of local transmitting boundaries. Contours of displacements presented in Fig. 8 show the effect of applying the transmitting boundary conditions on eliminating the peak in the crest displacement. The zone with displacements greater than 0.55 m located at a distance of about 20–40 m from the z ¼ 0 face in the model with rigid boundaries (Fig. 8a) is eliminated in the model with transmitting boundaries (Fig. 8b), which predicts maximum displacements on the plane of symmetry. Fig. 10. Contours of total displacement (m) in 3D ðB=H ¼ 3:2Þ and 2D, at t ¼ 15 s, for amax ¼ 0:3 g. Deformation magnification scale is 2.

6. Three-Dimensional Effects 6.1. Analysis results Figs. 9–17 present contours of predicted displacements, maximum shear strains, and excess pore pressure ratio (ru) in 3D and 2D models, at t ¼ 15 s. For 3D, the model with B=H ¼ 3:2 is selected. The following points are notable: 1. Displacement contours (Figs. 9–11): Patterns of displacement contours do not vary significantly with varying amax, and they are similar in shape for 2D and 3D models (on the plane of symmetry). For the type of soil material considered here, in 3D models, maximum values occur on the plane of symmetry at a point between the toe and crest of the slope as indicated by the location of the contour with the highest value. In 2D, as well, maximum values occur between the toe and crest of the slope. This is in accordance with some centrifuge experimental observations (which are closer to 2D configuration) reported by Taboada-Urtuzuastegui et al. [50].

Fig. 11. Contours of total displacement (m) in 3D ðB=H ¼ 3:2Þ and 2D, at t ¼ 15 s, for amax ¼ 0:5 g. Deformation magnification scale is 2.

Fig. 9. Contours of total displacement (m) in 3D ðB=H ¼ 3:2Þ and 2D, at t ¼ 15 s, for amax ¼ 0:1 g. Deformation magnification scale is 2.

2. Shear strain contours (Figs. 12–14): Patterns of maximum shear strain contours vary with amax while they are similar in shape for 2D and 3D models (on the plane of symmetry) for the same amax. With the increase in amax, an area of high shear strain (between 0.04 and 0.06) initiates near the lateral boundaries, due to the assumed rigid boundaries and perfect stick conditions, and expands towards the base of the model. Additionally, as amax increases, the contours of maximum shear strain on the plane of symmetry in 3D models show lower values below the slope face compared to those in 2D (see Figs. 12–14), which indicates that the 3D effect of the lateral boundaries is more significant when amax is higher. This can be more clearly observed from the plots of variation of shear strain index ratio as a function of B/H and amax, presented and discussed later (see Fig. 24). This can be attributed to the fact that the predicted excess pore

ARTICLE IN PRESS 878

A. Azizian, R. Popescu / Soil Dynamics and Earthquake Engineering 26 (2006) 870–887

Fig. 12. Contours of maximum shear strain in 3D ðB=H ¼ 3:2Þ and 2D, at t ¼ 15 s, for amax ¼ 0:1 g. Deformation magnification scale is 2.

Fig. 13. Contours of maximum shear strain in 3D ðB=H ¼ 3:2Þ and 2D, at t ¼ 15 s, for amax ¼ 0:3 g. Deformation magnification scale is 2.

Fig. 14. Contours of maximum shear strain in 3D ðB=H ¼ 3:2Þ and 2D, at t ¼ 15 s, for amax ¼ 0:5 g. Deformation magnification scale is 2.

Fig. 15. Contours of excess pore pressure ratio (ru) in 3D ðB=H ¼ 3:2Þ and 2D, at t ¼ 15 s, for amax ¼ 0:1 g. Deformation magnification scale is 2.

Fig. 16. Contours of excess pore pressure ratio (ru) in 3D ðB=H ¼ 3:2Þ and 2D, at t ¼ 15 s, for amax ¼ 0:3 g. Deformation magnification scale is 2.

Fig. 17. Contours of excess pore pressure ratio (ru) in 3D ðB=H ¼ 3:2Þ and 2D, at t ¼ 15 s, for amax ¼ 0:5 g. Deformation magnification scale is 2.

ARTICLE IN PRESS pressures near the slope face are lower when amax is higher, because of more tendency of dilation (as explained later in this section). More dilation, leading to higher confining effective stresses, will result in a stiffer structure and therefore in more significant constraint effect of the lateral boundaries on the shear strains predicted on the plane of symmetry. Centrifuge experiments on clean sand at the same relative density as considered here (40%), show the same trend (see e.g. Taboada-Urtuzuastegui et al. [50]), where with the increase in acceleration magnitude more dilation results in lower deformations under the slope face than expected. 3. Excess pore pressure ratio (ru) contours (Figs. 15–17): Patterns of ru contours vary significantly with the increase in amax, although in all models the area in the vicinity of the slope has the lowest values of ru due to presence of static shear stress, while the free-field areas in the up- and downslope zones are predicted to liquefy. As indicated by the plots of shear stress (txy) vs. vertical effective stress (ty) of elements A and B (Fig. 18) located on the plane of symmetry of the model with B=H ¼ 3:2 (see Fig. 2), the tendency of plastic dilation in zones with static shear stress (element A) limits the build-up of excess pore pressure (Fig. 19), whereas the state of zero effective stress (and ru ¼ 1) is predicted in the free-field (element B). With the increase in amax, the predicted depth of the liquefied areas in the free-field increases and reaches the base of the model for amax greater than or equal to 0.3 g. Conversely, the larger the amax, the smaller is ru just below the slope crest. This is due to the fact that, for the soil properties considered in this study, the constitutive model predicts more dilation, and thus, less excess pore water pressure, for elements with static shear stress, when larger cyclic load is applied. This is illustrated in Figs. 20 and 21 by comparing the shear stress (txy) vs. vertical stress (sy) and ru vs. time plots of element A obtained from the analyses where amax ¼ 0:1, 0.3 and 0.5 g. The stress path of element A is shown in Fig. 22 in 3D principal stress space,

10

Shear Stress (τxy, kPa)

t = 20 s

t = 9.3 s

8

t=0s

6 4 2 0 -2 -4 t = 5-20 s -6 -8

0

5

t=0 s

A B 10 15 20 25 30 Effective Vertical Stress (σy, kPa)

35

40

Fig. 18. Predicted shear stress (txy) versus effective vertical stress (sy) for elements A and B, shown in Fig. 2b, in the model with B=H ¼ 3:2, for amax ¼ 0:1 g. Element A is near the slope and element B is in the free-field.

ru

A. Azizian, R. Popescu / Soil Dynamics and Earthquake Engineering 26 (2006) 870–887 1.2 1.1 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 -0.1

879

B

A

0

5

10 Time (s)

15

20

Fig. 19. Variation of ru with time for elements A and B shown in Fig. 2b, in the model with B=H ¼ 3:2, for amax ¼ 0:1 g. Element A is near the slope and element B is in the free-field.

illustrating the contraction/dilation behavior of this element when the slope is subjected to large seismic acceleration ðamax ¼ 0:5 gÞ. In this case, the pore pressure build-up is faster; however, once the phase transformation surface is reached, the larger stress ratios result in larger tendency of dilation and thus less pore pressure. Additionally, as indicated by Figs. 15–17, for amax ¼ 0:5 g, the zone with ruo0.25 just below the slope crest extends to the lateral boundary with a greater width near the boundary, whereas for amax ¼ 0:1 g, the predicted ru for the zone below the slope crest is greater than 0.5, both on the plane of symmetry and near the lateral boundary. 6.2. Regression models In this section, the analysis procedure implementing transmitting boundaries, described in the previous sections, is applied to obtain regression models for the curves representing the 3D effects of lateral boundaries on the selected responses (Rd, Rs, Rr) as a function of the model width and amax. The regression models are obtained by means of response surface methodology, using DesignExperts software [51]. Only the two significant factors, i.e. model width to height ratio ðB=HÞ and maximum acceleration (amax), are selected for this study. The crest displacement as well as the shear strain and ru indices on the symmetry face ðz ¼ B=2Þ of the model with B=2 ¼ 100 m ðB=H ¼ 20Þ are used for assessing the 3D effects instead of the corresponding 2D results. This is because the predicted crest displacements on the symmetry face of the 3D model with B=2 ¼ 100 m are expected to be very close to the predicted displacement in a 2D analysis yet small differences may exist due to minor numerical errors induced by using a relatively coarse mesh and by having imperfect transmitting boundaries. Similar numerical errors are expected to occur in the other 3D models with different B=2 values. Crest displacement ratio (Rd): the curve plotted in Fig. 23 is best estimated by the following second-order polynomial

ARTICLE IN PRESS A. Azizian, R. Popescu / Soil Dynamics and Earthquake Engineering 26 (2006) 870–887

880

10 8

t=0s

t = 9.3 s

2

0.4

0.3g

0.2

0

0

-2

-0.2

-6

0.3g

0.6

4

0.1g

0.5g

0.8

6

-4

0.1g 0

2.5

0.5g 5

7.5 Time (s)

10

12.5

15

t = 3.8 s 0

5

(a) 10

10 15 20 Effective Vertical Stress (σy, kPa)

25

30

Fig. 21. Variation of ru with time for element A (shown in Fig.2b), in the model with B=H ¼ 3:2, obtained from the analyses where amax ¼ 0:1, 0.3 and 0.5 g.

amax = 0.3g

8 Shear Stress (τxy, kPa)

1

t = 20 s

ru

Shear Stress (τxy, kPa)

t = 15 s

amax = 0.1g

6

t = 15-20 s

t = 15 s

4

t=0s

t=0s

2

t = 20 s

25

0 20

-2

Dilation

t = 3.9 s -6

0

5

10 15 20 25 Effective Vertical Stress (σy, kPa)

(b)

30

35

σ'2 (kPa)

-4 15

Contraction 10

10

Shear Stress (τxy, kPa)

5

amax = 0.5g

8

t = 3.8 s

t = 20 s

6

0

20

5 10 σ' 15 1 (k Pa )

2 0 -2

5

20

15 10 (kPa) σ' 3

25 0

t=0s

-4 -6

25

0

4

t = 15 s

t = 3.8 s 0

5

(c)

10 15 20 Effective Vertical Stress (σy, kPa)

25

30

Fig. 20. Shear stress (txy) versus effective vertical stress (sy) for element A (shown in Fig. 2b) in the model with B=H ¼ 3:2, obtained from the analyses where amax is: (a) 0.1 g; (b) 0.3 g; and (c) 0.5 g.

regression model: Rd ¼ 0:335 þ 0:205ðB=HÞ  0:016ðB=HÞ2 Rd ¼ 1:0; for B=H46.

for B=H  6, ð2Þ

The adjusted and predicted R-squared values are 0.89 and 0.81, which indicate a reasonable goodness-of-fit for the equation. (Note that the adjusted R-squared is the value of R-squared adjusted for the number of terms in the regression equation relative to the number of data points involved in obtaining it. The adjusted R-squared decreases as the number of terms in the equation increases if those

Fig. 22. Stress path in 3D principal stress space for element A (shown in Fig. 2b), in the model with B=H ¼ 3:2, obtained from the analysis where amax ¼ 0:5 g.

additional terms do not add value to the predictive power of the equation. The predicted R-squared is calculated based on the difference between actual and predicted values, by first predicting where each point should be from an equation that contains all other points except the one in question—see [51].) This regression equation (Eq. (2)) is in terms of B=H only, not amax, since RSM analysis shows that the factor amax is insignificant in this case. Also, the equation is only applicable to B=Ho6, beyond which the ratio is equal to one. Also, since in the analyses performed, the smallest value of B=H was 0.8 (corresponding to B=2 ¼ 4 m) extrapolation of the above equation to smaller values should be done with caution as it is apparent that this second-order polynomial does not go to zero with B=H ¼ 0. For the geometry and assumptions used in this

ARTICLE IN PRESS A. Azizian, R. Popescu / Soil Dynamics and Earthquake Engineering 26 (2006) 870–887

881

1.2

Crest Disp. Ratio

1

0.8 Predicted (0.1g) Predicted (0.3g)

0.6

Predicted (0.5g) RSM

0.4

0.2

0 0

5

10 B/H

15

20

Fig. 23. Crest displacement ratio (Rd): RSM model and analysis results for all values of amax.

1.2

Shear Strain Index Ratio

1 Predicted (0.1g) 0.8 Predicted (0.3g) Predicted (0.5g)

0.6

RSM (0.1g) RSM (0.3g)

0.4

RSM (0.5g) 0.2

0 0

5

10

15

20

B/H Fig. 24. Shear strain index ratio (Rs): RSM model and analysis results.

study, these results indicate that the 3D configuration does not affect crest displacements on the symmetry plane for B=H46. Shear strain ratio (Rs): both factors B=H and amax and their interaction resulted significant. The following secondorder polynomial regression model represents the predicted results. Both adjusted and predicted R-squared values are equal to 0.97 that shows a very high goodness-of-fit. Rs ¼ 0:605 þ 0:144ðB=HÞ  0:543ðamax Þ  0:012ðB=HÞ2 þ 0:062ðB=HÞ:ðamax Þ; for B=H  7, Rs ¼ 1:0; for B=H47.

ð3Þ

Depending on amax, the value of B=H at which Rs is very close to one ranges from 6.2 to 7.25; therefore, the B=H  7 limit is selected for the above equation. In other words, if

B=H47 the 3D effects in terms of predicted shear strains are negligible. The recorded results and the regression model are shown in Fig. 24. The analysis results show that as amax increases, values of Rs decrease. This decrease is greater when B=H is lower. As discussed earlier, this is because the numerical model predicts more dilation when a soil element with static shear stress is subjected to larger loading, which for the soil material analyzed here, results in more pronounced effect of the lateral boundaries. In other words, with increasing amax, 3D effects are predicted to become more significant. ru ratio (Rr): The predicted results shown in Fig. 25 indicate that, except for the smallest model width ðB=H ¼ 0:8Þ, values of Rr (predicted on the symmetry plane) increase with the decrease in model width to a maximum value of 1.1, due to proximity of the lateral

ARTICLE IN PRESS A. Azizian, R. Popescu / Soil Dynamics and Earthquake Engineering 26 (2006) 870–887

882

1.15 1.1 1.05 Predicted (0.1g) ru Index Ratio

1

Predicted (0.3g)

0.95

Predicted (0.5g) RSM (0.1g)

0.9

RSM (0.3g)

0.85

RSM (0.5g) 0.8 0.75 0.7 0

5

10

15

20

B/H Fig. 25. ru index ratio (Rr): RSM model and analysis results.

boundaries that induce local shearing on a vertical plane. For B=H ¼ 0:8 and amax ¼ 0:3 and 0.5 g, Rr resulted below 1, apparently due to numerical errors induced by use of too coarse a mesh. Setting aside the outlier values (for B=H ¼ 0:8), the following regression model can be obtained, as plotted in Fig. 25. Adjusted and predicted R-squared values are 0.96 and 0.83, respectively, which indicate a satisfactory goodness-of-fit.

Table 4 Effect of slope angle Slope

2H:1V 3H:1V 4H:1V

B=H ¼ 2

B=H ¼ 5

Rd

Rs

Rr

Rd

Rs

Rr

0.74 0.73 0.73

0.77 0.72 0.74

1.04 1.03 1.03

1.00 0.99 0.99

0.97 0.95 0.95

1.01 1.01 1.01

Rr ¼ 1:173  0:034ðB=HÞ  0:406ðamax Þ þ 1:75E  03ðB=HÞ2 þ 0:291ðamax Þ2 þ 0:034ðB=HÞðamax Þ; for B=H  6, Rr ¼ 1; for B=H46. ð4Þ Eqs. (2)–(4) can be used for estimating results of a 3D analysis from 2D analysis predictions. For example, if the 3D model has a B=H ratio of 2, and the earthquake peak acceleration is 0.2, Eqs. (2)–(4) indicate that, for the 3D model on the plane of symmetry: 1) crest displacement is about 0.68 of that predicted in 2D model; 2) shear strains predicted under the slope are on average about 0.76 of those predicted in 2D model; and (3) excess pore pressures are on average about 8% larger than those predicted in 2D model. Additionally, for B=H ratio as low as one, the 3D crest displacement and shear strain index are approximately 50% of those predicted in 2D analyses, whereas excess pore pressure is expected to result about 10% higher in a 3D than in a 2D analysis.

analyzed here, the slope angle does not significantly influence the results in terms of 3D vs. 2D results. 6.4. Updated material properties For the material properties calibrated based on cyclic undrained simple shear tests performed at UBC on Fraser River sand (http://www.civil.ubc.ca/liquefaction/) and labelled as set #2 in Table 2, the ratios Rd, Rs, and Rr for the model with B=H ¼ 3:2 and amax ¼ 0:3 g are 0.80, 0.80, and 1.02, respectively. The ratios obtained for the first set of material properties (set #1) were 0.87, 0.88, and 1.03, respectively. The regression models predict these ratios as 0.83, 0.84, and 1.02, respectively. The second set of material properties corresponds to a soil with slightly lower liquefaction strength. The difference in 3D effects on predicted strains and displacement indicate that liquefaction strength can be an important factor.

6.3. Effect of slope angle 6.5. Effect of transverse loading The results presented so far were obtained for 3H:1V slope. A limited study [41] was performed to evaluate the effect of slope angle. The results obtained for 4H:1V and 2H:1V slopes as well as the 3H:1V slope (with B=H ¼ 2 and 5, and amax ¼ 0:3 g) are given in Table 4. The values in Table 4 indicate that for the type of soil and geometry

Effect of transverse loading, which cannot be accounted for in a 2D analysis, is partially examined here by performing two analyses. As opposed to the analyses described so far, where due to symmetric loading only half of the structure was considered, the whole structure is

ARTICLE IN PRESS A. Azizian, R. Popescu / Soil Dynamics and Earthquake Engineering 26 (2006) 870–887

10 0

Legend 0: Active degree-of-freedom 1: Prescribed displacement 3. Prescribed acceleration

1 13

z=0 x

z

110

: 113

Coners: 113010

000

3110 11

z=B

fa ce :

1130

1130 10

L= 1

40 m

00

Base : 113

010 1 nt: Fro

0 01

00

1 13

10 0 113

113

f ace

100

11 3

10 0

y

113

ble r ia Va : B 00 001 k: 1 c a B

883

110

Fig. 26. Boundary conditions of the 3D model subject to transverse loading. The boundary condition codes are according to the definitions given in Table 1, with the following order: x, y, z components of solid phase motion, and x, y, z components of fluid phase motion.

analyzed here. Details of boundary conditions are shown in Fig. 26. The main differences compared to the previously used boundary conditions (described earlier) are:

Fig. 27. Contours of total displacement (m) induced by transverse loading for B=H ¼ 2 and 8, and amax ¼ 0:3 g at t ¼ 15 s. Deformation magnification scale is 2.

1. The earthquake acceleration is applied along the zdirection instead of the x-direction. 2. The planes with transmitting boundary conditions are switched, i.e. in this case, x ¼ 0 and L planes have transmitting boundaries. In other words, the planes perpendicular to the direction of loading do not have transmitting boundaries. 3. Nodes with slaved degrees-of-freedom for horizontal fluid motion are now located on the z ¼ 0 and B planes and the very next planes of nodes. In this case, the sixth degree-of-freedom (horizontal component of fluid motion in the z-direction) of each node located on those faces is slaved to that of the corresponding node located on the very next vertical plane of nodes. Two structures with B=H ¼ 2 and 8, subject to amax ¼ 0:3 g, are analyzed. Figs. 27–29 show contours of results, with only one half of the mesh presented (from z ¼ 0 to z ¼ B=2) so that the results obtained for the symmetry plane can be shown. At the selected time instant of t ¼ 15 s that is after the end of shaking, the contours are almost symmetrical with respect to the z ¼ B=2 plane. 3D/ 2D ratios of responses (Rd, Rs, and Rr), as given in Table 5, can be calculated assuming that the mesh with B=H ¼ 8 is equivalent to a 2D configuration (i.e. the results on the plane at z ¼ B=2 are not influenced by the lateral boundaries). Also given in Table 5 are the values of these ratios calculated from the regression models presented before for seismic loading parallel to the slope. The following points are notable: 1. The large difference between the two groups of ratios obtained from the analyses with transverse and parallel

Fig. 28. Contours of maximum shear strain induced by transverse loading for B=H ¼ 2 and 8, and amax ¼ 0:3 g at t ¼ 15 s. Deformation magnification scale is 2.

loading, and the fact that the ratios obtained for transverse loading are much lower, indicate that the 3D effects may be more significant in the case of transverse loading. In other words, 3D configuration or proximity of the lateral boundaries may have greater importance for transverse loading than for parallel loading. Also, this may raise the question if B=H ¼ 8 can be considered sufficient for assuming 3D effects of lateral boundaries as insignificant. 2. The predicted displacements are much lower in this case compared to the case of parallel loading, by about 50% (also compare Fig. 27 with Fig. 10), which indicates that parallel loading is the critical mode of loading.

ARTICLE IN PRESS A. Azizian, R. Popescu / Soil Dynamics and Earthquake Engineering 26 (2006) 870–887

884

Fig. 29. Contours of excess pore pressure ratio (ru) induced by transverse loading for B=H ¼ 2 and 8, and amax ¼ 0:3 g at t ¼ 15 s. Deformation magnification scale is 2.

Table 5 Comparison between 3D/2D ratios of responses (for B=H ¼ 2, and amax ¼ 0:3 g) calculated from cases of parallel and transverse loading Loading

Rd

Rs

Rr

Transverse Parallel

0.49 0.68

0.59 0.72

0.75 1.04

3. Same as for displacements, shear strains, are also lower in this case. Zones with higher shear strains are located near the base of the analysis domain in the case of transverse loading, as opposed to the case of parallel loading in which those zones stretch from the base to the slope face. 4. In both cases, the predicted ru below the slope is lower than in other locations in the analysis domain (see Figs 29 and 16). However, another significant effect of transverse loading is that, in contrast to the case of parallel loading, the predicted ru in the mesh with lower B/H is lower than that in the one with higher B=H. In other words, the decreasing trend of Rr with increasing B=H, as observed in the case of parallel loading (Fig. 25), is not observed here for the transverse loading case. All above conclusions have been made based on the results of a very limited study. More investigations are necessary in this respect. 6.6. Comparison with previous studies As noted earlier, according to Duncan [1], with a few exceptions, most of the 3D slope stability analyses have

resulted in higher factors of safety compared to those resulting from 2D analyses, and 2D deformation analyses have given higher values of displacements, strains and stresses compared to 3D analysis results. None of the studies cited by Duncan [1] refer to dynamic, effectivestress, non-linear analyses, such as those performed here. The results presented here are in general agreement with the previous studies in the sense that the crest displacement and shear strain index predicted in 3D analyses are lower than those predicted in 2D. Excess pore pressure ratios predicted in 3D analyses, in the vicinity of local boundaries, are higher than those predicted in 2D analyses apparently due to local shearing. Ratios of 3D/2D factors of safety given by previous studies are herein compared to the 3D/2D ratio of the inverse of shear strain index, as explained hereafter. The factor of safety is normally defined as the ratio of the total available shear resistance to the total mobilized shear stress along the critical failure surface. Because with the increase in the shear stress mobilized on the failure surface the shear strains increase, the factor of safety can be qualitatively related to the inverse of shear strains occurring in the vicinity of the failure surface. In other words, if shear strains are higher, the mobilized shear stresses are also higher and therefore the factor of safety can be assumed to be lower. As the factor of safety is calculated by dividing the existing shear stresses by a limit shear stress (or shear strength), one may divide the existing shear strains by a limit value to normalize the mobilized shear strains with respect to a value corresponding to failure. However, since in this study only the 3D/2D ratio of response is of primary interest, there is no need for such normalization. It should also be mentioned that the shear strain index considered here is calculated over an arbitrarily selected zone close to the slope face (on the plane of symmetry) and not on a possible failure surface. Moreover, the selection of this index may not be appropriate when large shear strains occur in both 2D and 3D analyses. An index involving shear stresses would give a better image in this situation; however, because in this type of analysis the shear stresses are highly oscillatory, use of such an index was deemed impractical. As for quantifying the 3D effects as a function of slope width/height ratio, Fig. 30 gives a comparison between the results of the present study and a previous study by Arellano and Stark [10], based on the factor of safety approach to static stability analysis. The area enclosed by dashed lines in Fig. 30b is obtained from the present study from the inverse shear strain index ratio for values of amax ¼ 0:120:5 g. Although the two studies are essentially different (static vs. dynamic, and factor of safety vs. deformation analysis), the results are in the same range for the 3H:1V slope (Fig. 30b). The strong effect of the slope angle shown by the Arellano–Stark curves (Fig. 30a) results from an increase in the area of the vertical sides of the slide mass due to the way the slope geometry was defined in their study. In the study presented here, however, since the

ARTICLE IN PRESS A. Azizian, R. Popescu / Soil Dynamics and Earthquake Engineering 26 (2006) 870–887

Ratio of 3D/2D Factors of Safety

3.5

Table 6 Comparison with some previous studies (for 3H:1V slope)

1H:1V 3

3H:1V

3D/2D Factor of safety

3D/2D Factor of safety [10]

Present study Inverse of 3D/2D shear strain index ðamax ¼ 0:3 gÞ

[8] [52] [53]

1.3 1.2 1.2

1.36 2.15 (Extrapolated) 1.30

1.4 1.7 1.3

5H:1V

2

2.00 0.89 2.45

1.5

1

2

4

(a)

6 8 Width/Height Ratio

Hungr et al.[53]

that for width/height ratios greater than 5, the differences between F2D and F3D tend to lose significance. This limit corresponds to about 15% difference between 2D and 3D results. Fig. 30b shows that the curves obtained in the present study tend to decrease faster than those given by the former study. For a similar tolerance (i.e. about 15%) the present study indicates that 2D dynamic slope analyses are sufficiently accurate for B=H4325, with larger values corresponding to larger seismic accelerations.

Leshchinsky and Huang [8]

7. Conclusions

10

12

2.5 Arellano and Stark [10] Dennhardt and Forster [52]

2.25 2 1.75

increasing amax

Ratio of 3D/2D Factors of Safety

Study B/H

2.5

1

This study

1.5 1.25 1 0

(b)

885

2

4

6 8 Width/Height Ratio

10

12

Fig. 30. Comparison with some previous studies: (a) variation of 3D/2D factors of safety with width/height ratio redrawn after Arellano and Stark [10] for different slope angles; (b) results of the present study in terms of the inverse of 3D/2D shear strain index (dashed line) compared to other studies (for 3H:1V slope).

change in slope angle does not significantly change the area of the lateral boundaries, such an effect resulted less significant. Assuming that the ratio of 3D vs. 2D factors of safety is proportional to 1/Rs, Eq. (3) can be used to evaluate the 3D factor of safety based on 2D analysis results. For example, for a slope with B=H ¼ 2 subjected to amax ¼ 0:2g, one may infer from Eq. (3) and Fig. 30b that the 3D factor of safety is about 30–35% larger than the 2D factor of safety. Table 6 presents selected data from some other studies. The values from the present study included in Table 6 correspond to amax ¼ 0:3 g. These values are larger than the values from previous factor-of-safety studies, which may be an indication that, for small values of B/H, the 3D effects may be more significant in a dynamic analysis with large maximum acceleration than in a static analysis. These data points are also shown in Fig. 30b. It is worth noting that Chugh [54] presents static, finite element-based analysis of the same slope geometry analyzed by Arellano and Stark [10]. He has concluded

A comparison is made between the results of three- and two-dimensional seismic analyses of submarine slopes, using dynamic, effective stress, fully coupled, non-linear finite element analysis. Outcomes of this study provide some quantitative guidelines for geotechnical practice in this particular area of soil dynamics to extrapolate results of less expensive 2D, plane strain analyses to more realistic 3D behavior. The comparison is in terms of several factors and responses. The factors are model width (or distance between lateral boundaries) normalized as width/height ratio ðB=HÞ, maximum acceleration of earthquake (amax), and mesh size (fine/coarse). The responses are ratios of 3D to 2D results, namely, crest displacement, a shear strain index, and an excess pore water pressure ratio (ru) index. These ratios are referred to as Rd, Rs and Rr, respectively. Regression models are provided (see Eqs. (2)–(4) plotted in Figs. 23–25) to relate the responses (Rd, Rs and Rr) and the factors (B=H and amax) for practical extrapolation of 2D analysis results to 3D situation. These equations provide the ratio of each response in 3D analysis to that in 2D analysis and can be used to extrapolate 2D analysis results for situations similar to those considered in this study. The following conclusions can be drawn based on the results of this comparison:



For a tolerance of about 15%, the analysis results show that a 2D analysis is sufficiently accurate if the width/ height ratio ðB=HÞ of the slope is greater than 3–5. This value generally increases with the increase in maximum acceleration of the applied earthquake motion. Moreover, for B=H4627, the differences between 3D

ARTICLE IN PRESS A. Azizian, R. Popescu / Soil Dynamics and Earthquake Engineering 26 (2006) 870–887

886



 



analysis prediction on the symmetry plane of the slope and 2D analysis predictions resulted insignificant. It is found that the peak acceleration (amax) has a significant effect on the shear strain index ratio (Rs): for a constant B=H, Rs is lower for higher amax. This indicates that the 3D effects are predicted to be more significant when amax is larger, which can be attributed to the fact that, for the soil properties assumed here, the constitutive model predicts larger dilation for areas below the slope face (with static shear stress) when larger cyclic load is applied. More dilation results in stiffer soil, and thus, more pronounced constraint effects of the lateral boundaries. For the particular case analyzed here, with vertical lateral boundaries, it also resulted that the change in slope angle does not affect the 3D/2D ratios. It is also found that a soil with slightly lower liquefaction strength can result in slightly more significant 3D effects, apparently due to the impact the dilatant behavior under the slope has on the constraint effect of the boundaries. A limited study on the effects of transverse loading showed that the 3D effects of lateral boundaries are more significant when the slope is subjected to transverse loading than when it is subjected to seismic loading along the slope direction.

Acknowledgement The financial support provided by the National Sciences and Engineering Research Council of Canada (NSERC) through the COSTA-Canada project, the NSERC Research Grant No. RG203795-02, and the Graduate Fellowship provided by Memorial University of Newfoundland is acknowledged. The authors are also grateful to Professor L. Lye, Memorial University, for his valuable suggestions and help on using Response Surface Methodology.

References [1] Duncan JM. State-of-the-art: limit equilibrium and finite element analysis of slopes. ASCE J Geotech Geoenviron Eng 1996; 122(7):577–96. [2] Prevost JH. A simple plasticity theory for frictional cohesionless soils. Soil Dyn Earthquake Eng 1985;4(1):9–17. [3] Prevost JH. Dynaflow—A nonlinear transient finite element analysis program. Version 2002 Release 01.A. (first release: 1981). Princeton, NJ: Department of Civil Engineeering & Operation Research, Princeton Universtiy; 2002 http://www.princeton.edu/dynaflow/ description_df.htm. [4] Justo JL, Saura J. Three-dimensional analysis of Infiernillo dam during construction and filling of the reservoir. Int J Numer Anal Methods Geomech 1983;7:225–43. [5] Byrne RJ, Kendall J, Brown S. Cause and mechanism of failure, Kettleman Hills landfill B-19, Phase IA. In: Proceedings of the ASCE specialty conference on stability and performance of slopes and embankments—II, vol. 2, 1992; p. 1188–215.

[6] Stark TD, Eid HT. Performance of three-dimensional slope stability methods in practice. ASCE J Geotech Geoenviron Eng 1998;124(11): 1049–60. [7] Jeremic B. Finite element methods for 3D slope stability analysis. Slope Stability 2000. ASCE 2000;224–38. [8] Leshchinsky D, Huang CC. Generalized three-dimensional slopestability analysis. ASCE J Geotech Eng 1992;118(11):1748–64. [9] Duncan JM, Stark TD. Soil strength from back analysis of slope failures. In: Proceedings of the specialty conference on stability and performance of slopes and embankments-II. vol. 1. ASCE, Berkeley, CA, 1992; p. 890–904. [10] Arellano D, Stark TD. Importance of three-dimensional slope stability analyses in practice. Slope stability 2000, ASCE Geotech Special Pub 2000; (101): 18–32. [11] Lefebvre G, Duncan JM, Wilson EL. Three-dimensional finite element analysis of dams. ASCE J Soil Mech Found Div 1973;99(7): 495–507. [12] Martin HL. A three-dimensional analysis of the Storvass dam. Int J Numer Analyti Methods Geomech 1978;2:3–17. [13] Prevost JH, Abdel-Ghaffar AM, Lacy SJ. Nonlinear dynamic analysis of an earth dam. ASCE J Geotech Eng 1985;111(7):882–97. [14] Zienkiewics OC, Taylor RL. The finite element method. vol. I: Basic formulation and linear problems, vol. II: solid and fluid mechanics, dynamics and non-linearity. UK: McGraw-Hill; 1989. [15] Youd TL, Idriss IM, Andrus RD, Arango I, Castro G, Christian JT, et al. Liquefaction resistance of soils. Summary report from the 1996 NCEER and 1998 NCEER/NSF workshops on evaluation of liquefaction resistance of soils. ASCE J Geotech Geoenviron Eng 2001;127(10):817–33. [16] Biot MA. Mechanics of deformation and acoustic propagation in porous media. J Appl Phys 1962;33(4):1482–98. [17] Prevost JH. Nonlinear dynamic response analysis of soil and soil–structure interacting systems. In: Seco e Pinto P, editor. Soil dynamics and geotechnical earthquake engineering. Rotterdam: Balkema; 1993. [18] Popescu R, Prevost JH. Centrifuge validation of a numerical model for dynamic soil liquefaction. Soil Dyn Earthquake Eng 1993;12: 73–90. [19] Popescu R, Prevost JH. Numerical class ‘A’ predictions for models no. 1, 2, 3, 4a, 4b, 7, 11 and 12. In: Arulanandan K, Scott RF, editors. Proceedings international conference on verification of numerical procedures for the analysis of soil liquefaction, vol. 1. Rotterdam: Balkema; 1993. p. 1105–207. [20] Popescu R. Stochastic variability of soil properties: Data Analysis, Digital Simulation, Effects on System Behavior. Ph.D. Thesis, Princeton, NJ: Princeton University; 1995. [21] Keane CM, Prevost JH. An analysis of earthquake data observed at the wildlife liquefaction array site, imperial county, California. In: Proceedings of the 2nd US–Japan workshop on liquefaction, large ground deformations and effects on lifelines. New York, 1989. p. 39–53. [22] Popescu R, Prevost JH. Comparison between VELACS numerical ‘class A’ predictions and centrifuge experimental soil test results. Soil Dyn Earthquake Eng 1995;14(2):79–92. [23] Popescu R, Prevost JH, Ohbo N, Hayashi K. Numerical simulations of soil liquefaction. In: Proceedings of the 4th US–Japan workshop on liquefaction, large ground deformations, and effects on lifelines. 1992. p. 269–82. [24] Popescu R, Prevost JH, Deodatis G. Seismic analysis of a Wharf Seawall. In: Proceedings of the structural engineers world congress (SEWC), San Francisco, CA, 1998. [25] COSTA – Canada; 2000; http://www.costa-canada.ggl.ulaval.ca/english.html. [26] Locat J, Bornhold B, Byrne PM, Hart B, Hughes-Clarke J, Konrad JM, et al. COSTA-Canada, A Canadian contribution to the study of continental slope stability: an overview. In: Proceedings of the 54th Canadian geotechnical conference, Calgary. vol. 2. 2001, p. 730–37.

ARTICLE IN PRESS A. Azizian, R. Popescu / Soil Dynamics and Earthquake Engineering 26 (2006) 870–887 [27] UBC (Uniform Building Code). International conference of building officials, ICBO, vol. 2. CA: Whittier; 1994. [28] Howie JA, Shozen T, Vaid YP. Effect of ageing on stiffness of very loose sand. Can Geotech J 2002;39:149–56. [29] Vaid YP, Stedman JD, Sivathayalan S. Confining stress and static shear effects in cyclic liquefaction. Can Geotech J 2001;38:580–91. [30] Sivathayalan S, Vaid YP. Influence of generalized initial state and principal stress rotation on the undrained response of sands. Can Geotech J 2002;39:63–76. [31] Hardin BO, Richard FE. Elastic wave velocities in granular soils. ASCE J Soil Mech Found Div 1963;89(SM1):33–65. [32] Belloti R, Ghionna V, Jamiolkowski M, Lancellotta R, Manfredini G. Deformation characteristics of cohesionless soils from in-situ tests. In: Clemence SP, editor. Use of in-situ tests in geotechnical engineering. 1986. p. 47–73. [33] Vaid YP, Eliadorani A. Undrained and drained (?) stress strain response. Can Geotech J 2000;37:1126–30. [34] Das BM. Principles of geotechnical engineering. PWS-KENT, 1998. [35] Richart FE, Hall JR, Woods RD. Vibrations of soils and foundations. Englewood Cliffs, NJ: Prentice-Hall; 1970. [36] Chillarige ARV, Robertson PK, Morgenstern NR, Christian HA. Evaluation of the in situ state of Fraser River sand. Can Geotech J 1997;34:510–9. [37] Vaid YP, Chern JC. Cyclic and monotonic undrained response of saturated sands. In: Advances in the art of testing soils under cyclic conditions, ASCE, 1985. p. 120–47. [38] Vaid YP, Sivathayalan S, Stedman D. Influence of specimenreconstituting method on the undrained response of sand. ASCE Geotech Testing J GTJ9909 1999;22(3):187–96. [39] Vaid YP, Sivathayalan S. Static and cyclic liquefaction potential of Fraser Delta sand in simple shear and triaxial tests. Can Geotech J 1996;33:281–9. [40] Shahabi AA, Das BM, Tarquin AJ. An empirical relation for coefficient of permeability of sands. In: Proceedings of the 4th Australia–New Zealand conference on geomechanics, vol. 1. 1984, p. 54–57.

887

[41] Azizian A. Seismic analysis of submarine slopes: retrogressive and three-dimensional effects. Ph.D. thesis, Memorial University of Newfoundland, St. John’s, NL, Canada, January 2004. [42] Myers RH, Montgomery DC. Response surface methodology: process and product optimization using design experiments. New York: John Wiley & Sons; 1995. [43] Lye LM. Design of Experiments in Civil Engineering: Are we still in the 1920s? In: Annual conference of the Canadian society for civil engineering, GE069, 2002. p. 1–10. [44] Zangeneh N, Azizian A, Lye L, Popescu R. Application of response surface methodology in numerical geotechnical analysis. In: Proceedings of the 55th Canadian geotechnical conference, 2002. p. 321–29. [45] Nave CR. HyperPhysics. Georgia State University, 2000; http://hyperphysics.phy-astr.gsu.edu/hbase/hframe.html. [46] Seed HB, Peacock WH. Test procedures for measuring soil liquefaction characteristics. ASCE J Soil Mech Found Div 1971; 97(SM8):1099–119. [47] Castro G. Liquefaction and cyclic mobility of saturated sands. J Geotech Eng ASCE 1975;101(GT6):551–69. [48] Lysmer J, Kuhlemeyer RL. Finite dynamic model for infinite media. ASCE J Eng Mech Div 1969;95(EM4):859–77. [49] Engquist B, Majda A. Absorbing boundary conditions for the numerical simulation of waves. Math Comput 1977;31(139):629–51. [50] Taboada-Urtuzuastegui VM, Martinez-Ramirez G, Abdoun T. Centrifuge modeling of seismic behavior of a slope in liquefiable soil. Soil Dyn Earthquake Eng 2002;22:1043–9. [51] Design-Experts, http://www.statease.com. [52] Dennhardt M, Forster W. Problems of three-dimensional slope stability. In: Proceedings of the 11th international conference of soil mechanics and foundation engineering, 1985, p. 427–31. [53] Hungr O, Salgado FM, Byrne PM. Evaluation of a three-dimensional method of slope stability analysis. Can Geotech J 1989;26(4): 679–86. [54] Chugh AK. On the boundary conditions in slope stability analysis. Int J Numer Analyt Methods Geomech 2003;27:905–26.