The influence of the soil constitutive models on the seismic analysis of pile-supported wharf structures with batter piles in cut-slope rock dike

In the past decades, seaports suffered large structural disorders, economic and social losses, mostly due to earthquakes. Among the damaged port structures are pile-supported wharves founded in sloping rockfill or through the rock dikes. These structures are commonly used worldwide, and have been the subject of several experimental and numerical studies. Yang^[1] has studied the seismic response of the Seventh Street Terminal Port of Oakland using FLAC software,^[2] and used some modifications to take into account the effect of sloping ground. Several dynamic centrifuge tests were conducted at the University of California, Davis (UCD).^[3–10] All centrifuge models represented case histories, consisting of rock dike configurations and wharf deck supported by piles. These tests were carried out in order to better understand the seismic performance of wharf piles. Hwang et al.^[11] performed shaking table model tests, to reproduce the dynamic behavior of a gravity quay wall and a pile-supported wharf damaged during the 1995 Kobe earthquake. The results of shaking table model tests were compared with the field measurements and with the results of the previous model testing. McCullough^[3] and Dickenson and McCullough^[12] presented centrifuge modeling results of pile supported wharves and pile in sloping rockfill. These results were used to validate two geotechnical seismic analysis methods: Newmark sliding block method,^[13] and numerical modeling using the finite difference code FLAC.^[2] The authors also performed a numerical analysis of the Seventh Street Terminal (Berths 35 through 38), at the Port of Oakland subjected to the 1989 Loma Prieta earthquake and compared the displacement results to in-situ observations.^[14,15] Takahashi and Takemura^[16] carried out a centrifuge model test of the Takahama Wharf damaged during the 1995 Hyogo-ken Nambu earthquake, and performed a parametric study to establish the effect of the thickness of the liquefiable sand layer on the results, and to study the piled deck-caisson dynamic interaction. Dickenson et al.^[17] provided a synthesis of modeling considerations and summary of the 2D FLAC numerical modeling results of the Berth 404 at Pier 400 of Port of Los Angeles, and considered pile-deck connections and pile-rockfill interaction. The results were validated with the case history data and the results from large-scale tests. Huertas and Romanel^[18] investigated the coupled hydro-mechanical behavior of soils and studied the influence of the liquefaction on the stability of pile supported wharf in rockfill dike, considering the nonlinear elasto-plastic constitutive model UBCSand. Dickenson et al.^[19] investigated the seismic performance of pile-supported wharves subjected to long-duration ground motions and considered the evaluation of kinematic loads on pile and soil-pile interaction in rockfill dikes. Su et al.^[20] conducted 3D finite element analysis of pile-supported wharf and considered the effect of the ground slope and pile connectivity along the wharf deck. Souri et al.^[21] conducted physical modeling of wharf piles in sloping rockfill subjected to liquefaction-induced lateral spreading, and studied kinematic demands and full superstructure inertia in order to estimate bending moments of piles at different depths. Vytiniotis et al.^[22] studied the effect of prefabricated vertical drains installed behind the crest of the slope to control the lateral spreading in the fill and to reduce seismic damage of the wharf structure.

Recently, a large volume of research works used seismic fragility curves as a useful tool to evaluate the seismic performance of pile-supported wharves.^[23–30]

Several port structures use batter piles to provide lateral resistance and to limit deflections from lateral loads, such as ship berthing and seismic loads. In the past, engineers contest the use of batter piles due to the poor performance under earthquake, caused by the lack of proper design of the batter piles head connections. Actually, by resolving this problem, batter piles have smaller deformations than vertical ones. Roth et al.^[31] studied the role of batter piles in the seismic performance of wharves. Schlechter et al.^[32] carried out a series of large-scale centrifuge models, with and without batter piles of typical pile-supported wharf configurations. Razavi et al.^[33] presented some batter pile damages during severe earthquakes, and performed seismic analysis of a typical pier with and without batter piles. Heidary-Torkamani et al.^{[28, 29]} developed the seismic fragility curves of wharves with batter piles, considering different engineering demand parameters (EDPs) for each damage state, and also performed a sensitivity analysis in order to evaluate the influence of geotechnical parameter uncertainties in the seismic performance of the wharves. Li et al.^[34,35]) studied dynamic response characteristics and failure process of vertical and batter pile-supported wharf structure, established on the same conditions, and subjected to two different earthquake levels.

The soil is a complex material and exhibits nonlinearly behavior when it is subjected to very low strain levels Atkinson and Sallfors,^[36] as for example, far-field deformations around deep excavation and tunnel, and soils subjected to seismic ground motion. Researchers developed different soil constitutive models with various degrees of complexity, and it is important to choose the most representative soil behavior model, as it has a significant impact on the results.

A great amount of investigation was conducted to determine the influence of soil constitutive models at small-strain levels for soils around deep excavation and tunnel; however, less researches were done in the case of soils under seismic loading. Truty^[37] studied the deformations near the excavations in Berlin sand, using Mohr-Coulomb (MC), Hardening Soil (HS) and Hardening Soil with Small-strain stiffness (HSS) models, and presented the soil deformation, the bending moment and the deflection of the wall. Hamrouni et al.^[38] highlighted the impact of soil constitutive models on the prediction of the displacements of the surrounding soil in the execution of a shallow tunnel. Benz et al.^[39] described soil stiffness at small-strains and presented the available small-strain stiffness models. The authors also performed four numerical applications on: triaxial test, spread footing problem and excavation problem 2D and 3D, using the finite element method with and without small-strain stiffness formulation. Obrzud^[40] presented Hardening Soil standard and Hardening Soil Small-strain models in finite element program, and exposed the comparison result with Mohr-Coulomb modeling of excavation in Berlin sand, tunneling and shallow footing in Texas sand. Besseling^[41] evaluated the sensitivity of the site response to variation of the soil constitutive models. Four soil constitutive models implemented in PLAXIS software were considered: the linear elastic (LE), the Mohr-Coulomb (MC), the Hardening Soil (HS) and the Hardening Soil with Small-strain stiffness (HSS) models. Based on the comparison of site responses, the HSS model was adopted for modeling of the soil deposit in the study of the dynamic jetty analysis. Vakili et al.^[42] presented the numerical approach for adaptive constitutive modeling of soils using the finite element method, and gave results using this approach in numerical soil element tests (Oedometer and Triaxial tests) and staged tunneling construction. Vakili et al.^[43] evaluated the results for the vertical displacements of four selected nodes in shallow and deep tunnel projects using three soil constitutive models. Op de Kelder^[44] analyzed the horizontal deformations on the sheet pile wall as a result of excavation of a building pit, using the Hardening Soil and the Hardening Soil Small strain stiffness models implemented in the FE code PLAXIS 2D. For the calculation of shear modulus at very small strains (G₀) and the threshold shear strain (γ_0,7) of HSS model, the author employed the correlation of Alpan^[45] and the correlation of Benz and Vermeer.^[46] Thereafter, the obtained numerical results were compared to the inclinometer measurements. Alpan^[45] gave curves that describe the relationship between static E_stat (large-strain) and dynamic E_dyn (small-strain) stiffness moduli (E_stat with E_dyn= E₀). E₀ is the Young’s modulus at very small strains. According to Wichtmann and Triantafyllidis,^[47] it can be assumed that the correlation of Alpan estimated E_stat by the stiffness during first loading E₅₀ (E_stat = E₅₀). For the correlation of Benz and Vermeer, it is recommended to interpret E_stat as the secant modulus during large unloading and reloading cycle E_ur (E_stat = E_ur = 3E₅₀). Szerző and Batali^[48] used different constitutive laws for the soil and interface elements to analyze the numerical settlement results of piled raft foundations supporting a circular gas tank, and compared these to field measurements and two simplified methods: the Poulos-Davis-Randolph method^[49] and the Piglet software.^[50] Akbari Hamed[51] compared the experimental results of Triaxial and Oedometer tests to the numerical results in two finite element software, using various soil constitutive models, and also performed a comparison between the measured data and the numerical analysis results of the asphaltic core rockfill dam (Dam-X, Quebec). Yeganeh and Fatahi^[52] examined the effects of soil plasticity type and soil-structure interaction in the seismic analyses of the typical high rise buildings, constituted by 20-story reinforced concrete moment-resisting building with a mat foundation, using the explicit finite difference program FLAC3D.

Considering the few researches to evaluate the impact of choosing soil constitutive model on the seismic results and, in addition to the complexity of seismic behavior of wharf structures with batter piles in sloping rockfill, this study aimed to capture the more realistic seismic behavior of pile-supported wharf with batter piles in cut-slope rock dike. For this, three soil constitutive models with increasing levels of complexity were considered using finite element method. The techniques of the finite element modeling will be given in the following sections.

A two-dimensional finite element modeling was generated using the software PLAXIS 2D.^[53] The details of the modeling adopted in the dynamic analysis are listed below:

The geometry model boundaries are taken relatively far away to avoid direct influence of the boundary conditions Hamrouni et al.^[54] Horizontal fixity has been introduced at the vertical boundaries, while the bottom boundary of the model has been subjected to total fixity (both horizontal and vertical fixities).

A prescribed displacement is selected at the bottom boundary of the model in order to simulate the dynamic loading. Standard earthquake boundaries are specified at the right and left side of the model geometry, to avoid the reflection of the seismic waves on the model boundaries and perturbations.^[55] The relaxation coefficients C₁ and C₂were kept in default values 1 and 0.25, respectively.

The mesh coarseness was fine and refined around the dike with 15 node triangular elements. The influence of the mesh must be tested for each analysis. In Plaxis,^[56] the Average Element Size (AES) is calculated from the outer geometry dimensions (x_min, x_max, y_min, y_max) and the number of generated triangular elements (n_c), as presented in equation (1): (1)AES=[(xmax−xmin)(ymax−ymin)/nc]0.5AES = {\left[ {\left( {{x_{\user1{max}}} - {x_{\user1{min}}}} \right)\left( {{y_{\user1{max}}} - {y_{\user1{min}}}} \right)/{n_c}} \right]^{0.5}} Kuhlemeyer and Lysmer^[57] provided a guideline to use AES, which must be smaller than approximately one-tenth to one-eighth of the wavelength (λ) associated with the highest frequency of the input wave, equation (2): (2)AES≤λ10toλ8AES \le \frac{\lambda }{{10}}to\frac{\lambda }{8} The equation for the time-dependent movement of a system is shown in equation (3): (3)[M]U¨+[C]U˙+[K]U=F(t)\left[ M \right]\ddot U + \left[ C \right]\dot U + \left[ K \right]U = F(t) Here, [M], [C] and [N] are respectively the mass, the damping and the stiffness matrices. F(t) is the external force vector applied on the system. The displacement (U), the velocity (U̇) and the acceleration (Ü) can vary with time.

The damping matrix [C] of the system is proportional to the mass [M] and stiffness [K] matrices by means of Rayleigh coefficients, α_R and β_R,^[58,59] as presented in equation (4): (4)[C]=αR[M]+βR[K]\left[ C \right] = {\alpha _R}\left[ M \right] + {\beta _R}\left[ K \right] PLAXIS 2D uses implicit time integration based on the Newmark scheme.^[60] The displacement and the velocity at the point in time t + dt are expressed respectively in equations (5) and (6), where dt is the time step, α_N and β_N are the Newmark coefficients that determine the accuracy of the numerical time integration: (5)Ut+dt=Ut+U˙tdt+((12−αN)U¨t+αNU¨t+dt)dt2{U^{t + dt}} = {U^t} + {\dot U^t}dt + \left( {\left( {\frac{1}{2} - {\alpha _N}} \right){{\ddot U}^t} + {\alpha _N}{{\ddot U}^{t + dt}}} \right)d{t^2}(6)U˙t+dt=U˙t+((1−βN)U¨t+βNU¨t+dt)dt{\dot U^{t + dt}} = {\dot U^t} + \left( {\left( {1 - {\beta _N}} \right){{\ddot U}^t} + {\beta _N}{{\ddot U}^{t + dt}}} \right)dt In order to obtain a stable solution, the following condition in equations (7) and (8) must apply: (7)βN≥0,5{\beta _N} \ge 0,5(8)αN≥14(12+βN)2{\alpha _N} \ge \frac{1}{4}{\left( {\frac{1}{2} + {\beta _N}} \right)^2} However, during the last two decades, new and more accurate numerical methods have been developed. For example, the α-method (or HHT method),^[61] that uses the Newmark scheme to solve the modified equations of motion shown in equation (9): (9)MU¨t+(1−α)CU˙t+(1−α)KUt=(1−α)Ft++αFt−αCU˙t−dt−αKUt−dt\begin{array}{*{20}{c}}{M{{\ddot U}^t} + (1 - \alpha )C{{\dot U}^t} + (1 - \alpha )K{U^t} = (1 - \alpha ){F^t} + }\\{ + \alpha {F^t} - \alpha C{{\dot U}^{t - dt}} - \alpha K{U^{t - dt}}}\end{array} The numerical parameters α, α_N and β_N are selected such as in the following equations (10), (11) and (12) in order to obtain stable and accurate results: (10)0≤α≤130 \le \alpha \le \frac{1}{3}(11)αN=(1+α)24{\alpha _N} = \frac{{{{(1 + \alpha )}^2}}}{4}(12)βN=12+α{\beta _N} = \frac{1}{2} + \alpha For a damped Newmark integration scheme (α>0), the Rayleigh coefficient β_R should be reduced by the quantity α.dt.^[62]

The standard setting of the PLAXIS code for time integration is the damped Newmark scheme with α_N = 0.3025 and β_N = 0.6, that corresponds to α = 0.1.

The time step used in a calculation is constant and it is presented in equation (13): (13)dt=Δtn⋅mdt = \frac{{\Delta t}}{{n \cdot m}} Where ∆t is the duration of the dynamic loading, m is the max steps number (number of additional steps) and n is the number of sub steps (dynamic sub step). The upper limit of the max steps number is given equal to 10,000, and PLAXIS automatically calculates the number of sub steps. In order to visualize the complete signal, it is advised to use m equals the number of multipliers that defines the input signal and n equals 1.^[63,64]

To avoid unreliable numerical response in a finite element model, a careful evaluation of the critical time step dt_crit for the system should be executed.^[58,62]

The configuration of the wharf chosen for this study is the JCB01 model,^[10] selected from a series of centrifuge models carried out at the Center for Geotechnical Modeling of California University, Davis Campus (UC Davis).

This model consists of a wharf deck supported by seven rows of vertical piles included two pairs of batter piles, embedded in cut slope rock dike configuration. The dike is constructed by a thin layer of rockfill overlaid by a slope of loose sand (40% relative density). The side slope was 2 horizontal to 1 vertical (2:1). The foundation soils and backfill soil behind the wharf were all dense (74% relative density). The soil and structural geometry are shown in Fig. 1.

Figure 1

In numerical analysis, material behavior of the soils and structural elements are represented by material constitutive models.

4.1

Material model for soils and rockfill

The sandy soils used in the centrifuge model were Nevada sand.^[65] The rock was imported from a quarry on Catalina Island, off the coast of southern California, and crushed to an appropriate size in order to represent the rock-pile interaction.

The geotechnical parameter values used in numerical modeling are presented below according to the three soil constitutive models applied in this study:

4.1.1

Mohr-Coulomb model (MC model)

The most common material model used is the linear-elastic perfectly plastic model (Mohr-Coulomb model). It requires only five input parameters: Young’s modulus E, Poisson’s ratio ν, friction angle φ, cohesion C and dilatancy angle ψ.^[66] The stress-strain relationship is considered as elastic-perfectly plastic (Fig. 2(1)).

Figure 2

The soil is composed of granular particles, and is generally modeled in numerical analysis as a continuum. The individual particle size of the soil is in the order of 30 to 40 times less than the pile diameter or the dimensions of the foundation. In the case of rockfill used in the marine structures, the median rockfill size was 2.5 to 3 times smaller than the pile diameter. Then, the rock particles act as individual particles and less as a continuum. Diaz et al.^[67] introduced a concept of pseudo cohesion in rockfill so as to modify the modeled continuum behavior to match the measured behavior. McCullough^[3] added 15 kPa pseudo-cohesion for the rockfill in order to consider the individual rock particle interaction with the pile and took into account the difference between upslope and downslope spring stiffness. This concept was also adopted by Martin,^[68] and McCullough and Dickenson.[69] Kawamata^[70] considered the lateral reaction of the rockfill as a combination of reactions due to both friction and interlocking between large rock particles.

The FE code PLAXIS can handle cohesionless soils with cohesion equal to zero, but some options do not perform well. The material models manual of PLAXIS^[71] advises to enter at least a small value of cohesion (use C > 0.2 kPa). So, in this work, cohesion of the dense and loose sand is chosen equal to 0.5. For the rockfill, an artificial cohesion of 15 kPa and a dilation angle of 0° were applied.^[3]

The input parameters of soil used for the Mohr-Coulomb model are summarized in Table 1, and for more details refer to Boland et al.,^[10] McCullough,[3] Heidary-Torkamani et al.^[28,29] and Amirabadi et al.^[26]

Table 1

Soil parameters used for the MC model.

Parameters	Unit	Dense sand	Loose sand	Rockfill
φ	°	37	33.2	45
ψ	°	21	7	0
C	kN/m2	0.5	0.5	15
ν	-	0.3	0.3	0.2
E	kN/m2	6.032 × 104	1.638 × 104	2.232 × 104

The dry mass density was 1633 kg/m³, 1522 kg/m³ and 1611 kg/m³ respectively for dense sand, loose sand and rockfill. The material damping of soils was simulated with 5% of damping ratio (Rayleigh damping).

The permeability was 0.043 × 10^–3 m/s, 0.066 × 10^–3 m/s and 0.01 m/s respectively for dense sand, loose sand and rockfill. For water conditions, the ground was saturated. A water depth of 12.4 m, and volumetric weight of 10.3 kN/m³ (salted water) were defined.

4.1.2

Hardening Soil model (HS model)

The Hardening Soil model, developed by Schanz et al.,^[72] describes a hyperbolic stress – strain curve (Fig. 2(2)). It has the same failure criterion as the Mohr-Coulomb model and uses three types of soil stiffness for three loading conditions, namely triaxial loading, oedometer loading and unloading/reloading, as presented in equations (14), (15) and (16): (14)E50=E50ref(ccosφ−σ3'sinφccosφ+prefsinφ)m{E_{50}} = E_{50}^{ref}{\left( {\frac{{c\,\user1{cos}\,\varphi - \sigma _3^'\,\user1{sin}\varphi }}{{c\,\user1{cos}\,\varphi + {p^{ref}}\,\user1{sin}\varphi }}} \right)^m}(15)Eoed=Eoedref(ccosφ−σ3'sinφccosφ+prefsinφ)m{E_{oed}} = E_{oed}^{ref}{\left( {\frac{{c\,\user1{cos}\,\varphi - \sigma _3^'\,\user1{sin}\varphi }}{{c\,\user1{cos}\,\varphi + {p^{ref}}\,\user1{sin}\varphi }}} \right)^m}(16)Eur=Eurref(ccosφ−σ3'sinφccosφ+prefsinφ)m{E_{ur}} = E_{ur}^{ref}{\left( {\frac{{c\,\user1{cos}\,\varphi - \sigma _3^'\,\user1{sin}\varphi }}{{c\,\user1{cos}\,\varphi + {p^{ref}}\,\user1{sin}\varphi }}} \right)^m} where the E50refE_{50}^{ref} , EoedrefE_{oed}^{ref} , and EurrefE_{ur}^{ref} represent the stiffness at the reference pressure P^ref (at default setting P^ref equals 100 kN/m²). σ3′\sigma _3^' is the effective confining pressure in Triaxial test, it is negative for compression.^[73] The amount of stress dependency is given by the power m.

The HS model was implemented into PLAXIS finite element code. In many practical cases, it is appropriate to determine the stiffness moduli by using the default setting of PLAXIS: E50refE_{50}^{ref} is equal to the Young’s modulus defined in the Mohr Coulomb model, E50ref=1.25Eoedref$E_{50}^{ref} = 1.25 E_{oed}^{ref}$ and Eurref=3E50ref$E_{ur}^{ref} = 3 E_{50}^{ref}$ .^[71]

The soil parameters required in the Hardening Soil model for the materials used in this study, in addition to the soil parameters of MC model, are presented in Table 2.

Table 2

Soil parameters used for HS model

Parameters	Unit	Dense sand	Loose sand	Rockfill
E50refE_{50}^{ref}	kN/m2	6.032 × 104	1.638 × 104	2.232 × 104
EoedrefE_{oed}^{ref}	kN/m2	4.826 × 104	1.310 × 104	1.786 × 104
EurrefE_{ur}^{ref}	kN/m2	1.810 × 104	4.914 × 104	6.696 × 104
Pref	kN/m2	100	100	100
M	-	0.5	0.5	0.5
K0nc\user2{K}_0^{nc}	-	0.3982	0.4524	0.2929
Rf	-	0.9	0.9	0.9

4.1.3

Hardening Soil Small-strain model (HSS model)

The Hardening Soil model with small-strain stiffness is a modification of the Hardening Soil model. It takes account the degradation of shear stiffness with shear strain during loading (developed by Benz, 2007).^[74]

The model’s stiffness parameters in a triaxial test were presented in Fig. 2(3). The stress-strain curve for small-strains was described by a simple hyperbolic law, analogous to the hyperbolic law for larger strains.^[75]

The HSS model uses almost entirely the same parameters as HS model, except two additional parameters: The initial or very small-strain shear modulus G₀, and the shear strain level γ_0,7 at which the secant shear modulus G_S degrades to about 72% of initial shear modulus G₀.

Several works have demonstrated that when increasing strain amplitude, soil stiffness decays nonlinearly.^[36] The stiffness reduction curves exhibit a characteristic S-shape in logarithmic scale (Fig. 3), only at very small-strain levels (less than 10^–6) where the soil describes linear elastic behavior.

Figure 3

The stress dependency is given by the following equation (17). G0refG_0^{ref} is the reference shear modulus at very small-strains (ε < 10^–6). (17)G0=G0ref(ccosφ−σ3'sinφccosφ+prefsinφ)m{G_0} = G_0^{ref}{\left( {\frac{{c\,\user1{cos}\,\varphi - \sigma _3^'\,\user1{sin}\varphi }}{{c\,\user1{cos}\,\varphi + {p^{ref}}\,\user1{sin}\varphi }}} \right)^m}

The very small-strain shear modulus G₀ of non-cohesive soils can be determined using different methods, with different efforts and accuracy.^[47] In the absence of experimental data for the determination of G₀ parameter, approximations from empirical equations can be appropriate using the relation given by Hardin and Black^[78] (equation (18). Similarly, in the lack of experimental data for the determination of the threshold shear strain γ_0,7, correlations are also available, for example, a correlation between the threshold shear strain γ_0,7 and the Plasticity Index PI,^[79] or using the relationship of Hardin-Drnevich[80] as shown in equation (19): (18)G0ref=33(2.97−e)21+e(MPa)G_0^{ref} = \frac{{33\,{{(2.97 - e)}^2}}}{{1 + e}}(MPa)(19)γ0,7≈19G0[2c'(1+cos(2φ'))−σ1'(1+K0)sin(2φ')]{\gamma _{0,7}} \approx \frac{1}{{9{G_0}}}\left[ {2c'\left( {1 + \cos (2\varphi ')} \right) - \sigma _1^'(1 + {K_0})\sin (2\varphi ')} \right] where K₀ is the earth pressure coefficient at rest, and σ₁^‘ is the effective vertical stress (negative pressure).

In dynamic applications, the HSS model introduces hysteric material damping and shows a typical hysteretic behavior (Ishihara, 1996),^[81] in contrast to the Hardening Soil model.

The values of the two additional soil parameters G0refG_0^{ref} and γ_0,7 used in this analysis are presented in Table 3:

Table 3

Additional input parameters of HSS model.

Parameters	Unit	Dense sand	Loose sand	Rockӿll
γ0.7	-	4.918 × 10–4	3.629 × 10–4	4.246 × 10–4
G0refG_0^{ref}	kN/m2	1.051 × 105	9.108 × 104	1.038 × 105

4.2

Calibration of soil constitutive parameters in numerical soil element tests

The constitutive parameters of the MC model have to be calibrated to advanced models, using two soil element tests available in the two-dimensional finite element code PLAXIS. These tests are Oedometer test (one-dimensional compression) and Triaxial test. The following conditions are adopted: –

In the Oedometer test, the soil sample is subjected to four phases: two loading phases from 0 to 100 kPa and

–

200 kPa, followed by two unloading phases to 100 kPa and 50 kPa.

–

In Triaxial test, a confining pressure σ₃ of 100 kPa is applied.

The calibration is limited to the variation of elasticity modulus (E) for the MC model, which takes the values of stress dependent moduli of elasticity in the HSS model ( E50ref=1.25Eoedref$E_{50}^{ref} = 1.25 E_{oed}^{ref}$ , EurrefE_{ur}^{ref} and E0refE_{0}^{ref} , where E0refE_{0}^{ref} is the Young’s modulus at very small-strain at the reference pressure P^ref). The results of soil parameter calibrations in Oedometer and Triaxial tests are shown in Fig. 4.

Figure 4

According to the comparison between curves for both Oedometer and Triaxial tests, and for all the different soil materials, the best consistency between MC model and advanced models is observed when the value of the Young’s modulus in the MC model is assumed to be equal to E50refE_{50}^{ref} . Consequently, all the geotechnical characteristics of the soil layers presented previously will be used in the numerical analysis.

The numerical model is subjected to the 1989 Loma Prieta earthquake time history recorded at the Oakland Outer Harbor Station (Fig. 5):

Figure 5

In order to verify the numerical model using the finite element method of the PLAXIS 2D, the displacement of soils around the rock dike and the pore pressure are determined and compared to the numerical results predicted using the two-dimensional finite difference of the code FLAC.^[28]

The results of the total displacement of soils around the cut-slope rock dike for the three soil constitutive models are presented in Fig. 6, 7 and 8, respectively, for MC, HS and HSS models.

Figure 6

Figure 7

Figure 8

The maximum total displacement of soil observed for all three soil used constitutive models is located in the upper part of the loose sand embankment. The value of this maximum total displacement is obtained between 0.14 m and 0.15 m with MC and HS models, and between 0.08 m and 0.09 m with HSS model.

Heidary-Torkamani et al.^[28] have obtained 0.08973 m for the maximum vector displacement (Fig. 9). This displacement is included in the interval of the displacement determined with HSS model. It can be concluded that the results obtained from a 2D plane strain model (as seen in the paragraph 4.3 above), by using the HSS model give adequate results for soil and wharf displacements.

Figure 9

Concerning the pore pressure contour results obtained in the numerical model with PLAXIS 2D, the contours have the same shape for all three soil used constitutive models. So, it was plotted here with just one pore pressure contours (Fig. 10). The pore pressure around the zone of embedded pile is approximately 0.2 MPa. This result coincides with the total pore pressure contours presented by Heidary-Torkamani et al.^[28] in Fig. 11.

Figure 10

Figure 11

It can be justified that when the soils are subjected to dynamic loading, the MC, HS and HSS models of the finite element code PLAXIS 2D and the MC model of the finite difference code FLAC 2D^[28] give the same pore pressures, due to the fact that, all these models do not take into account the accumulation of pore pressures or liquefaction.^[71]

Fig. 12 presents the total pile displacements under the 1989 Loma Prieta earthquake for the three soil constitutive models.

Figure 12

According to the pile deflection shown in Fig 12, the total displacements at the base of the piles are small and are at about the same for all the models. Thereafter, these total displacements gradually increase until maximum values are reached in the upper part of the piles, and become different for the three models.

This is due to the fact that at depth, the piles are embedded in good soil (dense sand layer), able to take up the efforts of wharf and undergo small deformations during the earthquake. While in the dike embankment, liquefaction occurs and induces large deformations especially at the embankment top.

Analyzing the maximum total pile displacements, it can be noted that for the first three piles of the wharf and also for the batter piles, the maximum total displacements are reached at the piles head. For the other piles, the maximum total displacements are particularly localized at a height in the range of 20.2 m to 22.2 m.

The largest of the maximum total displacements for the different piles of the wharf is observed for all models at the last pile (pile 7), caused by the significant deformation of loose sand at the top of the embankment, as shown in Fig. 6, 7 and 8. The values obtained are 0.139 m, 0.131 m and 0.077 m for the MC, HS and HSS models, respectively.

Comparing the results of this study with those of the numerical analysis using the finite difference method (Fig. 9, where the maximum pile displacement equals 0.08606 m), the MC and HS models give a maximum total pile displacement overestimated by 61.51% and 52.2% respectively, while for the HSS model, the maximum total displacement is underestimated by 10.52%.

It can be concluded that the results obtained with the HSS model are close to the results of numerical analysis conducted by Heidary-Torkamani et al.,^[28] and the MC model results are larger than the results of advanced models (HS and HSS models). This is due to the fact that the MC model uses a constant Young modulus for soil and only Rayleigh damping, in contrast to the advanced models, which use three different types of stiffness modulus ( E50refE_{50}^{ref} , EoedrefE_{oed}^{ref} and EurrefE_{ur}^{ref} ). However, the HSS model is able to capture more realistic behavior of pile-supported wharf under earthquake, because it considers the reduction of soil stiffness with increasing shear strain amplitudes, and introduces hysteric material damping during seismic loading.^[73]

Table 5 shows the displacement time histories at the top and at the toe of the cut-slope rock dike.

The two-dimensional finite element method was used to study the sensitivity of the seismic response of pile-supported wharf to the variation of the soil constitutive models. Indeed, three constitutive models with different accuracies were adopted for this analysis, namely the MC model, the HS and HSS advanced models.

The JCB01 model selected from a series of centrifuge models carried out at the Center for Geotechnical Modeling of California University (UC Davis) was chosen for the study.^{[3,10,26,28,29]}

This model has a configuration of a typical wharf of Western United States ports. It consists of a deck supported by seven rows of vertical piles including two pairs of batter piles, embedded in cut slope rock dike overlaid an embankment of loose sand. The model was subjected to the 1989 Loma Prieta earthquake.

The wharf was modeled using the plane strain model. The concept of pseudo cohesion in rockfill was introduced. The analysis was carried out in time-domain and modification of the initial Newmark scheme defined with α-method (or HHT method)^[61] was used. The constitutive parameters of soil were calibrated through numerical soil element tests.

In this work, the performance of numerical modeling in the estimation of wharf structural displacement and pore pressure were compared to the numerical results of Mohr Coulomb model in FLAC2D, carried out by Heidary-Torkamani et al.^[28] The verification shows that the model based on higher stiffness under small-strain (HSS model) comes closest to the results of Heidary-Torkamani et al.^[28] contrary to the MC and HS models.

This is due to the fact that the MC model considers a constant Young’s modulus for soil and uses only Rayleigh damping, contrary to the advanced models (HS and HSS), which uses three different types of stiffness modulus ( E50refE_{50}^{ref} , EoedrefE_{oed}^{ref} and EurrefE_{ur}^{ref} ), in order to accurately predict the soil behavior. In addition, the HSS model takes into consideration the rigidity in small ground deformations. It includes two parameters G0refG_0^{ref} and γ_0,7, and incorporates the hysteretic damping of materials in addition to Rayleigh damping.

The use of approximations from empirical equations developed by Hardin and Black^[78] in order to determinate G₀ parameter for HSS model, gives appropriate result.

The displacement of piles in sloping rockfill and the displacement time histories of the rock dike at the top and at the toe were also highlighted. It can be observed that the earthquake induced significant lateral movement at the top of the cut-slope rock dike, due to deformations at the upper part of the embankment during the liquefaction of loose sand. This lateral ground displacement affects the behavior of the wharf and the displacement of piles, specially the piles closest to the top of the rock dike.

At the base of the piles, the total displacements are small, and then gradually increase until reaching maximum values at the top of the first three piles (seaward) as well as for the batter piles. But for the other piles, the maximum values of the total displacement are reached at the upper part of the piles when deformation of loose sand is important. It can be also seen that the batter piles have performed adequately.

According to the analysis results, the MC and HS models overestimate the maximum total displacement in piles by 61.51% and 52.2%, respectively, whereas for the HSS model, the maximum total displacement is underestimated by 10.52%, relative to the result of the numerical analysis established with the finite difference method.^[28]

The displacement time history at the dike top is greater than its displacement at the toe due to lateral movement at the embankment top. A difference of about 67%, 40.21% and 33.33% is recorded for the maximum lateral displacement for MC, HS and HSS models, respectively, and a difference of about 45.45%, 40.74% and 38% is estimated regarding its permanent lateral displacement for MC, HS and HSS models, respectively. Concerning predicted vertical displacement, the top of the rock dike has a tendency to settle contrary to its toe, which heaves.

The influence of the soil constitutive models used in this study is demonstrated on the results of piles and cut-slope rock dike displacements. It is concluded that the MC model is not appropriate for modeling seismic behavior of pile-supported wharf structures, and it is recommended to use an advanced model with nonlinear behavior under small-strain and hysteretic damping.

The response of the piles embedded into the rock dike is dictated by the resulting soil deformation mechanisms. The choice of the appropriate soil constitutive model has significant effects on the numerical results, and depends on many factors related to the type of analysis to perform, expected precision of predictions and available soil parameters.

The use of 2D finite element analysis with plane-strain and HSS model is an efficient tool for predicting the seismic behavior of pile-supported wharf. However, some suggestions and perspectives have emerged from this work. First, it would be interesting to use a model better suited to liquefaction analysis and allowing to determine the generation of excessive pore pressure during an earthquake (such as the UBCsand model), and to compare its results with those obtained in this paper.

Secondly, it can be noticed that the use of complex soil constitutive model (as advanced models) is computationally expensive, especially in dynamic calculation. And in order to remedy this, it is advisable to use adaptive constitutive modeling. That can considerably reduce the computational effort. This approach consists of using different soil constitutive models with different levels of complexity, but from the same hierarchical soil model family. Then, more sophisticated models are assigned to the regions that are directly influenced by the geotechnical problem. On the other hand, simpler models are assigned to other regions with low concentration of stress and strain. Lastly, it will be interesting to evaluate and to compare the results of this work with advanced 3D finite element analysis, focusing on the differences between the two methods.