Influence of spatial variability of soil shear-wave velocity considering intralayer correlation on seismic response of engineering site

Shear-wave velocity uncertainty significantly impacts seismic response analysis at engineering sites. Previous studies focused on the correlations of shear-wave velocities within soil layers. However, the influence of the vertical spatial variability of shear-wave velocity within individual soil layers was not considered. This study selected a typical Site Class II engineering site to investigate this influence based on the equivalent linear seismic response analysis method of layered soil deposits. Existing field test data on shear-wave velocity were used to discretize the shear-wave velocity within soil layers into a one-dimensional (1D) random field through covariance matrix decomposition and local average process theory. Monte Carlo simulations generated random shear-wave velocity profiles with different vertical correlation distances. Three artificial records with different earthquake return periods were used as input motions at the engineering bedrock for the 1D free-field model. Considering soil shear-wave velocity uncertainty, the peak shear strain increased by 19–27% with the input ground motion amplitude. Moreover, the site peak shear strain variability increased by 25–34% with the vertical correlation distance. The proposed 1D random field model effectively simulated the interlayer correlation and spatial variability of shear-wave velocity within soil layers, demonstrating its significance in seismic response analysis.


Introduction
The seismic response analysis of engineering sites is crucial for the seismic design of engineering structures.Previous studies reported that the seismic response of engineering sites is affected by many uncertainties [1][2][3], including earthquake ground motions, shear-wave velocities of soil layers (Vs), 1330 (2024) 012027 IOP Publishing doi:10.1088/1755-1315/1330/1/012027 2 mechanical parameters of soil constitutive models, and analysis methods.Vs is closely related to the shear modulus of the soil layer and a significant parameter affecting the stiffness of the site soil layer.Therefore, quantifying the impact of Vs uncertainty on site response analysis holds practical engineering significance.
Elkateb et al. [4] reported that soil parameter uncertainty manifests in two aspects: variability of soil properties between different layers and spatial variability within the same layer.Most studies on the impact of Vs uncertainty considered Vs as a mutually independent stochastic variable.For instance, using a site in Nanjing as an example, Chen et al. [5] analyzed the impact of Vs on the shape of ground acceleration response spectrum with a coefficient of variation of 15%, considering fluctuations within one standard deviation from the mean.Chen et al. [6] took a typical layered valley site as an example and obtained the first two moments of Vs based on actual drilling data.Assuming that Vs follows a normal distribution, they employed Monte Carlo methods to analyze the impact of Vs uncertainty on the site transfer function and the amplification function of the site acceleration response spectrum.Borja et al. [7] and Guzel et al. [8] used different dynamic elastoplastic soil constitutive models for the Luodong site in Taiwan.Assuming that Vs follows a log-normal distribution, they quantified the impact of Vs uncertainty on site responses through the Monte Carlo method.
Existing research indicates that Vs uncertainty significantly affects the results of site response analysis.However, in site modeling considering soil parameter uncertainty, the variability in layer thickness and bedrock depth, as well as the correlation in Vs between adjacent soil layers, were not considered.For example, building upon Toro's model [9], which accounts for the interlayer correlation of Vs, Rathje and Kottke [10][11] developed the Strata program for site response analysis based on the equivalent linear analysis method.In addition, Hashash et al. [12] added the function of randomly generating Vs profiles and conducting frequency-domain equivalent linear and time-domain nonlinear one-dimensional (1D) site response analysis in Deepsoil based on the Toro model.Griffith et al. [13][14] and Teague et al. [15] introduced a method for inverting Vs profiles based on multiple surface wave tests.To determine the validity of current models for Vs uncertainty, they compared this approach with the method of generating random samples based on the Toro model.Passeri et al. [16] further considered the uncertainty of the cumulative propagation time of shear waves from the calculated depth to the surface based on the Toro model, effectively reducing the discreteness of random samples when considering layer thickness and Vs uncertainties.Based on more field Vs measurement data, Hallal et al. [17] modified the range of random samples, such as the cumulative propagation time of Vs and soil layer thickness, based on the improved model proposed by Passeri et al. [16].The generated soil layer thickness and Vs profile samples were more consistent with the statistical characteristics of the actual site.
Overall, most previous studies have emphasized the importance of considering Vs uncertainty in seismic site response analysis using statistical models or the surface wave inversion method [5][6][7][8][9][10][11][12][13][14][15][16][17].However, only the influence of the correlations and Vs uncertainties among the adjacent soil layers on the site response were investigated with deterministic soil layer thicknesses in the existing publications.The spatial variability of soil properties within a single soil layer can also be significant, particularly in a thick soil layer.The influence of uncertainties in Vs distribution within the soil layer remains unexplored.
To address this gap, this study improved the existing empirical statistical models of Vs profiles developed by Toro [9] by considering the uncertainties of soil properties within each soil layer.Monte Carlo simulation was performed to generate a series of spatially correlated random Vs profiles.The equivalent linear approach was used to quantify the influence of Vs variability within soil layers on the seismic site response.Three artificially synthesized ground motions of different intensity levels were used as inputs to compute the response of the random field.Finally, the influence of the vertical spatial variability of Vs on the site response was quantified.

Traditional interlayer correlation model of Vs profiles
Random Vs profiles are typically generated according to a baseline Vs profile, a measure of variability, and an interlayer correlation in Monte Carlo simulations.Toro [9] statistically investigated the Vs profiles of 557 engineering sites and indicated that the Vs of the ith soil layer, Vs(i), can be described as a nonstationary lognormally distributed random field as follows:   ,0 ln ( ) exp ln ( ) where Vs,o(i) is the baseline value for the ith soil layer, Zi represents the number of standard deviations from the mean value of ln[Vs,o(i)], and ln Vs is the standard deviation of the natural logarithm of Vs.
To account for the interlayer correlation of Vs profiles, Toro [9] developed an empirical model of Zi based on the statistical analysis of the ensemble of shear-wave velocity profile, correlated with the layer above as follows: where Zi−1 is the standard normal variable of the (i-1)th layer above it, i  is the normal random variable with zero mean and unit standard deviation, and is the interlayer correlation coefficient, a function of burial depth, d, and thickness of the soi layer, t, expressed as follows: where are the depth-and thickness-dependent correlation coefficients, respectively.
where ρ200, d0, b, ρ0, and  are model parameters related to the site class.In summary, the critical parameters required to generate randomized Vs profiles of engineering sites using Toro's model primarily include the baseline shear-wave velocity profile (Vs,o), standard deviation of the natural logarithm of Vs (lnVs), and model parameters used to estimate the interlayer correlation coefficients (ρ200, d0, b, ρ0, and  ).Generally, Vs,o profiles can be determined from the results of site-specific boreholes or in situ testing, and the other parameters can be estimated based on Table 1 for different site classes according to the average 30 m shear-wave velocity Vs30.

Improved model considering Vs distribution uncertainties in single soil layer
The traditional random field model of randomized Vs profiles, as discussed in Section 2.1, typically assumes that the Vs of an individual soil layer remains constant and uses the average Vs of the soil layer as the representative value.However, this model neglects the variability of the velocity profile within the soil layer [17].Discretizing the random fields is necessary to consider the influence of the intralayer variability of Vs on seismic site response.Common methods include the centroid, local averaging, and Karhunen-Loève expansion methods.Notably, compared with other methods, the local averaging method has the advantages of fewer discrete elements, higher computational efficiency, and better suitability for stationary random field simulations [18].Based on the interlayer correlation model introduced in the previous section, this study adopted the local averaging method to achieve Vs discretization within the soil layer, and the spatial correlation matrix was decomposed using the covariance matrix decomposition method.Figure 1 illustrates the detailed process of generating a 1D random field considering the vertical spatial variability of Vs. Figure 2 shows the discrete schematic of the individual soil layer.This study assumed that Vs is lognormally distributed at any given depth of the kth soil layer.The vertical variability of Vs can be derived as follows [19]: where VSk is the average shear velocity of the k th soil layer considering the interlayer correlation of different soil layers, zk j is the burial depth of the j th sublayer in the k th soil layer, H and H are the mean and standard deviation of the natural logarithm of the average shear velocity of the sublayers in k th soil layer, respectively, Z(zk j ) is the standard normal variable corresponding to the j th sublayer in the k th soil layer, and γH is the variance reduction coefficient, reflecting the variance attenuation of the spatial averaging characteristics of the random field.In addition, 0<γH≤1.0indicates that the overall spatial variability of the soil parameters is consistently smaller than that of local points.The matrix of the standard normal variable Z(zk j ), Z, for each soil layer can be computed using Equations ( 7) and (8).
where U is the matrix of normally distributed random variables, R represents the spatial correlation coefficient matrix, comprising the correlation coefficients of Vs among sublayers with different burial depths within a single soil layer, and L is the lower triangle matrix obtained from the Cholesky decomposition of R. Using L to linearly transform U, the independent standard normal random variables can be converted into correlated standard normal random variables.In addition, the spatial correlation coefficient, ρ(r), can be computed using the widely used exponential spatial correlation coefficient function, as shown in Equation (9).
where r is the vertical distance between two sublayers, and  is the vertical correlation distance, a critical parameter reflecting the spatial correlations of soil mechanical parameters.Parameters between two points within the correlation distance are considered to be strongly correlated; beyond this distance, the correlation is negligible.Equation (10) presents the variance reduction coefficient corresponding to the exponential spatial correlation coefficient function.
Therefore, for the improved model, two additional critical parameters (r, ), in addition to those used in Toro's model, are considered in developing randomized Vs profiles considering intralayer correlations.

Realization of randomized Vs profile
This study obtained the baseline shear-wave velocity profile (Vs,o) from the in situ test results of boreholes from an engineering site in Beijing.The engineering site primarily comprised five different soil layers: backfill soil, silty clay, fine to medium sand, fine silt, and gravelly sand.For this engineering site, the interlayer correlation coefficients (ρ200 = 0.98, d0 = 0, b = 0.344, ρ0 = 0.99, and  = 0.98) were determined based on Table 1 through linear interpolation.Toro [9] proposed lnVs values ranging from 0.27 to 0.37 for different site types.However, several studies indicated that those values generally represent the variability of Vs profiles across the entire site and may overestimate the variability for most site-specific studies.For instance, Chen et al. [20] indicated that the standard deviation of the natural logarithm of Vs,  lnVs, are generally less than 0.15 for most engineering sites with layered soil deposit.Stewart et al. [21] also reported  lnVs values of 0.15 and 0.22 for engineering sites with a thickness of the overlayer above the base rock smaller and larger than 50 m, respectively.Therefore, this study generated an ensemble of the Vs profiles for lnVs = 0.15.interlayer correlations can effectively reduce the variability of the generated Vs profiles.In the case where interlayer correlation was not considered, the two typical profile samples, identified by the red and blue curves, fluctuated around the baseline profile (dashed curve) in Figure 3(a).This indicates that the realization of Vs for each soil layer is completely independent.The generated velocity profiles featured sandwiched soft soil layers with shear-wave velocities significantly smaller than those of the adjacent soil layers, inconsistent with the overall baseline velocity profile trend.These unrealistic realizations of randomized Vs profiles result from not considering the interlayer correlations.As shown in Figure 3(b), the two typical profile samples (red and blue curves) exhibited the same increasing trend as the burial depth increased, without the presence of underlying softer soil layers, consistent with the baseline profile.Therefore, considering interlayer correlations in the simulation effectively reduces the likelihood of generating unrealistic velocity profiles with significant velocity differences among adjacent soil layers.Based on the 400 randomized realizations of the Vs profiles, Figure 4 compares the correlation coefficient matrix of Vs between different soil layers before and after considering the interlayer correlation.Figure 4(a) shows that the correlation coefficients between different soil layers were close to zero, resulting in velocity profiles with unrealistically large velocity differences among different soil layers.After adopting Toro's model with  lnVs = 0.15, the shear velocities of different soil layers exhibited positive interlayer correlations.In addition, the correlation coefficients gradually decreased with increasing vertical distances between the two soil layers.For instance, the correlation coefficient of Vs between the two adjacent soil layers, backfill soil, and silty clay layers was 0.41, whereas the correlation coefficient between the backfill soil and gravelly sand layers decreased to 0.11.[18], discretizing a single soil layer into five to seven sublayers along the vertical direction effectively captures the inherent variability of soil properties within the same soil layer.Therefore, for the specific site presented in Table 2, the thickness of the soil sublayer was set to 1 m, further dividing the soil layers with different thicknesses into five to 20 sublayers.Moreover, the exponential function, shown in Equation ( 9), was adopted for the spatial correlation coefficient matrix of Vs among different soil sublayers.Based on the statistical characterization of geotechnical variability, Phoon [22] indicated that the vertical correlation distances  of the soil properties vary from 1.0 to 3.0 m. Figure 5   This study generated 2000 Vs profiles considering vertical correlation, with correlation distances  varying from 1.0 m to 3.0 m. Figure 6 shows the 400 realizations of the random Vs profile, average profile, and corresponding standard deviations obtained through Monte Carlo simulations with a  of 1.0 m.A comparison with the traditional Toro's method revealed the Vs within the single soil layer achieved vertical spatial variability.Moreover, with the given baseline Vs profile and variation coefficients, the mean and standard deviations of the profiles generated using the proposed method are consistent with those generated using the traditional Toro's method.

Seismic input motions
Three synthetic ground motions were generated based on the target acceleration response spectra with different return periods in the I0 site according to the Chinese seismic design code (GB50909-2014) [23].Figure 7(a) shows the compatibility of the three target and synthetic acceleration response spectra.Notably, E1, E2, and E3 represent the seismic ground motions with return periods of 100, 475, and 2450 years, respectively.Figure 7(b) presents the three synthetic ground motion acceleration time histories generated using the time-domain intensity envelope method.The time-domain intensity envelope method comprised the following three steps.First, the Kaul model [24] was used to transform the target response spectrum into a power spectrum.Subsequently, a stationary random process was synthesized through the triangular series superposition method.The regulated power spectrum was transformed into nonstationary random seismic motion based on the three-segment time-domain intensity envelope function model proposed by Amin and Ang [25].Finally, synthetic ground motions that closely match the target spectrum were obtained by iteratively calculating the Fourier spectrum.

Frequency-domain linear equivalent method
One-dimensional seismic response analysis of engineering sites can be performed using different soil constitutive models through time-domain or frequency-domain analysis methods.Among these, the frequency-domain equivalent linear analysis method [26][27][28] has been widely applied owing to its simple principles and fast calculation speed.The layer properties in the equivalent linear method are adjusted through an iterative process involving a series of linear analyses.The values of shear modulus and damping ratio are used to perform iteration of linear analysis, and the iterations are continued until the shear strains from consecutive analyses match within a predefined tolerance.Therefore, this study conducted further research based on this method.The dynamic shear modulus ratio and damping ratio curves of sand and clay were adopted in the seismic site response analysis following the statistical results of Seed and Idriss [29] and Vucetic and Dobry [30]. Figure 8 shows the mean dynamic shear modulus ratio and damping ratio curves of sand and clay based on the existing dynamic experimental results on sand and clay [29,30].

Influence of Vs variability on seismic site response
This study performed 7200 frequency-domain linear equivalent analyses to investigate the effects of uncertainties in intralayer Vs profiles on the seismic response of the engineering site under different ground motion intensities.This section presents the seismic site responses in terms of peak response distributions for different vertical correlation distances.Figures 9 and 10 illustrate the effects of Vs variability on the peak shear strain and peak ground acceleration (PGA) distributions along the burial depths under three different ground motion intensities considering interlayer correlations (Toro's model).The variability of the computed site responses in terms of peak shear strain and PGA were quantified based on the standard deviations, as shown in Figures 9(d) and 10(d).The dispersion of the peak shear strain and PGA significantly increased with increasing input motion intensity.The mean values for the computed maximum peak shear strain of the given site were 0.020, 0.039, and 0.099% under the three ground motion intensities, namely E1, E2, and E3, respectively.The standard deviations for the computed maximum peak shear strain were 0.004, 0.010, and 0.027% under E1, E2, and E3, respectively.Moreover, the mean values for the computed maximum PGA were 0.141, 0.231, and 0.382g under E1, E2, and E3, respectively.The standard deviations for the computed maximum PGA were 0.025, 0.040, and 0.082g under E1, E2, and E3, respectively.
For comparison,  show the peak shear strain and PGA distributions of the site after considering the intralayer correlations of Vs with correlation distances (  ) of 1.0 and 3.0 m.With smaller correlation distances (i.e.,  = 1.0 m), the overall shear strain and PGA distributions were similar to those obtained using the traditional model considering only the interlayer correlations of Vs profiles under E1 and E2 ground motions.With an increase in the bedrock motion intensity to E3, the standard deviation for the computed maximum peak shear strain increased to 0.029%.Moreover, at larger correlation distances (i.e.,  = 3.0 m), the average remained the same.However, the standard deviations for the computed maximum peak shear strain of the target site increased to 0.005，0.012,and 0.033% under E1, E2, and E3, respectively.15 illustrates the means and standard deviations of the maximum peak shear strain and maximum PGA for different seismic intensity levels and vertical correlation distances.The mean and standard deviation of the maximum PGA remained relatively stable regardless of changes in the correlation distance.Considering the vertical spatial variability of Vs within layers, the mean of maximum peak shear strain on the site remained unchanged.However, the standard deviation increased with the vertical correlation distance, exhibiting more pronounced behavior under E3.As the vertical correlation distance increased from 1 m to 3 m, the standard deviation of maximum peak shear strain on the site increased from 0.029% to 0.033%.
Figure 16 shows the coefficient of variation for the maximum peak shear strain and maximum PGA on the site under various seismic intensity levels and correlation distances.For maximum PGA, the coefficient of variation averaged 17, 18, and 22% under E1, E2, and E3, respectively.Therefore, considering spatial variability has little impact on the PGA at the site.For the maximum peak shear strain at the site, the average coefficients of variation under E1, E2, and E3 were 22, 29, and 31%, respectively.Without considering spatial correlation, the coefficients of variation for different seismic intensity levels were 19, 25, and 27%, respectively.The uncertainty in the maximum peak shear strain at the site increased with an increase in the vertical correlation distance and the strengthening of the

Conclusions
This study quantified the influence of the vertical spatial variability of Vs within layers on site response analysis under different design ground motion intensity levels.Taking a typical layered engineering site of Site Class II in Beijing as an example, 1D random fields of soil Vs were generated using the Monte Carlo method for various vertical correlation distances.In addition, site response analysis was performed using the frequency-domain equivalent linear method.The major conclusions are as follows: (1) Building upon the existing model that accounts for interlayer Vs correlation, this study introduced a 1D random field model that considers the vertical spatial variability of Vs within soil layers.This model better reflects the statistical characteristics of actual sites and can more effectively simulate the distribution characteristics of Vs along depth.(2) Under different seismic design levels, the coefficients of variation for the maximum peak shear strain and maximum PGA on the site increased with the seismic intensity.The coefficients of variation for the maximum peak shear strain and maximum PGA averaged 22, 29, and 31% and 17, 18, and 22%, under E1, E2, and E3, respectively.(3) As the vertical correlation distance increased, the maximum PGA on the site remained unaffected.
However, the coefficient of variation for the maximum peak shear strain on the site gradually increased, reaching a maximum of 25, 32, and 34% under E1, E2, and E3, respectively.This indicated that the vertical spatial variability of Vs within soil layers significantly influenced the seismic response of the site.Several areas for future research have been identified, including finite element modeling for actual sites and investigating the influence of the number and thickness of sublayers on the results to obtain more universal findings.

Figure 1 .
Figure 1.Generation process of one-dimensional random field considering the vertical spatial variability of Vs.Figure 2. Discretization of individual soil layers.

Figure 2 .
Figure 1.Generation process of one-dimensional random field considering the vertical spatial variability of Vs.Figure 2. Discretization of individual soil layers.

Figure 3
shows examples of 400 randomized Vs profiles for the engineering site based on Toro's model with and without considering interlayer correlations.The figure indicates that considering 1330 (2024) 012027 IOP Publishing doi:10.1088/1755-1315/1330/1/0120277

Figure 3 .
Example of randomized Vs profiles based on Toro's model: (a) without interlayer correlation and (b) with interlayer correlation (lnVs = 0.15)

Figure 4 .
Correlation coefficient matrix of Vs among soil layers: (a) without interlayer correlation and (b) considering interlayer correlation following Toro [9].To further consider the intralayer variability of Vs, each individual soil layer was discretized into several sublayers.As recommended by Chen et al.
illustrates the changes in the exponential functions as the vertical correlation distance  varied from 1.0 m to 3.0 m.

Figure 6 .
One-dimensional Vs random field considering interlayer correlation: (a) 400 realizations (=1 m) and (b) comparison of mean and standard variation.

Figure 7 .
Acceleration response spectra corresponding to different return periods in the Ⅰ0 site and synthetic acceleration time histories compatible with the target spectra: (a) compatibility of target and synthetic acceleration response spectra and (b) synthetic acceleration time histories.

Figure 9 .Figure 10 .Figure 11 .Figure 12 .Figure 13 .Figure 14 .
Peak shear strain distribution considering the correlation among Vs layers: (a) under E1, (b) under E2, (c) under E3, and (d) standard deviation of peak shear strain.PGA distribution considering the correlation among Vs layers: (a) under E1, (b) under E2, (c) under E3, and (d) standard deviation of PGA.Peak shear strain distribution considering the vertical spatial variability of Vs with =1.0 m: (a) under E1, (b) under E2, (c) under E3, and (d) standard deviation of peak shear strain.PGA distribution considering the vertical spatial variability of Vs with =1.0 m: (a) under E1, (b) under E2, (c) under E3, and (d) standard deviation of PGA.1330 (2024Peak shear strain distribution considering the vertical spatial variability of Vs with =3.0 m: (a) under E1, (b) under E2, (c) under E3, and (d) standard deviation of peak shear strain.PGA distribution considering vertical spatial variability of Vs with =3.0 m: (a) under E1, (b) under E2, (c) under E3, and (d) standard deviation of PGA.The box plot shown in Figure

Figure 15 .Figure 16 .
.1088/1755-1315/1330/1/012027 14 spatial correlation.Specifically, it increased to a maximum of 25, 32, and 34% under E1, E2, and E3, respectively.Box plot of site seismic response under different ground motion intensity levels and correlation distances: (a) maximum PGA and maximum peak shear strain under E1, (b) maximum PGA and maximum peak shear strain under E2, and (c) maximum PGA and maximum peak shear strain under E3.Coefficient of variation of site seismic response under different ground motion intensity levels and correlation distances: (a) maximum peak shear strain and (b) maximum PGA.

Table 2 .
Table 2 lists the geotechnical properties of each soil layer.Vs monotonically increased with the burial depth of the soil layer.The average 30 m shear-wave velocity Vs30 was 257 m/s.Geotechnical properties of soil layers.