Effect of Soil Anisotropy on Ground Motion Characteristics

1. Introduction

To improve the seismic design of structures, it is necessary to conduct site-specific seismic motion analysis to understand the characteristics of ground motion at the site. Earthquakes have resulted in significant economic losses and casualties. Earthquakes have the potential to inflict significant destruction upon various types of structures, encompassing inland buildings [1], hydraulic structures [2], and subterranean constructions [3]. Numerous scholars have extensively researched the seismic resilience of architectural structures [4,5,6,7]. The analysis of structures for seismic resistance utilizes a range of methodologies, encompassing response spectrum analysis, the base shear method, and time history analysis. The response spectrum method is extensively utilized as the predominant seismic analysis approach in engineering practice. The response spectrum is obtained by calibrating the ground motion characteristics at the specific site. Consequently, numerous scholars have researched the seismic design of structures based on the seismic motion characteristics at the site [7,8,9].

Soil dynamic constitutive models can be categorized into total stress models and effective stress models, considering the aspects of stress transmission and inter-particle contact. The effective stress model offers superior capability in addressing the issue of seismic liquefaction [10,11,12]. Regarding the issue of seismic ground motion response, researchers commonly employ the total stress model for analysis. The total stress models include elastic-plastic model [13,14,15,16,17,18,19], nonlinear model [20,21,22] and equivalent linear model [23]. Traditional elastoplastic models such as the Mohr-Coulomb model [24] and the Drucker-Prager model [25], can be integrated with boundary interface theory [15,16,17], kinematic hardening theory[13,14], and nested yield surface theory [18,19] to account for the effects of cyclic loading. Although elastoplastic model theory effectively describes the hysteresis characteristics and nonlinearity of soils, its computational complexity and significant workload make it inconvenient for engineering applications. Researchers employ nonlinear models [21,22] based on the Masing criteria to depict the nonlinearity and hysteresis characteristics of soils. While these nonlinear models simplify the workload associated with elastoplastic models, they still entail a certain level of complexity. The linear model offers the advantages of low workload and low complexity, allowing for the consideration of soil nonlinearity and hysteresis characteristics through the application of an equivalent concept. Hence, the most widely utilized approach at present is the equivalent linear model, which has evolved from linear models. Schnabel [18] initially pioneered the development of the frequency-domain equivalent linear model within the SHAKE program. Idriss et al. [26] modified the expression of damping in the SHAKE program and developed the SHAKE91 program. With the development of commercial software, the SHAKE and SHAKE91 programs have been integrated into software such as EERA and DEEPSOIL. The SAHKE program demonstrates limited proficiency in effectively addressing anisotropic challenges.

Soil anisotropy has aroused great interest in recent years. Zhang et al. [27,28] and Teng [29] investigated the impact of soil anisotropy on excavation-induced effects in excavations. Wei et al. [30], Soe et al. [31], and Zhang et al. [32] investigated the impact of soil anisotropy on tunnel design and construction. Peric et al. [33], and Ai et al. [34,35] investigated the impact of anisotropy on the design of pile foundations. It has been found that there is anisotropy in the soil under small strain conditions. Bentil [36] has delved into the anisotropy of the small-strain shear modulus by conducting bending element experiments. The consideration of soil anisotropy in seismic response analysis of soil layers has been scarce among scholars, primarily due to the complexities involved in studying anisotropy, seismic loads, and their associated intricacies. Considering the lack of relevant researches, and to maintain the coherence and rigor of scientific inquiry, this study investigates the impact of anisotropy on the seismic response of the site using a finite element method with an anisotropic time-domain equivalent linear model.

2. Soil dynamic characteristics

2.1. Linear viscoelastic model

The viscoelastic Kelvin model (a spring connected in parallel with a sticky pot) is used to reflect the hysteresis of the soil under cyclic loading. The stress-strain relationship is described by Equation (1):

$$\tau =G\gamma +{\eta }_G\dot{\gamma }$$

(1)

where $G$ is the shear modulus; $\tau$ is the shear stress; $\gamma$ is the shear strain; ${\eta }_G$ is the shear viscosity coefficient, as given in Equation (2):

$${\eta }_G=\frac{2GD}{\omega }$$

(2)

where $D$ is the damping ratio; $\omega$ is the circular frequency.

2.2. Modulus and damping models

The key to the effective linearization method is to determine the relationship between the shear modulus ratio and the damping ratio with the shear strain. Hyperbolic models are widely used to describe nonlinear soil behavior under cyclic loading, such as the Pyke model and the Stokoe model [37]. In this study, the improved Stokoe mode is selected to fit the relationship between the shear modulus ratio and shear strain, as given in Equation (3):

$$G_{norm}=\frac{G}{G_{\max }}=\frac{1}{1+{\left(r/r_r\right)}^{\alpha }}$$

(3)

where $G_{norm}$ is the normalized shear modulus, $r_r$ is the reference strain; $\alpha$ is the fitting parameter; the definition of $r_r$is different from the Harden-Drnevich model ($r_r={\tau }_{max}/G_{max}$).

Zhang et al’s. [38] formula is adopted to describe the relationship between $D$ and $G_{norm}$ as follow:

$$D=K_1{G}_{norm}^2+K_2{G}_{norm}+K_3$$

(4)

where $K_1\mathrm{、}{K}_2\mathrm{、}{K}_3$ are the model fitting parameters.

2.3. Effects of anisotropy

The notation for an anisotropic material used herein is: the $y$-axis (the vertical direction) represents the direction of the anisotropy, and the $x,z$-plane is the plane of isotropy. The stress-strain increment equation for an anisotropic material can be written as follows (5) [39].

$$\left\{\begin{array}{l}\delta {\epsilon }_x\\ \delta {\epsilon }_y\\ \delta {\epsilon }_z\\ \delta {\gamma }_{xy}\\ \delta {\gamma }_{yz}\\ \delta {\gamma }_{zx}\end{array}\right\}=\left[\begin{array}{cccccc}\frac{1}{E_h}& -\frac{v_{vh}}{E_v}& -\frac{v_{hh}}{E_h}& 0& 0& 0\\ -\frac{v_{hv}}{E_h}& \frac{1}{E_v}& -\frac{v_{hv}}{E_h}& 0& 0& 0\\ -\frac{v_{hh}}{E_h}& -\frac{v_{vh}}{E_v}& \frac{1}{E_h}& 0& 0& 0\\ 0& 0& 0& \frac{1}{G_{hv}}& 0& 0\\ 0& 0& 0& 0& \frac{1}{G_{vh}}& 0\\ 0& 0& 0& 0& 0& \frac{1}{G_{hh}}\end{array}\right]\left\{\begin{array}{l}\delta {\sigma }_x\\ \delta {\sigma }_y\\ \delta {\sigma }_z\\ \delta {\sigma }_{xy}\\ \delta {\sigma }_{yz}\\ \delta {\sigma }_{zx}\end{array}\right\}$$

(5)

where $E_v$ and $E_h$ Young's moduli in the vertical and horizontal directions, respectively; $V_{hh}$ and $V_{vh}$ Poisson's ratios for horizontal strains from a horizontal and vertical strain, respectively; $V_{hv}$ Poisson's ratio for vertical strains from a horizontal strain; $G_{vh}$ and $G_{hv}$ shear moduli in the vertical plane; and $G_{hh}$ shear modulus in the horizontal plane.

The anisotropy ratios ${AR}_G$ and ${AR}_E$ for the shear modulus and Young's modulus are respectively defined as:

$$A{R}_G=\frac{G_{hh}}{G_{vh}}$$

(6)

$$A{R}_E=\frac{E_h}{E_v}$$

(7)

Under the undrained condition, these Poisson's ratios [40] needs to satisfy the additional relationships ($v_{vh}=0.5,v_{hh}+v_{hv}=1$). Therefore, one can further get the following equations [40]:

$$v_{hh}=1-v_{hv}=1-A{R}_E\cdot v_{vh}$$

(8)

$$v_{hv}=\frac{E_h}{E_v}\cdot v_{vh}=A{R}_E\cdot v_{vh}$$

(9)

$$G_{hh}=\frac{A{R}_E\cdot E_v}{2(2-A{R}_E\cdot v_{vh})}$$

(10)

$$G_{vh}=G_{h\mathrm{v}}=\frac{1}{A{R}_G}\frac{A{R}_E\cdot E_v}{2(2-A{R}_E\cdot v_{vh})}$$

(11)

3. Simulation of time-domain equivalent linear model for anisotropic soil layers

3.1. Input Ground Motion

In the Eurocode 8 [41], sites are classified according to the average shear wave velocity in upper 30m thick soil profile. Therefore, the soil layer thickness is selected as 30m. Also, three seismic waves are chosen: Ei-Centro Wave, Shanghai Wave, and Kobe Wave. The amplitudes of these seismic waves were adjusted to have a peak acceleration of 0.1g. Figure 1 illustrates the three seismic waves.

3.2. Isotropic soil layer

This subsection is used to verify the validity of the established finite element method with an anisotropic time-domain equivalent linear model in the isotropic case studies. The shear modulus model chosen here is the Seed-Idriss mode [42], and the damping curve model used is the Idriss mode [26]. Table 1 presents the model parameters. The properties of the soil layers are shown in Table 2.

Figure 2 illustrates the comparison of the peak ground acceleration (PGA) with depth, considering three seismic records, namely Ei-Centro Wave, Shanghai Wave, and Kobe Wave, along with four soil conditions. They exhibit good agreements with the results obtained from the classical seismic analysis code EERA.

3.3. The Influence of Anisotropic Parameters

It should be noted that this study assumes a consistent relationship between $G_{vh}-\gamma$ and $G_{hh}-\gamma$. Different anisotropic parameter (${AR}_E$=1.00, 1.30, 1.60; ${AR}_G$=1.05, 1.30, 1.55, 1.80, 2.05) were selected for a 30-meter-thick layer of soft soil with a shear wave velocity of 150m/s. Detailed parameters are listed in Table 3.

${PGA}_{anios}/{PGA}_{ios}$is the ratio between the peak ground acceleration obtained at an anisotropic site and that obtained at an isotropic site. Figure 3 shows that, with ${AR}_E$ increasing, ${PGA}_{anios}/{PGA}_{ios}$ gradually decreases. Conversely, as ${AR}_G$ increases, PGA has an obvious increase. Especially when inputting Kobe wave, the PGA of anisotropic conditions exceeds that of isotropic conditions up to 14% with ${AR}_G$ increasing. Consequently, considering soil anisotropy is essential in analyzing site conditions for seismic response.

To assess the impact of anisotropy on the seismic response, the dominant frequency of the seismic waves and the fundamental frequency of the site are further discussed as listed in Table 4, where the fundamental frequency of the site is calculated using Equation (12) [43]:

$$f_g=\frac{1}{T_g}=V_{Si}/4{\displaystyle {\sum }_{i=1}^n{H}_i}$$

(12)

where the symbol "$f_g$" represents the fundamental frequency of the site; "$T_g$" refers to the site's characteristic period; the variable "$H_i$" denotes the thickness of individual soil layers; whereas "$V_{Si}$" represents the wave velocity specific to each respective soil layer.

As shown in Figure 3, the Ei-Centro and Kobe waves have significant influences on the ${PGA}_{anios}/{PGA}_{ios}$, while Shanghai wave has a relatively slight effect. From Table 4, compared with the dominant frequency of the Shanghai wave, those of the Kobe wave and Ei-Centro waves are closer to the fundamental frequency of the site. The waves with the dominant frequencies closer to the fundamental frequencies of sites may result in more significant effect of anisotropy. Therefore, detailed investigations on the effects of the soil anisotropy on the seismic response of Shanghai Yangshan Port site are conducted as an example.

4. Seismic response of Shanghai Yangshan Port

4.1. Input seismic wave

The seismic characteristics of Yangshan Port are examined in this study through the analysis of artificial seismic wave data from Shanghai. According to the site classification method specified in the "Code for Seismic Design of Buildings" (GB 50011-2010) [44], the Yangshan Port area belongs to the fourth category of site. Based on the linear elastic soil layer, the ground acceleration time history is inverted to bedrock to obtain the bedrock acceleration time history [45].

4.2. Calculation model and parameters of soil layer

Most of the topsoil in Yangshan port is soft clay, and there is sand soil in the lower layer. Through the geological investigation, the site of Yangshan Port consists of four typical soil layers: clay layer, silty clay layer, muddy silty clay, and sand soil layer. According to Hou [46], the ${AR}_E$ values for this area range from 1.6 to 2.40. Li [47] reported ${AR}_G$ values ranging from 1.08 to 1.39 for Shanghai soil, while Ng [48] stated the values were between 1.07 and 1.38. Considering the earlier discussion on the suppression effect of ${AR}_E$ and the promotion effect of ${AR}_G$ on PGA, an ${AR}_E$ value of 1.6 and an ${AR}_G$ value of 1.4 were chosen to model the unfavorable conditions. Detailed borehole data for the soil layers are provided in Table 5.

Resonance column and cyclic triaxial tests were conducted to obtain $G_{norm}-\gamma$ and $D-\gamma$ curves for clay, silty clay layer and muddy silty clay. The design confining pressures are 150 kPa, 200 kPa, and 250 kPa. The experimental results pertaining to silty clay are specifically analyzed in this study. Figure 4 illustrates the S-wave output signals of the bending element in the silty clay specimen subjected to a confining pressure of 150 kPa across different input frequencies. Furthermore, the $G_0$ value derived from the bending element tests is approximately 1.1 to 1.2 times greater than the results obtained from the resonant column tests, which aligns with the outcomes reported by Yang et al. [49] and Gu et al. [50,51]. The reason behind this is that the results of the bending element tests specifically relate to the localized stiffness of the shear wave propagation path, while the resonant column offer insights into the overall stiffness characteristics of the specimens [49,50,51]. This confirms the accuracy of the resonant column test. Simultaneously, Yang et al. [49], Gu et al. [50,51], and Youn et al. [52] propose to determine $G_0$ based on the resonant column test.

The cyclic triaxial experiments were performed employing strain control. To mitigate the influence of cyclic loading on subsequent tests, the number of cycles for each loading stage was limited to six. Figure 5 depicts the outcomes of the cyclic triaxial tests conducted on the silty clay. The modulus degrades as the strain increases, while it amplifies with the elevated confining pressure. Through fitting analysis, we obtained the $G_{norm}-\gamma$ and $D-\gamma$ curves for clay in the Yangshan Port site. Figure 6 depicts the curve of best fit for the silty clay specimen. For the other soil layers, the model parameters calibrated from the resonance column and the cyclic triaxial tests are also provided in Table 6.

4.3. Results

In this case study, the Shanghai wave is adjusted to have the peak ground acceleration of to 0.1g. Figure 7 and Figure 8 show the changes in peak acceleration along the depth and the ground peak acceleration spectrum, respectively. It is evident that, when the proposed calculation model reduces to the isotropic scenario, it aligns well with the results obtained from EERA. However, when considering the soil anisotropy, the PGA increases by 19.70%; the peak value of the ground response spectrum undergoes a significant increase of 28.82% at 1 Hz.

To further explore the design response spectrum of Yangshan Port, the ground motion responses of the site under different seismic intensities (0.1g, 0.15g, 0.2g, 0.25g, 0.3g) are studied. Based on the 5% damped ground acceleration response spectrum of different intensities, the least squares method by Andreotti [53] et al. was used to calibrate the design response spectrum. The shape function of the design response spectrum is as follows [54]:

$$\beta (T)=\left\{\begin{array}{ll}1+({\beta }_{\max }-1)\frac{T}{T_0}& 0\le T\le T_0\\ {\beta }_{\max }& T_0\le T\le T_g\\ {\beta }_{\max }{\left(\frac{T_g}{T}\right)}^{\chi }& T_g\le T\le T_m\end{array}\right.$$

(13)

where $T_0$ is the period associated with the initial inflection point, established at 0.1s. $T_g$ signifies the characteristic period. Since this study focuses on a certain soil layer in Yangshan Port, the characteristic period ($T_g$) can be calculated according to Equation (12), which is 0.92s. $T_m$ is the cutoff period, $\chi$ represents the attenuation index, and ${\beta }_{max}$ denotes the plateau value.

The gray lines in Figure 9 are the ground response spectrums at different intensities with a consideration of soil anisotropy. The upper and lower bounds of the gray area are the 16th to 84th percentil, respectively, showing the region with a probability range of 16% to 84% for the occurrence of seismic events. Red line represents the average of the response spectrum curve. Using the shape function given in Equation (13), the 5% damped design response spectrum for the Yangshan Port is calibrated as the blue line.

In addition, based on the same process, the design response spectrum for the isotropic site is calibrated as the green line, of which plateau value ($\beta$) is 2.55, the characteristic period ($T_g$) is 0.92s, and the attenuation index ($\chi$) is 0.69.

Compared with the isotropic one, the anisotropic one has a higher design response spectrum, with a plateau value ($\beta$) of 2.97, a characteristic period ($T_g$) of 0.92s, and an attenuation index ($\chi$) of 0.50. After accounting for anisotropy, the plateau value ($\beta$) of the standard design response spectrum increased by 16.47%. Therefore, when using the isotropic design response spectrum for the seismic design of structures at Shanghai Yangshan port, it may lead to dangerous results.

5. Conclusions

This study establishes a model for analyzing the dynamic behavior of anisotropic soil layers with equivalent linear viscoelasticity. The effects of different anisotropy values of three typical seismic waves on site ground motions are investigated. It is found that the anisotropy ratios for the shear modulus (${AR}_G$) has a promoting effect on the peak ground acceleration (PGA), while the anisotropic ratio of undrained Young's modulus(${AR}_E$) has an inhibitory effect on peak ground acceleration (PGA) in the site earthquake response problem. The impact of anisotropy on the ground motion of the site becomes more significant when seismic waves have the dominant frequencies closer to the fundamental frequencies of sites.

The study futher takes the Yangshan Port site as an example to calibrated the seismic motion parameters of the Yangshan Port soil layer based on cyclic triaxial tests and resonant column tests, and studied the seismic motion characteristics of the anisotropic site of Yangshan Port.

It is found that the consideration of site anisotropy leads to significant increases in both the peak ground acceleration (PGA) and the peak of ground acceleration spectrum for the Yangshan Port site. The design response spectrum calibrated for an anisotropic site is higher than the design response spectrum calibrated for an isotropic site. When conducting seismic design of buildings, the design response spectrum obtained by using isotropic site ground motion response calibration may be more dangerous.