Employing a variable permeability model in numerical simulation of saturated sand behavior under earthquake loading

Computers and Geotechnics 55 (2014) 211–223

Contents lists available at ScienceDirect

Computers and Geotechnics journal homepage: www.elsevier.com/locate/compgeo

Employing a variable permeability model in numerical simulation of saturated sand behavior under earthquake loading H. Shahir a,⇑, B. Mohammadi-Haji a, A. Ghassemi b a b

Kharazmi University, P.O. Box 31979-37551, Tehran, Iran Qazvin Islamic Azad University, P.O. Box 34185-1416, Qazvin, Iran

a r t i c l e

i n f o

Article history: Received 25 February 2013 Received in revised form 3 August 2013 Accepted 16 September 2013 Available online 6 October 2013 Keywords: Liquefaction Coupled analysis Permeability Settlement

a b s t r a c t Despite advances in the numerical analysis of saturated sand behavior under earthquake loading, accurate prediction of liquefaction-related phenomena by numerical simulation remains a challenge. Variation of the coefficient of permeability is a key issue which has not obtained due attention in most previous modeling. In this study, a revised form of a recently proposed variable permeability function was implemented in a fully coupled dynamic model adopting modern two-surface plasticity constitutive law to evaluate the effects of permeability variations on the results of numerical modeling. The variable permeability model is comprised of a simple function relating the permeability coefficient of soil mass to the excess pore water ratio. In this study, the constants of the variable permeability function were attained based mainly on theoretical evidence and experimental observation. Well-documented centrifuge experiments were examined to evaluate how well the proposed model captures the main features of soil response to earthquake loading. The results indicate that the proposed function greatly enhanced the capability of numerical modeling to predict the behavior of saturated sand under cyclic loading. Particularly, the variable permeability model with proposed constants significantly improved the amount of liquefaction-induced settlement predicted by numerical modeling. Ó 2013 Elsevier Ltd. All rights reserved.

1. Introduction Liquefaction is an important earthquake-induced phenomena that may result in significant damage to engineered structures such as quay walls, bridges, tunnels, and lifeline facilities. Liquefaction is generally attributed to the loss of strength of saturated loose sandy deposits as a consequence of the simultaneous generation of excess pore water pressure and reduction of effective stress through the liquefied layers subjected to seismic excitation. Despite attempts to clarify the mechanisms involved [1–4], a comprehensive theory describing the overall aspects of liquefaction and post-liquefaction behavior of sandy soils has not yet been developed. Accurate prediction of land subsidence in liquefied layers is important to many civil engineering applications. In recent decades, research has focused on factors affecting settlement of saturated sandy layers subjected to earthquake loading using different methodologies, including laboratory experiments and numerical simulations. For example, Jafarzadeh and Yanagisawa [5] used a small-scale shaking table and Dobry and Taboada [6]

⇑ Corresponding author. E-mail addresses: [email protected] (H. Shahir), bahareh.mohammadi.ac@gmail. com (B. Mohammadi-Haji), [email protected] (A. Ghassemi). 0266-352X/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compgeo.2013.09.007

used centrifuge geotechnical modeling tests to experimentally investigate the behavior of saturated sand under cyclic loading. Data from physical models, especially from centrifuge tests, generally provide a credible basis for calibration of model constants in numerical simulations. The results of some of these experiments were used in the current study to validate the numerical model. Generally, numerical simulation of liquefaction-induced settlement is difficult. Some researchers have not reported the settlement obtained by their numerical models [7,8] and other models gave lower settlement values compared to what occurred under real conditions [9]. The common weakness of the numerical models in simulating the surface settlement of liquefied soils can be mainly attributed to the superimposed effects of improper constitutive law for saturated sandy soils under seismic loading and/or unrealistic modeling of water discharge from the soil mass during earthquake loading, which is responsible for the reduction of saturated mass volume. Among the remarkable constitutive models developed to simulate the response of sandy soils subjected to cyclic loading, multi-surface plasticity and bounding surface plasticity models have been successfully employed to numerically model liquefaction-related phenomena [10–12]. In this study, a version of the bounding surface plasticity model modified by Dafalias and Manzari [13] was applied for numerical analysis. Bounding surface

212

H. Shahir et al. / Computers and Geotechnics 55 (2014) 211–223

Nomenclature A0 Ad B b c ch cz D d e e0 ec ff fs fv f(p) f(s) G G0 H h h0 k K0 kb kd ki kl Kp M m

dilatancy function of fabric change strain-displacement matrix stress vector critical state parameter plastic modulus parameter fabric dilatancy parameter dilatancy coefficient stress vector current void ratio critical state constant critical void ratio fluid property factor fluid property factor fluid property factor vector in u–p formulation vector in u–p formulation elastic shear modulus model parameter permeability matrix positive scalar-valued function plastic modulus parameter hydraulic conductivity pore shape factor build up phase coefficient of permeability dissipation phase coefficient of permeability initial coefficient of permeability liquefaction state coefficient of permeability plastic modulus mass matrix yield surface parameter

models, such as that proposed by Dafalias and Manzari, use the same set of model parameters to simulate soil behavior for a range of densities and confining pressures using a unified plasticity framework based on critical state concepts. This model has shown satisfactory results for simulating the monotonic and cyclic behavior of dry and saturated sands, especially to capture the volumetric response of soil responsible for surface settlement [12,14,15]. The numerical efficiency of the model is good because only the yield surface must be updated at each increment in response to kinematic and isotropic hardening [12]. Conventional belief about the behavior of undrained saturated sands during earthquake loading contradicts reported observations indicating a gradual increase of settlement from initial stages of loading [16–18]. In other words, significant water drainage from the soil mass takes place during shaking while excess pore water is being generated. Consequently, a considerable amount of settlement occurs during the shaking. This settlement is caused by high values of outward water flux from the soil mass caused by the increase of permeability in the soil while approaching liquefaction [5,19–21]. After shaking ends, the soil continues to settle at a decreasing rate from drainage caused by the hydraulic gradient remaining in the soil mass. Neglecting variations in the coefficient of permeability is a main reason for the under-prediction of settlement in most previous numerical models. Microscopic changes in the characteristics of saturated pores between sand grains are responsible for the variation of soil permeability. A detailed description of the changes in permeability in saturated sand under cyclic loading will be presented later in this paper. Some researchers account for the variation in permeability by assuming a constant increased permeability value [12,21,22]. This

Mb Md nb nd nb nd P p Pc Q ru S S0 T U Zmax

a b1 b2

cf n kc

lf g m r0 w

bounding stress ratio dilatancy stress ratio plastic modulus parameter dilatancy parameter positive material constant positive material constant pore pressure vector mean effective stress existing confining stress discrete gradient operator coupling the motion and flow equations excess pore pressure ratio compressibility matrix specific surface tortuosity solid displacement fabric dilatancy parameter constant constant constant unit weight of fluid critical state constant critical state constant fluid viscosity stress ratio poisson’s ratio effective stress tensor state parameter

is not a realistic model of the mechanisms involved in the process. Recently, Shahir et al. [15] suggested a variable permeability function with respect to excess pore pressure ratio. They determined the constants of this function by adjusting the results of the numerical model using the results of centrifuge measurements. Although this variable permeability model established an appropriate framework to account for permeability variation, there is no theoretical basis for the proposed constants of this function. In addition, these constants were found in the course of model calibration for one test only and may not be appropriate for numerical modeling of centrifuge tests under different conditions. The present study uses a revised form of the variable permeability model function proposed by Shahir et al. [15] implemented in a fully coupled dynamic analysis adopting a two-surface plasticity constitutive law. The model investigates the effects of permeability changes in saturated sand under earthquake loading on the results of numerical modeling. The original variable permeability function was modified by selecting different values for the function constants. Theoretical evidence and experimental observations were utilized to discover the rough values for these parameters for sand. These values were then calibrated such that the results of numerical simulations of benchmark problems such as settlements and excess pore pressure are a better match to experimental records. The benchmark problems were four well-documented centrifuge experiments covering a range of values of the following parameters: relative density, initial permeability, maximum base acceleration, and height of the sand layer. The focus of the investigation was the degree of improvement in predicting liquefaction-induced settlement using the variable permeability function with the proposed constants in the numerical model.

H. Shahir et al. / Computers and Geotechnics 55 (2014) 211–223

2. Permeability characteristics of liquefied sand An appropriate variable permeability function for saturated sand under cyclic loading as a function of excess pore pressure ratio (ru) calls for accurate estimation of the permeability coefficient during extreme liquefaction (ru = 1). Although some approaches estimate the level of permeability increase during liquefaction, direct measurement of the permeability coefficient from laboratory tests is difficult. In this study, theoretical evidence is used with experimental observation to determine the range of increase in permeability during liquefaction. This range is the basis for selecting the constant for the variable permeability function that controls the value of the permeability coefficient during liquefaction applied for all simulations in this study. 2.1. Theoretical evidence Darcy’s law for the rate of creep flow through porous media is valid for a wide range of soil types for the normal range of hydraulic gradients in geotechnical problems. Despite its simplicity, this law calls for suitable estimation of hydraulic conductivity (k). The fundamental equation for k combines three influence factors as follows [23]:

k ¼ ff f v fs

ð1Þ

where ff is the effect of fluid properties (density and viscosity of the fluid); fv, is the void space (interconnected effective porosity of the mass and tortuosity of the flow paths); and fs is solid grain surface characteristics (specific surface area, roughness, and roundness of particles). For the permeability of a specific saturated sand deposit, only changes in fv are responsible for the variation of permeability. Porosity of the soil mass and tortuosity of the flow paths through the porous medium are the main factors affecting the void space properties. Permeability is related to the effective porosity of a soil mass, which is defined as the fraction of the volume of interconnected voids to the total volume [24]. Interconnected voids in a given soil mass can be visualized as a number of passages through which fluid flows. It is easy to understand why soil permeability increases as effective porosity of the mass increases. The value for effective porosity is always less than the total porosity, which is the ratio of the volume of all pores in the mass to the volume of the mass. Another important feature of flow through porous media is tortuosity of the flow. A common definition of tortuosity (T) is the ratio of the length of the actual path of the fluid particles to the shortest path length in the direction of the flow. Tortuosity is always greater than one and decreases as porosity increases [24]. A lower tortuosity factor results in higher permeability of the soil mass. The permeability coefficient (k) can be related to the void ratio and tortuosity factor through the Kozeny–Carman (K–C) equation. This equation has been proven to predict fairly good values for the permeability coefficient of sandy soils and incorporates both the void ratio (e) and T:

k¼

cf 1 e3 2 2 lf K o S0 T 1 þ e

ð2Þ

where cf is the unit weight of the fluid; lf is the fluid viscosity; Ko is the pore shape factor; and S0 is the specific surface. Variation in soil permeability over the course of liquefaction can be justified by tracking changes in void space properties during earthquake loading. Shaking a saturated sand deposit eventually results in a decrease of the whole void ratio of the soil mass. On the other hand, the spreading of the liquefied region throughout the soil mass during shaking causes a loss of contact between the

213

sand grains and creates additional pathways for water flow. In other words, the effective (interconnected) porosity of the soil mass temporarily increases during liquefaction of the soil, especially in the upper layers of the sand deposit, through which high hydraulic gradients take place. Quicksand conditions are probable when the granular layer under gravity loads is subjected to a critical upward pore fluid flow during liquefaction. Therefore, an increase in soil permeability during liquefaction can be expected, considering the increase in effective porosity of the soil mass. Measurement of effective (or interconnected) porosity is not routine in geotechnical engineering. Several hydrological studies have evaluated the effective porosity of soils based on field and laboratory techniques [25]. McWorter and Sunada [26] is a widely accepted study that compared total and effective porosity values for different geomaterials. According to this reference, the mean values of the total and effective porosity for fine sand were about 0.43 and 0.33, respectively. Dead-end pores gradually collapse during shaking and, as a consequence, the effective porosity of the soil increases from its initial value (mean value of 0.33). Assuming that all dead-end pores are removed from the soil mass during liquefaction, the effective porosity of the soil mass reaches the total porosity value (mean value of 0.43). By assuming a level of change in effective porosity and accepting the validity of the K–C equation, the permeability coefficient of sand increases more than three-fold during liquefaction. Another effect of liquefaction on soil permeability is alteration of the tortuosity factor through the liquefied soil. Direct measurement of T is not possible using conventional laboratory tests. A simplified mathematical model, however, shows that the maximum tortuosity for a mono-sized assembly of spherical particles is p/2 [24]. In practical problems where the K–C pffiffiffi equation for granular soils is used, T is usually taken to equal 2 [27]. During the liquefaction state, T decreases significantly as the sand grains float in the surrounding pore water. The K–C equation assumes that the permeability coefficient increases two- or three-fold if T approaches unity. It can be concluded that the permeability coefficient increases significantly during earthquake loading of saturated sand deposits. Although the authors found no straightforward method to estimate the level of permeability change during liquefaction, the issues discussed herein indicate that an increase in the permeability coefficient of saturated sand during liquefaction is in the order of 6 to 9 times the initial (static) value, considering the superimposed effects of effective porosity and changes in T. 2.2. Experimental observations Several experimental studies quantitatively evaluated the increase in permeability of saturated soil during cyclic loading. Arulanandan and Sybico [21] studied variations in the permeability coefficient by measuring changes in the electrical resistance of saturated sand deposits during centrifuge tests. They introduced pore shape factor and tortuosity as the main causes of the growth of the permeability coefficient during liquefaction, as shown in the K–C equation. According to their study, the average increase in the permeability coefficient of sand is 6–7 times its initial value. Jafarzadeh and Yanagisawa [5] conducted shaking table tests to examine settlement of saturated sand columns under cyclic loading. They assumed laminar flow through the soil mass during the application of input motion and employed Darcy’s law to calculate the increased permeability coefficient by measuring the cumulative water flowing upward through the sand column. They concluded that the average increase in the permeability coefficient was 5-6 times greater than that for the static condition. Ha et al. [28] conducted shaking table tests on five kinds of sand with different effective grain sizes and coefficients of uniformity.

214

H. Shahir et al. / Computers and Geotechnics 55 (2014) 211–223

They examined the dissipation pattern of excess pore pressure after liquefaction using analysis of the rate of dissipation. Permeability during dissipation was calculated using measured settlement and dissipation velocity. They assumed linear variation for the permeability coefficient over time during the buildup and dissipation phases and indicated that the increase in the permeability coefficient during liquefaction was 1.4–5 times higher than the initial permeability for the sands tested. Recently, Su et al. [29] estimated the increase in the permeability coefficient of saturated sand subjected to earthquake loading using the results of a centrifuge test. Using the law of conservation of mass, they obtained the average discharge velocity from the rate of surface settlement. By introducing this discharge velocity and the hydraulic gradient of the flow calculated from recorded excess pore pressures into Darcy’s law, the increased permeability coefficient was back-calculated and found to be as much as 6 times greater than its static value. They then conducted numerical simulations using the constant increased permeability coefficient that showed good agreement with experimental measurements. Except for the results of Ha et al. [28] that imply that in-flight permeability of saturated sand at the time of liquefaction state is a maximum of 5 times its static value, the results indicate that the average increase in the permeability coefficient during testing is 5–7 times the initial value. During each experiment, the permeability coefficient gradually increased from its initial static value to a maximum in-flight permeability value that corresponded to the liquefaction condition. At the end of shaking, excess pore pressure dissipated and the permeability coefficient decreased until it reached its initial value again. It is evident that in-flight permeability is significantly higher than the average permeability. To obtain a rough estimate of inflight permeability, two simple assumptions were made: (a) linear variation for the increase and decrease of the permeability coefficient over time, and (b) the duration of liquefaction state was much less than the duration of the buildup and dissipation period. Based on the ideal assumptions of triangular variation of the permeability coefficient over time, it can be shown that the ratio of in-flight permeability to average permeability is about 2. Consequently, it is possible to estimate the in-flight permeability of saturated sand at the time of liquefaction to be 10–14 times that of its initial static value. Both experimental observation and theoretical evidence confirm that the ratio of the permeability coefficient during liquefaction to its initial value is in the order of 10. The variable permeability function applied to the proposed numerical model takes into account variations in permeability (see Section 4). 3. General aspects of numerical modeling 3.1. Finite element formulation for porous media

where M is the mass matrix; U is the solid displacement vector; B is the strain-displacement matrix; rh is the effective stress tensor; Q is the discrete gradient operator coupling the motion and flow equations; P is the pore pressure vector; S is the compressibility matrix; and H is the permeability matrix. The vectors f(s) and f(p) include the effects of body forces, external loads, and fluid fluxes. Eq. (3) is the equation of continuity of motion; the first term is the inertia force of the mixture followed by the internal force caused by soil skeleton deformation and by the internal force caused by pore fluid pressure. Eq. (4) is the equation for continuity of fluid flow; the first and third terms represent the rate of volume change for the soil skeleton and the fluid phase, respectively, and the second term is the rate of pore fluid seepage. In this study, numerical simulation of saturated sand behavior under seismic loading was done by numerical integration of Eqs. (3) and (4) and their finite element formulation using the Open System for Earthquake Engineering Simulation (OpenSees), an open-source finite element framework developed at the Pacific Earthquake Engineering Research Center. This evolving software is a reliable tool for simulating the seismic response of structural and geotechnical systems. The version of OpenSees used in this study used the variable permeability model [15]. 3.2. Constitutive law A critical state two-surface plasticity model was employed to model sand behavior. This model was originally developed by Manzari and Dafalias [31] and refined by Dafalias and Manzari [13]. The formulation is based on the bounding surface plasticity theory [32] within a critical state soil mechanics framework [33] that yields a comprehensive multi-axial constitutive model for simulating the monotonic and cyclic behavior of sand. The formulation of the model is defined by bounding (peak), dilatancy, and critical surfaces. A schematic representation of these surfaces in the stress ratio p-plane is shown in Fig. 1. In this model, the isotropic hypoelasticity assumption was adopted using elastic shear modulus as a function of current pressure and void ratio [34]:

G ¼ G0 Pat

 0:5 ð2:97  eÞ2 P Pat 1þe

where G0 is the model parameter. The yield surface is a circular cone with its apex at the origin. The size of the yield surface is normally considered a constant (no isotropic hardening) and is a small value for most applications. The critical surface directly corresponds to the critical stress ratio in the triaxial space. The critical state of a soil is attained when the stress ratio g = q/p equals the critical stress ratio (M), which is a material constant. Similarly, the bounding and dilatancy lines are introduced by defining the peak (bounding)

In this study, a fully coupled u–P formulation was applied to model the two-phase porous medium of saturated sand. The primary unknowns in this formulation are displacement of the solid phase (u) and pore fluid pressure (P). The fully coupled effective stress u–P formulation is a simplified case for the general set of equations governing the behavior of saturated porous media. This formulation is applicable to dynamic problems where high-frequency oscillations, such as soil deposits under earthquake loading, are not important [30]. Using the finite element (FE) method for spatial discretization, the u–P formulation is:

€þ MU

Z

BT r0 dV  QP  f ðsÞ ¼ 0

ð3Þ

V

Q T U_ þ HP þ SP_  f ðpÞ ¼ 0

ð4Þ

ð5Þ

Fig. 1. Schematic representation of the two-surface model in the p-plane.

H. Shahir et al. / Computers and Geotechnics 55 (2014) 211–223

215

stress ratio (Mb) and dilatancy stress ratio (Md). In the current model, Mb and Md are related to the critical stress ratio using the state parameter as follows:

of liquefaction, when the permeability coefficient equals a times the initial permeability (ki):

M b ¼ M expðnb wÞ;

kl ¼ aki

M d ¼ M expðnd wÞ

ð6Þ

where nb and nd are positive material constants. The state parameter proposed by Been and Jefferies [35] is w = e  ec where e is the current void ratio of the soil element and ec is the critical void ratio corresponding to the existing confining stress (Pc). In the current model, the power relation was used to define the critical void ratio as:

 ec ¼ e0  kc

Pc Pat

n ð7Þ

where e0, kc, and n are critical state constants. In Fig. 1, the bounding and dilatancy surfaces are delineated by dashed lines indicating their change with w, and the critical surface by a solid line. All surfaces are fully determined by the value of w and this increases the numerical efficiency of the model. The state parameter also includes the combined effect of density (void ratio) and confining stress. One main feature of the current constitutive model is its applicability to all densities and confining pressures using the same set of material constants. The plastic modulus (Kp) and dilatancy coefficient (D) are related to the distance from the bounding and dilatancy surfaces as follows:

Kp ¼

2 ph b : n 3

D ¼ Ad d : n

ð8Þ ð9Þ

The vectors b and d in Fig. 1 are the vectors connecting the current stress state to its image on the bounding and dilatancy surfaces, respectively. The mean effective stress is p and h is a positive scalar-valued function. Ad is a function of the effects of the fabric change phenomenon that arises during stress increment reversal after dilative plastic volumetric strain. The distance-dependent plastic modulus is the main feature of the classical bounding surface model [32]. In the current model, the dilatancy coefficient is defined based on distance dependency in the bounding surface model. Dafalias and Manzari [13] compared the values of their model constants for Toyoura sand with laboratory data for a wide range of pressures and densities and presented the proper values. Shahir et al. [15] and Shahir and Pak [36] calibrated the model constants for Nevada sand using monotonic and cyclic triaxial test data performed by Arulmoli et al. [37]. Full mathematical formulation of the model can be found in Dafalias and Manzari [13]. 4. Variable permeability model Changes in effective (interconnected) porosity of the soil mass and tortuosity of the flow paths through the porous medium play a major role in altering the permeability of the saturated sand during earthquake loading. The rate of change depends on the type of contact between the soil particles (macroscopically on the effective stress state). Connecting the permeability coefficient to the excess pore pressure ratio (ratio of excess pore pressure to initial vertical effective stress of soil) can help clarify microscopic events in macroscopic analyses. Shahir et al. [15] suggested a general variable permeability function for excess pore pressure ratio (ru). Their equation shows a gradual increase in permeability up to the onset

kb ¼ ki ð1 þ ða  1Þrbu1 Þ r u < 1 during pore water pressure build up phase r u ¼ 1 during liquefaction state

kd ¼ ki ð1 þ ða  1Þrbu2 Þ ru < 1 during pore water pressure dissipation phase

ð10Þ where ki is the initial (static) coefficient of permeability and a, b1 and b2 are constants. Shahir et al. [15] implemented the above formulation using OpenSees to update the coefficient of permeability at the end of each time step during seismic analysis and applied it to simulate the behavior of a saturated sand layer subjected to earthquake loading in a centrifuge experiment. By comparing the numerical results with centrifuge experiment records, they calibrated the constants as follows:

a ¼ 20; b1 ¼ 1; b2 ¼ 8:9 Using Eq. (10) with above values as constants has two main drawbacks. First, using (a = 20) in this equation causes the permeability coefficient to equal 20 times the initial permeability during liquefaction (ru < 1), which is much higher than that reported in some experimental observations and higher than that attained from theoretical evidence (see Section 2). Second, using the value of 1.0 for b1 allows a linear relationship between the permeability coefficient and the excess pore water pressure ratio, which results in a considerable increase of soil permeability before liquefaction (ru = 1). According to the established theory for variation of permeability during earthquake loading of saturated sand (see Section 2), however, the maximum rate of increase for permeability reasonably occurs at ru values near unity once the separation of contacting sand particles begins. Values greater than unity for b1 seem to express the phenomenon more realistically. The applied variable permeability model follows the general formulation developed by Shahir et al. [15] with some modifications. As stated previously, a = 10 was assumed in numerical analysis of the permeability coefficient during complete liquefaction. An appropriate value for b1 is also critical for accurate simulation of the settlement and rate of excess pore pressure buildup. Although it is reasonable to assume the value of b1 to be greater than unity, there is no convincing logic to determine its exact value. Rahmani et al. [38] utilized the original formulation proposed by Shahir et al. [15] (where b1 = 6.0 was assumed) for numerical simulation of a centrifuge experiment and reported good agreement between simulation results and experiment measurements. The value of b2 is a key parameter for accurate simulation of excess pore pressure variation during the dissipation period. Based on a series of shaking table tests on five type of sand, Ha et al. [28] indicated that the average rate of permeability decrease during dissipation (end of shaking to completion of consolidation) is almost one third of the rate of permeability increase during the buildup phase. Numerical simulation of this requires applying a value for b2 that is greater than that for b1 in the variable permeability function. In this study, parameters b1 and b2 were calibrated such that the results of numerical simulations give a better match to the experimental records (see Section 5.3.1). Based on the results, appropriate values for the constants in the general variable permeability function proposed by Shahir et al. [15] are b1 = 2, b2 = 10, and a = 10. The following equation is achieved for numerical simulation of saturated sand behavior under cyclic loading (Fig. 2) by substituting the above values into Eq. (10):

216

H. Shahir et al. / Computers and Geotechnics 55 (2014) 211–223

Fig. 2. Schematic graph of proposed permeability function.

kb ¼ ki ð1 þ 9r 2u Þ

r u < 1 during pore water pressure build up phase

kl ¼ 10ki

r u ¼ 1 during liquefaction state

kd ¼ ki ð1 þ 9r 10 u Þ

r u < 1 during pore water pressure dissipation phase

ð11Þ 5. Validation of the proposed model 5.1. Benchmark centrifuge tests In this study, the results of four well-documented centrifuge experiments were applied as benchmarks for validation of the numerical model. The benchmark tests were selected to represent different conditions for relative density, initial permeability, maximum base acceleration, and height of sand layer. 5.1.1. Benchmark experiment 1 Centrifuge model test No. 1 was selected from the VELACS project conducted by Taboada and Dobry [16] at Rensselaer Polytechnic Institute. In this test, a laminar box containing a uniformly leveled layer of Nevada no.120 sand (Gs = 2.67, emax = 0.887, emin = 0.511, D10 = 0.08 mm) with a relative density of approximately 40% that is fully saturated with water was subjected to a centrifugal acceleration of 50 g. This centrifugal acceleration resulted in prototype soil permeability 50 times greater than the permeability of the sand specimen. The height of the sand layer was 10 m in the prototype scale. The laminar box was excited horizontally at the base using the target prototype accelerogram with a peak acceleration of 0.23 g. 5.1.2. Benchmark experiments 2 and 3 A series of four centrifuge experiments on buildings situated at the top of a liquefiable layer in a rectangular flexible shear beam container (FSB3) were performed using the large centrifuge facility at the Center for Geotechnical Modeling at the University of California, Davis by Dashti [17]. Free field results for two of these four models were used in the present study to validate the numerical model (denoted here as experiments 2 and 3). Experiment 2 was for a liquefiable soil layer (Nevada sand; Dr = 30%; Cu  2; Gs = 2.65; D50 = 0.14 mm) with a prototype thickness of 3 m located beneath a 2 m thick Monterey sand layer placed to minimize capillary rise during liquefaction. The centrifugal acceleration was 55 g and the models were saturated with fluid having a viscosity 22 times greater than the viscosity of water so that the permeability at prototype scale would be twice real soil permeability. An actual earthquake ground motion record from the Kobe Port Island station during the 1995 Kobe earthquake with peak base accelerations of about 0.62 g (large Kobe) was applied to the base of the model. The water level was 1.1 m below the surface. In experiment 3, the thickness of the liquefiable layer was the same as for experiment 2, but the relative density of sand was 50% and the specimen was

subjected to the moderate Kobe Port Island event with a peak acceleration of 0.17 g. 5.1.3. Benchmark experiment 4 This centrifuge experiment was done by Gonzalez et al. [18] at Rensselaer Polytechnic Institute. This centrifuge model consisted of uniform saturated Nevada no. 120 sand placed in a laminar box at a relative density of 55% and topped by a steel plate. The soil deposit height was 24 m in the prototype scale. The steel plate applied a uniform prototype surcharge of 140 kPa at the soil surface. Several holes were made through the steel plate to provide drainage from the soil surface. The pore fluid used in this experiment had a viscosity 40 times greater than that of water and the model was spun up to a centrifuge acceleration of 80 g. The soil column was excited horizontally at the base by 50 cycles of a sinusoidal acceleration with prototype frequency of 1.5 Hz and peak acceleration of 0.25 g. Table 1 provides a summary of the important aspects of the benchmark experiments. Input motions for all experiments are presented in Fig. 3. 5.2. Description of numerical models In all benchmark experiments considered in the current study, a leveled uniform sand layer placed in a laminar box was subjected to one-directional earthquake excitation for which one-dimensional behavior is expected. Experimental records of pore pressure and/or deformations at similar elevations have been reported [18] and confirm that numerical modeling of a single soil column is sufficient to simulate the main features of soil behavior under horizontal base excitation. In this study, a column of 8-noded cubic elements with u–p formulation was applied for all numerical analyses. Each node had 3 displacement and one pore water pressure degree(s) of freedom; however, lateral nodes at equal depths were tied together in the direction of excitation and all nodes were fixed in the opposite horizontal direction. The mesh base nodes were fixed in all directions. Pore water pressures were free to develop for all nodes except those at the ground surface. The values of the constants for the constitutive model for Nevada sand are shown in Table 2. 5.3. Results and discussion 5.3.1. Evaluation of parameters b1 and b2 As stated in Section 4, the value of b1 (the power of ru in the variable permeability function for the buildup phase) should be greater than unity to ensure that the maximum rate of increase in permeability occurs at excess pore pressure ratios close to unity. Since there is no explicit approach to determine the exact value of this constant experimentally or theoretically, numerical simulation was employed in this study to estimate the proper value of b1.

217

H. Shahir et al. / Computers and Geotechnics 55 (2014) 211–223 Table 1 Summary of benchmark experiments used for numerical modeling. Researcher(s)

Tabaoda and Dobry Gonzalez et al. Dashti Dashti

Model characteristics

Input motion

Thickness of liquefiable layer* (m)

Spinning acceleration (g)

Fluid viscosity**

Permeability coefficient*(m/s)

Relative density (%)

Type

amax (g)

Duration*** (s)

10

50

1  lwater

3.30  103

40

Earthquake

0.23

11

4

24

80

40  lwater

1.17  10

55

Sinusoidal

0.25

34

3

55

22  lwater

1.42  104

30

0.62

15.5

3

55

22  lwater

1.22  104

50

Large Kobe 1995 Earthquake Moderate Kobe 1995 Earthquake

0.17

8.5

*

Values in prototype scale. Values in model scale. *** Time between first and last exceedance of acceleration from 0.05 g. **

Fig. 3. Input motions applied in benchmark experiments.

For this purpose, analyses for different values of b1 were done to simulate the benchmark centrifuge experiments. Figs. 4 and 5 indicate the evolution of settlement and excess pore pressure (depth 7.5 m) over time for values of b1 obtained from numerical simulation of benchmark experiment 1. Fig. 4 shows that using lower values of b1 in numerical simulations increased the surface settlement of the sand layer. Fig. 5 shows that using a higher value of b1 in numerical simulation resulted in higher values for peak excess pore pressure. Fig. 6 compares the results of numerical simulations of benchmark experiment 1 with records of excess pore pressure and surface settlement to obtain a suitable value for b1 to find the best match to the experimental records. In this figure, the settlement factor (ratio of final settlement from numerical modeling

to experimental record) and maximum excess pore pressure factor (ratio of maximum excess pore pressure from numerical modeling to experimental value) are presented for different values of b1. The maximum excess pore pressure ratio was based on the maximum excess pore pressure at the bottom of the liquefiable sand layer. The results indicate that b1 = 1.0 for the variable permeability function gave the best result for predicted settlement (settlement ratio of about 0.8). On the other hand, numerically predicted excess pore pressure improved by selecting higher values for b1. It was important in this study to achieve a satisfactory level of accuracy when simulating both settlement and excess pore pressure using the same set of function constants. From this figure, it can be resolved that b1  1.5 gives the best results; however, similar sensitive analyses conducted for other benchmark problems indicate

218

H. Shahir et al. / Computers and Geotechnics 55 (2014) 211–223

Table 2 Dafalias-Manzari material parameters for Nevada sand [36]. Parameter function

Parameter index

Value

Elasticity

G0

m

150 0.05

M c kc e0 n

1.14 0.78 0.027 0.83 0.45

Critical state

Yield surface

m

0.02

Plastic modulus

h0 ch nb

9.7 1.02 2.56

Dilatancy

A0 nd

0.81 1.05

Fabric dilatancy

Zmax cz

5 800

Fig. 5. Evolution of excess pore pressure over time at depth of 7.5 m for b1 (benchmark experiment 1).

Fig. 4. Evolution of settlement over time for b1 (benchmark experiment 1).

that the optimum value of b1 is generally in the range of 1 to 3. An average value of b1 = 2, thus, proposed for the variable permeability function. A similar approach was adopted for the calibration of b2. As mentioned, the value of b2 plays the key role in proper simulation of excess pore pressure dissipation after the end of shaking. To evaluate the effects of different values for b2 on the variation in excess pore pressure over time, numerical modeling of benchmark experiment 1 was done. Fig. 7 shows that higher values of b2 in the numerical simulation increased the excess pore pressure during the dissipation phase. Fig. 8 shows the results of numerical modeling of benchmark experiment 1 for different values of b2 compared to experimental records. Excess pore pressure factor (ratio of maximum excess pore pressure from numerical modeling to experimental values) at two moments during the dissipation period is depicted against b2. The results confirm that using a greater value for b2 in the variable permeability function increased the accuracy of the predicted excess pore pressures during the dissipation period. The results, however, also indicate that values greater than 10 do not alter the numerical prediction significantly. Based on these findings, b2 = 10 was used in the proposed variable permeability model.

Fig. 6. Variation of settlement and excess pore pressure factor vs. b1 (benchmark experiment 1).

Numerical modeling of all four benchmark experiments using the variable permeability function with the proposed constants validated the proposed numerical model. 5.3.2. Simulation of centrifuge tests using Nevada sand The results of numerical analysis simulating experiment 1 are presented in Figs. 9 and 10. Fig. 9 compares the time history of surface settlement from the numerical simulation using the proposed variable permeability equation with the experimental time history and numerical simulation using constant (static) permeability. The results show that the numerical ultimate settlement (t = 50 s) is about 27% less than the recorded value, but the proposed model significantly improved the results of conventional numerical modeling using the constant permeability coefficient (settlement with about 64% error). The end of strong shaking, defined as the last exceedance of the acceleration from 0.05 g, is depicted on the settlement time history. It can be seen that the amount of settlement at the end of shaking is about 81% and 87% of the value of final settlement for numerical and experimental results, respectively. Fig. 10 shows variation in generated excess pore pressure by depth at the end of strong shaking. The excess pore pressure was under-predicted with a maximum difference of about 15% from

H. Shahir et al. / Computers and Geotechnics 55 (2014) 211–223

219

Fig. 7. Evolution of excess pore pressure over time at a depth of 7.5 m for different values of b2 (benchmark experiment 1).

Fig. 10. Comparison of numerical and experimental excess pore pressure development by depth at the end of shaking (benchmark experiment 1).

Fig. 8. Variation of excess pore pressure factor versus b2 (benchmark experiment 1).

Fig. 11. Measured and computed peak horizontal displacement by depth (benchmark experiment 1).

Fig. 9. Measured and computed time histories for settlement (benchmark experiment 1).

the experimental results. Although the focus of this study is improvement of numerical models for predicting liquefaction-induced surface settlement, which is important for leveled ground under earthquake loading, the model also simulated other features of the dynamic response with a good degree of accuracy. Fig. 11 shows the profile of measured and computed peak horizontal displacement by depth. The maximum discrepancy between numerical and experimental results is about 18% at the ground surface.

220

H. Shahir et al. / Computers and Geotechnics 55 (2014) 211–223

Fig. 12. Measured and computed time histories for settlement (benchmark experiment 2).

Fig. 14. Measured and computed time histories for settlement (benchmark experiment 3).

Fig. 13. Measured and computed time histories for excess pore pressure at middle of liquefiable sand layer (benchmark experiment 2). Fig. 15. Measured and computed time histories for excess pore pressure at middle of liquefiable sand layer (benchmark experiment 3).

The results of numerical simulation for benchmark experiment 2 are depicted in Figs. 12 and 13 for variation of settlement and excess pore pressure over time. Although the numerical simulation does not show abnormal fluctuations for settlement for the centrifuge test, the surface settlement of the model at the end of excitation was predicted using the proposed variable permeability model. Fig. 13 shows that the model maximum excess pore pressure for the middle of the liquefiable sand layer using the proposed model is in good agreement with the experimental record. Some researchers have assumed a constant increase in permeability coefficient to account for permeability variation in liquefaction modeling. For example, Arulanandan and Sybico [21] used an increased permeability coefficient of 3.7 ki (initial value for permeability coefficient) to simulate a centrifuge experiment where the variation of permeability was recorded. They obtained good agreement between the numerical results and the experimental measurement for settlement. Balakrishnan [22] employed an increase factor of 10 for permeability to adjust the results of the numerical modeling for centrifuge measurements. This factor brought the settlement value closer to the observed value, but caused a significant reduction in the peak and residual pore pressures, causing

large variation from the experimental measurements. Taiebat et al. [12] used a fully coupled effective stress numerical model with a critical state two-surface plasticity model to simulate the VELACS model no. 1. They showed that increasing permeability to 4 times the initial value gave the best fit to the centrifuge data. They concluded that the coefficient of permeability should not be a constant parameter during shaking and drainage processes when they compared the predicted and measured short-term and longterm pore pressures. Figs. 12 and 13 show the results of numerical simulation assuming a constant increase in the permeability coefficient and the results of the proposed variable permeability model. A constant increase in permeability of 5 ki was used as the average value for initial and in-flight permeability (10 ki according to Eq. (11)). The results show that using a constant increase in the permeability coefficient in liquefaction simulation improves the prediction of settlement over using a constant initial permeability coefficient. It cannot properly simulate the response of pore water pressure, especially during the dissipation phase, which is in agreement with

H. Shahir et al. / Computers and Geotechnics 55 (2014) 211–223

221

Fig. 18. Comparison of volumetric strain from numerical simulation for variable permeability function and constant permeability. Fig. 16. Measured and computed time histories for settlement (experimental benchmark 4).

for surface settlement and pore pressure variation during buildup, liquefaction, and dissipation phases. Figs. 14 and 15 present similar results for benchmark experiment 3. For this model, the rate of settlement was identical for the numerical and experimental results. The predicted final settlement and maximum excess pore pressure for the developed model was about 10% less than that measured in the lab. Fig. 16 compares the time history of surface settlement observed in experiment 4 and the numerical results from simulations using the variable permeability model. The numerical curve of the model traces the trend of the experimental time history. The results show that the numerical model predicts the settlement of the sand layer at the end of the excitation phase with less than 15% error. Shahir et al. [15] increased the factor of permeability only by adjusting the results of the numerical modeling (such as settlement) according to the measurements of centrifuge experiment 4. Using this method, they proposed that the permeability coefficient should be 20 times the initial permeability during liquefaction, which contradicts available experimental observations and cannot be justified by existing theories. The settlement in experiment 4 predicted by the modified variable permeability model is slightly less than that from the model proposed by Shahir et al. [15]. The revised function has the advantage of predicting settlement and pore pressures to an accuracy suitable for experiments under different conditions using reasonable values for the variable permeability function. Fig. 17(a) and (b) shows the excess pore pressure time histories for numerical and experimental models at two depths (13.0 and 3.9 m, respectively). The excess pore pressure time histories at different depths indicate that liquefaction spreads through the soil from top to bottom of the specimen and requires about 7 s to propagate from a depth of 3.9 to 13.0 m. A similar pattern having nearly the same propagation velocity was observed in the numerical results with a few seconds delay over the experimental records. At both depths, the experimental and numerical time histories intersect after 25 to 30 s. Fig. 17. Measured and computed time histories for excess pore pressure for experimental benchmark 4 at depths of (a) 13 m; and (b) 3.9 m.

results reported by other researchers. On the other hand, the results obtained by the numerical model equipped with a variable permeability function with the proposed constants demonstrated good capability in modeling the main features of soil response

5.3.3. Predicted settlement for variable permeability function vs. constant permeability Fig. 18 compares average volumetric strain predicted by the numerical simulation using the variable permeability function and constant permeability. The figure shows that, although numerical simulations under-predict the average volumetric strains and, consequently, settlement, the proposed model greatly enhanced

222

H. Shahir et al. / Computers and Geotechnics 55 (2014) 211–223

the accuracy surface settlement predictions over the results of numerical simulations using constant permeability. It should be emphasized that proper modeling of permeability changes during cyclic loading of saturated sand removes a major cause of deviations between numerical simulation and experimental results. Other reasons can be cited for deviation observed between experimental data and the results of numerical simulations obtained using the variable permeability function. For example, since no vertical acceleration input was applied in the numerical models in this study, the discrepancy between numerical and experimental results can be attributed to the assumption made in the numerical model that does not completely correspond to the real test conditions (e.g., the measured vertical acceleration in benchmark experiment 1 is about 0.04 g [16]). Employing the same set of constitutive model constants for different Nevada sands in the experiments and the same constants for variable permeability functions for all different conditions may not be in complete accordance with the real physics of the problem. Moreover, values of some modeling parameters, such as the fluid compressibility coefficient, that have not been measured in the experiments may significantly affect the results. It is notable that all models include free field conditions (without slope) under horizontal excitation, which does not cause large shear strains. Fig. 11 shows lateral displacement of the benchmark number 1 for a surface of about 7 cm, which is equivalent to an average shear strain of 0.7%. This level of shear strain is on the order of mediumstrains and less than the level of failure strains. The level of shear strain in the other benchmarks are on the same order (benchmark 2 is about 1%; benchmark 3 is 0.2%; benchmark 4 is 0.5%). Therefore, the small deformation theory applied in this study obtains appropriate results. However, employing the large deformation theory may enhance numerical modeling results to some extent.

6. Conclusions In this study, a revised form of a recently proposed variable permeability function was implemented in a fully coupled dynamic numerical model to simulate the behavior of saturated sand under earthquake loading. For this purpose, a well-known critical state two-surface plasticity constitutive law adequately verified in simulating cyclic behavior of sands was employed in the numerical analyses. The variable permeability function which relates the permeability coefficient to the excess pore pressure ratio permits the gradual increase/decrease of permeability coefficient with generation/dissipation of excess pore water pressure. The constants of the variable permeability function were calibrated such that the results of numerical simulation give a better match to the experimental records. In this relationship, the maximum permeability coefficient taking place during the liquefaction of sand is considered to be equal 10 times of the initial (static) permeability coefficient which is in agreement with several experimental observations and also can be justified by theoretical evidence. For validation purposes, four well-documented centrifuge experiments using a uniform level layer of Nevada sand under different conditions were simulated for factors such as relative density, initial permeability, maximum base acceleration, and height of sand layer using the proposed model. Generally, the results of numerical simulations, including excess pore pressure time histories at different depths and settlement of ground surface, are in good agreement with the experimental measurements. This implies that the coupled processes of soil mass deformation and fluid flow through the mass can be appropriately simulated using the proposed numerical model. Based on the findings of this study, the proposed numerical model is capable of simulating the behavior of saturated sand

under earthquake loading with a suitable degree of accuracy. In particular, surface settlement of the sand layer can be predicted by the numerical model involving the variable permeability function with a significant improvement relative to most previous numerical models using the constant permeably coefficient. To prove the comprehensiveness of the variable permeability function with proposed constants, experiments under different conditions such as sand type and geometry of the liquefiable layer, should be examined using the established model. References [1] Castro G. Liquefaction and cyclic mobility of saturated sands. Journal of the geotechnical engineering division. ASCE, vol. 101, No. GT6; 1975. p. 551–69. [2] Seed HB. Soil liquefaction and cyclic mobility evaluation for level ground during earthquakes. J Geotech Eng Div 1979;105(GT2):201–55. [3] Poulos SJ, Castro G, France JW. Liquefaction evaluation procedure. J Geotech Eng, ASCE 1985;111(6):772–92. [4] Ishihara K. Liquefaction and flow failure during earthquakes. Geotechnique 1993;43(3):351–451. [5] Jafarzadeh F, Yanagisawa E. Settlement of sand models under unidirectional shaking. In: Ishihara K, editor. First international conference on earthquake geotechnical engineering, IS-Tokyo; 1996. p. 693–8. [6] Dobry R, Taboada VM. Possible lessons from VELACS No. 2 results. In: Arulanandan K, Scott RF, editors. Proceedings of the international conference on verification of numerical procedures for the analysis of soil liquefaction problems, vol. 2, Rotterdam: Balkema; 1994. p. 1341–52. [7] Byrne PM, Park SS, Beaty M, Sharp M, Gonzalez L, Abdoun T. Numerical modeling of liquefaction and comparison with centrifuge tests. Can Geotech J 2004;41:193–211. [8] Andrianopoulos K, Papadimitriou A, Bouckovalas G. Use of a new bounding surface model for the analysis of earthquake induced liquefaction phenomena. In: 4th International conference on earthquake geotechnical engineering. Paper No. 1443; 2007. [9] Yang Z, Elgamal A, Parra E. A computational model for liquefaction and associated shear deformation. J Geotech Geoenviron Eng, ASCE 2003;129(12):1119–27. [10] Elgamal A, Yang Z, Parra E. Computational modeling of cyclic mobility and post liquefaction site response. Soil dyn Earthquake Eng 2002;22(4):259–71. [11] Elgamal A, Yang Z, Parra E, Ragheb A. Modeling of cyclic mobility in saturated cohesionless soils. Int J Plast 2003;19:883–905. [12] Taiebat M, Shahir H, Pak A. Study of pore pressure variation during liquefaction using two constitutive models for sand. Soil Dyn Earthq Eng 2007;27:60–72. [13] Dafalias YF, Manzari MT. Simple plasticity sand model accounting for fabric change effects. J Eng Mech 2004;130(6):622–34. [14] Jeremic B, Cheng Z, Taiebat M, Dafalias YF. Numerical simulation of fully saturated porous materials. Int J Numer Anal Methods Geomech 2008;32(13):1635–60. [15] Shahir H, Pak A, Taiebat M, Boris Jeremic. Evaluation of variation of permeability in liquefiable soil under earthquake loading. Comput Geotech 2012;40:74–88. [16] Taboada VM, Dobry R. Experimental results of model no. p. 1 at RPI’’, In: Arulanandan K, Scott RF, editors. Verification of numerical procedures for the analysis of soil liquefaction problems. Rotterdam: A.A. Balkema; 1993. p. 3–18. [17] Dashti Sh. Toward developing an engineering procedure for evaluating building performance on softened ground, PhD Thesis. Berkeley: University of California; 2009. [18] Gonzalez L, Abdoun T, Sharp MK. Modeling of seismically induced liquefaction under high confining stress. Int J Phys Model Geotech 2002;2(3):1–15. [19] Ishihara K. Review of the predictions for model 1 in the VELACS program. In: Arulanandan K, Scott RF, editors. Verification of numerical procedures for the analysis of soil liquefaction problems. Rotterdam: A.A. Balkema; 1994. p. 1353–1368. [20] Schofield AN. Dynamic and earthquake centrifuge geotechnical modeling. In: Proceeding of international conference on recent advances in geotechnical earthquake engineering and soil dynamics. University of Missouri-Rolla, Mo.; 1981. p. 1081–1100. [21] Arulanandan K, Sybico J. Post-liquefaction settlement of sand. In: Proceeding of the wroth memorial symposium. England: Oxford University; 1992. [22] Balakrishnan A. Liquefaction remediation at a bridge site, PhD dissertation. Davis University of California; 2000. [23] Aubertin M, Bussière B, Chapuis RP. Hydraulic conductivity of homogenized tailings from hard rock mines. Can Geotech J 1996;33:470–82. [24] Ghassemi A, Pak A. Pore scale study of permeability and tortuosity for flow through particulate media using Lattice Boltzmann method. Int J Numer Anal Meth Geomech 2011;35(8):886–901. [25] Brady MM, Kunkel LA. A practical technique for quantifying drainage porosity. In: Proceedings of 2003 petroleum hydrocarbons and organic chemicals in ground water: prevention, assessment, and remediation, Costa Mesa, USA; 20– 22 August, 2003. [26] McWorter DB, Sunada DK. Groundwater hydrology and hydraulics. Ft. Collins: Water Resources Publications; 1977.

H. Shahir et al. / Computers and Geotechnics 55 (2014) 211–223 [27] Das BM. Advanced soil mechanics. 3rd ed. Taylor & Francis; 2008. [28] Ha I.S., Park Y.H., Kim M.M. (2003). Dissipation Pattern of Excess Pore Pressure after Liquefaction in Saturated Sand Deposits, Transportation Research Record 1821, paper No. 03–3655. [29] Su D, Li X, Xing F. Estimation of the apparent permeability in the dynamic centrifuge tests. Geotech Test J 2009;32(1). Paper ID GTJ100972. [30] Zienkiewicz OC, Shiomi T. Dynamic behavior of saturated porous media; the generalized Biot formulation and its numerical solution. Int J Numer Meth Eng 1984;8:71–96. [31] Manzari MT, Dafalias YF. A critical state two-surface plasticity model for sands. Geotechnique 1997;47(2):255–72. [32] Dafalias YF. Bounding surface plasticity. I: Mathematical foundation and hypoplasticity. J Eng 1986;112(9):966–87.

223

[33] Schofield AN, Wroth CP. Critical state soil mechanics. New York: McGraw-Hill; 1968. [34] Richart FE, Hall JR, Woods RD. Vibrations of soils and foundations. Englewood cliffs, New Jersy: Prentice Hall; 1970. 414 pp. [35] Been K, Jefferies MG. A state parameter for sands. Geotechnique 1985;35(2):99–112. [36] Shahir H, Pak A. Estimating liquefaction-induced settlement of shallow foundations by numerical approach. Comput Geotech 2010;37:267–79. [37] Arulmoli K, Muraleetharan KK, Hossain MM, Fruth LS. VELACS Laboratory testing program: soil data report. The Earth Technology Corporation; 1992. [38] Rahmani A, Ghasemi Fare O, Pak A. Investigation of the influence of permeability coefficient on the numerical modeling of the liquefaction phenomenon. ScientiaIranica 2012;19(2):179–87.