Three-Dimensional Site Response Analysis of Clay Soil Considering the Effects of Soil Behavior and Type

Abstract

To understand changes in bedrock motion at the ground surface, frequency effects, and spatial distribution within the soil, it is important to look at how a site responds to earthquakes. This is important for soil–structure interaction in structural and geotechnical earthquake engineering. This study deals with the effect of classifying clays according to shear wave velocity (stiff/medium/soft) and nonlinearity in behavior (linear/nonlinear) on the analysis of the site response. A 3D soil model with a combination of free fields and quiet boundaries and advanced constitutive models for soil to obtain accurate results was used to conduct this study. A strong TABAS earthquake was used to excite the compliant base of the model after converting the velocity record of TABAS to an equivalent surface traction force using a horizontal force–time history proportional to the velocity–time history. This study reveals that the site response analysis is affected by the type of clay soil and the soil material behavior, with soft clay soil causing higher PGV and PGV values in the linear case and lower values in the nonlinear case due to soil yielding, which causes soil response attenuation. This results in extremely conservative and expensive building designs when linear soil behavior is adopted. On the other hand, the applied earthquake exhibits greater attenuation at longer frequencies and greater amplification at mid and short frequencies. However, at frequencies near the applied earthquake frequency, neither attenuation nor amplification occurs. Furthermore, nonlinear soil behavior is crucial for soil evaluation and foundation design due to higher octahedral shear strain and settlement values, especially in softer soils, resulting from extensive plastic deformation.

1. Introduction

Site response analysis, also known as free-field reaction, is the process of estimating how soil-sediment interactions interact with ground motion in the absence of structures. Similarly, it serves as an introduction to soil–structure interaction (SSI) analysis and how to calculate free-field response [1].

A great deal of experimental and analytical work has been done to characterize the cyclic shear behavior of soils because it is well known that during an earthquake, soils close to the surface and far enough away from built structures will repeatedly deform in shear, resulting in a variety of stresses and strains on the soil [2]. The clay soil types and their structural variations influence buildings and infrastructure. Many variables affect how much swell and shrinkage occurs, including moisture content, confining pressure, dry density, the kind and quantity of clay minerals, soil structure, and climate change. The different structures eventually sustain significant damage as a result of these volume changes. Loss of strength due to moisture causes significant damage to buildings and foundations [3]. Soft clay soils can cause several technical problems, including bearing capacity failure, excessive settlement, and slope instability, both during and after construction. This is due to their low shear strength combined with their high compressibility. This is why geotechnical engineers have traditionally found it difficult to build structures on soft clay soils [4].

Soil collapses under the influence of strong earthquakes are one of the important reasons for the collapse of buildings above them, as the soil’s characteristics impact seismic vibrations that increase pressures arising in the soil. The response spectra calculated on-site can also exceed the design seismic levels given in the codes, such as American Codes, and the maximum considered earthquakes can exceed their specified levels in some locations. This was demonstrated by reviewing the causes of the destruction and collapse of buildings as a result of the 6 February 2023, earthquake that struck Turkey and Syria [5]. For all these reasons, it became evident that several studies are required to analyze the site response using advanced soil-material models and 3D-analytical models with suitable boundary conditions to get general recommendations and insights to prevent similar damage and collapse in the future.

Many researchers have investigated how soil deposits react and how the conditions affect the largest expected earthquakes, which is a crucial part of the procedure of seismic design. One-dimensional ground response analysis has been adopted by many researchers, such as [6,7,8,9,10,11]. The study of Garala and Madabhushi (2019) [6], which investigated the significance of the nonlinear 1D ground response analysis of the soft clay and the friction piles using sinusoidal excitations, pointed out the importance of using more realistic motions in the future.

The site response is more or less amplified due to local heterogeneities and nonlinear behavior of the subsoil, as is shown for Vietnam by Nguyen et al. (2020) [7], for Turkey by Silahtar A (2023) [10], and for Japan by Zhang et al. (2024) [12]. Such a result is also found in the study by Guzel et al. (2023) [11].

The study by Che et al. (2024) [13] also shows the impact of the long-term and short-term entrance movement on the retrospective site responses to class soil deposits, where significant differences of frequencies appear > 0.5 Hz, especially in the more softening sites concerning the largest ${PGA}_{input}$ levels.

The inadequacy of the one-dimensional site response analyses has been mentioned by many researchers to consider the heterogeneity in the strata below the earth’s surface [14,15,16]. To account for the impact of spatial differences between these layers’ thickness and stiffness, two-dimensional site response assessments were performed [17,18,19,20]. However, these analyses are still in the initial phase for small sites in meters, despite their significance to the design of the shallow bedrock region. Most of the research has focused on large-scale basins that were several kilometers in length.

The types of soil and intensities of motions affected the 2D response of both the basin sites and small-scale subsurface profiles, as found by Khanbabazadeh and Iyisan (2014) [21] and Chandran and Anbazhagan (2020) [22], respectively. The study by Chala and Ray (2024) [23] demonstrated the role of local soil conditions and types in shaping ground-surface responses and influencing the intensity of ground shaking.

The 2D case does not apply to all geotechnical study scenarios and is only appropriate for analysis under plain strain conditions. Its use is limited to situations with narrow widths and long lengths, including slope stability, continuous foundations, retaining walls, etc. Even in these situations, if there is a change in the soil along the longitudinal direction, the 2D solutions become an inaccurate generalization, and the use of 3D becomes necessary. Although 3D models are preferable to 2D models, limited 3D models have recently been used to analyze large sites’ responses, such as [24,25], without paying attention to small sites’ responses for a particular building construction, which is an essential step in analyzing soil–structure interaction.

It is evident from the previous research that:

• The simplest approach is 1D with linear behavior of a simple soil structure [6,7,8,9,10,11].
• Several publications show that non-linear behavior and more complex soil structure have led to amplification and deamplification [7,10,11,12,13].
• 2D simulation showed that the geometry of the basin (e.g., edge inclinations) has an effect [14,15,16,17,18,19,20].

Shortcomings are as follows:

• No 3D simulation combines advanced nonlinear soil material and suitable boundaries.
• No building–soil interaction with 3D advanced site response analysis.
• No study on the effect of clay soil types on 3D site response analysis.

The key objectives and questions of this paper are as follows:

• How big is the range of deviation of the ground motion due to soil behavior at the foundation edge?
• How is the 3D field of the ground motion within a 3D basin? How big are the spatial deviations?
• Is there a frequency dependency in the site response? For which frequency is the site response amplified or deamplified?

This research is limited to using one earthquake, specifying one type of soil that is clay, and using one layer of soil over the entire height of the model to achieve the desired objective of this study.

2. Materials and Methods

To address the research questions mentioned above, we ran a 3D simulation of the Opensees FEM [26] and STKO [27] as a pre- and post-processor, observing that to differentiate the ‘commands’ of Opensees and STKO, they are written in quotation marks.

2.1. Finite 3D Model

2.1.1. Geometry

Utilizing eight-node ‘SSPbrick’ hexahedral elements, which are effective against shear/volumetric lock issues in finite element computations, 3D homogenous soil sediments with bedrock 50 m below the ground surface were modeled [28]. The studied soil’s self-weight is considered a body force acting on each element and is set as the weight unit of the soil. Figure 1 presents the dimensions and mesh of the used FE model of the studied soils. Although the input wave’s frequency content is the same for all three types of soil, the stiff soil’s vertical mesh is the coarsest because it has the highest shear wave velocity. The number of mesh entities and partitions used to divide the model for parallel analysis (using Openseemp solver) are displayed respectively in

Table A1. The number of mesh entities that comprise the three studied clay soil models.

Entity Name	Clay Soil Type
	Soft Clay	Medium Clay	Stiff Clay
nodes	69,762	40,746	31,074
elements	173,269	103,993	80,901
edges	27,016	18,184	15,240
faces	75,465	45,357	35,321
solids	59,072	33,584	25,088
interactions	11,716	6868	5252
partitions	8	8	8

and Figure A1, which are in Appendix A at the end of the paper. For more details, see Reference [29].

In dynamic studies, the features of the shear wave velocity “$V_s$” (assuming vertical wave propagation) and the frequency content “$f$” of the input wave can have an impact on the accuracy of wave transmission. The size of the element “$\Delta l$” should be less than one-tenth to one-eighth of the wavelength associated with the highest frequency component of the input wave, according to Kuhlmeyer and Lysmer (1973) [30]; that is,

$$\Delta l\le {\displaystyle \frac{V_s/f}{(8\to 10)}}$$

(1)

The study of the fast Fourier transform (FFT) is used to determine the earthquake’s frequency content. The mesh size selection is influenced by the frequency content of the input motion to achieve precise wave propagation. To have an accurate representation, the maximum region size should be selected using Equation (1). Eliminating the high-frequency components may enable a coarser mesh if the highest input motion frequencies demand an extremely fine mesh (and a similarly small-time step). If the most input power is detected in low-frequency components (e.g., 80% to 90%), the history can be filtered to exclude higher frequencies without appreciably altering the results [31].

In dynamic analysis, it is crucial to replicate the appropriate degree of frequency-independent damping of materials. A certain amount of energy must be dispersed via damping because of the material’s elastic behavior. Damping for geological materials usually lies between 2 and 5% of the critical. However, substantial energy dissipation can happen during plastic flow if one of the material’s plasticity-forming models is employed in the analysis. Consequently, at low deformation levels, small Rayleigh damping is required to prevent system resonance. As a result, in linear analysis, as opposed to nonlinear analysis, the damping selection is more crucial [32].

The Kuhlemeyer and Lysmer criteria may necessitate an extremely fine spatial grid and a correspondingly small temporal step for dynamic inputs with high peak velocity and small-time step. As a result, doing a reasonable analysis can take a lot of time and memory. When this happens, it could be possible to adjust the input by seeing that the majority of the input log’s energy is contained in lower-frequency components. By using a Fourier transform (FFT) and filtering the log to exclude high-frequency components that only marginally contribute to the load, a coarser grid can be used without substantially altering the outcomes [32,33].

Thus, the chosen f in Equation (1) eliminates the high frequencies that add very little to the seismic loading and covers the band that includes the most significant energy. It also includes the basic frequency of structure and soil. This shortens the computation time and produces a suitable mesh without appreciably altering the outcome.

2.1.2. Boundary Conditions

For soil domain borders, free and quiet field boundaries are utilized in addition to standard static boundary conditions. A quiet boundary simulation employing viscous dashpots is used in Figure 2 to show how the lateral borders of the internal soil core grid and the free-field grid are coupled. The primary grid border receives the unbalanced pressures from the free field.

Figure 3 shows the Lysmer–Kuhlemeyer dashpots that are incorporated, in the shear and normal directions, into (a) the lateral edges of the interior soil domain, using (distributed) zero-length elements and viscous uniaxial materials, to connect it to the adjacent free fields and also to connect between the adjacent free fields and (b) the base of the model using a ‘DistributedLK3D’, which is a boundary condition (mp constraint with viscous materials inserted in the three dimensions x, y, z) available in STKO, not in Opensees. The viscous material of these dashpots required a damping coefficient of “$\mathit{C}$”, which is defined as:

$$C=\rho ·V·A$$

(2)

where $\rho$ is the mass density of the (soil or bedrock); $V$ is the (soil or bedrock) p- and s-wave velocities, which result in three components of damping coefficient $\mathit{C}$: two in the tangential or shear directions ${\mathit{C}}_{\mathit{v}}$ and one in the normal direction ${\mathit{C}}_{\mathit{p}}$; and $\mathit{A}$ is the area of the (soil or bedrock), noting that unit area is used here due to activation of the ‘distributed’ option, which is an STKO option that does not exist in Opensees. This option distributes the material values assigned to the element by area unit and distributes them to the area by calculating the attributed area automatically. These values were calculated using the characteristics of the soils for LK dashpots at the lateral boundaries and of the bedrock for LK dashpots at the base.

The free field is represented by using a solid element with equal DOF conditions (displacement constraints for ${\mathit{U}}_{\mathit{x}},{\mathit{U}}_{\mathit{y}},{\mathit{U}}_{\mathit{z}}$) between its faces in the same direction for both X and Y short directions. The Node-to-Node Links option of interaction type is used for tying the nodes of those faces together, and thus the simple shear deformation pattern is achieved. The two kinds of Lysmer–Kuhlemeyer dashpots, in addition to equal DOF conditions within free fields are presented in Figure 3 for the studied 3D model.

A brief illustration of the boundary conditions used is introduced in this section as follows:

• Compliant Base

The study conducts a site response analysis of three-dimensional soil deposits under an elastic half-space. To simulate the finite stiffness base half-space, Lysmer–Kuhlemeyer (1969) [34] dashpots are incorporated into the 3D soil model base using ‘DistributedLK3D’. The mass density, shear wave velocity, and “unit” area of the base are multiplied to get the unit factor. This method takes into account the base layer’s finite stiffness by allowing energy to radiate back into the foundation material. The base of the soil model is excited with a horizontal force–time history proportionate to a time history of the ground motion velocity under investigation. For further details on this modeling method, see Joyner and Chen (1975) [35] and Lysmer (1978) [36], and other sources.

• Quiet Boundaries

Modeling of geomechanics problems includes media that are better represented, on the scale of analysis, as infinite. Surface and near-surface constructions are considered to be situated in a semi-space, and deep subterranean excavations are assumed to be surrounded by a limited medium. For static studies, boundary elements could be placed at some distance from the area under study to achieve appropriate boundary conditions. Lysmer and Kuhlemeyer (1969) [34] have developed viscous boundaries that rely on the application of independent dashpots at the model boundaries in both the shear and normal directions. In dynamic studies, these quiet limits serve as absorbing barriers that let outgoing propagating waves pass through while blocking reflections back into the model. The computational burden of reducing the problem by implementing a larger model to absorb most energy in waves reflected from far boundaries will be significant.

• Free-Field Boundaries

Boundary conditions on the model’s sides must account for the free-field motion that can arise if the structure is not present. To achieve free-field conditions and lessen wave reflection, these boundaries should be placed at suitable distances. For soils with a large amount of material damping, relatively small distances could be set for this condition [37]. On the contrary, the large distances required in the case of low material damping could result in an impractical large model. By using the free-field calculations in tandem with the main grid analysis, a different process for free-field situations has been established, preserving the non-reflective characteristics of the borders [33].

2.1.3. Ground Motion

The model base is excited in the X-direction using the horizontal component “RSN139_DAY-L1” of the velocity–time series of “TABAS” [38]. The value of the velocity history of the Tabas earthquake multiplied by the bedrock’s ${\mathit{C}}_{\mathit{v}}$ yields the real force that is applied to the base’s face at each time step. Figure 4a,b displays acceleration and velocity time histories of the input wave “TABAS”. Figure 4c presents the corresponding Fourier spectrum of the velocity of “TABAS” obtained using the GiD+OpenSees Interface [39]. It can be seen that the frequencies involved are lower than 4.8 Hz, with a maximum peak of $f_{seismic load} = 0.12 \mathrm{H}\mathrm{z},$ which is the same as the first peak of the seismic load. The value of $f = 4.8 \mathrm{H}\mathrm{z}$ is used in Equation (1) to determine the suitable mesh of the soil model. This value covers the fundamental frequencies ${\mathit{f}}_{\mathit{s}\mathit{o}\mathit{i}\mathit{l}}$ of the stiff, medium, and soft soils, which are equal to 1.44, 1, and 0.5, respectively, which is determined using the following formula:

$$f_{soil}={\displaystyle \frac{V_s}{4H}}$$

(3)

where $H$ and $V_s$ stand for the soil layer’s thickness and shear velocity, respectively.

Refer to [29] for further details and other information.

2.2. Material Parameters

The suggested parameter values for clay soils in the references [40,41] are adopted in this research.

Table 1. Characteristics of the studied clay soils and bedrock; given and calculated values for soil properties.

	Soil Properties	Soft Clay	Medium Clay	Stiff Clay	Bedrock
Given Values	Mass density $\rho$ (kg/m3)	1300	1500	1800	2400
	Ref. Shear Modul. $G_r$ (MPa)	13	60	150	4000
	Ref. Bulk Modul. $B_r$ (MPa)	65	300	750	6667
	Cohesion (MPa)	0.018	0.037	0.075	-
	Peak Shear Strain (at $p_r^{\prime }=0.08 \mathrm{M}\mathrm{P}\mathrm{a}$)	0.0001	0.0001	0.0001	-
	Friction Angle	0	0	0	-
	PressDependCoe.	0	0	0	-
Calculated Values	Weight (kN)	12.8	14.7	17.7	23.5
	VP (m/s)	252	505	726	2236
	VS (m/s)	100	200	289	1291
	Poisson’s Ratio ν	0.41	0.41	0.41	0.25
	Elastic Modulus E (MPa)	36.6	168.8	421.8	10000
	Cp (Normal) Damping Coefficient (kN·s/m)	328	758	1306	5367
	Cv (Shear) Damping Coefficient (kN·s/m)	130	300	520	3098

The first section of the table represents the soil’s characteristic features (Given Values), from which the values of other properties (Calculated Values) in the second section are derived.

presents the characteristics of the studied clay soils and the bedrock layer. Automatic surface generation is used for PIMY material, with the number of yield surfaces NYS = 20.

2.3. Material Behavior

Both linear elastic behavior and nonlinear elastic-plastic behavior of clay material are used in this research to explore the diverse responses of site soil when subjected to a strong earthquake. ‘PressureIndependMultiYield’ nD material (PIMY) [40,41] is employed to depict clayey soils’ nonlinear behavior, while ‘ElasticIsotropic’ nD material [42] is used for the linear elastic case. ‘InitialStateAnalysisWrapper’ nD material is used with the previous two materials for setting the initial conditions of the problem and analyses [43,44].

2.3.1. ‘PressureIndependMultiYield’ Material

‘PressureIndependMultiYield’ nD is an elastoplastic material that mimics the monotonic and cyclic behavior of materials like clay and organic soils, exhibiting an insensitive shear response to confinement. The elasticity lies in the volumetric response, while the plasticity exhibits only in the deviatoric response. The concept of multiple Von Mises-type surfaces, with an associative flow rule, formulates this plasticity. The response of this material is elastic and linear when static gravity loads are applied. In contrast, during the next dynamic (fast) loading phase, the stress–strain response is elastic-plastic [40,41].

The soil’s nonlinear (shear–strain) hyperbolic backbone curve and a piecewise linear representation of the plasticity with multiple surfaces are presented in Figure 5. In stress space, yield surfaces—a set of many plasticity surfaces—define regions with constant shear modulus and provide a linear representation of the hyperbolic backbone curve. The yield surface ${\mathit{f}}_{\mathit{m}}$ with shear modulus ${\mathit{H}}_{\mathit{m}}$ for $\mathit{m}=1,2,\dots \mathit{N}\mathit{Y}\mathit{S}$ is represented by each line segment in Figure 5, where NYS is the total number of yield surfaces [45]. The failure surface, or zero shear modulus ${\mathit{H}}_{\mathit{N}\mathit{Y}\mathit{S}}$, is represented by the outer yield surface ${\mathit{f}}_{\mathit{N}\mathit{Y}\mathit{S}}$, which is equivalent to the ultimate shear strength ${\mathit{\tau }}_{\mathit{m}\mathit{a}\mathit{x}}$.

The maximal octahedral shear force ${\mathit{\tau }}_{\mathit{f}}$, which is a function of the effective current confinement ${\mathit{p}}_{\mathit{i}}^{\prime }$, is determined by the friction angle $\mathit{\varphi }$ and cohesion $\mathit{c}$:

$${\tau }_f={\displaystyle \frac{2\sqrt{2}\sin \varphi }{3-\sin \varphi }}{p}_i^{\prime }+{\displaystyle \frac{2\sqrt{2}}{3}}c$$

(4)

The hyperbolic curve for octahedral shear (stress–strain) to generate automatic surfaces at a fixed confinement ${\mathit{p}}^{\prime }$ is determined as follows:

$$\tau ={\displaystyle \frac{G\gamma }{1+\frac{\gamma }{{\gamma }_r}{\left(\frac{p_r^{\prime }}{p^{\prime }}\right)}^d}}$$

(5)

as ${\mathit{\gamma }}_{\mathit{r}}$ satisfies the following equation at ${\mathit{p}}_{\mathit{r}}^{\prime }$:

$${\tau }_f={\displaystyle \frac{2\sqrt{2}\sin \varphi }{3-\sin \varphi }}{p}_r^{\prime }+{\displaystyle \frac{2\sqrt{2}}{3}}c={\displaystyle \frac{G_r{\gamma }_{max}}{1+{\gamma }_{max}/{\gamma }_r}}$$

(6)

The key elements of the multi-yield surface plasticity model, which is independent of the applied stress, are the yield surface, hardening law, and the flow rule [40,47,48,49].

This command is used to create a PressureIndependMultiYield nDMaterial object:

nDmaterial PressureIndependMultiYield $tag $nd $rho $refShearModul $refBulkModul $cohesi $peakShearStra <$frictionAng=0. $refPress=100. $pressDependCoe=0. $noYieldSurf=20 <$r1 $Gs1 [40,41]

For more details about this material and its parameters, refer to [40,41].

2.3.2. ‘InitialStateAnalysisWrapper’ Material

The ‘InitialStateAnalysis’ command can be used to specify initial conditions by using the ‘InitialStateAnalysisWrapper’ nD material. The field of the initial stress is developed while the original problem geometry is maintained when this material is used [43,44]. The ‘InitialStateAnalysis’ command can be used with the ‘InitialStateAnalysisWrapper’ to generate a gravitational state of stress in the soil elements. At the end of the analysis, there should be non-zero stress and strain in the soil model with zero displacement [44]. In other words, the soil is stable under the influence of its weight; that is, its vertical displacement is zero with initial stresses and deformations in it. The ‘InitialStateAnalysisWrapper’ is a bit tricky to use for dynamic analysis. Sometimes, setting the displacement to zero appears to be interpreted as an initial displacement in later steps, resulting in undesirable vibrations [43,44]. To get rid of these undesirable vibrations caused by the initial state analysis, a static analysis of the vertical loads from the soil self-weights should be carried out for the resting phase after the initial state analysis has been turned off.

This command is used to construct an InitialStateAnalysisWrapper nDMaterial object:

nDMaterial InitialStateAnalysisWrapper $matTag $nDMatTag $nDim [43,44]
$matTag	unique integer tag identifying nDMaterial object
$nDMatTag	the tag of the associated nDMaterial object. (Here, it is the number of “$\mathit{\$}\mathit{t}\mathit{a}\mathit{g}$” that is used to define the material of ‘ElasticIsotropic’ or ‘PressureIndependMultiYield’.)
$nDim	number of dimensions (2 for 2D, 3 for 3D)

2.3.3. Update Material Stage

As previously mentioned, there are two stages to gravity loading (from soil self-weight only): elastic loading and plastic loading, which may be modified using updateMaterialStage. This material update phase is only applied when using nonlinear soil behavior (where ‘PressureIndependMultiYield’ is included in the ‘InitialStateAnalysisWrapper’ definition), and there is no need to apply this phase when studying linear behavior only (where ‘ElasticIsotropic’ is included in the InitialStateAnalysisWrapper definition). The script for defining the material update phase is as follows:

updateMaterialStage -material $tag -stage $sNum [40,50]
$tag	Material number. (Here, it is the number of “$\mathit{\$}\mathit{m}\mathit{a}\mathit{t}\mathit{T}\mathit{a}\mathit{g}$” that is used to define ‘InitialStateAnalysisWrapper’ material.)
$sNum	desired stage: 0—linear elastic 1—plastic

Two stages can be followed to perform the seismic analysis: First, the material stage is set to 0, with elastic coefficients B_r and G_r when fixed gravity loads are applied. After the application of these gravity loads, the material stage is adjusted to 1. For Stage 1: the (stress–strain) deviatoric response, in the subsequent fast dynamic analysis, turns to elastic-plastic, while the volumetric behavior is still linear and elastic [40,50].

In this research, these two material stages are applied when the ‘PressureIndependMultiYield’ material is used. However, only one stage is applied for ‘ElasticIsotropic’ material.

2.4. Analysis Steps

Analyses were performed with Penalty constraints, Parallel Reverse Cuthill–McKee Numberer, Mumps System, “Linear/Krylov–Newton” algorithms for linear/nonlinear cases of soil material, respectively, Norm Displacement Increment Test, and “Load Control/TRBDF2” for static and transient integrators, respectively.

The stages of construction and analyzing the model are divided into two stages:

 

STAGE 1: Construction of the soil and static analyses.

• Add the first-stage model, consisting of soil, by the ‘modelSubset’ command.
• Add constraints to the model, which are fixing the soil and free-field columns by fixing Ux, Uy, and Uz for the vertical surfaces perpendicular to the X and Y axes and the bottom of the soil, respectively.
• Run the initial state analysis (TCL script: InitialStateAnalysis on).
• Update the soil material to the elastic state (stage 0 for updateMaterialStage) (applied only when nonlinear material is used).
• Run static analysis of the vertical or gravity loads for the elastic state.
• Update the soil material to the plastic state (stage 1 for updateMaterialStage) (applied only when nonlinear material is used).
• Run static analysis of the vertical loads of the plastic state, (applied only when nonlinear material is used).
• Turn off the initial state analysis (TCL script: InitialStateAnalysis off).
• Wipe analysis (TCL script: wipeAnalysis).
• Run static analysis of the vertical loads of the rest phase.
• Record the results of the soil by MPCORecorder.
• Add Rayleigh damping of the soil.
• Remove fixities (defined in analysis step 2) by the ‘removeSpConstraints’.
• Add forces ${\mathit{F}}_{\mathit{X}}$, ${\mathit{F}}_{\mathit{Y}}$, and ${\mathit{F}}_{\mathit{Z}}$ from removed reactions, using a constant time series to apply them.

 

STAGE 2: Application of the viscous boundaries LK at the sides and base of the soil and equal DOF for the free-field columns.

• Add the second-stage model, consisting of viscous boundary “$\mathit{L}\mathit{K}1$” at the sides of the soil and equal DOF for free-field columns, by the ‘modelSubset’ command.
• Add constraints to the model, which are the viscous boundary “$\mathit{L}\mathit{K}2$” at the base of the soil and equal DOF for the free-field columns by the ‘constraintPattern’.
• Add the seismic load after introducing a dynamic force down the soil in the X-direction.
• And then perform nonlinear dynamic analysis.

3. Results

In this section, the site response analysis’s dynamic analysis results for three different soil types—stiff, medium, and soft clay—are shown. These results include some essential parameters for soil–structure interaction studies that structural and geotechnical engineers need to know, like ground surface acceleration–time histories, PGA, PGV, settlement, and soil strain, and the site response’s frequency effect on amplification and deamplification (or attenuation). This will be done through three subsections, namely: Section 3.1. Scale of Deviations of the Ground Motion (Ground Shaking) at the Mark Point due to Soil Type and Soil Behavior, Section 3.2. Frequency Effect of the Site Response, and Section 3.3. Spatial Distribution of the Ground Motion within the Soil.

3.1. Scale of Deviations of the Ground Motion (Ground Shaking) at the Mark Point Due to Soil Type and Soil Behavior

As shown in Figure 6, the corner of the foundation cavity base is taken as a mark point to record the output at the top of the studied soil.

Figure 7 shows the changing surface velocity–time histories Vx with the types of clay soils. It can be observed that all of these velocities transmitted from the base of the clay soil to its surface differ from the TABAS earthquake that excited the base. The linear case of soil material gives higher values for the velocity–time histories at the ground surface than the nonlinear case, whose values do not differ much from TABAS values except for soft clay, for which a significant degree of nonlinearity is responsible.

Figure A2 in Appendix A presents the variation of the surface acceleration–time histories A_X with clay soil types. It can be observed that all of these accelerations, transmitted from the base of the clay soil to its surface, differ from the basic TABAS earthquake that excited the base. The linear soil material exhibits higher surface acceleration–time history values, as the nonlinear behavior of soil results in larger hysteretic energy dissipation compared to linear soil’s small Rayleigh damping. Therefore, it was neither practical nor cost-effective to construct buildings to withstand severe earthquakes by relying on the elastic behavior of the soil.

Figure A3 in Appendix A shows the variation of the surface horizontal displacement-time histories U_X with the types of clay. It can be observed that the soft clay gives the largest difference in the horizontal displacement from the U_X of TABAS, especially for the linear case of soil material. The horizontal displacement time records for nonlinear soils do not reach zero at the end due to the anticipated large plastic deformation during ground motion.

Table 2. Differences in (PGV) peak ground velocity, (PGA) peak ground acceleration, and (PGD) peak ground displacement at mark point for all studied clay soils in both linear and nonlinear cases.

	TABAS Earthquake	Stiff Clay LS	Medium Clay LS	Soft Clay LS	Stiff Clay NLS	Medium Clay NLS	Soft Clay NLS
PGV (m/s)	0.223	0.392	0.515	0.541	0.246	0.232	0.198
Diff. (%) from TABAS	0	75.74	130.88	142.64	10.2	4.26	−10.99
Diff. (%) from LS	–	–	–	–	−37.29	−54.84	−63.31
Diff. (%) from Stiff	–	–	31.38	38.07	–	−5.39	−19.23
PGA (g)	0.324	0.59	0.737	0.826	0.214	0.154	0.144
Diff. (%) from TABAS	0	82.07	127.69	155.13	−33.82	−52.41	−55.61
Diff. (%) from LS	–	–	–	–	−63.65	−79.1	−82.6
Diff. (%) from Stiff	–	–	25.05	40.13	–	−28.1	−32.92
PGD (m)	0.152	0.159	0.173	0.218	0.143	0.163	0.156
Diff. (%) from TABAS	0	4.35	14	43.29	−5.83	7.25	2.71
Diff. (%) from LS	–	–	–	–	−9.76	−5.92	−28.33
Diff. (%) from Stiff	–	–	9.24	37.31	–	13.9	9.07

shows the peak ground velocity, acceleration, and horizontal displacement values at the soil surface and their differences from the TABAS earthquake, soil with linear behavior, and stiff soil, for all soil types.

The PGV of ground surface velocity increases with decreasing shear velocity in linear soil, i.e., its value is higher for linear soft soil. Conversely, PGV decreases with the decrease of shear velocity in nonlinear soil, i.e., its value is higher for nonlinear stiff soil. This confirms that selecting the linear elastic case of soil material for analysis and design of all buildings is very conservative, far from realistic, and not economical. Additionally, the maximum difference between nonlinear and linear clay is about 63.31% for soft clay, while the maximum difference between stiff and soft clay for the linear case of material is about 38.07%. Consequently, the behavior of the soil has the most influence on PGV compared to its type.

It is also observed that the PGA of the ground-surface acceleration increases when the linear soil’s shear velocity decreases, meaning that soft clay soil has a greater PGA value. In contrast, for the nonlinear situation of soil material, PGA has the largest value for stiff clay soil. Consequently, the soil behavior has the most influence on PGA values compared to the types of clay, as the maximum difference between nonlinear and linear clay is about 82.60% for soft clay, while the maximum difference between soft and stiff clay for the linear case of material is about 40.13%. Conversely, the maximum difference between soft and stiff clay for the nonlinear case of material is only about 32.92%. Consequently, considering the nonlinearity of soil behavior leads to more economic results compared to the linear case, which is very conservative, especially for soft clay. The nonlinear stiff clay soil is more critical for determining the seismic hazard of short buildings than other types.

Furthermore, it can be found that the PGD for surface displacement does not differ significantly from the PGD value for the TABAS earthquake, except for the linear soft clay, where the difference value is about 43.29% from TABAS PGD.

Figure 8 shows the changing of surface settlement–time history at the corner of the foundation cavity base with the types of clay. It can be observed that nonlinear clay soils have higher values for settlement Uz (or vertical displacement), where this settlement is approximately small in the first three seconds of the earthquake and then begins to increase extremely with time until reaching its highest value at the earthquake end. The maximum settlement for all studied soils at the mark point and the differences from stiff soil are presented in

Table 3. The maximum settlement for all studied soils and its difference from stiff soil in both linear and nonlinear cases.

Soil Type	Max UZ (m)	Diff. (%) from Stiff
Stiff Clay LS	0.003	–
Medium Clay LS	0.006	124.67
Soft Clay LS	0.016	477.98
Stiff Clay NLS	−0.246	–
Medium Clay NLS	−0.266	8.07
Soft Clay NLS	−0.318	29.08

. It is noted that the highest value of the settlement is (0.318) m for the nonlinear soft clay, which differs about 29.08% from the nonlinear stiff clay, while this value is very negligible for linear clay and gets the highest value of 0.016 m for soft soil.

3.2. Frequency Effect of the Site Response

This section investigates the impact of bedrock motion frequency on the site response by calculating the simulated amplitude spectra recorded at the mark point on the site surface using the fast Fourier transform and then dividing them by the TABAS amplitude spectra. Figure 9 shows the log–log presentation of this ratio vs. the frequency according to velocity amplitude, while Figure A4 in Appendix A presents the log–log presentation of this ratio vs. the frequency according to acceleration amplitude.

For the Fourier spectrum amplification of the velocity, it can be observed that the amplification appears in the linear soils for the frequencies (0.3–10), (0.24–10), and (0.14–0.99) for the stiff, medium, and soft types, respectively. However, the deamplification appears for frequencies (2.2–11.4), (1.5–10.6), and (0.8–16.7) for stiff, medium, and soft linear soils, respectively, with increasing values and thresholds of deamplification as soil stiffness decreases. For nonlinear soils, these amplifications extend from (0.05–0.07, 0.2–9.3), (0.05–0.1, 0.2–7.42), and (0.05–0.07, 0.15–7.42) for stiff, medium, and soft nonlinear soils, respectively. The maximum amplification values of the Fourier amplitude ratio are 4.47, 6.89, and 11 for stiff, medium, and soft linear soil, respectively, and equal to 8, 8.22, and 6.89 for stiff, medium, and soft nonlinear soil, respectively.

However, for deamplification for the Fourier spectrum of velocity, it appears for frequencies (0.42–10.43), (0.42–10), and (0.44–14.82) for stiff, medium, and soft nonlinear soils, respectively, with increasing values and thresholds of deamplification as soil stiffness decreases. Furthermore, there is no effect of the soil behavior nor soil type near the frequency, which is equal to the first maximum peak of seismic load frequency $f_{seismicload} = 0.12 \mathrm{H}\mathrm{z}$, as the amplification ratio equals 1 in frequency ranges (0.05–0.3), (0.05–0.24), and (0.05–0.14) for stiff, medium, and soft linear soils, respectively, and in frequency ranges (0.07–0.2), (0.1–0.2), and (0.07–0.15), respectively, for stiff, medium, and soft nonlinear soils.

The Fourier amplitude ratios of acceleration’s maximum amplification values are 4.66, 6.75, and 11.87 for stiff, medium, and soft linear soil, respectively, and equal to 2.09, 8.85, and 9.29 for stiff, medium, and soft nonlinear soil, respectively.

Consequently, the amplification and deamplification of the site response according to soil type are greater in linear and nonlinear soft soil, respectively, as the amplification is greater for mid frequencies, while the deamplification is greater for longer frequencies.

Notably, the Fourier velocity spectrum amplification’s initial peak was amplified at frequencies that were near the linear soil’s fundamental frequencies. However, the first peaks for nonlinear soil occur at very short frequencies, which is similar to the results of the Fourier spectrum for acceleration in both linear and nonlinear soils. For the comparison between linear and nonlinear soils, linear soil leads to greater amplification at small and medium frequencies, and deamplification (or attenuation) is limited to large periods. While nonlinear soils tend to attenuate more at medium and large frequencies.

3.3. Spatial Distribution of the Ground Motion within the Soil

Figure 10 shows the final settlement UZ for all studied soil at the end of the TABAS earthquake (t = 21 s). These results are presented for half the section of the interior soil profile. The figures in this table confirm that the linear soils give negligible settlements distributed at the lateral sides of the soil for stiff and medium soil clay, in addition to their distribution in the middle of the soft soil, with almost non-existent settlements at the base of all linear soils. Figure A5, which is in Appendix A at the end of the paper, presents the octahedral shear strains at t = 6.7 s (corresponding to the maximum acceleration of TABAS).

It can also be observed from Figure 10 that nonlinear soils have higher settlements near the region of the foundation cavity for medium and soft clay, as their maximum values are 0.303 and 0.467 m for medium and soft clay, respectively. It was also noted that there are small upward vertical displacements at the compatible base for nonlinear soil, as its maximum value increases from 0.0103 m for soft soil to 0.092 m for stiff soil.

At t = 21 s, after the earthquake ended, Figure 11 depicts the octahedral shear strains. The octahedral shear strain values for nonlinear soils are found to be significantly higher than those for linear soils, as their maximum values are about 0.0228 and 0.00378 for nonlinear and linear soft clay, respectively. It is also noted that the octahedral shear strains increase with the depth of soil, as the distribution of these strains in the linear soils is more uniform than in the nonlinear soils. However, these strains are almost nonexistent on all the surfaces of linear soils and just on the lateral sides of the ground surface of nonlinear soils.

4. Discussion

As can be observed, this paper’s findings agree with those of earlier studies by Chala and Ray (2024) [23], Garala and Madabhushi (2019) [6], Nguyen et al. (2020) [7], Khanbabazadeh and Iyisan (2014) [21], and other researchers. According to their research, the clay layer’s stiffness and shear strength affect the clay soil’s response and its ability to amplify or attenuate. Furthermore, because of the significant amount of nonlinearity produced by the nonlinear method, equivalent linear analyses yield larger values of the peak ground accelerations at the surfaces than do nonlinear methods. Nguyen et al. (2020) [7] suggested that the attenuation in Hanoi’s soft soil Class D site, which corresponds to the medium and stiff clay soils in this research, should be taken into consideration when performing the seismic design and response of upper structures in Site Class D. This is because the NL results revealed this attenuation. The surface PGA of the EQL analysis is higher than the NL analysis through a series of one-dimensional (1D) site response assessments. Since the NL approach is considered to better capture the real reaction, the EQL approach should only be utilized if the shear stress level is less than 0.4% [51,52]. This finding aligns with the current study’s findings, which demonstrate that the stiff, medium, and soft nonlinear clays had, at the end of the earthquake, maximum values of octahedral shear strain, respectively, of 0.498%, 0.711%, and 2.28%. In the current investigation, it has been demonstrated that, in contrast to adopting linear behavior, the nonlinear dynamic analysis of clay soils, by embracing their nonlinear behavior, results in lower values of time histories of PGV, PGA, velocity, and acceleration at the ground surface. Additionally, the current research indicates that material behavior has a greater influence on the soil response than clay soil stiffness and shear strength, pointing out that stiff soil has greater PGA and PGV values than softer soils for the nonlinear behavior of clay. By contrast, soft clay exhibits higher PGA and PGV values for linear behavior. This is consistent with the results of the frequency effect on increasing the amplification of site response through increasing the ratio of Fourier amplitude for soft linear soils and increasing the attenuation (deamplification) of the same soils in the nonlinear case through decreasing this ratio, and this is consistent with the result of Che et al. (2024) [13]. These results are in good agreement with those of Garala and Madabhushi (2019) [6], who discovered that the strength and stiffness of the clay, along with the intensity of the input motion, determine whether the bedrock acceleration is amplified or attenuated as it propagates through the soft clay. Clay can increase the bedrock movement if the shear stresses generated by the propagation of the shear wave exceed its shear strength. Because of the considerable plastic deformation that is anticipated to happen during ground shaking, it is evident from the current study’s findings that the horizontal displacement–time histories for nonlinear soil do not equal zero after the record. However, except for the linear soft clay, peak ground displacement (PGD) does not significantly differ from base excitation.

In conclusion, this study’s findings highlight the critical role that clay soil type and its behavior play in determining the three-dimensional site behavior of clay soils.

5. Conclusions

An important part of earthquake engineering is site response analysis. The main focus of this research, which has enormous practical implications, is how clay soil behavior (linear/nonlinear) and type affect 3D site response analysis. The compliant base of a 3D soil model is excited by the strong TABAS earthquake to show the impact of the mentioned factors. This study investigates ground motion deviations, their frequency-dependent impact on the site response, and their spatial distribution within the soil, with the results providing the following summary:

• The site response analysis is significantly impacted by both the type of behavior (nonlinear or linear) and clay soil categorization based on the shear wave velocity, with behavior having a greater influence than type.
• For structural analysis and design, utilizing the elastic or linear case of clay soil material is extremely conservative, unrealistic, and expensive. On the other hand, soil response measures, including PGV, PGA, and acceleration records at the surface of the soil, are amplified by linear soil behavior. On the other hand, when clay soil is subjected to intense excitation, soil response is attenuated by nonlinear behavior.
• Clayey soils exhibit varying shear wave velocity classifications, with soft soils responding more to linear behavior and less to nonlinear behavior.
• Depending on soil type, site response is amplified and deamplified more in linear and nonlinear soft soils, respectively. In the Fourier amplitude spectrum, middle and short frequencies show more amplification, while longer frequencies show more attenuation. Moreover, site response is barely affected by attenuation (deamplification) and amplification at frequencies close to the applied seismic load frequency.
• The crucial state for foundation design and soil status evaluation is the nonlinear case of soil behavior. Due to high nonlinearity and extensive plastic deformation, the values of octahedral shear strain and settlement are higher in nonlinear soils than in linear soils, and their values increase with softer soils.