Introduction

Soil liquefaction caused by earthquake loading is a major threat to infrastructure. Loose saturated sands when subjected to seismic loading develop high excess pore water pressures that may lead to soil liquefaction. In recent years effort has been made to understand the liquefaction phenomenon on the basis of large-scale experiments and field measurements. Owing to the high cost of large-scale tests and lack of availability of field data coupled with recent advancement in the field of numerical modelling, a shift has been observed. Numerical modelling has been identified as an effective means of understanding the liquefaction phenomena and its consequences. Researchers have tried to model the free-field 1-D seismic response of saturated soils, and predictions of 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 pressure-dependent-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.

Centrifuge tests on layered sands

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 (Dr) of 90% at bottom, overlaid by a 6 m thick loose Ottawa F65 sand layer with an initial Dr of 40%. Finally, a 2-m-thick coarse Monterey 0/30 sand crust with an initial Dr 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-shear-beam 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.

Table 1 Initial soil properties in centrifuge
Full size table
Fig. 1
figure 1

(adapted from [12])

The soil profile and instrument layout in centrifuge test

Full size image

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.

Fig. 2
figure 2

The acceleration (left) and Arias Intensity time histories (middle) as well as the 5%-damped response spectrum of the Kobe-L motion (right)

Full size image

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.

Fig. 3
figure 3

The measured acceleration and excess pore pressure time histories at different depths during the Kobe-L motion

Full size image

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 ru = \(\Delta u\)/\({\sigma }_{zo}^{\text{'}}\) = 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.

Fig. 4
figure 4

The measured settlement histories at different depths during the Kobe-L motion

Full size image

Numerical model

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 Dr = 90%) at bottom, overlaid by 8 m thick loose Ottawa F65 sand layer (Dr = 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 sand. A time scaling factor (which determines the number of time steps in ABAQUS) of 0.4 was chosen after multiple iterations in order to ensure convergence.

Fig. 5
figure 5

Numerical model (left) and applied acceleration time history along with rest period (right)

Full size image
Table 2 Hypoplastic material properties of Ottawa F65 sand
Full size table

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.

$$-{\rho }_{r}{a}_{s}+\text{div}\left({\sigma }_{tot}\right)+{\rho }_{r}g=0$$
(1)

with \({\rho }_{r}=\left(1-n\right){\rho }_{s}+{n \rho }_{w}\) where, \({\rho }_{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, \({\rho }_{w}\) is the density of water, \(g\) is an external applied acceleration and \({\rho }_{s}\) is the density of the solid grains.

The development of the total stress state \({\sigma }_{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.

$$\left(\frac{\alpha -n}{{K}_{s}}+\frac{n}{{K}_{w}}\right)\frac{\partial {p}_{w}}{\partial t}+\alpha \text{div}\left({\text{v}}_{\text{s}}\right)+\frac{1}{{\rho }_{w}}\text{div}\left\{{\rho }_{w}\frac{k}{{\mu }_{w}}\left[-\text{grad}\left({p}_{w}\right)+{\rho }_{w}\left(g-{a}_{s}\right)\right]\right\}=0$$
(2)

with the Biot’s constant \(\alpha =1-\frac{{K}_{T}}{{K}_{s}}=1-\frac{\frac{E}{3-6\nu }}{{K}_{s}}\) where, \({K}_{s}\) is the bulk modulus of the solid grains, \({K}_{w}\) is the bulk modulus of water, \({\mathrm{v}}_{s}\) is the current velocity of the solid phase, \(k\) is the permeability of the solid skeleton, \({\mu }_{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 \(\nu\) 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:

$$M{a}_{s}+C{v}_{s}+K{u}_{s}={F}_{ext}$$
(3)

with:

$${F_{int}}=K u$$
(4)

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 \({\gamma }_{sat}\) is used and the nodal internal force vector consists of an integration of total stresses. The explicit central-difference integration rule satisfies the equilibrium of the equation of motion at the time \({t}_{i}\) at the beginning of a time-increment. With known nodal accelerations \({a}_{s}^{{t}_{i}}\) and deformations \({u}_{s}^{{t}_{i}}\) at the time \({t}_{i}\) and known nodal velocities \({v}_{s}^{{t}_{i-0.5}}\) time \({t}_{i-0.5}\) at the midpoint of the previous time-increment the velocity solution \({v}_{s}^{{t}_{i+0.5}}\) at the time \({t}_{i+0.5}\) at the midpoint of the current time-increment and the nodal deformations \({u}_{s}^{{t}_{i+1}}\) at the time \({t}_{i+1}\) at the end of the current time increment are calculated as per the following:

$${v}_{s}^{{t}_{i+0.5}}={v}_{s}^{{t}_{i-\mathrm{0,5}}}+\frac{\mathrm{\Delta }{t}_{i-1}+\mathrm{\Delta }{t}_{i}}{2}{a}_{s}^{{t}_{i}}$$
(5)
$${u}_{s}^{{t}_{i+1}}={u}_{s}^{{t}_{i}}+\mathrm{\Delta }{t}_{i}{v}_{s}^{{t}_{i+0.5}}$$
(6)

The nodal accelerations \({a}_{s}^{{t}_{i}+1}\) at the time \({t}_{i+1}\) at the end of the current time-increment are obtained by solving the equation of motion as shown:

$${a}_{s}^{{t}_{i+1}}={M}^{-1}\left({F}_{ext}^{{t}_{i+1}}-{F}_{int}^{{t}_{i+1}}\right)$$
(7)

where, \({F}_{ext}^{{t}_{i+1}}\) is the applied external load vector and \({F}_{int}^{{t}_{i+1}}\) is the nodal internal force vector at the time \({t}_{i+1}\).

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 \(\mathrm{\Delta }{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 ID 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 (Dr = 90%) and loose (Dr = 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.

Fig. 6
figure 6

Prediction of undrained triaxial experimental results by calibrated hypoplastic model with intergranular strain parameters

Full size image
Fig. 7
figure 7

Shear modulus reduction curves and damping ratio curves for Ottawa F65 sand as obtained from hypoplastic constitutive model parameters with intergranular strain

Full size image

Validation of numerical model

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 measurements (Fig. 8). At the top of the liquefiable layer, accelerations showed similarity for the first few cycles, but with progression of time and building up of pore pressures, the numerical model started to overpredict accelerations. This overprediction could be due partly to the underestimation of damping and overestimation of soil’s dilative behaviour. The numerical model considers no additional damping other than the inherently occurring material damping in the constitutive model and one occurring due to a pre-defined linear bulk viscosity (value = 0.42) of the explicit numerical model in ABAQUS. 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 fully-coupled boundary value problem.

Fig. 8
figure 8

Comparison of numerical and experimental results for acceleration at bottom of model (left) and -2 m from top (right)

Full size image

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.

Fig. 9
figure 9

Comparison of numerical and experimental results for excess pore water pressure at -2 m from top (left) and bottom of model (right)

Full size image

Total settlements

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 free-field. 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.

Fig. 10
figure 10

Comparison of numerical and experimental results for settlement at top of the model

Full size image

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 settlement predictions. 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.