Seismic site response of layered saturated sand: comparison of finite element simulations with centrifuge test results

A numerical model based on the finite element framework was developed to predict the seismic response of saturated sand under free-field conditions. The finite element framework used a non-linear coupled hypoplastic model based on the u-p formulation to simulate the behaviour of the saturated sand. The u-p coupled constitutive model was implemented as a user-defined routine in commercial ABAQUS explicit 6.14. Results of centrifuge experiments simulating seismic site response of a layered saturated sand system were used to validate the numerical results. The centrifuge test consisted of a three-layered saturated sand system subjected to one-dimensional seismic shaking at the base. The test set-up was equipped with accelerometers, pore pressure transducers, and LVDTs at various levels. Most of the constitutive models used to date for predicting the seismic response of saturated sands have underestimated volumetric strains even after choosing material parameters subjected to rigorous calibration measures. The hypoplastic model with intergranular strains calibrated against monotonic triaxial test results was able to effectively capture the volumetric strains, reasons for which are discussed in this paper. The comparison of the numerical results to centrifuge test data illustrates the capabilities of the developed u-p hypoplastic formulation to perform pore fluid analysis of saturated sand in ABAQUS explicit, which inherently lacks this feature.

accelerations, pore pressure generation, and settlements have been made. The numerical results have often been validated against centrifuge test results. The prediction capability depends mainly on the chosen constitutive model in addition to the nature of the numerical framework. The prediction of liquefaction triggering requires the use of a coupled stress-field analysis with a suitable elastoplastic constitutive model [1].
Over the years a large number of elastoplastic models has been developed to specifically model the behaviour of saturated sands subjected to seismic loading [2][3][4][5]. The liquefaction prediction capability of these constitutive models has been evaluated by comparing the simulation results to element level tests or centrifuge experimental results for single saturated soils [6][7][8]. Jeremic et al. [9] and Tasiopoulou et al. [10] developed the element capable of modelling pore fluid and soil skeleton responses based on the u-p-U formulation. This formulation allows water accelerations to be considered which are critical to model soil-structure interactions. Jeremic et al. [9] used the constitutive model proposed by Dafalias and Manzari [2,4] to capture the liquefaction behaviour of level and sloping grounds. Tasiopoulou et al. [10] used the framework developed by Jeremic et al. [9] to capture the behaviour of saturated sands under seismic loading and compared the results to centrifuge test results. Though the numerical model was able to capture the acceleration and excess pore pressure responses, the model highly underestimated vertical displacements on top of the sand layer under free-field conditions (primarily due to volumetric strains). The authors along with other probable causes, recognized use of constant permeability as one the reasons leading to such underprediction of settlements during liquefaction. Variation of permeability affects prediction of settlement due to liquefaction, which was recognised by Shahiri et al. [11]. Ramirez et al. [12] compared seismic response of saturated sands with two different numerical frameworks namely OpenSees and FLAC with coupled BrickUP elements. The pressuredependent-multiyield version 02 (PDMY02) [13] and bounding surface SANISAND constitutive models [4,14] were used. Ramirez et al. [12] noted that calibration of material parameters was essential and determines the efficiency of predictions for boundary value problems. Two sets of calibrations were performed, first using results of element tests (Class C) alone and second also using centrifuge experimental results (Class C1). It was concluded that model calibration based only on triaxial tests was not sufficient and that calibration based on free-field centrifuge test results was helpful in making better predictions of key response parameters. Ramirez et al. [12] also noted that all the models, regardless of the calibration methodology, notably underestimated volumetric settlements in the free-field. Failure to capture changes in porosity and permeability by the constitutive models was suggested as one of the probable causes.
Research to date has noted that the choice of the constitutive model and calibration of parameters play a critical role in capturing the soil response under seismic loading [9]. The efficiency of the chosen model to capture changes in soil porosity during seismic shaking holds the key. Hence, in this work, a coupled u-p formulation based on the hypoplastic model with intergranular strains was used in ABAQUS 6.14 explicit to predict soil response of layered saturated sand under seismic excitation. Coupled u-p hypoplastic formulation was implemented as a VUMAT (User defined routine) in commercial software (ABAQUS explicit), which inherently is incapable of performing pore fluid coupled analysis. The hypoplastic model with intergranular strains incorporates void ratio as a state variable which is critical in capturing volumetric changes during seismic loading. Centrifuge tests were performed at the University of Colorado Boulder on layered saturated sand deposits subjected to 1D seismic loading. The results of one test under free-field conditions were used to access the performance of the hypoplastic model for predicting the liquefaction phenomenon in terms of accelerations, pore pressures, and importantly, settlements.

Test set-up and materials
A series of centrifuge experiments were conducted using the 5.5 m-radius, 400g-ton centrifuge facility at the University of Colorado Boulder (in the U.S.) to evaluate seismic site response and soil-structure interaction in a layered liquefiable soil deposits in the free-field as well as in presence of one or multiple structures [12,[15][16][17]. In this paper, we focus on one centrifuge experiment that modelled free-field conditions, without any superstructures present.
The total thickness of the soil deposit was 18 m in prototype scale, which composed of a 10 m-thick dense Ottawa F65 sand layer with an initial relative density (D r ) of 90% at bottom, overlaid by a 6 m thick loose Ottawa F65 sand layer with an initial D r of 40%. Finally, a 2-m-thick coarse Monterey 0/30 sand crust with an initial D r of 90% was placed at the top. To enhance the uniformity of the soil properties across the layer, an automated sand pluviator [16] was used to prepare the soil deposit inside a flexible-shearbeam container [18]. The initial soil properties for each layer are summarized in Table 1, and the soil profile with the layout of instruments is shown in Fig. 1. The three instrumented arrays were placed to measure the settlement, acceleration, and pore water pressure at the base of the container, middle of each layer, and interface between various layers.
Following the preparation of the model, a solution of hydroxy propyl methyl cellulose (HPMC) with a viscosity 70 times greater than water was used to saturate the soil deposit under the centrifuge scaling laws [16,18]. The model was then spun to a nominal centrifugal acceleration of 70g estimated at the centre of loose Ottawa F65 sand layer prior to the application of one-dimensional shaking.
A scaled version of the north-south fault normal horizontal component of the 1995 Kobe, Japan earthquake recorded at the Takatori station was applied at the base of the container, which was referred to as Kobe-L motion in this paper. Figure 2 shows the acceleration, Arias Intensity time histories, and 5% damped response spectra of the motion recorded at the base of the container.

Experimental results
The response in free-field at the central instrumentation array during the Kobe-L motion is presented to show the overall dynamic behaviour of the liquefiable deposit. The transverse acceleration and excess pore pressure time histories measured at top, bottom of liquefiable layer, and the bottom of container are shown in Fig. 3. The time window of the acceleration records is restricted to the period of strong shaking to show the details of seismic response in a porous medium. During the Kobe-L motion, the acceleration time histories displayed a de-amplification of shear waves from the bottom up, particularly at higher frequencies as excess pore pressures increased within the soil column. Initial liquefaction (defined as r u = u /σ ' zo = 1.0) was observed in all layers during the motion, but at a slightly slower rate in the lower dense layer of Ottawa F65 sand. The dilation of soil under large shear strains and low confining stresses resulted in the large amplitude spikes in acceleration response, which were visible particularly at z = 2 m and 8 m.
As shown in Fig. 4, the settlement time histories are presented at different depths within the liquefiable layer. The volumetric deformation was most significant at the ground surface and accumulated during 0-50 s as excess pore pressures built up shown in Fig. 3.

Plane strain model
A plane strain model of 18 m height and 1.5 m width with 25 Lagrangian elements was created in ABAQUS 6.14 [19], as shown in Fig. 5a. The dimensions of the model were to represent free-field conditions as in the centrifuge prototype scale model. The maximum element size at each depth was determined as per Ramirez et al. [12]. The bottom of the model was completely fixed in the vertical direction, and a lateral periodic boundary condition was imposed by constraining the opposite side nodes on the lateral boundaries to undergo the same displacement. This boundary condition is reasonable for strong earthquakes since the energy dissipation in the soil due to hysteretic damping supersedes the energy radiation from the boundaries. It is to be noted that the periodic boundary condition provides exact results in the case of level ground with free surface subjected to base shaking without any structure [20]. The assigned fixities or boundary conditions forced the soil column to follow a shear beam mechanism and roughly reproduce the conditions in a flexible-shear-beam container. The recorded base acceleration in the centrifuge was applied as input to the base of the numerical model hence, assuming a rigid rock condition at base (Fig. 5b). The long rest period was chosen in order to be also able to capture the dissipation of developed excess pore pressures with time and compare the same to as observed in the centrifuge test. The total thickness of the soil model was composed of a 10 m dense Ottawa F65 sand layer (initial D r = 90%) at bottom, overlaid by 8 m thick loose Ottawa F65 sand layer (D r = 40%). The 2 m thick Monterey 0/30 sand layer as used in the centrifuge tests was not modelled due to lack of triaxial data to calibrate the constitutive material parameters. Instead it was modelled as 2 m thick dense Ottawa F65 sand layer due to its similarity to it with hypoplastic material parameters as per Table 2 and basic material properties such as density and permeability of original Monterey sand was considered as per Table 1. The numerical calculation steps considered a linear bulk viscosity of 0.42 and quadratic bulk viscosity of 1.2 as per the work of Kelm [21] and Hamann [22] which successfully modelled 1D wave propagation in saturated

Coupled u-p based hypoplastic model
The material behavior of Ottawa F65 sand was modelled by the hypoplastic model developed by Kolymbas [23,24] and von Wolffersdorff [25], and further extended for intergranular strain by Niemunis and Herle [26]. The model is suitable for modeling the behavior of granular materials as it can capture phenomenon such as dilatancy, contractancy, material softening, the dependency of stiffness and strength on the pressure and the void ratio, as well as a different stiffness for loading, unloading, and reloading. The intergranular strain concept helps the model overcome the so-called ratcheting during cyclic loading [26], which is of high relevance as it controls the accumulation of strains over loading cycles. The state of the granular material is characterised by the effective stress and void ratio. The coupled two phase hypoplastic model based on the u-p formulation was used in order to simulate the behaviour of saturated sand [27,28]. The u-p formulation neglects the acceleration of water phase with respect to the solid phase unlike the u-p-U formulation which includes it. The solid phase is considered to be compressible. Zienkiewicz et al. [29] suggested that the u-p formulation can work well for cases considering quasi static and low frequency dynamic forces. The formulation provides sufficient mathematical and physical accuracy for modeling of earthquake loads [29]. Implementation of the u-p-U formulation where acceleration of fluid phase is also considered, would only lead to a marginal improvement [29]. The linear momentum balance equation of a mixture of the solid and water phase is given by the following. The relative acceleration of the water phase with respect to the solid phase is neglected.
with ρ r = (1 − n)ρ s + nρ w where, ρ r is the density of the mixture of the solid and the water phase, a s is the current acceleration of the solid phase, n is the porosity, ρ w is the density of water, g is an external applied acceleration and ρ s is the density of the solid grains.
The development of the total stress state σ tot with time is obtained by the definition of a hypoplastic constitutive model in an incremental formulation for the mixture of the solid and the water phase with consideration of the concept of effective stress.
with the Biot's constant α = 1 − K T K s = 1 − E 3−6ν K s where, K s is the bulk modulus of the solid grains, K w is the bulk modulus of water, v s is the current velocity of the solid phase, k is the permeability of the solid skeleton, µ w is the viscosity of fluid, K T is the bulk modulus of the solid skeleton, E is the Young's modulus of the solid skeleton and ν is the Poisson's ratio of the solid skeleton.
A non-linear dynamic analyses procedure with an explicit central-difference time integration rule was used to solve the governing equations described above. For the explicit dynamic analyses the explicit central-difference integration rule together with a diagonal or "lumped" mass matrix M is used: (1) −ρ r a s + div(σ tot ) + ρ r g = 0 where, C is the damping matrix, F ext is a nodal force vector due to an applied external load, F int is the nodal internal force vector and rest all parameters are as described before. For the computation of the mass matrix the density of the saturated soil γ sat is used and the nodal internal force vector consists of an integration of total stresses. The The coupled hypoplastic model was developed as a user defined single VUMAT capable of performing coupled pore fluid analysis, which is inherently not possible in ABAQUS explicit. During analyses every time the user subroutine is called, the current strain-increments and the old stress-state at the integration points of a finite element for a given time-increment t i are passed in at the interface. Within the above defined VUMAT, computation of the new total stress-state is split into two parts: one for the computation of the new effective stress-state and another for the computation of the new pore water pressure. Hamann and Grabe [27] have verified the developed coupled u-p formation for simple boundary value problems under both static and dynamic loading. The numerical solutions were compared to analytical solutions or to solutions obtained by the finite element method using an implicit equation solver. The benchmark problems showed a good agreement to analytical solutions. The coupled u-p formulations was found to have higher numerical stability compared to other implicit formulations especially in scenarios involving large deformations due to its explicit formulation.

Calibration of hypoplastic material parameters
The primary 8 parameters of the hypoplastic model excluding the intergranular strain parameters were determined on the basis of index and triaxial tests performed in the laboratory as per [30]. A series of strain-controlled drained monotonic triaxial compression tests were conducted on Ottawa F65. Undrained monotonic triaxial compression tests for I D 40 % (loose) and 90 % (dense) were isotropically consolidated to three levels of effective stress 50, 100 and 200 kPa. The confining stresses were chosen according to stresses experienced in the centrifuge model. All the monotonic tests were performed under a strain-controlled axial compression rate of 0.35 mm/min. The purpose of these tests was to provide data for calibrating the constitutive model at an element level before analyzing the response of the larger layered soil system. Ramirez et al. [12] demonstrated that elaborate calibration of certain sophisticated constitutive models also did not lead to satisfactory centrifuge predictions. It is to be noted that the hypoplastic constitutive routine was specifically calibrated solely against undrained monotonic triaxial tests. This was done in order to highlight the effectiveness of model parameters calibrated against simple monotonic undrained triaxial tests to predict complicated undrained seismic response. The hypoplastic material parameters of Ottawa F65 sand where calibrated against triaxial results as shown in Fig. 6. Calibration was carried out for Ottawa F65 sand for both dense (D r = 90%) and loose (D r = 40 %) conditions as shown in Fig. 6. The decrease in deviatoric stresses following the peak (i.e., softening) in dense specimens was captured well. One may note that the peak strength is fairly well reproduced by the hypoplastic for dense Ottawa F65 sand but not so well for the loose Ottawa F65 sand. This could be expected since a single value of parameter α is shown to reproduce the peak strength of these tests under both loose and dense conditions. The hypoplastic model underestimates ε v (ε 1 ) as it does include any independent parameter to control ε v (ε 1 ) behavior under drained triaxial conditions [31]. It only includes α which controls peak stress and as it better predicts peak stress under dense condition, the same is reflected in case of volumetric strains too. In general, predicting the stress-strain response of soil under undrained loading, particularly for denser samples is difficult [12]. But the hypoplastic model with intergranular strains could well predict Ottawa F65 sand undrained behavior. In the absence of such experimental tests for Ottawa F65 sand, the intergranular parameters as determined by Niemunis and Herle [26] for Hochstetten sand was adopted. The parameters have been seen to work reasonably well for finite element simulations simulating other similar sands without much variation in results [26]. The reduction of the shear modulus with increasing shear strain is controlled by the intergranular strain parameters R, βr, and χ. Low-strain maximum-shear modulus value of Ottawa F65 sand was determined as per Ramirez et al. [12]. Simulations of element level cyclic undrained triaxial test were performed. The normalized shear modulus reduction and damping v/s strain curves were obtained by simulating the stress-strain loops for different strain amplitudes. In each loop the secant shear modulus and the damping ratio were measured. Predicted variation of shear modulus and damping ratio with increasing strain by the hypoplastic model with intergranular strain were compared to values expected for Ottawa F65 sand as detailed in Ramirez et al. [12]. Parameters R, βr, and χ for Ottawa F65 sand were calibrated to best match shear modulus and damping ratio curves detailed by Ramirez et al. [12]. The chosen intergranular strain parameters were able to predict shear modulus and damping ratio curves in the range suggested by Ramirez et al. [12] as depicted in Fig. 7. The corresponding calibrated material parameters are as tabulated in Table 2.

Acceleration time history
The acceleration was compared at the bottom of the soil column (base) and top of the liquefiable layer (2 m below) between the numerical model and the centrifuge It is also noteworthy that the constitutive material parameters have been calibrated according to drained monotonic triaxial tests. The calibrated material parameters under an element test conditions seem to overestimate stiffness and accelerations in a fullycoupled boundary value problem.
Excess pore pressures Figure 9 compares the development of excess pore pressures at 2 m below the surface and at the base of model. The numerical and experimental maximum value and dissipation of excess pore pressures over time matched well for the loose Ottawa F65 sand layer (− 2 m below). The numerical model slightly underestimated the magnitude and overpredicted pore pressure dissipation for the dense Ottawa F65 sand layer. The dense sand layer showed highly dilative behaviour in triaxial tests at lower stress levels (Fig. 5), against which the material parameters were calibrated. This can lead to an underprediction of excess pore pressures. Whereas the numerical model was able to capture excess pore pressures comparatively well in the looser sand layer. The calibrated material parameters also predicted less dilatancy for looser sand at the element (triaxial) level, likely aiding the centrifuge predictions. Figure 10 describes the settlement on surface of the soil model over time. The numerical simulations captured the volumetric strain behaviour of Ottawa F65 sand in the freefield. Ramirez et al. [12] and Wichtmann et al. [31] discussed that the PDMY02 and SANISAND models underestimated volumetric settlements governed by consolidation and sedimentation mechanisms as well as those due to partial drainage. The better performance of the coupled u-p hypoplastic model can be attributed to its ability to better capture changes in soil properties during shaking such as void ratio. The implementation of the u-p formulation by the user also assures the proper functioning of the multiphase formulation, which may not be the case when one uses built-in features in the available commercial softwares.

Conclusions
A u-p based formulation coupled with hypoplastic constitutive model was developed as a single VUMAT to be used with commercial ABAQUS Explicit software, which does not possess in-built capability to model pore fluid pressures. The formulation considered the hypoplastic model with intergranular strains to model the multiphase behaviour of saturated sand. Centrifuge tests were conducted on layered saturated sands, where the model was subjected to 1D horizontal excitation at the base. Centrifuge results in terms of acceleration, excess pore pressures, and settlement were compared to the simulation results. The coupled formulation was successfully able to capture the seismic behaviour of saturated loose and dense Ottawa F65 sand based on single set of calibrated parameters from monotonic triaxial tests. The chosen hypoplastic model was able to capture the volumetric changes during seismic shaking, unlike many other constitutive models used to date. The inclusion of void ratio as a state variable in the constitutive model led to better estimation of porosity changes during seismic shaking, leading to more accurate The work demonstrated the efficiency of material parameters calibrated against monotonic triaxial tests in predicting seismic behaviour of saturated sand. This paves a way for development of a simple yet efficient framework to help practitioners perform realistic seismic simulations with commercially available software with easily calibrated material parameters. The improvisation in prediction capabilities of the routine with material parameters calibrated against advanced cyclic laboratory tests is being pursued. The developed coupled formulation is planned to be further developed to incorporate changes in permeability during seismic shaking and made compatible to perform couple u-p simulations for varied constitutive material models.