One dimensional time-domain non-linear site seismic response analysis program integrating two hysteresis models of soil

In the numerical simulation of site seismic responses, traditional equivalent linearization methods typically realized in the frequency domain cannot satisfactorily analyze the high-degree non-linearity of soil under strong input motions. Therefore, the “true” non-linear numerical methods performed in the time domain are often utilized in such cases. However, a crucial element of the time-domain non-linear method, which is the hysteresis model of soil that describes the rule controlling the loading–unloading behavior of soil, has no established guidelines for earthquake engineering. Different researchers presented different models, revealing the epistemic uncertainty related to the dynamic properties of soil. Thus, the time-domain non-linear method should consider this uncertainty in practice. Therefore, in this study, a one-dimensional (1D) time-domain non-linear site seismic response analysis program was developed. The developed program was coded using Fortran95 and integrates two kinds of soil hysteresis models (i.e., extended Masing model and dynamic skeleton curve model). In both models, the damping correction was introduced to calibrate the hysteresis loop area toward the damping ratio measured in the dynamic triaxial test or resonant column test. Moreover, the temporospatial finite difference algorithm was used to resolve the 1D non-linear wave equation, and its precision was demonstrated in comparison with the results of the frequency-domain program for the linear case. Finally, the non-linear seismic response of a specific site was calculated by the proposed program. The findings of the fitting were compared to those of the two popular time-domain non-linear programs DEEPSOIL (Hashash, V6.1) and CHARSOIL (Streeter et al., CHARSOIL, Characteristics Method Applied to Soils, 1974 March 25). Simultaneously, the Japanese KIK-net strong motion observation station data were applied to validate the reliability of this program.

1 Introduction

Many data on earthquake investigation and observation show that (Seal, 1988; Gao et al., 1996; Beresnev, 2002; Kokusho et al., 2008; Kim et al., 2013) local site conditions will directly affect the ground motion amplitude, spectral characteristics, and the distribution of earthquake disaster degree. The amplification effect of ground-to-ground motion must be considered in seismic zoning, seismic safety evaluation, and seismic design analysis of engineering sites. When the ground motion parameters of engineering site design are determined, it has to consider the influence of site conditions on ground motion amplification by using corresponding adjustment methods (Li, 2013). The adjustment relationship of site ground vibration parameter should be established based on many obtained strong vibration observation data and numerical calculation results of site model. In addition, the statistical relationship of site characteristic data and ground motion property spectrum has to be analysis to determine the adjustment model of site ground vibration parameters. This requires collecting sufficient basic data for research in addition to ensuring reasonable and reliable results. Site model with numerical calculation by means of the study of field adjustment should involve the soil seismic response calculation, the construction of the theoretical model, and the reasonable descriptions on the non-linear change of stress and strain of soil constitutive relationship. In this case, it can support the large-scale site numerical model for computing platform and improve the field response calculation efficiency.

At present, the frequency-domain equivalent linearization method and time-domain non-linear (TNL) method are the most widely applied in soil seismic response analysis. The influence of site soil conditions is analyzed by using a frequency-independent equivalent linearization method in China (Bo, 1998; Li et al., 2005; Li, 2010), which is widely used in engineering due to simple concept and small calculation. However, the equivalent linearization method is often no longer applicable in weak sites with large strain and strong non-linearity (Qi et al., 2010). To ensure the calculated results being consistent with the observed results, the TNL stepwise integration rule is proposed, which is a seismic response analysis method with clear physical significance and can more truly reflect the non-linear physical process of soil under stress state (Li, 1992). The true non-linear analysis method can describe the dynamic stress-strain non-linear hysteresis model of soil (Luan et al., 1994).

There are three aspects for the one-dimensional (1D) time-domain calculation method of the site non-linear seismic response, including the model of the soil skeleton curve, the soil loading-unloading-reloading criterion, and the time-domain integration method. The soil loading-unloading-reloading criterion was firstly established by Masing in 1926 (Masing, 1926). Later, it is revised continually and the empirical function form of the skeleton curve is proposed, which is a constitutive model that can better fit the test results (Pyke, 1979; Matasovic et al., 1993; Wang et al., 1981; Li, 1992; Luan, 1992; Chen, 2006; Qi et al., 2010; Martin et al., 1982; Yee et al., 2013). At present, the TNL constitutive models are widely applied in China, including Wang Zhiliang’s extended Masing criterion (Wang et al., 1981), Li Xiaojun’s dynamic skeleton curve (Li, 1992), Luan Maotian’s true non-linear model (Luan et al., 1992), and Qi Wenhao’s UE model (Qi et al., 2010). With the development of computer technology and engineering fluctuation, various site seismic reaction analysis programs have been developed based on soil constitutive model and time-domain numerical integration methods, which are widely applied in the international geotechnical seismic engineering. At present, the internationally well-known CHARSOIL program proposed by Streeter in 1974 (Streeter et al., 1974) is based on the finite difference method and uses the Ramberg-Osgood constitutive model to reflect the dynamic non-linear characteristics of soil. It is the first TNL soil seismic response analysis program in the world. However, it uses rigid boundary conditions, so sometimes the calculation results are not consistent with the actual situation. DEEPSOIL (Hashash, 2009) is a non-linear seismic response analysis program of 1D soil layer widely used abroad. It adopts a variety of numerical analysis methods (frequency-domain equivalent linearization methods and TNL methods), and introduces several constitutive models to describe the non-linear changes of soil in the TNL method. It can be applied for different soil layer earthquake responses. In addition to that, similar programs include DesRA-2C (Idriss et al., 1968), MASH (Martin et al., 1978), D-MOD (Kavazanjian and Matasovic, 1995), and ONDA (Diego et al., 2006).

The current dilemma of the TNL methods is analyzed as follows. On the one hand, the key elements of the current TNL method, namely the soil constitutive model describing the behavior of soil addition and unloading hysteresis, is not a universally accepted rule in seismic engineering. Therefore, most research focuses on developing a more scientific and reasonable method which can actually describe the non-linear dynamic changes of soil constitutive relationship. Based on previous studies and combination of the experiment data, a simplified constitutive model for complex mathematical correction with a variety of test soil parameters can be more reasonable to describe the dynamic skeleton curve model of the relationship. However, it ignores the analysis of cognitive uncertainty related to soil dynamic properties revealed by different constitutive models. In practical engineering applications, the TNL method considers this uncertainty and the good combination of developing constitutive models and engineering applications. On the other hand, more mature international geotechnical TNL method generally exist interface interaction with poor intuition and efficiency and more complex modelling process. In addition, it gives an upper limit on the calculation scale of site model and seismic input. The research trend of the field adjustment model based on many statistical results is to adopt a variety of computational methods and constitutive models for large-scale site seismic response, so as to eliminate the uncertain impact caused by the single frequency-domain equivalent linearization method and the small calculation site type. Therefore, it is imperative to develop efficient and diversified geotechnical non-linear calculation methods with independent property rights.

This work develops a new non-linear seismic response analysis program based on Fortran95, which integrates two hysteresis models of soil, namely extended Masing criterion and dynamic skeleton curve. Meanwhile, it introduces the damping correction coefficient to calibrate the lag ring area to match the damping ratio measured in dynamic triaxial test or resonance column test. It can store data by using dynamic array technology. In principle, there is no limit on the number of site models and the calculation scale of time and space points of input and output data. The results of this work can be used for calculation of large-scale site seismic response and uncertainty analysis of site calculation method, based on which the influence characteristics of site conditions on ground motion can be studied.

2 Soil skeleton curve model

The hyperbolic skeleton model is adopted to describe the soil skeleton curve, namely, a function of strain-stress. The soil skeleton curve model is described below by taking the shear deformation as an example. The vertices of stress-strain hysteresis cycles corresponding to different strain amplitudes are connected to each other to form the soil skeleton curve, as shown in Figure 1. If any point on the skeleton curve is connected to the origin, the slope of the resulting line is the cut modulus G corresponding to the point. In addition, due to the non-linear nature of the soil, the modulus is a function of the strain, namely $G=G\left(\gamma \right)$. The stress-strain relationship can be expressed as follows:

$\mathit{\tau }\left(\mathit{\gamma }\right)=\mathit{G}\left(\mathit{\gamma }\right)\mathit{\gamma }(1)$

FIGURE 1

FIGURE 1. Soil skeleton curve.

If a hyperbolic function is employed to describe a 1D soil constitutive model, Eq. 1 above can be expressed as follows:

$\mathit{\tau }\left(\mathit{\gamma }\right)={{\mathit{G}}_{}}_{\mathbf{max}}\frac{{\mathit{\gamma }}_{\mathbf{r}}\mathit{\gamma }}{\left|\mathit{\gamma }\right|+{\mathit{\gamma }}_{\mathbf{r}}}(2)$

In the equation above, ${G_{}}_{\max }$ is the maximum shear modulus, and the tangent shear modulus corresponding to the zero strain was shown in Figure 1. Then, if the shear wave velocity is c and the mass density is ρ, then below relation can be obtained:

${{\mathit{G}}_{}}_{\mathbf{max}}=\mathit{\rho }{\mathit{c}}^{\mathbf{2}}(3)$

In Eq. 2, ${\gamma }_{\mathrm{r}}$ is the reference strain that can be obtained from the test data regression. The Eq. 4 can be acquired based on Eqs 1, 2:

$\mathit{G}\left(\mathit{\gamma }\right)=\frac{{\mathit{\gamma }}_{\mathbf{r}}}{\left|\mathit{\gamma }\right|+{\mathit{\gamma }}_{\mathbf{r}}}{{\mathit{G}}_{}}_{\mathbf{max}}(4)$

Thus:

$\frac{{{\mathit{G}}_{}}_{\mathbf{max}}}{\mathit{G}\left(\mathit{\gamma }\right)}=\mathbf{1}+\frac{\left|\mathit{\gamma }\right|}{{\mathit{\gamma }}_{\mathbf{r}}}(5)$

Typically, dynamic triaxial or resonant column tests can provide discrete data between the shear modulus ratio $G\left(\gamma \right)/{G_{}}_{\max }$ and the shear strain γ. Based on the model given in Eq. 5, the reference strain ${\gamma }_{\mathrm{r}}$ can be regressed by the least square method. For example, the discrete data points in Figure 2 can be regressed and substituted into Eq. 5, based on which the shear modulus ratio curve can be obtained (the solid line in the figure). The discrete points in Figure 2 are not the original test data, but the modulus ratios under different strains are obtained by fitting the test data using Eq. 5. Therefore, the curve by fitting the hyperbolic model shown in Eq. 5 naturally coincides with the discrete data points. Those processed by the actual calculation are the discrete data points shown in Figure 2.

FIGURE 2

FIGURE 2. The shear modulus ratio curve.

In this work, the hyperbolic model can describe the soil skeleton curve.

3 Soil loading and unloading criterion

After the soil skeleton curve is defined, the corresponding loading and unloading criteria should be set to fully describe the loading-unloading-reloading dynamic process of the soil under the complex cyclic load. This work provides two addition and unloading criteria by previous study: extended Masing criteria proposed by Wang Zhiliang (Wang et al., 1981) and dynamic skeleton curve proposed by Li Xiaojun (Li, 1992).

3.1 Extended Masing criterion

The extended Masing criterion is widely used in seismic engineering. It is developed on the basis of Masing Criterion (Masing, 1926) and proposed by Wang Zhiliang. Its basic rules are as follows:

1) During the initial loading, the stress-strain relationship follows the skeleton curve $\tau \left(\gamma \right)=f\left(\gamma \right)$.

2) During the unloading and reverse loading, the tangent modulus at the initial unloading is equal to the maximum shear modulus ${G_{}}_{\max }$ of soil, and the stress-strain relation curve shows the “double” relationship with the original skeleton curve:

$\frac{\mathit{\tau }-{\mathit{\tau }}_{\mathit{c}}}{\mathbf{2}}=\mathit{f}\left(\frac{\mathit{\gamma }-{\mathit{\gamma }}_{\mathit{c}}}{\mathbf{2}}\right)(6)$

If the hyperbolic model shown in Eq. 2 is used, Eq. 6 can be expressed as follows:

$\mathit{\tau }={\mathit{\tau }}_{\mathit{c}}+\mathbf{2}{{\mathit{G}}_{}}_{\mathbf{max}}\frac{{\mathit{\gamma }}_{\mathbf{r}}\left(\mathit{\gamma }-{\mathit{\gamma }}_{\mathit{c}}\right)}{\left|\mathit{\gamma }-{\mathit{\gamma }}_{\mathit{c}}\right|+\mathbf{2}{\mathit{\gamma }}_{\mathbf{r}}}(7)$

In the above equations, ${\tau }_c$ and ${\gamma }_c$ are the stress and strain values corresponding to the nearest inflection point of the stress-strain curve, respectively.

3) If the stress-strain curve of unloading and reverse loading intersects the skeleton curve, the subsequent stress-strain relation curve follows the skeleton curve, satisfying the “upper skeleton curve” rule;

4) If the unloaded and reverse loaded stress-strain curve intersects the previous load-unloading curve, the subsequent stress-strain curve follows the previous stress-strain curve, which is the upper circle rule.

According to the above extended Masing criterion, the 1D non-linear dynamic stress-strain relation of soil can be given for any dynamic loading process.

Figure 3 shows the hysteretic process of soil loading, unloading, and reloading under the extended Masing criterion. During the initial loading, the stress-strain relationship follows the skeleton curve, namely the 0→2 curve. If unloading is performed at point 2 and loading is realized in the direction after unloading to zero stress point 2′, the unloading-reverse loading curve follows curve 2→2′→1 in Figure 3. The unloading-reverse loading curve intersects the skeleton curve at point 1, which is symmetric about the origin. When the stress-strain relationship reaches point 1 along the curve 2→2′→1, the subsequent reverse loading curve follows the skeleton curve 1→1′. If the stress-strain relationship is reverse unloaded and reloaded along the curve 2→2′→1, the stress-strain curve follows the line 1→2. The expression of the curve 2→2′→1 can be determined by Eq. 7 according to the above rule (2). In the equation, ${\tau }_c$ and ${\gamma }_c$ are corresponding stress and strain at point 2, respectively. The slope of tangent line of curve 2→2′→1 at point 2 is equal to that of the skeleton curve at point 0, which is the maximum shear modulus G_max.

FIGURE 3

FIGURE 3. The extended Masing criterion (without considering damping correction).

In the above unloading-reverse loading process 2→2′→1, if reverse loading is carried out at point 3, the reverse unloading-reverse loading curve still can be determined according to Eq. 7. It illustrates that the curve determined in such way necessarily intersects the skeleton curve at point 2 and the unload-reverse loading curve 2→2′→1. Therefore, this reverse unloading-reloading curve will follow the 3→4→2 curve. If this reverse unloading-reloading process reaches point 2, then the skeleton curve is followed if loading continues, and the above unloading-reloading curve 2→2′→1 is followed if unloading.

If the reverse unloading and reloading process 3→4→2 is unloaded and reverse loaded after reaching point 4, the unloading and reverse loading curve 4→5 can be obtained by substituting the stress and strain corresponding to point 4 into Eq. 7. In this case, the extension line of the curve should meet the corresponding “big circle,” that is, the unloading and reverse loading curve 2→2′→1 of the upper layer reaches point 3. If the unloading and reverse loading process 4→5→3 reaches point 3, the subsequent reverse loading curve follows the curve 3→1 with the principle of “upper big circle” if the revserse loading continues. If this unload-reverse loading process 4→5→3 reaches point 3, the reverse unloading and reverse loading process 3→4→2 above is repeated.

The figures above illustrate the rules that should be followed in the stress-strain relationship of soil under the extended Masing criterion.

The soil-like non-linear dynamic parameter test demonstrates the relationship between shear modulus ratio and damping ratio-shear strain. An analytical representation of the soil skeleton curve can be obtained based on the relationship between the shear modulus ratio and shear strain shown in Figure 3. In combination with the above expansion Masing criteria, the soil retardation curve under equal amplitude cyclic load is obtained, including the hysteresis loops 2→2′→3→1→2 and 4→5→3→4. The area of the resulting lag loop is different from that of the damping ratio corresponding to the maximum strain (i.e., half of the difference between the maximum strain and the arrested loop). To correct this deviation, the damping correction coefficient $K\left({\gamma }_0\right)$ should be included in the stress-strain skeleton curve.

After that, the stress-strain relationship of soil is obtained based on the extended Masing criterion, which can be expressed as Eq. 8:

$\mathit{\tau }\left(\mathit{\gamma }\right)=\left\{\begin{array}{l}\begin{array}{cc}{\mathit{\tau }}_{\mathbf{c},\mathit{i}}+\mathit{K}\left(\left|\frac{{\mathit{\gamma }}_{\mathbf{0}}}{\mathbf{2}}\right|\right)\left[\frac{{\mathit{G}}_{\mathbf{0}}\left(\mathit{\gamma }-{\mathit{\gamma }}_{\mathbf{c},\mathit{i}}\right)}{\mathbf{1}+\left|\frac{\left(\mathit{\gamma }-{\mathit{\gamma }}_{\mathbf{c},\mathit{i}}\right)}{\mathbf{2}{\mathit{\gamma }}_{\mathbf{r}}^{\prime }}\right|}-\frac{{\mathit{\tau }}_{\mathbf{0}}}{{\mathit{\gamma }}_{\mathbf{0}}}\left(\mathit{\gamma }-{\mathit{\gamma }}_{\mathbf{C}}\right)\right]+\frac{{\mathit{\tau }}_{\mathbf{0}}}{{\mathit{\gamma }}_{\mathbf{0}}}\left(\mathit{\gamma }-{\mathit{\gamma }}_{\mathbf{c},\mathit{i}}\right),& \mathrm{N}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{k}\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{t}\mathrm{o}\mathrm{n}\mathrm{c}\mathrm{u}\mathrm{r}\mathrm{v}\mathrm{e}\\ \frac{{\mathit{G}}_{\mathbf{0}}\mathit{\gamma }}{\mathbf{1}+\frac{\left|\mathit{\gamma }\right|}{{\mathit{\gamma }}_{\mathbf{r}}}},& \mathrm{S}keletoncurve\end{array}\\ \end{array}\right.(8)$

If the stress-strain state point is under the skeleton curve, the above equation is applicable for calculation; otherwise, the following equations are applicable:

${\mathit{\gamma }}_{\mathbf{0}}=\left|{\mathit{\gamma }}_{\mathbf{c},\mathit{i}-\mathbf{1}}-{\mathit{\gamma }}_{\mathbf{c},\mathit{i}}\right|(9)$

${\mathit{\tau }}_{\mathbf{0}}=\left|{\mathit{\tau }}_{\mathbf{c},\mathit{i}-\mathbf{1}}-{\mathit{\tau }}_{\mathbf{c},\mathit{i}}\right|(10)$

${\mathit{\gamma }}_{\mathbf{r}}^{\prime }={\left(\frac{\mathbf{2}{\mathit{G}}_{\mathbf{0}}}{\left|{\mathit{\tau }}_{\mathbf{c},\mathit{i}-\mathbf{1}}-{\mathit{\tau }}_{\mathbf{c},\mathit{i}}\right|}-\frac{\mathbf{2}}{\left|{\mathit{\gamma }}_{\mathbf{c},\mathit{i}-\mathbf{1}}-{\mathit{\gamma }}_{\mathbf{c},\mathit{i}}\right|}\right)}^{-\mathbf{1}}(11)$

$\mathit{K}\left({\mathit{\gamma }}_{\mathbf{0}}\right)=\frac{\mathbf{\pi }{\mathit{\gamma }}_{\mathbf{0}}^{\mathbf{2}}\mathit{\lambda }\left({\mathit{\gamma }}_{\mathbf{0}}\right)}{\mathbf{2}{\mathit{\gamma }}_{\mathbf{0}}\left(\mathbf{2}{\mathit{\gamma }}_{\mathbf{r}}^{\prime }+{\mathit{\gamma }}_{\mathbf{0}}\right)-\mathbf{4}{\mathit{\gamma }}_{\mathbf{r}}^{\prime }\left({\mathit{\gamma }}_{\mathbf{r}}^{\prime }+{\mathit{\gamma }}_{\mathbf{0}}\right)\mathbf{ln}\left(\mathbf{1}+{\mathit{\gamma }}_{\mathbf{0}}/{\mathit{\gamma }}_{\mathbf{r}}^{\prime }\right)}(12)$

In the above equations, ${\tau }_{\mathrm{c},i}$ and ${\gamma }_{\mathrm{c},i}$ are the shear stress and shear strain for the current moment before the last inflection point, respectively; ${\tau }_{\mathrm{c},i-1}$ and ${\gamma }_{\mathrm{c},i-1}$ are the turning points before $\left({\gamma }_{\mathrm{c},i},{\tau }_{\mathrm{c},i}\right)$, $G_0$ is the elastic shear modulus, namely the skeleton curve at the origin of tangent slope ($G_{\max }$), and $\lambda \left({\gamma }_0\right)$ is the damping ratio corresponding to the strain amplitude ${\gamma }_0$, which can be calculated by fitting the function relation to the discrete data or directly interpolated to the discrete data.

Equation 12 defines the correction coefficient, which is introduced for better consideration of the damping fit. The soil stress-strain change corresponding to Figure 3 is given in Figure 4, which considers the damping correction according to the hysteresis curve. It suggests that after the damping correction, the “double” relationship between the unloading-reverse loading curve and the skeleton curve is no longer satisfied, and the area of the damping correction is the same as the damping circle corresponding to the test damping ratio. Moreover, for the correction case shown in the figure, the two branches of the hysteresis curve intersect in the skeleton curve. In this case, the unloading-reverse loading curve determined in Eq. 8 intersects the skeleton curve in advance without adopting the “upper skeleton curve” criterion. In Figure 4, the unloading-reverse loading curve 2→2′→3→1 intersects the skeleton curve even at the point 3′, but the upper skeleton curve criterion is followed at point 1. Therefore, both the “upper skeleton curve” criterion and the “upper large circle” criterion are based on whether the current stress-strain curve reaches and exceeds the previous inflection point. For the first loading, the first inflection point appears (point 2 in Figure 4), then the previous inflection point is symmetrical with the origin (point 1). Therefore, the unloading-reverse loading process from point 2 follows the curve 2→2′→3→1.

FIGURE 4

FIGURE 4. Extended Masing criterion (damping correction not considered).

3.2 Dynamic skeleton curve criteria

When the soil stress-strain arrest behavior is described by using the extended Masing criteria, it should repeatedly remember and identify the current inflection point and its previous inflection point. It will be more complex and detrimental to programs when the cyclic load is complex (such as ground motion input). Due to the complexity and uncertainty of the dynamic characteristics of the soil, the hysteresis criterion based on the dynamic skeleton curve proposed by Li (1992) remembers and identifies the current inflection point and the maximum inflection point of stress-strain history while showing the basic characteristics of the Masing criterion, so it is easier to implement.

The dynamic skeleton curve describes the loading-unloading hysteresis criterion for 1D soil based on the following three basic hypotheses:

1) During the initial loading, the stress-strain relationship curve coincides with the soil skeleton curve;

2) During the unloading and reverse loading, the stress-strain relationship curve directly points to the absolute maximum stress-strain point of the stress-strain history or its reverse symmetry point, enabling that the stress and strain relationship curve and the skeleton curve meet the “double” relationship for the equal amplitude cyclic load process;

3) The stress and strain relationship curve of the subsequent loading process after meeting with the skeleton curve will follow the skeleton curve.

The above hypotheses are not the exactly same as those of the extended Masing criterion, where it assumes that hypothesis (3) is the “upper skeleton curve” criterion. The biggest difference between the two is the above hypothesis (2). In the extended Masing criterion, the unloading-reverse loading stress-strain relationship curve points to the previous inflection point of the current inflection point. While in the dynamic skeleton curve criterion, it points directly to the absolute maximum stress-strain point of the stress-strain history or its reverse symmetry point.

According to the above hypotheses, the relationship between soil shear stress and strain corresponding to the hyperbolic dynamic skeleton curve criterion can be obtained as follows:

$\mathit{\tau }\left(\mathit{\gamma }\right)=\left\{\begin{array}{cc}{\mathit{\tau }}_{\mathbf{C}}+\mathit{K}\left(\left|\frac{{\mathit{\gamma }}_{\mathbf{0}}}{\mathbf{2}}\right|\right)\left[\frac{{\mathit{G}}_{\mathbf{0}}\left(\mathit{\gamma }-{\mathit{\gamma }}_{\mathbf{C}}\right)}{\mathbf{1}+\left|\frac{\left(\mathit{\gamma }-{\mathit{\gamma }}_{\mathbf{C}}\right)}{\mathbf{2}{\mathit{\gamma }}_{\mathbf{r}}^{\prime }}\right|}-\frac{{\mathit{\tau }}_{\mathbf{0}}}{{\mathit{\gamma }}_{\mathbf{0}}}\left(\mathit{\gamma }-{\mathit{\gamma }}_{\mathbf{C}}\right)\right]+\frac{{\mathit{\tau }}_{\mathbf{0}}}{{\mathit{\gamma }}_{\mathbf{0}}}\left(\mathit{\gamma }-{\mathit{\gamma }}_{\mathbf{C}}\right),& \left|\mathit{\gamma }\right|\le {\mathit{\gamma }}_{\mathbf{M}}\\ \frac{{\mathit{G}}_{\mathbf{0}}\mathit{\gamma }}{\mathbf{1}+\frac{\left|\mathit{\gamma }\right|}{{\mathit{\gamma }}_{\mathbf{r}}}},& \left|\mathit{\gamma }\right|>{\mathit{\gamma }}_{\mathbf{M}}\end{array}\right.(13)$

In the above relationship,

${\mathit{\gamma }}_{\mathbf{0}}=\pm {\mathit{\gamma }}_{\mathbf{M}}-{\mathit{\gamma }}_{\mathbf{C}}(14)$

${\mathit{\tau }}_{\mathbf{0}}=\pm {\mathit{\tau }}_{\mathbf{M}}-{\mathit{\tau }}_{\mathbf{C}}(15)$

${\mathit{\gamma }}_{\mathbf{r}}^{\prime }={\left(\frac{\pm \mathbf{2}{\mathit{G}}_{\mathbf{0}}}{\pm {\mathit{\tau }}_{\mathbf{M}}-{\mathit{\tau }}_{\mathbf{C}}}-\frac{\pm \mathbf{2}}{\pm {\mathit{\gamma }}_{\mathbf{M}}-{\mathit{\gamma }}_{\mathbf{C}}}\right)}^{-\mathbf{1}}(16)$

$\mathit{K}\left({\mathit{\gamma }}_{\mathbf{0}}\right)=\frac{\mathbf{\pi }{\mathit{\gamma }}_{\mathbf{0}}^{\mathbf{2}}\mathit{\lambda }\left({\mathit{\gamma }}_{\mathbf{0}}\right)}{\mathbf{2}{\mathit{\gamma }}_{\mathbf{0}}\left(\mathbf{2}{\mathit{\gamma }}_{\mathbf{r}}^{\prime }+{\mathit{\gamma }}_{\mathbf{0}}\right)-\mathbf{4}{\mathit{\gamma }}_{\mathbf{r}}^{\prime }\left({\mathit{\gamma }}_{\mathbf{r}}^{\prime }+{\mathit{\gamma }}_{\mathbf{0}}\right)\mathbf{ln}\left(\mathbf{1}+{\mathit{\gamma }}_{\mathbf{0}}/{\mathit{\gamma }}_{\mathbf{r}}^{\prime }\right)}(17)$

${\gamma }_{\mathrm{M}}and{\tau }_{\mathrm{M}}$ are the absolute values of the maximum strain and the strain and stress history before the current moment, or the absolute values the maximum inflection point of the strain and stress before the current moment, respectively; ${\gamma }_{\mathrm{C}}$ and ${\tau }_{\mathrm{C}}$ are the strain and stress corresponding to the last inflection point before the current moment, respectively; $G_0$ refers to the elastic shear modulus, namely the slope of the tangent line of the skeleton curve at the origin; $\lambda \left({\gamma }_0\right)$ is the damping ratio related to the strain amplitude; and $K\left({\gamma }_0\right)$ is the correction coefficient introduced by considering the damping fitting. In the above format, the selection criterion of ± sign is as follows: if the strain increment is positive, the + sign is selected; otherwise, the - sign is selected.

Equation 17 is the same form as Eq. 12, but the definition ${\gamma }_{\mathrm{r}}^{\prime }$ is different.

Without the damping correction, $K$ =1.0. Figure 5 shows the soil stress-strain relationship curve based on the dynamic skeleton curve criterion. For the initial loading, point 2 and point 1 are the first inflection point and the corresponding previous one, respectively. At this point, point 1 coincides with the historical maximum stress-strain point, so the unloading-reverse loading curve 2→2′→3→1 from point 2 is exactly the same curve under the extended Masing criterion shown in Figure 3. For subsequent unloading inflection points 4 and 6 in the figure, the previous inflection points are inflection points 3 and 5 (defined here as the reverse unloading points). When the historical maximum stress-strain point is point 2, and the reverse point is point 1. Therefore, following the dynamic skeleton curve criterion, the unloading-reverse loading process starting from unloading points 4 and 6 will follow the curves 4→1 and 6→1 in Figures 2.2 and 2.3, respectively, in which the reverse point directly points to the historical maximum stress-strain point. According to the extended Masing criteria, the curves of the two unloading-reverse loading processes will point to the previous inflection point, namely inflection points 3 and 5 (Figure 3), which is the difference between the two criteria. Figure 6 illustrates the stress-strain relation curve obtained by considering the dynamic skeleton curve criterion after damping correction.

FIGURE 5

FIGURE 5. Dynamic skeleton curve criterion (damping correction not considered).

FIGURE 6

FIGURE 6. Dynamic skeleton curve criterion (damping correction considered).

3.3 Time-space finite difference method

In this work, the finite difference method is employed to discretize the time-domain and the space-domain. The dynamic equilibrium equation of the 1D wave can be written as follows:

$\frac{\partial \mathit{\tau }}{\partial \mathit{z}}=\mathit{\rho }\frac{\partial \mathit{v}}{\partial \mathit{t}}(18)$

The deformation coordination condition is expressed as follows:

$\frac{\partial \mathit{\gamma }}{\partial \mathit{t}}=\frac{\partial \mathit{v}}{\partial \mathit{z}}(19)$

The non-linear constitutive equation between stress and strain is expressed as follows:

$\mathit{\tau }=\mathit{\tau }\left(\mathit{\gamma }\right)(20)$

In the several equations above, $\tau =\tau \left(z,t\right)$ is the shear stress, $v=v\left(z,t\right)$ refers to the particle motion velocity, and $\gamma =\gamma \left(z,t\right)$ is the medium shear strain. Eqs 18–20 constitute a complete non-linear 1D wave equation.

Before the dynamic equation difference is solved differentially, each soil layer should be subdivided according to corresponding stability conditions of the differential format. Because the space-time center difference scheme is adopted to solve the dynamic equation in this work, the thickness of each soil layer after the subdivision layer should satisfy Eq. 21 below:

$∆{\mathit{z}}_{\mathit{l}}\ge {\mathit{c}}_{\mathbf{s},\mathit{l}}∆\mathit{t}(21)$

$∆z_l$ is the thickness of the first soil layer (sequence numbered from free surface to substrate), $c_{\mathrm{s},l}$ is the shear wave speed of the first soil layer, and $∆t$ refers to the interval of calculation time.

To ensure smaller calculation error caused by spatial discretization, the thickness of each calculated soil should satisfy the below relationship as much as possible:

$\frac{∆{\mathit{z}}_{\mathit{l}}}{∆{\mathit{z}}_{\mathit{m}}}=\frac{C_{s,l}}{C_{s,m}}(22)$

After spatial-temporal discretization of Eqs 18–20 by using the central difference method, the following explicit difference stepwise integral scheme can be acquired:

$\left.\begin{array}{c}{\mathit{\rho }}_{\mathit{l}}^{\prime }\frac{{\mathit{v}}_{\mathit{l}}^{\mathit{p}}-{\mathit{v}}_{\mathit{l}}^{\mathit{p}-\mathbf{1}}}{∆\mathit{t}}=\frac{{\mathit{\tau }}_{\mathit{l}+\mathbf{0.5}}^{\mathit{p}-\mathbf{0.5}}-{\mathit{\tau }}_{\mathit{l}-\mathbf{0.5}}^{\mathit{p}-\mathbf{0.5}}}{∆{\mathit{z}}_{\mathit{l}}^{\prime }}\\ \frac{{\mathit{\gamma }}_{\mathit{l}+\mathbf{0.5}}^{\mathit{p}+\mathbf{0.5}}-{\mathit{\gamma }}_{\mathit{l}+\mathbf{0.5}}^{\mathit{p}-\mathbf{0.5}}}{∆\mathit{t}}=\frac{{\mathit{v}}_{\mathit{l}+\mathbf{1}}^{\mathit{p}}-{\mathit{v}}_{\mathit{l}}^{\mathit{p}}}{∆{\mathit{z}}_{\mathit{l}}}\\ {\mathit{\tau }}_{{\mathit{l}}^{\prime }}^{{\mathit{p}}^{\prime }}=\mathit{\tau }\left({\mathit{\gamma }}_{{\mathit{l}}^{\prime }}^{{\mathit{p}}^{\prime }}\right)\end{array}\right\}(23)$

They can be rewritten as follows:

$\left.\begin{array}{l}{\mathit{v}}_{\mathit{l}}^{\mathit{p}}={\mathit{v}}_{\mathit{l}}^{\mathit{p}-\mathbf{1}}+\frac{∆\mathit{t}}{{\mathit{\rho }}_{\mathit{l}}^{\prime }∆{\mathit{z}}_{\mathit{l}}^{\prime }}\left({\mathit{\tau }}_{\mathit{l}+\mathbf{0.5}}^{\mathit{p}-\mathbf{0.5}}-{\mathit{\tau }}_{\mathit{l}-\mathbf{0.5}}^{\mathit{p}-\mathbf{0.5}}\right)\\ {\mathit{\gamma }}_{\mathit{l}+\mathbf{0.5}}^{\mathit{p}+\mathbf{0.5}}={\mathit{\gamma }}_{\mathit{l}+\mathbf{0.5}}^{\mathit{p}-\mathbf{0.5}}+\frac{∆\mathit{t}}{∆{\mathit{z}}_{\mathit{l}}}\left({\mathit{v}}_{\mathit{l}+\mathbf{1}}^{\mathit{p}}-{\mathit{v}}_{\mathit{l}}^{\mathit{p}}\right)\\ {\mathit{\tau }}_{{\mathit{l}}^{\prime }}^{{\mathit{p}}^{\prime }}=\mathit{\tau }\left({\mathit{\gamma }}_{{\mathit{l}}^{\prime }}^{{\mathit{p}}^{\prime }}\right)\end{array}\right\}{\mathit{\tau }}_{{\mathit{l}}^{\prime }}^{{\mathit{p}}^{\prime }}=\mathit{\tau }\left({\mathit{\gamma }}_{{\mathit{l}}^{\prime }}^{{\mathit{p}}^{\prime }}\right)(24)$

In the equations above.

$∆{\mathit{z}}_{\mathit{l}}^{\prime }=\frac{∆{\mathit{z}}_{\mathit{l}}+∆{\mathit{z}}_{\mathit{l}-\mathbf{1}}}{\mathbf{2}},\mathit{l}=\mathbf{1,2},\dots ,\mathit{L}$

${\mathit{\rho }}_{\mathit{l}}^{\prime }=\frac{{\mathit{\rho }}_{\mathit{l}}∆{\mathit{z}}_{\mathit{l}}+{\mathit{\rho }}_{\mathit{l}-\mathbf{1}}∆{\mathit{z}}_{\mathit{l}-\mathbf{1}}}{\mathbf{2}∆{\mathit{z}}_{\mathit{l}}^{\prime }},\mathit{l}=\mathbf{1,2},\dots ,\mathit{L}$

${\mathit{\rho }}_{\mathbf{0}}=\mathbf{0},∆{\mathit{z}}_{\mathbf{0}}=\mathbf{0}$

$v_l^p$ represents the velocity of the medium at the top of the soil layer calculated at the lth time $p∆t$, ${\tau }_{l+0.5}^{p-0.5}$ and ${\gamma }_{l+0.5}^{p+0.5}$ represent the shear stress and strain of the medium at the midpoint of the soil layer calculated at the lth time $p∆t+0.5∆t$, respectively; ${\rho }_l$ is the mass density of the medium calculated at the lth time, and L represents the total number of soil layers calculated, excluding the input substrate.

The Eq. 24 is adopted to combine the boundary conditions with initial conditions, it can gradually integrate the velocity of the medium at the top of each layer and the shear stress and shear strain response of the medium at the midpoint of the layer. The corresponding acceleration and displacement can be derived from the velocity response quantity. The initial conditions for calculation are as follows:

${\mathit{v}}_{\mathit{l}}^{\mathbf{0}}=\mathbf{0},\mathit{l}=\mathbf{1,2},\dots ,\mathit{L}(25)$

${\mathit{\tau }}_{\mathit{l}-\mathbf{0.5}}^{-\mathbf{0.5}}=\mathbf{0},\mathit{l}=\mathbf{1,2},\dots ,\mathit{L}(26)$

The boundary conditions are defined as follows:

${\mathit{\tau }}_{-\mathbf{0.5}}^{\mathit{p}-\mathbf{0.5}}=\mathbf{0}(27)$

${\mathit{v}}_{\mathit{L}+\mathbf{1}}^{\mathit{p}}={\mathit{v}}_{\mathbf{b}}^{\mathit{p}}(28)$

In Eq. 28, $v_{\mathrm{b}}^p$ refers to the motion speed of the medium at the base top surface is calculated at time $p∆t$. In practical engineering, only the incident wave field (ascending wave field) in the base half space is given. If the substrate is regarded as an elastic body, the radiation effect of energy from the soil layer to the substrate has to be considered in reaction calculation. In this way, it can estimate the real motion $v_{\mathrm{b}}^p$.

${\mathit{v}}_{\mathbf{b}}^{\mathit{p}}=\frac{\mathbf{2}∆\mathit{t}}{{\mathit{\rho }}_{\mathit{L}}∆{\mathit{z}}_{\mathit{L}}+{\mathit{\rho }}_{\mathbf{b}}{\mathit{c}}_{\mathbf{s}\mathbf{b}}∆\mathit{t}}\left[\frac{{\mathit{\rho }}_{\mathit{L}}∆{\mathit{z}}_{\mathit{L}}-{\mathit{\rho }}_{\mathbf{b}}{\mathit{c}}_{\mathbf{s}\mathbf{b}}∆\mathit{t}}{\mathbf{2}∆\mathit{t}}{\mathit{v}}_{\mathbf{b}}^{\mathit{p}-\mathbf{1}}+{\mathit{\rho }}_{\mathbf{b}}{\mathit{c}}_{\mathbf{s}\mathbf{b}}\left({\mathit{v}}_{\mathbf{I}\mathbf{b}}^{\mathit{p}}+{\mathit{v}}_{\mathbf{I}\mathbf{b}}^{\mathit{p}-\mathbf{1}}\right)-{\mathit{\tau }}_{\mathit{L}+\mathbf{0.5}}^{\mathit{p}-\mathbf{0.5}}\right](29)$

In the above equation, ${\rho }_{\mathrm{b}}$ and $c_{\mathrm{s}\mathrm{b}}$ are the density and shear wave velocity of the substrate medium, respectively; $v_{\mathrm{I}\mathrm{b}}^p$ is the velocity of the incident wave field at the top of the substrate at time $p∆t$, and ${\tau }_{L+0.5}^{p-0.5}$ represents the shear stress of the medium at the middle point of the Lth soil layer calculated at time $p∆t+0.5∆t$, and it satisfies the below relationship:

${\mathit{\tau }}_{\mathit{L}+\mathbf{0.5}}^{\mathit{p}-\mathbf{0.5}}=\frac{{\mathit{\tau }}_{\mathit{L}+\mathbf{0.5}}^{\mathit{p}}-{\mathit{\tau }}_{\mathit{L}+\mathbf{0.5}}^{\mathit{p}-\mathbf{1}}}{\mathbf{2}}(30)$

Figure 7 presents the explicit space-time finite difference grids. Among them, solid points indicate the known initial conditions or boundary conditions, hollow points represent the spacetime points to be solved, square points mark the stress (or strain) spacetime points, and circular points represent the velocity spacetime points. According to the grid, the speed and stress at each space-time point can be gradually found.

FIGURE 7

FIGURE 7. Space-time explicit finite difference grid.

In the above time-domain numerical integration method, the discretized time interval ∆t should be small enough to ensure computational stability due to the application of central difference algorithm. Thickness of soil calculation stratification should meet the requirements of Eq. 21.

4 Numerical examples

This section employs several basic numerical examples to validate the reliability of the time-domain non-linear calculation program developed for this paper.

4.1 Example 1

This example verifies the accuracy of the time-space finite difference method, which is introduced in the field seismic response by comparing with an ideal two-layer field and a complex multi-layer field of an actual major project.

As for the mechanical parameters of the two models (Tables 1, 2), site 1 is an ideal two-storey site, and site 2 is a complex multi-storey site for a practical major project, which contains a hard interlayer with thickness of 15.9 m. In this example, all media are assumed to be linear elastic media, and the damping ratio is set to 0.05. Meanwhile, site-2 was chosen to verify the accuracy of the time-space finite difference method when the difference between the upper and lower soil layers is significant. Whereupon, a comparison is performed between the numerical solution and the analytical solution.

TABLE 1

TABLE 1. Computational model of site 1.

TABLE 2

TABLE 2. Computational model of site 2.

The time history curves of upgoing waves at the top of input half space in terms of accelerated speed, velocity and displacement (incident wave) are indicated in Figure 8. Under the ground motion input as indicated in Figure 8, the time-history curve of ground motion on site 1 surface obtained with the time-space finite difference method as described in Section 3.3 is compared with the calculation result obtained with the frequency domain method, as shown in Figure 9, with the comparison result of response spectrum curve presented in Figure 10. Figures 11, 12 indicate the corresponding results of site 2.

FIGURE 8

FIGURE 8. Time history curve of upgoing wave at the top of input half space (A) Accelerated speed. (B) Velocity. (C) Displacement.

FIGURE 9

FIGURE 9. Comparison between ground motion time history curves of site 1. (A) Accelerated speed. (B) Velocity. (C) Displacement.

FIGURE 10

FIGURE 10. Comparison of ground motion response spectrum curves of site 1.

FIGURE 11

FIGURE 11. Comparison between ground motion time history curves of site 2. (A) Accelerated speed. (B) Velocity. (C) Displacement.

FIGURE 12

FIGURE 12. Comparison of ground motion response spectrum curves of site 2.

It can be observed that the calculation results of the two methods are basically consistent. It needs to be noted that the time-space finite difference method adopts the central difference method to replace the time and spatial derivatives of the wave equation, and then it determines the numerical solution of the fluctuation problem by explicit numerical integration, while the frequency domain method determines the analytical expression of the steady-state solution followed by the transient reaction by Fourier transformation. Moreover, when referring to the time-space finite difference method, the artificial boundary processing method is adopted to handle the infinite domain problem of the input half space, while the frequency domain method may eliminate the infinite domain truncation problem, without considering the exact accuracy of the numerical error of the discrete Fourier transformation itself. The consistency of the calculation results as presented in Figures 9–12 indicates that it is reliable to solve one-dimensional site seismic response with time-space finite difference method by the numerical program developed in this work.

4.2 Example 2

In this section, the data from a certain practical engineering site and the KiK-net station was introduced. With contrastive analysis of DEEPSOIL (Hashash, V6.1) and CHARSOIL (Streeter et al., 1974), the non-linear time-domain programs used worldwide, the reliability of the program developed herein was further demonstrated.

Presently, DEEPSOIL is the most extensively adopted time-domain non-linear soil layer seismic response analysis and calculation program in the world, representing the highest level of time-domain non-linear analysis. Therefore, in this section, the actual project site with silt soil layer is used as the calculation model (the mechanical parameters of the model are shown in Table 3). The Northridge strong seismic observation data is used as input for a comparative analysis of the proposed method (the constitutive model adopts the extended Masing criterion, same as below) and the DEEPSOIL time domain algorithm. Figures 13, 14 depict the calculated surface acceleration time history and response spectrum. In comparison, the proposed method is consistent with the time-history curve of surface acceleration calculated by DEEPSOIL, with a marginally lower amplitude than the DEEPSOIL results. According to the comparison curves of the reaction spectrum, the proposed method for fitting the long-period reaction spectrum is relatively consistent. While the high-frequency component is less than the DEEPSOIL calculation results. The reason for the difference is inferred to be the difference between the backbone curves of the non-linear constitutive model of rock and soil and the finite difference method utilized by the two programs.

TABLE 3

TABLE 3. Computational model of site 3.

FIGURE 13

FIGURE 13. Comparison between ground motion time history curves of site 3 (Accelerated speed).

FIGURE 14

FIGURE 14. Comparison of ground motion response spectrum curves of site 3.

Observation station data from the Japanese KiK-net network at a class II site are used concomitantly. The observation point is used as the input; and by comparing and analyzing the calculation results of various time-domain non-linear programs (the proposed method, DEEPSOIL, CHARSOIL) with the actual surface observation records corresponding to the observation points, the reliability of the non-linear program in this paper is verified.

Figure 15 depicts the comparison results of the seismic response spectrum curve at the surface. In accordance with the figure, the research methodology and DEEPSOIL calculation findings are comparable to the actual surface observation record. The overall reaction spectrum curve is essentially consistent with the actual surface observation record, especially the curves of the long period (T > 2 s), which completely overlap, with the exception of the large fitting error of individual reaction spectrum periods (such as T = 0.5 s ∼ 1 s). Whereas the GHARSOIL calculation result represents a significant error.

FIGURE 15

FIGURE 15. Comparison of ground motion response spectrum curves of different programs with actual observation records.

In terms of this section, the non-linear method is comparable to the internationally recognized DEEPSOIL program’s calculation results and is superior to the CHARSOIL when compared to actual observation records. In contrast to DEEPSOIL, the developed procedure is simple, quick, and has a high modeling efficiency, particularly when the site soil layer is thick and the soil is abundant.

5 Conclusion

In order to solve the problems such as the uncertainty of site non-linear seismic response results caused by a single calculation method and the restriction of traditional site non-linear calculation program on the calculation scale, a site non-linear seismic response analysis program is developed based on fortran95 software platform for the time-space finite difference method in this work. To facilitate the consideration of the uncertainty effect of various soil constitutive models under the TNL method, the program is embedded in two relatively mature soil constitutive models at the present stage, namely the extended Masing criterion model and the dynamic skeleton curve model. Simultaneously, port is reserved for the further introduction of various soil constitutive models, and it can be a tool platform for the uncertainty research work of soil non-linear dynamic change. Moreover, the program uses dynamic array technology to store relevant data, and it has no limitation on the computational scale of the site model, the time points and space points of the input and output data principally. In addition, the program can be used for large-scale site seismic response calculation work, and facilitate the research work of site adjustment coefficient in the new generation of zoning map in China. Finally, the reliability of the developed calculation program are verified based on the practical examples, and through a comparative analysis of a relatively mature time-domain non-linear program and testing with actual records from the KiK-net strong motion observation station in Japan, the time-nonlinear program developed in this paper demonstrates a high degree of calculation reliability and modeling efficiency. Consequently, it is applicable to the calculation of more complex site conditions and large-scale site seismic response. Due to the singularity and uniqueness of the example model, additional real site models should be included in the follow-up work to validate the procedure’s applicability. Moreover, an uncertainty analysis of various site calculation methods should be carried out.

Data availability statement

The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

Author contributions

All authors listed have made a contribution to conception and design of the study, and approved it for publication. Main program development: YZ. Collection and processing of data: YZ and JY. Writing-original draft: JY. Writing-reviewing and editing: YZ and JY. All authors contributed to manuscript revision, read, and approved the submitted version.

Funding

This study is supported by the Science for Earthquake Resilience of China Earthquake Administration (XH22013YA).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Beresnev, I. A. (2002). Nonlinearity at California generic soil sites from modeling recent strong-motion data. Bull. Seismol. Soc. Am. 92 (2), 863–870. doi:10.1785/0120000263

CrossRef Full Text | Google Scholar

Bo, J. S. (1998). Site Classification and design response spectrum adjustment method (postdoctoral research Report). China Earthquake Administration: Institute of Engineering Mechanics, 55–58.

Google Scholar

Chen, X. L. (2006). Study on soil dynamic properties, and nonlinear seismic response of complex site and its methods. China Earthquake Administration, Institute of Engineering Mechanics.

Google Scholar

Diego, C. F. L. P., Carlo, G. L., and Ignazio, P. (2006). ONDA: Computer code for nonlinear seismic response analyses of soil deposits. J. Geotech. Geoenviron. Eng. 132 (2), 223–236.

Google Scholar

Gao, S., Liu, H., and Davis, P. M., (1996). Localized amplification of seismic waves and correlation with damage due to the northridge earthquake: Evidence for focusing in santa monica. Bull. Seismol. Soc. Am. 869 (1B), 209–230.

Google Scholar

Hashash, Y. M. A. (2009). DEEPSOIL V6.1, 12, user manual and tutorial.

Google Scholar

Idriss, I. M., and Seed, H. B. (1968). Seismic response of horizontal soil layers. Soil Mechanics and Foundations Division. J. Soil Mech. Found. Div. 94 (4), 1003–1031. doi:10.1061/jsfeaq.0001163

CrossRef Full Text | Google Scholar

Kavazanjian, E., and Matasovic, N. (1995)., 46. New Orleans, Louisiana, USA, 1066–1080.Seismic analysis of solid waste landfillsGeoenvironmental 2000, Geotech. Spec. Publ.2

Google Scholar

Kim, B., and Hashash, Y. M. A. (2013). Site response analysis using downhole array recordings during the march 2011 tohoku-oki earthquake and the effect of long-duration ground motions. Earthq. Spectra 29 (S1), 37–54. doi:10.1193/1.4000114

CrossRef Full Text | Google Scholar

Kokusho, T., and Sato, K. (2008). Surface-to-base amplification evaluated from KiK-net vertical array strong motion records. Soil Dyn. Earthq. Eng. 28 (9), 707–716. doi:10.1016/j.soildyn.2007.10.016

CrossRef Full Text | Google Scholar

Li, P. (2010). The effect of site types on platform value of the design response spectrum. China Earthquake Administration: Harbin: Institute of Engineering Mechanics.

Google Scholar

Li, X. J. (1992). A method for analysing seismic response of nonlinear soil layers. South China J. Seismol. 12 (4), 1–8.

Google Scholar

Li, X. J. (1992). A simple function expression of the soil dynamic constitutive relationship. Chin. J. Geotechnical Eng. (05), 90–94.

Google Scholar

Li, X. J. (2013). Adjustment of seismic ground motion parameters considering site effects in seismic zonation map. Chin. J. Geotechnical Eng. 35 (S2), 21–29.

Google Scholar

Li, Y. Y., Xu, Y., and Li, D. M. (2005). Analysis of earthquake responses for different site categories. Earthq. Res. Shanxi 4, 27–33.

Google Scholar

Luan, M. T., and Lin, G. (1994). An effective time-domain computational method for nonlinear seismic response analysis of soil deposit. J. Dalian Univ. Technol. 34 (2), 228–234.

Google Scholar

Luan, M. T., and Lin, G. (1992). Computational model for nonlinear analysis of soil site seismic response. Eng. Mech. 9 (1), 94–103.

Google Scholar

Luan, M. T. (1992). Ramberg-osgood constitutive model with variable parameters for dynamic nonlinear analysis of soils. Earthq. Eng. Eng. Vib. (02), 69–78.

Google Scholar

Martin, P. P., and Seed, H. B. (1982). One-dimensional dynamic ground response analyses. J. Geotech. Engrg. Div. 108 (7), 935–952. doi:10.1061/ajgeb6.0001316

CrossRef Full Text | Google Scholar

Martin, P. P., and Seed, H. B. (1978). MASH-A computer program for the non-linear analysis of vertically propagating shear waves in horizontally layered deposits. Calif: Report EERC, University of California at Berkeley, 78–23.

Google Scholar

Masing, G. (1926). Eigenspannungeu and verfertigung beim Messing. Zurich: Proceedings of the 2nd International Congress on Applied Mechanics.

Google Scholar

Matasovic, N., and Vucetic, M. (1993). Cyclic characterization of liquefiable sands. J. Geotech. Engrg. 119 (11), 1805–1822. doi:10.1061/(asce)0733-9410(1993)119:11(1805)

CrossRef Full Text | Google Scholar

Pyke, R. M. (1979). Nonlinear soil models for irregular cyclic loadings. J. Geotech. Engrg. Div. 105 (6), 715–726. doi:10.1061/ajgeb6.0000820

CrossRef Full Text | Google Scholar

Qi, W. H., Wang, Z. Q., and Bo, J. S. (2010). Development and verification of a method for analyzing the nonlinear seismic response of soil layers. J. Harbin Eng. Univ. 31 (4), 444–450.

Google Scholar

Qi, W. H., and WangBo, Z. Q, J. S. (2010). Development and verification of a method for analyzing the nonlinear seismic response of soil layers. J. Harbin Eng. Univ. 31 (04), 444–450.

Google Scholar

Seed, H. B., Romo, M. P., Sun, J. I., Jaime, A., and Lysmer, J. (1988). The Mexico earthquake of september 19, 1985: Relationships between soil conditions and earthquake ground motions. Earthq. Spectra 4 (4), 687–729. doi:10.1193/1.1585498

CrossRef Full Text | Google Scholar

Streeter, V. L., Wylie, E. B., and Richart, F. E. (1974). CHARSOIL, 100(3), 247-263, characteristics method applied to soils.

Google Scholar

Wang, Z. L., and Han, Q. Y. (1981). Analysis of wave propagation for the site seismic response, using the visco-elastoplastic model. Earthq. Eng. Eng. Vibraion (01), 117–137.

Google Scholar

Yee, E., Stewart, J., and Tokimatsu, K. (2013). Elastic and large-strain nonlinear seismic site response from analysis of vertical array recordings. J. Geotech. Geoenviron. Eng. 139 (10), 1789–1801. doi:10.1061/(asce)gt.1943-5606.0000900

CrossRef Full Text | Google Scholar