Seismic response of multi-support structures to spatially varying ground motions considering uncertainty in poroviscoelastic soil media

Abstract

This study investigates the seismic response of multi-support structures, such as cable-stayed bridges, subjected to spatially and randomly varying ground motions, with a particular focus on uncertainties in soil properties, namely porosity, shear modulus, and layer depth. A refined Improved Descriptive Sampling (IDS) method, validated through 500 realizations, outperforms conventional Monte Carlo and Latin Hypercube techniques in representing soil uncertainty. A refined state space formulation combined with a ground motion simulation approach, benchmarked against Wavelet Packet Transform and Expansion Optimal Linear Estimation-based methods, effectively models non-stationary excitations. The adopted poroviscoelastic soil model incorporates fluid–solid interactions, revealing substantial reductions in the soil amplification function by 52.9%, 56.3%, and 85% for porosity, shear modulus, and layer depth (Coefficient of variation, equal to 30%, 20%, and 50%) along with frequency shifts of the soil layer of 6.7%, 3.5%, and 24.5%, respectively. Ground acceleration time histories at bridge supports indicate peak reductions from 3–3.5 to 2–2.5 m/s². The bridge’s displacement and acceleration responses in the horizontal direction (Degrees of Freedom (DOFs) 1 and 2, corresponding to the tops of the towers) and in the vertical direction (DOF 3, corresponding to the center of the bridge) exhibit reduced amplitudes, broader envelopes, and smoother oscillations, reflecting enhanced damping and period lengthening. These results underscore the importance of probabilistic modeling and soil randomness representation in improving the accuracy of seismic response analyses for multi-support structures founded on porous media.

Graphical abstract

Keywords

spatial variability / probabilistic modelling / IDS / amplification / power spectral density function

1 Introduction

In earthquake response analysis, it has long been assumed that all points on the ground surface beneath a structure's foundation experience identical ground motions. This simplification treats seismic excitation as uniform across the entire foundation area. However, real earthquake records clearly show that ground motions vary both temporally and spatially [1–3]. Recognition of these variations has underscored the importance of studying spatially varying seismic ground motions and their impact on structural response, particularly for multi-support structures such as bridges. The dynamic behavior of bridges and bridge-type structures can be strongly influenced by spatial variability, and numerous studies have demonstrated its critical effects on the performance of both long and short span bridges [4–8]. Spatial variability arises from three main sources: “the wave passage effect”, results from differences in arrival times of seismic waves at-different supports of the structure; “the incoherence effect”, is caused by the scattering of seismic waves due to reflections and refractions as they propagate through heterogeneous soil layers; finally “site-response effect”, reflects differences in local soil conditions, which can lead to variations in ground motion characteristics across support points.

The concept of “spatial variability” is widely recognized as one of the main sources of soil uncertainty [9–11]. It refers to the systematic variation of soil properties across a deposit due to natural processes such as deposition, erosion, and weathering. In contrast, the complex and heterogeneous nature of soil formation also introduces “inherent (or random) variability”, which results in unpredictable property fluctuations even within small, apparently uniform areas. Quantifying both forms of uncertainty is essential for improving the realism and reliability of dynamic soil analyses and for enabling more accurate predictions of soil behavior under seismic loading. Several studies have investigated the influence of soil uncertainties on the seismic response of soil deposits. Wang and Hao [12] examined the effects of random variations in soil properties on the-site amplification, modeling both the seismic input at the bedrock and-soil properties such as shear modulus, damping ratio, density, and groundwater level as stochastic variables. Rathje et al. [13] analyzed the impact of input motion selection and soil property uncertainties on equivalent-linear site response simulations, treating shear wave velocity and nonlinear soil properties as random processes. Sadouki et al. [14] studied the response of randomly inhomogeneous layered media to harmonic excitations, modeling density and shear modulus as spatial random fields. Also, Djilali Berkane et al. [15,16] investigated how natural spatial variability in soil properties affects response spectra spatial variation across sites with different conditions., considering parameters such as layer thicknesses, shear wave velocity, and damping ratios as Gaussian random variables. Hamrouni et al. [17] showed that soil spatial variability strongly influences the deformation and failure probability of reinforced earth walls, particularly under serviceability limit states using Monte Carlo simulations (MCS) and Karhunen–Loève expansion of anisotropic random fields.

The role of uncertainties in the seismic response of multi-support systems, such as bridges, has been extensively examined. Falamarz-Sheikhabadi and Zerva [18] highlighted the importance of uncertainties in selecting incoherency coefficient and apparent propagation velocities when simulating asynchronous seismic excitations in compliance with Eurocode 8. Chen et al. [19] highlighted the significant impact of uncertainties in spatially varying ground motions (SVGMs) on the seismic response and reliability of large-span bridges. Zhang et al. [20] incorporated uncertainties in seismic magnitude and source location to model the propagation of seismic waves from the source to the bridge site. Guajardo et al. [21] proposed a numerical framework integrating spatial soil variability and three-dimensional nonlinear dynamic finite element models into the seismic performance analysis of multi-span bridges. Gurjar and Basu [22] developed a rational framework combining the conditional simulation of multicomponent SVGM field with the seismic performance assessment of a medium-span reinforced concrete highway bridge. Ba et al. [23] introduced a hybrid frequency wavenumber domain method and finite element method approach that accounts for spatial variability in seismic ground motion and simulates the full propagation process from the source to the underground structure, demonstrating that variations in velocity structure and source parameters significantly affect structural deformation and failure direction.

It is well-established that natural soil deposits are never fully solid, as-they inherently contain pore spaces between-particles. The degree of porosity, however, can vary considerably depending on-particle size, sorting, and packing. These pore spaces permit the penetration of air, water, and roots, making porosity a defining characteristic of all-natural soils. Consequently, investigating the influence of uncertainties associated with porous soil deposits is essential for accurately evaluating the seismic response of long-span structures. Yet, to date, no study in the existing literature has specifically addressed this issue. The aim of this paper is therefore to investigate the influence of such uncertainties on the response of a multi-support structure founded on porous soil deposits. To this end, ground motions are first simulated efficiently in the frequency domain and then transformed into the time domain using the Fast Fourier Transform (FFT), within the framework of the Spectral Representation Method (SRM). Non-stationary seismic accelerations are subsequently generated through an appropriate shape function. The simulated ground motions are then applied to the supports of a multi-support structure to analyze the effects of uncertainties in the surrounding medium’s properties on its response and to derive practical insights.

2 State space formulation of multi-support structures

The equation of motion for a system subjected to multi-support excitations and modeled as a multi degree of freedom system is expressed as follows [24–26]

(1)\begin{array}{rl}& \left[\begin{array}{cc}\mathit{m}& {\mathit{m}}_{\mathrm{g}}\\ {{\mathit{m}}_{\mathrm{g}}}^{\mathrm{T}}& {\mathit{m}}_{\mathrm{g}\mathrm{g}}\end{array}\right]\left\{\begin{array}{c}{\ddot{\mathit{u}}}^{\mathrm{t}}\\ {\ddot{\mathit{u}}}_{\mathrm{g}}\end{array}\right\}+\left[\begin{array}{cc}\mathit{c}& {\mathit{c}}_{\mathrm{g}}\\ {{\mathit{c}}_{\mathrm{g}}}^{\mathrm{T}}& {\mathit{c}}_{\mathrm{g}\mathrm{g}}\end{array}\right]\left\{\begin{array}{c}{\dot{\mathit{u}}}^{\mathrm{t}}\\ {\dot{\mathit{u}}}_{\mathrm{g}}\end{array}\right\}\\ +& \left[\begin{array}{cc}\mathit{k}& {\mathit{k}}_{\mathrm{g}}\\ {{\mathit{k}}_{\mathrm{g}}}^{\mathrm{T}}& {\mathit{k}}_{\mathrm{g}\mathrm{g}}\end{array}\right]\left\{\begin{array}{c}{\mathit{u}}^{\mathrm{t}}\\ {\mathit{u}}_{\mathrm{g}}\end{array}\right\}=\left\{\begin{array}{c}0\\ {\mathit{p}}_{\mathrm{g}}(t)\end{array}\right\}.\end{array}

In Eq. (1), \mathit{m} denotes the mass matrix corresponding to the superstructure (or the non-support degrees of freedom (DOFs)), {\mathit{m}}_{\mathrm{g}\mathrm{g}} represents the mass matrix associated with the support DOFs, and {\mathit{m}}_{\mathrm{g}} refers to the coupling mass matrix that accounts for the inertial forces transmitted to the superstructure DOFs as a result of support motions. The damping matrices (\mathit{c},{\mathit{c}}_{\mathrm{g}\mathrm{g}},{\mathit{c}}_{\mathrm{g}}) and the stiffness matrices (\mathit{k},{\mathit{k}}_{\mathrm{g}\mathrm{g}},{\mathit{k}}_{\mathrm{g}}) are partitioned in a similar manner and have analogous definitions. The vector {\mathit{u}}^{\mathrm{t}} contains the total displacements of the superstructure DOFs, whereas u_{\mathrm{g}}represents the ground displacement vector at the support DOFs. Likewise, {\dot{\mathit{u}}}^{\mathrm{t}}, {\dot{\mathit{u}}}_{\mathrm{g}}, {\ddot{\mathit{u}}}^{\mathrm{t}}, and {\ddot{\mathit{u}}}_{\mathrm{g}} denote the corresponding velocity and acceleration vectors for the superstructure and support DOFs, respectively. Finally, {\mathit{p}}_{\mathrm{g}}(t) stands for the vector of forces developed at the support DOFs.

The total displacements can be decomposed into two components, as expressed below

(2)\left\{\begin{array}{c}{\mathit{u}}^{\mathrm{t}}(t)\\ {\mathit{u}}_{\mathrm{g}}(t)\end{array}\right\}=\left\{\begin{array}{c}{\mathit{u}}^{\mathrm{s}}(t)\\ {\mathit{u}}_{\mathrm{g}}(t)\end{array}\right\}+\left\{\begin{array}{c}\mathit{u}(t)\\ 0\end{array}\right\},

where {\mathit{u}}^{\mathrm{s}} denotes the vector of structural displacements produced by the static application of the prescribed support displacements {\mathit{u}}_{\mathrm{g}} at each time instant (commonly referred to as the quasi-static displacement), and \mathit{u}\left(t\right) is the vector of dynamic displacements such that [25,26]

(3)\mathit{k}{\mathit{u}}^{\mathrm{s}}+{\mathit{k}}_{\mathrm{g}}{\mathit{u}}_{\mathrm{g}}=0,

or, in a simplified form,

(4){\mathit{u}}^{\mathrm{s}}=\iota {\mathit{u}}_{\mathrm{g}},

with

(5)\iota {=-{\mathit{k}}^{-1}\mathit{k}}_{\mathrm{g}}.

The first of the two partitioned equations given in Eq. (1) can be expressed as

(6)\mathit{m}{\ddot{\mathit{u}}}^{\mathrm{t}}+{\mathit{m}}_{\mathrm{g}}{\ddot{\mathit{u}}}_{\mathrm{g}}+\mathit{c}{\dot{\mathit{u}}}^{\mathrm{t}}+{\mathit{c}}_{\mathrm{g}}{\dot{\mathit{u}}}_{\mathrm{g}}+\mathit{k}{\mathit{u}}^{\mathrm{t}}+{\mathit{k}}_{\mathrm{g}}{\mathit{u}}_{\mathrm{g}}=0.

By applying Eq. (2), Eq. (6) can be rewritten as follows

(7)\begin{array}{rl}\mathit{m}\ddot{\mathit{u}}+\mathit{c}\dot{\mathit{u}}+\mathit{k}\mathit{u}=& -\left(\mathit{m}{\ddot{\mathit{u}}}^{\mathrm{s}}+{\mathit{m}}_{\mathrm{g}}{\ddot{\mathit{u}}}_{\mathrm{g}}\right)-\left(\mathit{c}{\dot{\mathit{u}}}^{\mathrm{s}}+{\mathit{c}}_{\mathrm{g}}{\dot{\mathit{u}}}_{\mathrm{g}}\right)\\ & -\left(\mathit{k}{\mathit{u}}^{\mathrm{s}}+{\mathit{k}}_{\mathrm{g}}{\mathit{u}}_{\mathrm{g}}\right).\end{array}

Owing to Eq. (3), the last term in Eq. (7) vanishes. Using Eqs. (4) and (5), Eq. (7) becomes

(8)\mathit{m}\ddot{\mathit{u}}+\mathit{c}\dot{\mathit{u}}+\mathit{k}\mathit{u}=-\left(\mathit{m}\iota +{\mathit{m}}_{\mathrm{g}}\right){\ddot{\mathit{u}}}_{\mathrm{g}}-\left(\mathit{c}\iota +{\mathit{c}}_{\mathrm{g}}\right){\dot{\mathit{u}}}_{\mathrm{g}}.

Moreover, the damping term in Eq. (8) is typically negligible compared with the inertia term. On other hand, for structures modeled with a lumped-mass idealisation, the mass matrix is diagonal, resulting in a zero m_{\mathrm{g}} matrix. Under these assumptions, Eq. (8) reduces to the simplified form

(9)\mathit{m}\ddot{\mathit{u}}+\mathit{c}\dot{\mathit{u}}+\mathit{k}\mathit{u}=-\mathit{m}\iota {\ddot{\mathit{u}}}_{\mathrm{g}}.

The matrix \iota is an influence coefficient matrix of size (n\times r) which describes the displacements resulting from the unit support displacements, where n is the number of super structure non-supports DOFs and r is the number of components of the input ground motion. In other words, it quantifies the influence of support displacements on the structural displacements.

Using the state space method, the system response is evaluated by considering both displacement and velocity as independent variables referred as “states”. The two response variables are expressed as a state vector, \mathit{z} which can be written as

(10)\mathit{z}\left(t\right)=\left\{\begin{array}{c}\mathit{u}(t)\\ \dot{\mathit{u}}\left(t\right)\end{array}\right\}.

The state space velocity can be expressed as

(11)\dot{\mathit{z}}\left(t\right)=\mathit{A}\mathit{z}\left(t\right)+\mathit{f}\left(t\right),

where \mathit{A} is a matrix and \mathit{f}\left(t\right) a vector expressed as

(12)\begin{array}{c}\mathit{A}=\left[\begin{array}{cc}0& \mathit{I}\\ -{\mathit{m}}^{-1}\mathit{k}& -{\mathit{m}}^{-1}\mathit{c}\end{array}\right],\mathit{f}\left(t\right)=\left\{\begin{array}{c}0\\ -\iota {\ddot{\mathit{u}}}_{\mathrm{g}}\left(t\right)\end{array}\right\}\end{array}.

3 Ground motion simulation

3.1 Method for ground motion simulation

Seismic ground motions often vary spatially due to wave passage effects, incoherence, and site-specific conditions. Such variations are particularly important for multi-support structures, as they may induce differential motions between supports. In this subsection, a formulation is presented to simulate spatially and randomly varying ground motions, accounting for both temporal non-stationarity and spatial correlation. Figure 1 illustrates the physical configuration used to simulate spatially and randomly varying ground motions at the surface of a poroviscoelastic medium.

The model consists of n_{\mathrm{s}} surface stations, each underlain by a soil column of height h_i, resting on a horizontally extended bedrock. Vertically propagating incident SH waves excite the bedrock at multiple locations (1′,2′,...,n_{\mathrm{s}}^{\prime }), and the corresponding surface ground motions are simulated at stations 1,2,...,n_{\mathrm{s}}. The distances (d_{ij}) between stations exhibit spatial variability, while random variation in the soil layer properties capture the heterogeneity of the porous medium. This configuration forms the formulation, used to simulate spatially coherent, non-stationary ground motions across the surface profile.

Ground motions can be simulated in time domain, but this approach is computationally demanding. However, it has been shown that their simulation is more efficient in the frequency domain, where they can be generated using the FFT technique as follows [16,27,28]

(13)Y_j\left(\mathrm{i}{\omega }_n\right)=\sum _{m=1}^j{B}_{jm}({\omega }_n)\left[\cos {\alpha }_{jm}\left({\omega }_n\right)+\mathrm{i}\mathrm{s}\mathrm{i}\mathrm{n}{\alpha }_{jm}\left({\omega }_n\right)\right],

where B_{jm}\left({\omega }_n\right) and {\alpha }_{jm}\left({\omega }_n\right) are the amplitude and phase of the simulated ground motions at frequency {\omega }_n, respectively, and defined as

(14a)\begin{array}{c}{B}_{jm}\left({\omega }_n\right)=\sqrt{\mathrm{\Delta }\omega }\left|{\mathit{L}}_{jm}\left(\mathrm{i}{\omega }_n\right)\right|,0\le {\omega }_n\le {\omega }_N,\end{array}

(14b)\begin{array}{c}{\alpha }_{jm}\left({\omega }_n\right)={\beta }_{jm}\left({\omega }_n\right)+{\phi }_{mn}\left({\omega }_n\right),0\le {\omega }_n\le {\omega }_N,\end{array}

with

(14c)\begin{array}{c}{\beta }_{jm}\left({\omega }_n\right)={\tan }^{-1}\left(\frac{\mathrm{I}\mathrm{m}\left[L_{jm}\left(\mathrm{i}{\omega }_n\right)\right]}{\mathrm{R}\mathrm{e}\left[L_{jm}\left(\mathrm{i}{\omega }_n\right)\right]}\right),0\le {\omega }_n\le {\omega }_N,\end{array}

where ω is the excitation frequency, and {\phi }_{mn}\left({\omega }_n\right) is the random phase angles uniformly distributed over the range \left[0,2\pi \right]. N is the total number of discrete frequencies ({\omega }_n=i_{\mathrm{f}}\mathrm{\Delta }\omega, where i_{\mathrm{f}} is the frequency number sampling such that i_{\mathrm{f}}=1,N_{\mathrm{f}}, and \mathrm{\Delta }\omega ={\displaystyle \frac{{\omega }_{\mathrm{u}}}{N_{\mathrm{f}}}} with {\omega }_{\mathrm{u}} is an upper cut-off frequency).

The time histories y_j(t) are then obtained by applying the inverse FFT of Y_j(\mathrm{i}\omega ). Since earthquake ground motions are inherently nonstationary, they begin at zero, intensify to a steady phase, and eventually-decay back to zero. The nonstationary nature of seismic accelerations must be incorporated into analysis. Nonstationary seismic accelerations can be derived using a shape function, as described below [16,29]

(15){{\ddot{u}}_j\left(t\right)=E\left(t\right)y}_j(t).

The shape function E\left(t\right) used in this study is defined as

(16)

In Eq. (16), T denotes the total motion duration, while, t_0 and t_n correspond to specific time intervals within this duration. The coefficient c is an attenuation parameter assumed equal to 0.155.

In Eqs. (14a)–(14c), \mathit{L}\left(\mathrm{i}\omega \right) is the lower triangular matrix obtained by the Cholesky decomposition of the ground surface cross-power spectral density (PSD) matrix {\mathit{S}}_{jk}^{\mathrm{S}\mathrm{u}\mathrm{r}\mathrm{f}\mathrm{a}\mathrm{c}\mathrm{e}}\left(\mathrm{i}\omega \right) defined by Eq. (17)

(17){\mathit{S}}_{jk}^{\mathrm{S}\mathrm{u}\mathrm{r}\mathrm{f}\mathrm{a}\mathrm{c}\mathrm{e}}\left(\mathrm{i}\omega \right)={\mathit{H}}_j\left(\mathrm{i}\omega \right){\mathit{H}}_k^{\ast }\left(\mathrm{i}\omega \right){\mathit{S}}_{j^{\prime }{k}^{\prime }}^{\mathrm{b}\mathrm{e}\mathrm{d}\mathrm{r}\mathrm{o}\mathrm{c}\mathrm{k}}\left(\mathrm{i}\omega \right),j,k,j^{\prime },k^{\prime }=1,n_{\mathrm{s}}.

According to the following equation

(18){\mathit{S}}_{jk}^{\mathrm{S}\mathrm{u}\mathrm{r}\mathrm{f}\mathrm{a}\mathrm{c}\mathrm{e}}\left(\mathrm{i}\omega \right)=\mathit{L}(\mathrm{i}\omega ){\mathit{L}}^{\mathrm{H}}(\mathrm{i}\omega ),

where {\mathit{L}}^{\mathrm{H}}(\mathrm{i}\omega ) is the Hermitian matrix. The ground surface cross-PSD matrix (Eq. (17)) is formulated from the cross-PSD function between the bedrock spatial points (j^{\prime } and k^{\prime }) and the corresponding ground surface points (j and k) in Fig. 1 [16]. The cross-PSD function of ground motions at n_{\mathrm{s}} bedrock stations is expressed as

(19){\mathit{S}}_{j^{\prime }{k}^{\prime }}^{\mathrm{b}\mathrm{e}\mathrm{d}\mathrm{r}\mathrm{o}\mathrm{c}\mathrm{k}}\left(\mathrm{i}\omega \right)={\mathit{S}}_{\mathrm{g}}\left(\omega \right){\gamma }_{j^{\prime }{k}^{\prime }}\left(d_{j^{\prime }{k}^{\prime }},\mathrm{i}\omega \right),j^{\prime },k^{\prime }=1,n_{\mathrm{s}},j^{\prime }\ne k^{\prime }.

In Eq. (19), {\gamma }_{j^{\prime }{k}^{\prime }}(d_{j^{\prime }{k}^{\prime }},\mathrm{i}\omega ) is the spatial coherence function at the bedrock. For relatively long distances d_{j^{\prime }{k}^{\prime }}, the Harichandran and Vanmarcke [1] spatial coherence model may be more efficient. This model takes the following equation

(20)\begin{array}{rl}{\gamma }_{j^{\prime }{k}^{\prime }}\left(d_{j^{\prime }{k}^{\prime }},\mathrm{i}\omega \right)=& C\exp \left[\frac{-2d_{jk}}{a\theta \left(\omega \right)}\left(1-C+aC\right)\right]\\ & +\left(1-C\right)\exp \left[\frac{-2d_{jk}}{\theta \left(\omega \right)}\left(1-C+aC\right)\right],\end{array}

with \theta \left(\omega \right)is given by the following equation

(21)\theta \left(\omega \right)=k_0{\left[1+{\left(\frac{\omega }{2\pi f_0}\right)}^b\right]}^{-\frac{1}{2}},

where C, a, k_0, f_0, and b are the model parameters. In this study the values obtained by Harichandran and Vanmarcke [1] are used (C = 0.736, a = 0147, k_0 = 5210, f_0 = 1.09 Hz, and b = 2.78). In the same equation (Eq. (19)), d_{j^{\prime }{k}^{\prime }} is the distance between two support points j’ and k’.

In the definition of the ground surface cross-PSD function according to Eq. (17), {\mathit{H}}_j(\mathrm{i}\omega ) is the transfer (or amplification) function between the bedrock and surface for the station j while {\mathit{H}}_k^{\ast }\left(\mathrm{i}\omega \right) is the conjugate of the amplification function (AF) at station k. The bedrock may be subjected to different kind of excitations expressed in terms of the density function S_{\mathrm{g}}\left(\omega \right). Usually, the bedrock is assumed to be excited by the white noise filtered model suggested by Clough and Penzien [24]

(22)S_{\mathrm{g}}\left(\omega \right)=\frac{{\omega }^4}{{\left({\omega }_{\mathrm{f}}^2-{\omega }^2\right)}^2+{\left(2{\xi }_{\mathrm{f}}{\omega }_{\mathrm{f}}\omega \right)}^2}\times \frac{{\omega }_{\mathrm{g}}^4+{\left(2{\xi }_{\mathrm{g}}{\omega }_{\mathrm{g}}\omega \right)}^2}{{\left({\omega }_{\mathrm{g}}^2-{\omega }^2\right)}^2+{\left(2{\xi }_{\mathrm{g}}{\omega }_{\mathrm{g}}\omega \right)}^2}{S}_0,

where S_0 is a scale factor depending on the ground motion intensity, {\omega }_{\mathrm{g}} and {\xi }_{\mathrm{g}} are the circular frequency and damping ratio of the Clough and Penzien [24] filter, respectively, and {\omega }_{\mathrm{f}} and {\xi }_{\mathrm{f}} are those of the second filter in Eq. (22). For this study, the following values are used: {\omega }_{\mathrm{g}}=6\mathrm{p}, {\xi }_{\mathrm{g}}=0.2, {\omega }_{\mathrm{f}}=0.16\mathrm{p}, {\xi }_{\mathrm{f}}=0.6, S_0=0.0034{\mathrm{m}}^2/{\mathrm{s}}^3.

In soil dynamics and earthquake engineering ground motion amplification can be effectively assessed using the transfer (or amplification) function of a soil profile. This function is often modeled by a poroviscoelastic soil column model of depth h_j (Fig. 1), referred as Biot’s column, that characterizes the biphasic nature of soil, with both the solid skeleton and pore fluid contributing to its dynamic response. The coupling between these phases, manifested through fluid-solid interactions and pore pressure, fundamentally influences the seismic behavior of the soil [30–32]. For vertically propagating SH waves, the complex transfer function {\mathit{H}}_j\left(\mathrm{i}\omega \right) represents the ratio between the surface displacement amplitude at (point j) and the amplitude at the soil-bedrock interface (point j′) [15,16,33].

When a poroviscoelastic soil column overlies a poroviscoelastic half-space and is subjected at its base to vertically propagating harmonic SH waves (Fig. 1), the total out-of-plane displacement and shear stress in each layer are [32]

(23)\begin{array}{rl}{u}_y^j\left({\mathrm{z}}_j,t\right)=& \left[B_{\mathrm{S}\mathrm{H}}^j\exp \left(\mathrm{i}\frac{\omega }{v_{\mathrm{S},j}^{\ast }}{\mathrm{z}}_j\right)+B_{\mathrm{S}\mathrm{H}}^{\prime j}\exp \left(-\mathrm{i}\frac{\omega }{v_{\mathrm{S},j}^{\ast }}{\mathrm{z}}_j\right)\right]\exp \left(\mathrm{i}\omega t\right),\\ & (j=l,h),\end{array}

(24)\begin{array}{rl}{\sigma }_{yz}^j\left({\mathrm{z}}_j,t\right)=& \frac{\mathrm{i}\omega }{v_{\mathrm{S},j}^{\ast }}{G}_j^{\ast }\left[B_{\mathrm{S}\mathrm{H}}^j\exp \left(\mathrm{i}\frac{\omega }{v_{\mathrm{S},j}^{\ast }}{\mathrm{z}}_j\right)-B_{\mathrm{S}\mathrm{H}}^{\prime j}\exp \left(-\mathrm{i}\frac{\omega }{v_{\mathrm{S},j}^{\ast }}{\mathrm{z}}_j\right)\right]\\ & \cdot \exp \left(\mathrm{i}\omega t\right),(j=l,h).\end{array}

In these equations, the superscript \mathrm{l} corresponds to the soil layer (column) while superscript \mathrm{h} refers to the half-space. Here B_{\mathrm{S}\mathrm{H}}^j and B_{\mathrm{S}\mathrm{H}}^{\prime j} (for j = l, h) represent the amplitudes of the incident and reflected SH-wave in the layer and the half-space.

Applying the relevant boundary conditions, specifically: 1) zero shear stress at the free surface ({\sigma }_{y\mathrm{z}}^l\left(0,t\right)=0) leads to B_{\mathrm{S}\mathrm{H}}^{\prime \mathrm{l}}=B_{\mathrm{S}\mathrm{H}}^{\mathrm{l}}, and 2) continuity of displacements and stresses at the layer half-space interface () allows expressing the amplitudes of the incident and reflected waves in the half-space in terms of those in the layer [32].

According to the above definition and using Eq. (23), the AF due to SH wave propagating vertically from the bedrock (half-space) is

(25)H\left(\mathrm{i}\omega \right)=\frac{1}{\exp \left(\mathrm{i}{\displaystyle \frac{\omega }{v_{\mathrm{S},\mathrm{l}}^{\ast }}}{h}_{\mathrm{l}}\right)+\exp \left(-\mathrm{i}{\displaystyle \frac{\omega }{v_{\mathrm{S},\mathrm{l}}^{\ast }}}{h}_{\mathrm{l}}\right)},

where h_{\mathrm{l}} is the soil layer thickness (also referred to as the depth), and v_{\mathrm{S},\mathrm{l}}^{\ast } represents the complex shear wave velocity in the layer considering the viscous skeleton damping (ξ), such that

(26)v_{\mathrm{S},\mathrm{l}}^{\ast }=v_{\mathrm{S}}\sqrt{1+\mathrm{i}2\xi }.

Several researchers derived different expressions for the shear wave velocity for elastic porous media. As an example, but not exclusively, Sadouki et al. [30] gave the following expression

(27)v_{\mathrm{S}}=\sqrt{\frac{G\left({\rho }_{22}-\mathrm{i}{\displaystyle \frac{b}{\omega }}\right)}{\left({\rho }_{11}{\rho }_{22}-{{\rho }_{12}}^2\right)-\mathrm{i}{\displaystyle \frac{b}{\omega }}\left({\rho }_{11}+{\rho }_{22}+{\rho }_{12}\right)}},

where G is the solid skeleton elastic shear modulus in the layer, and {\rho }_{11} and {\rho }_{12} are the mass coefficients that account for the inertia effects of the moving fluid. These coefficients are associated with the densities of the solid ({\rho }_{\mathrm{s}}) and the fluid \left({\rho }_{\mathrm{f}}\right) through the following mixtures relationships {\rho }_{11}=\left(1-\varphi \right){\rho }_{\mathrm{s}}-{\rho }_{12}, {\rho }_{22}=\varphi {\rho }_{\varphi }-{\rho }_{12}, {\rho }_{12}=\varphi \left(1-a\right){\rho }_{\mathrm{f}} with a={\displaystyle \frac{1}{2}}\left(1-{\displaystyle \frac{1}{\varphi }}\right) is the dynamic tortuosity. The parameter b denotes the dissipation coefficient, defined as b={\varphi }^2{\displaystyle \frac{\eta }{\kappa }} [34,35], where \varphi is the porosity, \eta is the fluid viscosity, and \kappa is the coefficient of permeability (in m/s). Al Rjoub [36] derived a similar equation as Eq. (27) but defined the dissipation coefficient as b={\varphi }^2{\displaystyle \frac{\eta }{\kappa }} where k is the intrinsic permeability (in m²). There is a relationship between intrinsic permeability (k in m²) and coefficient of permeability (\kappa in m/s) [32,37] expressed as k={\displaystyle \frac{\kappa .\eta }{{\rho }_{\mathrm{f}}\cdot g}} with g = 9.81 m/s².

In the case of monophasic soil layer, a same equation as Eq. (27) is valid by replacing Eq. (27) by

(28)v_{\mathrm{S}}=\sqrt{\frac{G}{{\rho }_{\mathrm{s}}}}.

3.2 Validation of the simulation method

To evaluate the accuracy and performance of the proposed simulation approach, two comparative studies are conducted against methods reported in recent literature. The first benchmark compares the present results with those obtained using the Wavelet Packet Transform (WPT)-based method introduced by Ji et al. [38], which enhances the classical SRM by incorporating time–frequency localization. The second comparison focuses on the Expansion Optimal Linear Estimation (EOLE)-based procedure developed by Peng et al. [39], which introduces a virtual continuous process and achieves computational efficiency by requiring only a single decomposition of the expanded correlation matrix. Both studies are performed using the same input excitation and simulation settings, ensuring a consistent evaluation of the proposed method’s capability and accuracy to reproduce non-stationary ground motions.

3.2.1 Comparison with wavelet packet transform-based simulation method

As an illustrative example, the results of the proposed simulation are compared with those of Ji et al. [38]. The 1979 Imperial Valley earthquake acceleration time history is used as the seed ground motion, with its Pseudo-Spectral Acceleration (PSA) serving as the target spectrum. Ji et al. [38] employed the WPT algorithm to enhance the SRM for simulating spatially varying non-stationary ground motions while preserving the characteristics of the seed recording. In this study, the bedrock excitation is represented by the PSD function proposed by Kaul [40], which corresponds to the target PSA, rather than by the filtered white noise model (Eq. (22)) suggested by Clough and Penzien [24]

(29)S_{\mathrm{g}}\left(\omega \right)=\frac{-{\displaystyle \frac{\xi }{\pi \omega }}{PSA}^2\left(\omega \right)}{\mathrm{l}\mathrm{n}\left(-{\displaystyle \frac{\pi }{\omega T}}\ln p\right)},

where T is the time duration, \xi is the damping ratio, and p is the probability coefficient, taken as 0.85 following Kaul [40]. Equation (29) is incorporated into the cross-PSD function (Eq. (19)) for simulating ground motions at three stations as described in Subsection 3.1. The stations are located 200 m apart along a straight line. The seed recording and its pseudo response spectrum (target PSA) are shown in Figs. 2(a) and 2(b), respectively. To introduce non-stationarity into the simulated stationary ground motions, the envelope function (Eq. (16)) is applied with parameters t = 37.05 s, t₀ = 5.44 s, t_n = 10.76 s. The resulting non stationary ground motions are compared with those of Ji et al. [38] at the same stations (Figs. 2(c), 2(e), and 2(g)), along with their corresponding normalized energy build-up (Figs. 2(d), 2(f), and 2(h)). The simulated non-stationary ground motions exhibit strong visual similarity to both the seed recording and Ji et al. results [38], capturing comparable temporal non-stationary characteristics. Although exact coincide is not expected due to the inherent randomness of the simulation method, their overall behavior remains consistent and statistically representative. In addition, the normalized energy build-up of the simulated samples (Figs. 2(d), 2(f), and 2(h)) is nearly identical to that reported by Ji et al. [38], indicating a similar trend in energy accumulation.

3.2.2 Comparison with the Expansion Optimal Linear Estimation-based simulation method

In their study, Peng et al. [39] proposed an efficient approach for simulating multivariate non-stationary ground motions, based on a virtual continuous process combined with the EOLE technique. As an application example, they simulated ground motions at nine support points distributed along an eight-span bridge model. For each location, 10 samples were generated, and representative acceleration time histories were obtained using the Clough–Penzien spectrum (Eq. (22)) to characterize the stationary ground motions. Three typical soil conditions were used to compute the parameters of the PSD as shown in Table 1 [41–43].

The scale factor S_0 value is 200 cm/s². The time-dependent modulation function applied to the ground motion-follows Eq. (16). The model of Harichandran and Vanmarcke [1] (Eq. (20)) was employed to describe the spatial coherence of ground motions, incorporating the wave-passage effect. Typical seismic ground motion samples at points 1, 5, and 8, exhibiting characteristic features of non-stationarity in the time domain, are compared in Fig. 3 with those generated by the simulation method of Peng et al. [39], demonstrating excellent agreement between the two methods. While minor differences are observed, due to the inherent randomness of the simulation process, the generated motions remain statistically consistent and successfully reproduce the expected temporal patterns. Similarly to the first comparative study, the simulated samples exhibit a normalized energy build-up consistent with the findings of Peng et al. [39], suggesting that both follow a similar trend in energy accumulation.

4 Sampling methods for soil parameter uncertainty

Due to the inherent variability of soil properties, which are often addressed through stochastic or probabilistic approaches, their randomness cannot be neglected. A wide range of methods has been developed to account for this variability. In general, a deterministic problem is solved repeatedly using a large number of simulated random variables and the statistical properties of the resulting responses are then evaluated [30,44,45]. Commonly used methods include Monte Carlo Simulation (MCS), Latin Hypercube Sampling (LHS), Importance Sampling (IS), Rejection Sampling (RS), Improved Descriptive Sampling (IDS) among others [46–48]. Each method offers distinct strengths and limitations, making them suitable for different problems contexts. MCS is often considered the standard approach due to its simplicity and flexibility, though it can be computationally demanding. LHS provides a more efficient alternative, particularly when computational resources are limited. IS is highly effective for capturing rare events but requires careful implementation to avoid bias. RS, while flexible, can be inefficient and is generally applied in specialized scenarios where other methods may be unsuitable. IDS, an enhancement of traditional descriptive sampling methods, is designed to better capture the variability of input parameters in stochastic and probabilistic analyses. It combines elements of stratified sampling and LHS with additional measures to ensure that the sampled data more accurately reflect the statistical properties (e.g., mean, variance, skewness.) of the input distributions. This makes it particularly effective when dealing with non-normal distributions or cases where higher moments play a significant role in the output. Compared with simpler methods, IDS offers a more accurate representation of input parameters leading to more reliable results, especially in problems where accuracy is critical. However, its increased complexity and computational cost may restrict its use to cases where such rigor justified. In fact, as highlighted in the literature, soil properties frequently deviate from normality and are better described by lognormal or skewed distributions rather than Gaussian ones [9,49]. Monte Carlo and Latin Hypercube methods primarily reproduce mean and variance, but may fail to capture skewness and kurtosis [50,51], which are critical in geotechnical variability. The IDS method is particularly suitable in this context, as it preserves higher-order statistical moments and provides a more representative sampling of non-normal distributions with fewer realizations [51]. This makes IDS especially advantageous for probabilistic modeling of soil uncertainties.

Ultimately, the choice of method depends on the specific characteristics of the problem under consideration, including-model complexity, the nature of soil variability, and the available computational resources.

Several of these methods will be applied to a simple problem involving the calculation of the mean and standard deviation of the AF for a single-layer soil profile.

In this example, the soil layer thickness (depth, h_{\mathrm{l}}) and the shear wave velocity (v_{\mathrm{S}}) are treated as random variables, each following a lognormal distribution. By modeling these key parameters as random variables, the inherent variability and uncertainty in soil properties can be systematically incorporated into the analysis. The damping ratio (ξ) is considered deterministic and is fixed at 5%. The AF, which quantifies the amplification of seismic waves as they travel through the soil layer, is highly sensitive to these input parameters. Consequently, evaluating the mean and standard deviation of the AF is crucial for reliable seismic risks assessment.

For this simplified case, the AF is defined as a particular instance of Eq. (25), where the shear wave velocity is given by Eq. (28) and the complex shear modulus as {v_{\mathrm{S}}}^{\ast }=\sqrt{{\displaystyle \frac{G\left(1+2\mathrm{i}\xi \right)}{{\rho }_{\mathrm{s}}}}}.

It is important to note that the number of random parameters included in the computational model significantly affects both the complexity and duration of the calculations. As the number of random parameters rises, the computational effort and analysis time may increase considerably. Engineers should therefore aim to minimize the number of random parameters by identifying and prioritizing those variables that have the greatest influence on results. The soil layer thickness (h_{\mathrm{l}}) and shear wave velocity (v_{\mathrm{S}}) are both lognormally distributed.

The accuracy of the estimated statistics for the AF is highly dependent on the number of realizations performed for each parameter set. Estimates based on a small number of realizations are prone to large standard errors. To ensure reliable and accurate AF statistics, a sufficiently high number of realizations is often required. Moreover, the required number of realizations to achieve stable outcomes depends on the specific stochastic sampling method employed.

Figure 4 illustrates how the mean and standard deviation of the AF are affected by the number of realizations, for each of the three considered methods. In this study, MCS, LHS, and IDS are used to compare the efficiency of these methods in producing accurate statistics for a soil layer characterized by a given thickness h₁ = 30 m and a shear wave velocity v_{\mathrm{S}} = 300 m/s, each with a coefficient of variation (C_v) of 20%. The results demonstrate that increasing the number of realizations lead to a reduced dispersion in both the mean and standard deviation of the AF. When the number of realizations is low (for example, 10, as shown in Figs. 4(a) and 4(b)), all three stochastic sampling methods exhibit large standard errors in their statistical estimates. However, with 500 samples, the standard deviation is effectively eliminated for the IDS method (Fig. 4(f)), and with 1000 samples for the LHS method (Fig. 4(h)). In contrast, MCS requires more than 1000 samples to reach comparable accuracy. Based on this analysis, the IDS method was selected for all subsequent analyses, with the number of realizations set at 500 samples.

5 Multi-support structural response

5.1 Configuration of the studied cable-stayed bridge

In this Subsection, a simplified model of cable-stayed bridge is considered (Fig. 5), and its seismic behavior is analyzed under ground motions generated by the proposed method. The bridge consists of a main span of 250 m and two side spans of 80 m each. It is supported by two towers, each 50 m high, with the deck located 30 m above the ground level [25]. The towers are considered inextensible, and the geometric nonlinearity of the pre-tensioned cable members is neglected when deriving the stiffness matrix. The following relationships are assumed:

(EI{)}_{\mathrm{t}\mathrm{o}\mathrm{w}\mathrm{e}\mathrm{r}}=0.25(EI{)}_{\mathrm{d}\mathrm{e}\mathrm{c}\mathrm{k}}=0.25EI, {\displaystyle \frac{AE}{l_1}}={\displaystyle \frac{12EI}{s^3}}, {\left({\displaystyle \frac{AE}{l_{\mathrm{d}\mathrm{e}\mathrm{c}\mathrm{k}}}}\right)}_{\mathrm{d}\mathrm{e}\mathrm{c}\mathrm{k}}=0.8{\left({\displaystyle \frac{AE}{l_1}}\right)}_{\mathrm{c}\mathrm{a}\mathrm{b}\mathrm{l}\mathrm{e}}, {\displaystyle \frac{12EI}{s^{3}m}}=20{\left({\displaystyle \frac{\mathrm{r}\mathrm{a}\mathrm{d}}{s}}\right)}^2, with s=125.

In these expressions, {\left(EI\right)}_{\mathrm{d}\mathrm{e}\mathrm{c}\mathrm{k}} represents the bending stiffness (flexural stiffness) of the bridge deck, (EI{)}_{\mathrm{t}\mathrm{o}\mathrm{w}\mathrm{e}\mathrm{r}} represents the bending stiffness of the towers. The term AE refers to the axial stiffness of the stay cables, defined as the product of the cross-sectional area and the Young’s modulus of the cable material. l_{\mathrm{d}\mathrm{e}\mathrm{c}\mathrm{k}} denotes the total length of the bridge deck, l_1is the length of the cable connecting the top of the tower to the end of the bridge, and l_2 is the length of the cable joining the tower top to the center of the bridge.

The parameters m_1 and m_2 correspond to the masses lumped at the tops of the towers, whereas m_3 represents the mass lumped at the center of the bridge. The masses are related to a reference mass m by the expressions m_1=m_2=20\mathrm{m} and m_3=60\mathrm{m}.

5.2 Procedure for obtaining the bridge response

The bridge response is obtained by solving Eq. (9) using the state-space method, according to the following steps.

1) Generation of the influence coefficient matrix \mathit{\iota }, as described in the Appendix A in Electronic Supplementary Material (Eq. (A4)).

2) Generation of the mass submatrix corresponding to the superstructure DOFs (Eq. (A5)).

3) Calculation of the eigenvalues and natural frequencies.

4) Generation of the Rayleigh damping matrix \mathit{c} by assuming 5% critical damping (ξ) for all modes.

5) Simulation of the ground acceleration vector {\ddot{\mathit{u}}}_{\mathrm{g}} corresponding to the support DOFs, considering the effects of several parameters (Subsection 5.3).

6) Calculation of the \mathit{A} matrix and \mathit{f}(t) vector (Eq. (12)).

7) Calculation of the state vector \mathit{z} (Eq. (10)) at each time step.

8) Calculation of the derivative of the state vector \dot{\mathit{z}} (Eq. (11)) at each time step, which represents the state-space form of Eq. (9) in terms of the dynamic displacement of the mass.

These steps are implemented in a MATLAB code to obtain the bridge response, as in Subsection 5.4. The simulation of the ground acceleration vector {\ddot{\mathit{u}}}_{\mathrm{g}} corresponding to the support DOFs is carried following the algorithm illustrated in Fig. 6.

5.3 Effects of soil parameter uncertainty on surface ground motions

This section examines the influence of uncertainties in key parameters of the soil environment, associated with spatial variability, on the simulated ground motions at the bridge supports (Table 2). As summarized in Table 2, uncertain parameters include the depth of the sand layer, the shear modulus of the solid skeleton, and the porosity, while other parameters such as grain density, fluid density, permeability, viscosity, and Poisson’s ratio were treated as deterministic values. Before addressing this topic, the effects of these uncertainties on both the soil AF and the PSD are investigated. Figures 7 and 8 illustrate the impact of uncertainties in porous medium properties, specifically porosity and shear modulus, as well as the depth of the soil layer beneath each bridge support on soil amplification. The uncertain parameters are modeled as lognormally distributed random variables, characterized by their mean values and C_v.

Figures 7(a) and 7(b) illustrate the effect of porosity uncertainty on the AF at the bridge supports. As the C_v for porosity increases from 0% (deterministic case) to 30%, the mean amplitude of the AF at the fundamental frequency (1.42 Hz) decreases significantly. Specifically, at C_v = 10%, the AF amplitude is reduced by approximately 19.65% across all supports, with supports 1 and 4 (typically located at the bridge’s outer spans) showing slightly more pronounced reductions than supports 2 and 3 (closer to the towers). At C_v = 30%, the AF amplitude decreases by 52.9% across all supports. This reduction can be attributed to the damping effect of fluid–solid interactions in the porous medium, which attenuates the amplification of seismic energy. In fact, the present Biot-based formulation reveals a frequency-dependent behavior. Specifically, when porosity uncertainty is introduced, the shear wave velocity tends to increase, and consequently, the soil amplification magnitude decreases [52,53]. This behavior arises from the dynamic coupling between the pore fluid and the solid matrix described by Biot’s theory: at relatively higher frequencies, fluid flow within the pores becomes restricted, effectively stiffening the composite medium. However, porosity uncertainty may render the soil either softer or stiffer, that is, associated with a greater or lower fundamental period, respectively, depending on the specific variations in pore structure and distribution [53].

Additionally, a 6.7% leftward shift in the fundamental frequency is observed, indicating a transition to longer-period modes. This shift is more pronounced at supports 1 and 4, where the increasing dispersion of soil layer depth uncertainty, i.e., its spatial variability of soil stratigraphy (with increasing C_v), enhances the damping effect. In contrast, supports 2 and 3, located closer to the stiffer tower foundations, exhibit slightly less pronounced frequency shifts.

Figures 7(c) and 7(d) examine the impact of shear modulus uncertainty. As the C_v of shear modulus increases from 0% to 20%, the AF amplitude at the fundamental frequency decreases by 15.1% at C_v = 5%, 35.7% at C_v = 10% and by 56.3% at C_v = 20% across all supports. Furthermore, a 3.5% shift toward higher fundamental frequency is observed, suggesting a transition to shorter-period dynamic responses. This effect is characteristic of stiffer soils and is associated with reduced seismic amplification due to diminished impedance contrasts and energy-trapping effects [12,14,52], captured by the perturbation associated with increased uncertainty in the shear modulus.

Figures 7(e) and 7(f) address the effect of layer depth uncertainty. As the C_v for layer depth increases from 0% (deterministic case) to 10%, 20%, 30%, and 50%, the AF amplitude at the fundamental frequency decreases. For supports 1 and 4, an increase to C_v = 10%, corresponds to an approximate 55.4% reduction in AF amplitude, with the reduction becoming more pronounced at higher C_v values. At C_v = 50%, the AF amplitude-undergoes an extreme reduction of up to 85%, exhibiting a near-linear downward trend that approximates a straight-line decrease. This linear decay reflects a significant collapse in amplification caused by severe geometric variability within the soil profile, accompanied by a 24.6% leftward shift on the fundamental frequency, indicating a toward longer-period modes. Supports 2 and 3, which are located closer to the stiffer tower foundations, exhibit slightly less pronounced reductions at each C_v level but nevertheless, follow a similar near-linear decrease at C_v = 50%. These supports also experience notable frequency shifts and a reduction in high-frequency content, although mitigated somewhat by their proximity to the towers.

In summary, an increase in uncertainties around layer depth generally results in a reduction in seismic amplification. This behavior is attributed to enhanced scattering and attenuation of seismic waves due to irregular boundaries and increased energy dissipation within the saturated, viscous pore structure. As uncertainties increase, the constructive interference responsible for resonance becomes weaker, leading to reduced surface motion amplification. A similar trend was reported by Guellil et al. [44], who demonstrated that spatial variability in layer depth significantly affect the response of soi-foundation-structure system.

Figure 8 illustrates the influence of soil parameter uncertainties, porosity, shear modulus, and layer depth, on the PSD function, with distinct effects observed at supports 1 and 4. The PSD is shaped by both the AF and the filtered white noise model (Eq. (22)), which governs the spectral characteristics of ground motions.

Figures 8(a) and 8(b) show that increasing porosity uncertainty produces a marked reduction in peak PSD amplitude. Supports 1 and 4 experience more pronounced declines compared to supports 2 and 3, reflecting their higher sensitivity to fluid-solid interactions in porous media, as modulated by the AF and Eq. (22), which dissipate spectral energy. Supports 2 and 3, in contrast, retain slightly more spectral energy, showing milder reductions.

Figures 8(c) and 8(d) indicate that uncertainty in the shear modulus similarly reduces peak PSD amplitude, with supports 1 and 4 exhibiting greater attenuation. This effect arises from decreased soil stiffness, represented by the perturbation introduced by shear modulus variability, which enhances energy dissipation through the combined influence of the AF (Figs. 7(c) and 7(d)) and the filtered white noise-model (Eq. (22)). Supports 2 and 3 display milder reductions.

Figures 8(e) and 8(f) demonstrate that layer depth uncertainty results in the most dramatic decrease in peak PSD amplitude. Supports 1 and 4 undergo substantial attenuation, driven by geometric variability in the soil profile, which modifies both the AF and Eq. (22), leading to more diffuse spectral responses. Supports 2 and 3 experience smaller reductions, maintaining marginally higher spectral energy.

The consistent attenuation of PSD amplitudes and the greater sensitivity of supports 1 and 4, compared with supports 2 and 3, underscore the role of spatial variability in ground motions driven by inherent soil variability (soil heterogeneity). These findings, derived and validated using the IDS method with 500 realizations, highlight the critical influence of soil parameter uncertainties on the seismic response of multi-support structures. They emphasize the importance of probabilistic modeling and detailed geotechnical characterization to ensure reliable seismic performance assessments of bridges, particularly at supports founded on variable soil conditions.

The influence of uncertainties in the same key soil parameters on the simulated ground motions at the bridge supports is further illustrated in Figs. 9–11. These figures present the effects of random variations in porosity (Fig. 9), shear modulus (Fig. 10), and layer depth (Fig. 11) on the acceleration time histories simulated at the base of the cable-stayed bridge.

Figure 9 illustrates the influence of porosity variability, modeled as a random variable with C_v = 0%, 10%, 20%, and 30%, on the simulated acceleration time histories at the four-bridge supports. Increasing porosity variability leads to a systematic reduction in ground acceleration amplitudes across all supports. In the deterministic case (C_v = 0%), acceleration peaks range between 3.0 and 3.5 m/s². As C_v increases to 30%, these peaks attenuate significantly, dropping to around 2.0 to 2.5 m/s². This reduction in amplitude is accompanied by two notable dynamic features: 1) smoother acceleration time histories; where the simulated-signals become more regular and less oscillatory as porosity randomness increases. This behavior reflects a damping-like effect induced by pore distributions, which diffuse seismic energy and reduce sharp dynamic transitions; and 2) broader, less distinct peaks, indicating that the ground motion envelopes widen with increasing porosity uncertainty, resulting in less sharply defined acceleration peaks. This suggests a general loss of coherence in wave propagation [54], exacerbated by heterogeneous pore structures in porous media where fluid–solid interactions govern wave behavior.

Overall, the results highlight that greater porosity variability in the soil profile leads to notable reductions in ground motion intensity [30]. This attenuation reduces seismic inputs to the structure and confirms the necessity of incorporating porosity uncertainty into site response analysis to achieve more realistic seismic performance assessments.

Figure 10 illustrates the effect of shear modulus uncertainty on the simulated ground accelerations at the four support points of the bridge. The shear modulus is modeled as a lognormally distributed random variable with increasing C_v (C_v = 0%, 5%, 10%, and 20%). The results clearly indicate that greater variability in the shear modulus, arising from material uncertainties, leads to a significant reduction in ground motion amplitudes. This attenuation is primarily attributed to the combined effects of reduced AF and diminished PSD amplitudes.

In the deterministic case (C_v = 0%), the acceleration time histories exhibit sharp peaks, with maximum values between 3.0 and 3.5 m/s². As C_v increases to 10% and 20%, these peaks are progressively reduced, falling within the range of 2.0 to 2.5 m/s². The observed decline in amplitude is attributed to perturbation induced by soil stiffness uncertainty, which promotes wave scattering and energy dissipation, thereby attenuating the seismic waves reaching the surface. With increasing shear modulus variability, the acceleration signals become smoother and more regular, displaying broader waveforms and fewer abrupt transitions, indicative of damping-like behavior. This leads to reduced amplification and delayed oscillatory components in the ground motion, consistent with dispersive wave propagation in non-homogeneous media. These findings align with El Haber et al. [55], who demonstrated through two-dimensional stochastic modeling that spatial variability in shear-wave velocity induces wave scattering, coherence loss, and delayed seismic response, manifested as reduced amplification and phase-shifted oscillatory components consistent with dispersive propagation in heterogeneous media.

In summary, Fig. 10 shows that uncertainty in the shear modulus significantly affects both the spatial and temporal characteristics of ground motion. It results in lower peak accelerations, smoother and phase-shifted waveforms, and overall disruption of wave coherence due to local stiffness variability. These effects enhance energy dissipation across the soil profile and alter the seismic demand on structures. Such findings highlight the importance of incorporating shear modulus variability into probabilistic site response analyses to ensure realistic and reliable assessments of structural performance under SVGMs.

Figure 11 presents the simulated ground acceleration time histories at the four support points of the cable-stayed bridge, highlighting the influence of uncertainty in the soil layer depth. The variability in soil layer depth is modeled as a lognormal distribution with C_v ranging from 0% to 50%. Among the various sources of soil uncertainty considered in this study, randomness in layer thickness proves to have the most pronounced impact on the ground motion characteristics.

In the absence of uncertainty (C_v = 0%), the simulated accelerations exhibit well-defined, high-amplitude peaks ranging between 3.0 and 3.5 m/s². However, as the depth variability increases, these peaks are progressively dampened, and the overall amplitude of the acceleration signals is significantly reduced. For a C_v of 30%, the peak accelerations decrease to approximately 2.0−2.5 m/s², while at C_v = 50%, they fall even further, approaching 1.8 m/s² in some cases. This strong attenuation is primarily caused by the influence of depth variation on the travel time, impedance contrast, and resonance conditions of seismic waves propagating through the soil, as the thickness of the soft surface layer strongly affects both the frequency and amplitude of resonance peaks [56].

Beyond the amplitude reduction, the time histories exhibit a marked evolution in waveform shape. As depth uncertainty increases, the acceleration signals become broader and smoother, with sharp oscillatory components, typically associated with resonant ground motion, progressively diminished [16]. This transformation indicates a shift in the system’s dynamic response from well-defined resonance behavior toward a more diffuse and damped regime, driven by randomness in propagation path lengths and interference effects. This observation is consistent with foundational. 0site effect studies, which have shown that soil layering, depth variability, and near-surface wave interference significantly reshape seismic ground motion, reducing high-frequency content and producing prolonged, attenuated response characteristics [57].

Additionally, increased depth variability leads to greater dispersion and phase distortion across the support points. Subtle temporal shifts in the arrival time and peak timing of ground motion are observed, especially at higher C_v levels, indicating that wave passage effects and incoherence are amplified under uncertain soil conditions. Such phenomena may induce significant asynchronous excitation in long-span structures, amplifying internal force demands and deformation patterns. This observation is consistent with Konakli and Der Kiureghian [54], who noted that while low-frequency components of ground motion are typically assumed coherent under uniform soil conditions, this assumption breaks down when soil variability is present, requiring additional considerations to account for phase differences and temporal misalignments.

Finally, Figure 11 clearly demonstrates that layer depth uncertainty exerts a dominant influence on the amplitude, frequency content, and temporal coherence of seismic ground motions. It not only leads to considerable attenuation of acceleration peaks but also fundamentally alters waveforms characteristics, making them more dispersed and irregular. These findings are consistent with those of Roy and Jakka [58], who showed that combined uncertainties in shear wave velocity and layer thickness significantly increase the variability of amplification spectra, response spectra, and peak ground acceleration. This highlights the necessity of accurately characterizing stratigraphic variability in seismic hazard assessments, especially for infrastructure systems spanning multiple supports, where local soil conditions can vary significantly and drive complex dynamic interactions.

5.4 Effects of soil parameter uncertainty on bridge response

The effects of soil parameter uncertainties on bridge displacements and accelerations are shown in Figs. 12–17. Figures 12 and 13 present the impact of shear modulus variability on displacement and acceleration time histories at three DOFs: with DOFs 1 and 2-correspond to horizontal displacements in the direction of seismic excitation and DOF 3 represents vertical displacement.

In Fig. 12, displacement time histories show a clear reduction in peak amplitudes as the C_v of shear modulus increases from 0% to 20%. Horizontal DOFs 1 and 2 decrease from ±(0.07–0.08) m to ±(0.04–0.05) m, while DOF 3-shows a smaller reduction due to lower vertical dynamic demand. This reduction in amplitude, together with a broader response envelope and smoother oscillations, indicates increased damping within the porous medium, consistent with the AF reductions of 15.1%, 35.7%, and 56.3% reported by Figs. 7(c) and 7(d) at C_v values of 5%, 10%, and 20%.

Similarly, Figure 13 shows horizontal accelerations decreasing from ±(4–5) m/s² to ±(2–3) m/s², while vertical accelerations at DOF 3 are less affected. The smoother behavior and loss of high-frequency content, combined with a 3.5% rightward shift in fundamental frequency (Fig. 7(c)), indicate transitions to shorter-period modes caused by stiffness randomness. These findings are consistent with those of Tran et al. [59], who found that probabilistic soil modeling yields significant variations in spectral response, particularly smoothing and de-amplification effects.

Figures 14 and 15 demonstrate the impact of porosity uncertainty, with peak displacements and accelerations at horizontal (DOFs 1 and 2) decreasing from ±(0.07–0.08) to ±(0.04–0.05) m and ±(4–5) to ±(2–2.5) m/s², respectively, as C_v increases from 0% to 30%. Vertical (DOF 3) responses are less sensitive.

These reductions, accompanied by broader envelopes and smoother oscillations, are linked to fluid–solid interactions in porous media, reflected by a 6.7% leftward frequency shift (Fig. 7(a)). The dominance of horizontal DOFs underscores the influence of horizontal excitation, while variability across supports highlights site-specific amplification effects. These findings are consistent with previous studies [34,60,61] that showed porosity and fluid–solid coupling strongly damp seismic energy, shift dominant frequencies, and reduce peak responses, justifying stochastic modeling in seismic bridge design.

Figures 16 and 17 illustrate the effects of soil layer depth variability on the displacement and acceleration, evaluated at the same three DOFs. Figure 16 shows that as the C_v for the soil layer depth increases from 0% to 50%, the peak horizontal displacements at DOFs 1 and 2 progressively decrease from approximately ±(0.07–0.08) to ±(0.04–0.05) m. In contrast, the vertical displacements at DOF 3 exhibit a less marked reduction, reflecting the predominance of horizontal seismic input and the inherently lower vertical dynamic demand.

The response envelopes become increasingly broader and the oscillation patterns smoother with rising C_v, indicating a shift toward longer-period modes. This trend is consistent with the 24.6% leftward shift in the fundamental frequency observed at C_v = 50% in Fig. 7(e), which corresponds to an approximately 85% reduction in the AF amplitude. Figure 17 presents the corresponding acceleration time histories. Similar to the displacement results, peak horizontal accelerations at DOFs 1 and 2 decline from ±(4–5) m/s² in the deterministic case to ±(2–2.5) m/s² at C_v = 50%.

Vertical accelerations at DOF 3 also decrease but to a lesser extent, in line with the generally lower vertical energy transmission. The observed reduction in acceleration amplitudes, along with the broader envelopes and diminished high-frequency content, reflects the increased damping and altered resonance conditions resulting from geometric uncertainty in the porous medium, governed by Eq. (25). The greater sensitivity of the horizontal DOFs and the spatial variability across support points underline the critical role that layer depth uncertainty plays in site response and structural performance. These findings reinforce the importance of adopting probabilistic modeling approaches and ensuring accurate subsurface characterization, particularly in the seismic analysis of bridges founded on porous soils with uncertain parameters.

This interpretation is consistent with previous numerical investigations, including those by Allam and Datta [62], and Tonyalı et al. [63]. The former demonstrated that SVGMs, accounting for incoherence, wave-passage, and site-response effects, can significantly modify the seismic demands on cable-stayed bridges, especially in the horizontal direction. The latter further confirmed that variations in soil stiffness and wave velocities across bridge supports result in amplified responses, highlighting the necessity of incorporating spatial variability and soil variability into seismic design frameworks. Neglecting these factors may lead to unconservative predictions of structural demand, especially in long-span, multi-support systems.

The reliability of these results is further supported by the use of the IDS method with 500 realizations, which ensures an effective representation of the statistical variability in the response parameters.

6 Conclusions

This study advances the seismic analysis of multi-support structures by introducing and validating improved methodologies, including an enhanced IDS method, a robust state space formulation, and a refined simulation approach for SVGMs. The optimized IDS method, using 500 realizations, outperforms Monte Carlo and LHS by efficiently capturing soil property variability and providing stable statistical estimates of the AF with reduced computational cost. The proposed simulation approach, benchmarked against WPT and EOLE-based methods, reliably reproduces non-stationary ground motions.

The results demonstrate that uncertainties in porosity, shear modulus, and layer depth of a viscoelastic porous medium supporting a multi-support structure, significantly reduce the AF and PSD, reflecting a softening of the soil-structure system. Acceleration time histories at bridge supports (Figs. 9–11) show peak values decreasing from 3–3.5 to 2–2.5 m/s². Displacement and acceleration responses at horizontal (DOFs 1 and 2) and vertical (DOF 3) DOFs (Figs. 12–17) also exhibit reduced amplitudes (displacements from ±(0.07–0.08) to ±(0.04–0.05) m; accelerations from ±(4–5) to ±(2–3) m/s²), broader response envelopes, and smoother oscillations. These patterns indicate enhanced damping, longer-period modes, and site-specific amplification effects driven by fluid–solid interactions and random soil variability, with horizontal DOFs and supports showing pronounced sensitivity.

Overall, these methodological and physical insights underscore the necessity of probabilistic modeling and detailed geotechnical characterization to capture the complex dynamics of natural soil environments. From an earthquake engineering perspective, incorporating poroviscoelastic soil uncertainties into design is essential for developing resilient infrastructure, ensuring reliable performance-based seismic design, and safeguarding multi-support structures in seismically active regions with complex subsurface conditions.

The findings of this study have important implications for engineering practice and seismic design. The demonstrated reductions in peak accelerations and the broadening of response envelopes indicate that neglecting soil randomness may lead to conservative estimates of seismic demand, with direct consequences for the calibration of seismic design spectra and reduction factors prescribed in current codes such as Eurocode 8 and ASCE 7. Moreover, the smoother oscillations and enhanced damping observed in the structural response highlight the role of poroviscoelastic soil deposits in increasing effective damping, thereby suggesting the need for more refined damping provisions in design guidelines. It should be noted that soil-structure interaction effects were not included in this study; however, they may influence the actual response of bridges and could be investigated in future research. By explicitly quantifying the influence of uncertainties in porosity, shear modulus, and layer depth, the study supports the integration of probabilistic approaches into performance-based and reliability-based design frameworks. This is particularly relevant for long-span bridges where differential support motions govern seismic vulnerability, and resilience considerations are essential for ensuring both safety and serviceability under earthquake loading.

References

[1] Harichandran R S , Vanmarcke E H. . Stochastic variation of earthquake ground motion in space and time. Journal of Engineering Mechanics, 1986, 112(2): 154–174

[2] Zerva A Spatial Variation of Seismic Ground Motions: Modeling and Engineering Applications. 1st ed. Boca Raton: CRC Press, 2009

[3] Zerva A , Zervas V. . Spatial variation of seismic ground motions: An overview. Applied Mechanics Reviews, 2002, 55(3): 271–297

[4] Kiureghian A D , Neuenhofer A. . Response spectrum method for multiple support seismic excitations. Earthquake Engineering & Structural Dynamics, 1992, 21(8): 713–740

[5] Hao H , Chouw N. . Response of a RC bridge in WA to simulated spatially varying seismic ground motions. Australian Journal of Structural Engineering, 2008, 8(1): 85–98

[6] Zhang Y , Li Q , Lin J , Williams F. . Random vibration analysis of long-span structures subjected to spatially varying ground motions. Soil Dynamics and Earthquake Engineering, 2009, 29(4): 620–629

[7] Liu Z , Zhang S. . Influence of strong spatially varying near fault ground motion on steel box arch bridge. Bulletin of Earthquake Engineering, 2021, 19(13): 5439–5469

[8] Belkheiri N , Tiliouine B. . Quantification of the effects of the spatial variation of ground motions on the seismic response of highway bridges. Bridge Structures, 2023, 19(1–2): 3–14

[9] Phoon K K , Kulhawy F H. . Characterization of geotechnical variability. Canadian Geotechnical Journal, 1999, 36(4): 612–624

[10] Elkateb T , Chalaturnyk R , Robertson P K. . An overview of soil heterogeneity: Quantification and implications on geotechnical field problems. Canadian Geotechnical Journal, 2003, 40(1): 1–15

[11] Popescu R , Deodatis G , Nobahar A. . Effects of random heterogeneity of soil properties on bearing capacity. Probabilistic Engineering Mechanics, 2005, 20(4): 324–341

[12] Wang S , Hao H. . Effects of random variations of soil properties on site amplification of seismic ground motions. Soil Dynamics and Earthquake Engineering, 2002, 22(7): 551–564

[13] Rathje E M , Kottke A R , Trent W L. . Influence of input motion and site property variabilities on seismic site response analysis. Journal of Geotechnical and Geoenvironmental Engineering, 2010, 136(4): 607–619

[14] Sadouki A , Harichane Z , Chehat A. . Response of a randomly inhomogeneous layered media to harmonic excitations. Soil Dynamics and Earthquake Engineering, 2012, 36: 84–95

[15] Djilali Berkane H , Harichane Z , Guellil M E , Sadouki A. . Investigation of soil layers stochasticity effects on the spatially varying seismic response spectra. Indian Geotechnical Journal, 2019, 49(2): 151–160

[16] Djilali Berkane H , Harichane Z , Çelebi E , Elachachi S M. . Site dependent and spatially varying response spectra. Earthquake Engineering and Engineering Vibration, 2019, 18(3): 497–509

[17] Hamrouni A , Dias D , Sbartai B. . Soil spatial variability impact on the behavior of a reinforced earth wall. Frontiers of Structural and Civil Engineering, 2020, 14(2): 518–531

[18] Falamarz-Sheikhabadi M R , Zerva A. . Two uncertainties in simulating spatially varying seismic ground motions: incoherency coefficient and apparent propagation velocity. Bulletin of Earthquake Engineering, 2018, 16(10): 4427–4441

[19] Chen Z , Zheng S , Ding Z , Zhang J , Tai Y. . Seismic reliability evaluation of bridges under spatially varying ground motions using a four-parameter distribution. Engineering Structures, 2021, 247: 113157

[20] Zhang J , Wei K , Gao L , Qin S. . Effect of V-shape canyon topography on seismic response of deep-water rigid-frame bridge based on simulated ground motions. Structures, 2021, 33: 1077–1095

[21] Guajardo B , Pinto F , Astroza R. . Effects of soil spatial variability on the seismic response of multi-span simply-supported highway bridges. Bulletin of Earthquake Engineering, 2024, 22(5): 2643–2675

[22] Gurjar N , Basu D. . A framework to integrate conditional simulation of multicomponent spatially varying ground motion field with seismic performance assessment and its application to medium-span bridges. Soil Dynamics and Earthquake Engineering, 2025, 191: 109238

[23] Ba Z , Fu J , Zhu Z , Zhong H. . An end-to-end 3d seismic simulation of underground structures due to point dislocation source by using an FK-FEM hybrid approach. Frontiers of Structural and Civil Engineering, 2022, 16(12): 1515–1529

[24] Clough R WPenzien J Dynamics of Structures. 2nd ed. New York: McGraw Hill, 1995

[25] Datta T K Seismic Analysis of Structures. Singapore: John Wiley & Sons, 2010

[26] Chopra A K Dynamic of Structures, Theory and Applications to Earthquake Engineering. 5th ed. Pearson Education Limited, United Kingdom, 2020.

[27] Hao H , Oliveira C S , Penzien J. . Multiple-station ground motion processing and simulation based on SMART-1 Array Data. Nuclear Engineering and Design, 1989, 111(3): 293–310

[28] Bi K , Hao H. . Modelling and simulation of spatially varying earthquake ground motions at sites with varying conditions. Probabilistic Engineering Mechanics, 2012, 29: 92–104

[29] Zanardo G , Hao H , Modena C. . Seismic response of multi-span simply supported bridges to a spatially varying earthquake ground motion. Earthquake Engineering & Structural Dynamics, 2002, 31(6): 1325–1345

[30] Sadouki A , Harichane Z , Elachachi S M , Erken A. . Response of anisotropic porous layered media with uncertain soil parameters to shear body and love-waves. Earthquakes and Structures, 2018, 14(4): 313–322

[31] Sadouki A , Harichane Z , Erken A. . Effects of soil parameters uncertainties on the behaviour of anisotropic porous media to shear waves. International Journal of Geotechnical Earthquake Engineering, 2019, 10(2): 32–49

[32] Rahmani Kouadri R , Harichane Z , Bouaricha L. . Amplification of seismic wave motions in partially fluid-filled poroviscoelastic media. Soil Dynamics and Earthquake Engineering, 2022, 160: 107382

[33] Harichane Z , Afra H , Elachachi S M. . An identification procedure of soil profile characteristics from two free field accelerometer records. Soil Dynamics and Earthquake Engineering, 2005, 25(6): 431–438

[34] Biot M A. . Generalized theory of acoustic propagation in porous dissipative media. Journal of the Acoustical Society of America, 1962, 34(9A): 1254–1264

[35] Sharma M D. . Rayleigh waves in dissipative poro-viscoelastic media. Bulletin of the Seismological Society of America, 2012, 102(6): 2468–2483

[36] Al Rjoub Y S. . The reflection of P-waves in a poroelastic half-space saturated with viscous fluid. Soil Dynamics and Earthquake Engineering, 2013, 49: 218–230

[37] Bardet J P. . A viscoelastic model for the dynamic behavior of saturated poroelastic soils. Journal of Applied Mechanics, 1992, 59(1): 128–135

[38] Ji K , Cao X , Wang S , Wen R. . Efficient simulation of spatially correlated nonstationary ground motions by wavelet-packet algorithm and spectral representation method. Earthquake Engineering and Engineering Vibration, 2024, 23(4): 799–814

[39] Peng B , Xu J , Peng Y. . Efficient simulation of multivariate non-stationary ground motions based on a virtual continuous process and EOLE. Mechanical Systems and Signal Processing, 2023, 184: 109722

[40] Kaul M J. . Stochastic characterization of earthquakes through their response spectrum. Earthquake Engineering & Structural Dynamics, 1978, 6(5): 497–509

[41] Seya H , Talbott M E , Hwang H H. . Probabilistic seismic analysis of a steel frame structure. Probabilistic Engineering Mechanics, 1993, 8(2): 127–136

[42] Deodatis G. . Nonstationary stochastic vector processes: Seismic ground motion applications. Probabilistic Engineering Mechanics, 1996, 11(3): 149–167

[43] Liu Z , Liu Z , Ruan X , Zhang Q. . Spectral representation-based dimension reduction for simulating multivariate nonstationary ground motions. Soil Dynamics and Earthquake Engineering, 2018, 114: 313–325

[44] Guellil M E , Harichane Z , Berkane H D , Sadouki A. . Soil and structure uncertainty effects on the soil foundation structure dynamic response. Earthquakes and Structures, 2017, 12(2): 153–163

[45] Harichane Z , Guellil M E , Gadouri H. . Benefits of probabilistic soil-foundation-structure interaction analysis. International Journal of Geotechnical Earthquake Engineering, 2018, 9(1): 42–64

[46] McKay M D , Beckman R J , Conover W J. . A comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics, 1979, 21(2): 239–245

[47] Denny M. . Introduction to importance sampling in rare-event simulations. European Journal of Physics, 2001, 22(4): 403–411

[48] Tari M , Dahmani A. . Refined descriptive sampling: A better approach to Monte Carlo simulation. Simulation Modelling Practice and Theory, 2006, 14(2): 143–160

[49] Fenton G AGriffiths D V Risk Assessment in Geotechnical Engineering. Singapore: John Wiley & Sons, 2008

[50] Helton J C , Davis F J. . Latin hypercube sampling and the propagation of uncertainty in analyses of complex systems. Reliability Engineering & System Safety, 2003, 81(1): 23–69

[51] Hoshino N , Takemura A. . On reduction of finite-sample variance by extended Latin hypercube sampling. Bernoulli, 2000, 6(6): 1035–1050

[52] Poggi V , Edwards B , Fäh D. . Derivation of a reference shear-wave velocity model from empirical site amplification. Bulletin of the Seismological Society of America, 2011, 101(1): 258–274

[53] Sun M , Liu P , Chen Y , Bate B , Xue F. . Sand stiffness variability induced by stochastic distributions of calcite precipitates: A Monte Carlo-DEM study. Acta Geotechnica, 2025, 20(3): 1363–1377

[54] Konakli K , Der Kiureghian A. . Simulation of spatially varying ground motions including incoherence, wave-passage and differential site-response effects. Earthquake Engineering & Structural Dynamics, 2012, 41(3): 495–513

[55] El Haber E , Cornou C , Jongmans D , Lopez-Caballero F , Youssef Abdelmassih D , Al-Bittar T. . Impact of spatial variability of shear wave velocity on the lagged coherency of synthetic surface ground motions. Soil Dynamics and Earthquake Engineering, 2021, 145: 106689

[56] Stanko D , Markušić S. . An empirical relationship between resonance frequency, bedrock depth and VS30 for Croatia based on HVSR forward modelling. Natural Hazards, 2020, 103(3): 3715–3743

[57] Kawase H. . Site effects on strong ground motions. International Geophysics, 2003, 81(B): 1013–1030

[58] Roy N , Jakka R S. . A study on the effect of shear wave velocity (VS) and layer thickness uncertainty on seismic site response analysis considering statistical randomness. Journal of Earth System Science, 2025, 134(3): 133

[59] Tran T T , Salman K , Han S R , Kim D. . Probabilistic models for uncertainty quantification of soil properties on site response analysis. ASCE-ASME Journal of Risk and Uncertainty in Engineering Systems, Part A: Civil Engineering, 2020, 6(3): 04020030

[60] Zienkiewicz O C , Shiomi T. . Dynamic behaviour of saturated porous media; the generalised Biot formulation and its numerical solution. International Journal for Numerical and Analytical Methods in Geomechanics, 1984, 8(1): 71–96

[61] Shen YZhong ZLi LDu X Fluid-solid fully coupled seismic response analysis of layered liquefiable site with consideration of soil dynamic nonlinearity. In: Wang L, Zhang JM, Wang R, eds. Wang L, Zhang JM, Wang R, eds. PBD-IV Geotechnical, Geological and Earthquake Engineering, Vol 52. Cham: Springer, 2022, 1708–1716

[62] Allam S M , Datta T K. . Seismic behaviour of cable-stayed bridges under multi-component random ground motion. Engineering Structures, 1999, 21(1): 62–74

[63] Tonyali Z , Ates S , Adanur S. . Spatially variable effects on seismic response of the cable-stayed bridges considering local soil site conditions. Structural Engineering and Mechanics, 2019, 70(2): 143–152

RIGHTS & PERMISSIONS

Higher Education Press