Stochastic nonlinear ground response analysis considering existing boreholes locations by the geostatistical method

The field and laboratory evidence of nonlinear soil behaviour, even at small strains, emphasizes the importance of employing nonlinear methods in seismic ground response analysis. Additionally, determination of dynamic characteristics of soil layers always includes some degree of uncertainty. Most of previous stochastic studies of ground response analysis have focused only on variability of soil parameters, and the effect of soil sample location has been mostly ignored. This study attempts to couple nonlinear time-domain ground response analysis with variability of soil parameters considering existing boreholes’ location through a geostatistical method using a program written in MATLAB. To evaluate the efficiency of the proposed method, stochastic seismic ground responses at construction location were compared with those of the non-stationary random field method through real site data. The results demonstrate that applying the boreholes’ location significantly affects not only the ground responses but also their coefficient of variation (COV). Furthermore, the mean value of the seismic responses is affected more considerably by the values of soil parameters at the vicinity of the construction location. It is also inferred that considering boreholes’ location may reduce the COV of the seismic responses. Among the surface responses in the studied site, the values of peak ground displacement (PGD) and peak ground acceleration (PGA) reflect the highest and lowest dispersion due to variability of soil properties through both non-stationary random field and geostatistical methods.


Introduction
Prediction of ground response from an earthquake is a fundamental step to estimate the possible damages of structures. The amplitude of earthquake motions at the bedrock level can be drastically modified as seismic waves are transmitted through a soil deposit. Some historical earthquake events have demonstrated the effects of soil deposits on earthquake 1 3 ground motions. In recent years, several major studies have been performed to study the nature of earthquakes occurrence, the associated released energy, and the effects of site parameters.
Ground response analysis can be performed based on different considerations of problem geometry and seismic soil behaviour models. One-dimensional (1-D) ground response analyses are the most commonly used technique because of their simplicity and reasonable assumptions, including horizontal layering of soil deposits and SH-waves propagation (Sun and Chung 2008). Ground response analysis of horizontally-layered soil deposits can be conducted by employing various methods, which are categorized into Frequency-Domain (FD), including the equivalent-linear, (e.g., SHAKE (Schnabel et al. 1972)), Time-Domain (TD) (Hashash et al. 2020;Lo Presti et al. 2006;Matasovic and Vucetic 1995), Hybrid Frequency Time Domain (HFTD) (Asgari and Bagheripour 2010) methods.
Frequency-domain methods are the most widely used to estimate seismic site effects due to their simplicity. However, these methods are unable to consider true nonlinear soil models, but only linear elastic models. Since the nonlinearity of soil behaviour is well known, the linear approach was modified to the equivalent-linear method by Schnabel et al. (1972). Their efforts led to the development of an equivalent-linear program, SHAKE. Jishnu et al. (2013) used SHAKE to determine the surface response and amplification ground response analysis of a site in India. Garini et al. (2018) investigated soil amplification in three sites during the Tokachi-Oki earthquake using the SHAKE and two nonlinear programs. Torabi and Rayhani (2017) conducted a number of site response analyses in both nonlinear time-domain and equivalent-linear frequency-domain methods for a site in Canada. The results suggested that the nonlinear analysis would provide a more reliable prediction of ground motion at the natural site period. Kim and Hashash (2013) performed one-dimensional ground response analyses on Kiban-Kyoshin network (KiK-net) recorded data of the 2011 Tohoku earthquake to evaluate the applicability of equivalent-linear and nonlinear site response methods. The results revealed that the differences between the response spectra computed by equivalent-linear and nonlinear approaches for the stations where the estimated maximum strains are larger than 0.3% are more significant, and equivalent-linear site response models mispredicted the site amplifications. Yee et al. (2013) studied freefield downhole recorded ground motions from the 2007 Niigata-ken Chuetsu-oki earthquake at the site of a nuclear power plant. Their results indicated that the nonlinear site response provides a better representation of the high-frequency energy content in the initial parts of the ground motion relative to equivalent-linear analysis.
Although the equivalent-linear approach provides reasonable results when the ground motions are not significant, it may mispredict the ground responses due to severe seismic excitations (Tsai and Chen 2016). An alternative approach is to analyze the actual nonlinear ground response. For the first time, Masing (1926) presented nonlinear hysteretic behaviour by using stress-strain relations, which is named the Masing rules. Pyke (1980) and Vucetic (1990) proposed the extended Masing rules, which define unloading-reloading behaviour under general cyclic loading. As a result of this progress, a nonlinear solution of the shear wave propagation was carried out by a number of researchers. Using the hyperbolic model, Lee and Finn (1978) developed a 1-D site response analysis program. Matasovic and Vucetic (1995) modified the hyperbolic model and coupled it with extended Masing rules to better simulation of soil hysteretic behaviour. Cyclic soil behaviour also can be simulated by plasticity constitutive models. For example, Constantopoulos et al. (1973) obtained the ground response during an earthquake by Ramberg-Osgood soli model (Ramberg and Osgood 1943). Joyner and Chen (1975) used the Iwan model (Iwan 1967) for presenting a method for calculating the nonlinear seismic response of a profile with horizontal soil layers. Borja et al. (1999) utilized a bounding surface plasticity model to simulate cyclic soil behaviour response at a site in Taiwan. Hashash and Park (2001) developed a new nonlinear 1-D site response analysis model, which could consider the effect of large confining pressures. Phillips and Hashash (2009) presented two new formulations of soil damping for small and large strains. Assimaki et al. (2008) proposed a set of criteria for quantification of nonlinear behaviour susceptibility using data from three sites in the Los Angeles basin. The results indicated that nonlinear analyses are necessary at soft sites when the PGA of input acceleration is greater than 0.2 g. In a comprehensive study, Kaklamanos et al. (2015) selected a set of 191 ground motions records at six sites in the KiK-net. Then the approximate applicable thresholds of different site response analyses methods, including linear, equivalent-linear, and nonlinear analyses, were proposed. Angina et al. (2018) carried out a study to obtain the free-field seismic response of the site of Pisa tower. The results indicated that the response spectra obtained with the nonlinear code ONDA (Lo Presti et al. 2006) and corresponding amplification factors (spectrum ratio between the ground surface and base layer motions) were significantly lower than those computed using the equivalent-linear codes. It is worth mentioning that among the available 1-D ground response analysis codes, only ONDA (Lo Presti et al. 2006) and DEEPSOIL (Hashash et al. 2020) account for soil strength.
Due to the inherent variability in determining soil parameters, it is necessary to carry out non-deterministic analyses. During the past few years, there has been a considerable increase in recent researches dedicated to this topic. Rahman and Yeh (1999) coupled Monte Carlo Simulation (MCS) with the Finite Element Method (FEM) to study the variability of soil parameters in frequency-domain analysis. Wang and Hao (2002) assessed the influence of random variations of soil properties on the site amplification in the frequency-domain by employing Point Estimate Method (PEM) (Rosenblueth 1975). Nour et al. (2003) combined MCS with the FEM to simulate variability of the soil profile in seismic response analysis. Bazzurro and Cornell (2004) presented a statistical ground response study based on the nonlinear finite element program, which was presented by Li et al. (1992). Andrade and Borja (2006) performed a stochastic-deterministic analysis of site response using two computer codes (SHAKE, SPECTRA). Thompson et al. (2010) estimated the spectral amplifications from seismic properties using the square root of the impedance method. In their study, geostatistical approaches were used to increase the accuracy of the results. Rota et al. (2011) developed a procedure for estimating the stochastic site amplification of a case study in central Italy. Johari and Momeni (2015) proposed a procedure for estimating the stochastic site amplification of soil deposits with uncertain properties by combining the non-recursive algorithm of the HFTD site response approach and MCS. The influence of soil variability in the shear-wave velocity, density, thickness of profile, and soil properties was studied. The results of the sensitivity analysis showed that the damping ratio is the most influential parameter in enhancing the surface peak ground acceleration. Berkane et al. (2019) investigated the effects of the soil properties heterogeneity due to their natural randomness on the spatial variation of the response spectra at different sites. Their results showed that soil properties' randomness significantly affects the amplitudes of the response spectra. Johari et al. (2020) employed the sequential compounding method through the frequency-domain site response analysis to propose a method for conducting system reliability analysis of the surface PGA by considering soil layers crosscorrelation. Their findings in the studied site demonetrated the reasonable accuracy of the proposed method for calculating the reliability indices of ground surface PGAs.
The common feature of the mentioned stochastic ground response analysis studies are the use of linear or equivalent-linear methods and consideration of variability in soil 1 3 parameters using the random variables method, while it is expected that the use of nonlinear time-domain methods shall be more capable of considering the real behaviour of the soil. On the other hand, the random variables method does not have the ability to consider the existing boreholes' location. Therefore, this study aims to overcome the mentioned limitations and develop a method for combining the nonlinear behaviour of soil parameters and their variability by considering the borehole location effects through the geostatistical method for a real site. In order to evaluate the efficiency of the geostatistical method for considering the spatial variability of soil parameters, the results of this method were compared and discussed with that of the non-stationary random field method. The variability influence of such soil properties as Plastic Index (PI), shear wave velocity (V S ), and unit weight (γ) are assessed through MCS. The mean and COV of the ground motion parameters and the amplification factor for the non-stationary random field and geostatisticalbased stochastic analysis through the nonlinear method are determined and compared using a program prepared in MATLAB.

Nonlinear site response analysis
The main purpose of site response analysis is to obtain the wave amplitude at the top of soil layers based on bedrock with an incoming shear wave motion generated by an earthquake that propagates vertically upward, as shown in Fig. 1. The soil deposit is assumed to be made of several homogeneous horizontal layers with an infinite extent. Ground motions are considered to be generated horizontally at the soil-bedrock interface. In 1-D nonlinear site response analyses, the soil column is idealized as Multi Degrees of Freedom (MDOF) lumped mass system, in which each soil layer is depicted by a corresponding mass, nonlinear spring, and dashpot, as shown in Fig. 2.
The real nonlinear behaviour can be modeled via a time-domain analysis by the use of direct numerical integration (e.g., the central difference method, Newmark β method The mass matrix M is a diagonal matrix formed by lumping half of the mass of each two sequential layers at their common boundary, as illustrated in Fig. 2. Note that the mass of layer "i" is found by: where i and h i are the density and the height of layer i, respectively.
The stiffness matrix K should be updated at each time step to account for soil nonlinearity. According to the 1-D wave propagation assumption, the stiffness matrix can be defined using a simple shear model as follows: where G i is the shear modulus of layer i, that is obtained through nonlinear shear modulus models at each time step.
In a nonlinear analysis, soil damping is obtained from hysteretic loading-unloading cycles. The use of the viscous damping matrix [C] in the Eq. (1) is necessary to avoid unrealistic high frequencies resonance at very small strains. The viscous damping matrix  Rayleigh and Lindsay (1945), known as Rayleigh damping, is defined as follows: The scalar values of α and β are calculated from two natural frequencies corresponding to natural modes m and n, as follows: where ξ m and ξ n are the damping ratios for the target frequencies f m and f n .
Based on Stewart and Kwok (2008) recommendation, the f m and f n frequencies can be the natural frequencies corresponding to the first and third natural modes of the soil profile and the ξ m and ξ n can be set equal.

Nonlinear constitutive models based on hyperbolic relationships
Several empirical models have been proposed in the literature for estimating shear modulus reduction and damping curves. Among important contributions are the models of Vucetic and Dobry (1991), Darendeli (2001), Menq (2003), Kishida et al. (2009), and Yang and Woods (2014). The empirical models, such as Darendeli (2001), Menq (2003), and Kishida et al. (2009) are utilized to simulate the dynamic soil parameters in the site response analyses when no dynamic laboratory tests were conducted. A nonlinear simple-shear model is based on the following principles: • A basic formula for determining the stress-strain path at initial loading (backbone curve) • A series of criteria (rules) for determining the unloading and reloading soil behaviour.
In this paper, the Darendeli (2001) modified hyperbolic model was implemented in modeling the backbone curve and the extended Masing rule (Pyke 1980;Vucetic 1990) was used for governing the unloading and reloading behaviour of subsequent cycles. Additionally, in order to prevent the overestimation of damping during hysteretic cycles in the Darendeli model (Darendeli 2001), the second Masing rule was modified in accordance with Pyke (1980) proposed hypothesis. Darendeli (2001) proposed a model to estimate shear modulus reduction and damping curves from a soil test database using a first-order, second-moment Bayesian method. The model is based on the simple hyperbolic model proposed by Kondner and Zelasko (1963) and Hardin and Drnevich (1972), but added a curvature coefficient, 'a = 0.919', as follows:

Darendeli nonlinear model
where γ r , is the shear strain when G/G max equals 0.5, which can be obtained by the following equation: where PI is plasticity index, OCR is the over-consolidation ratio, and σ′ m is the mean effective confining pressure. Based on Darendeli assumption, total damping is formed from two main parts, namely, small strain damping (D min ) caused by internal friction and material viscosity and the damping corresponding to soil nonlinearity or hysteretic behaviour (D Masing ) (Darendeli 2001): where D is total damping, and F is a reduction factor. D Masing is proportional to the ratio of the dissipated energy to stored strain energy in one complete cycle of motion. D min is calculated by the following equation (Darendeli 2001): where f is loading frequency.

Extended Masing rules
The stress-strain relationship for initial loading is named the backbone curve. The Masing rules (Masing 1926) and extended Masing rules (Pyke 1980;Vucetic 1990) describe unloading-reloading behaviour under general cyclic loading conditions, as shown in Fig. 3. The Masing and extended Masing rules are explained as follows (Hashash and Park 2001): (1) The stress-strain curve follows the backbone curve for initial loading.
(2) The unloading and reloading curves have the same shape as the backbone curve but are enlarged by a factor of "C" with the origin shifted to the reversal point (γ rev , τ rev ).
where the "C" is the scale factor of the stress-strain relationship between the initial loading and the unloading and reloading cycles. The "C" factor can be considered as a constant coefficient equal to 2. An alternative way suggested by Pyke (1980) is to take it as a function of stress level: where the τ ult is the ultimate soil stress (τ ult = G max γ r ), the "±" symbol is equal to "+" when the shear strain increment (dγ) is positive and is equal to "-" when the shear strain increment is negative. (3) When the unloading or reloading curve exceeds the maximum past strain and intersects the backbone curve, the stress-strain path follows that of the backbone curve until the next reversal point.
(4) When the unloading or reloading curve intersects the curve from the previous cycle, then the stress-strain curve follows the path of the previous cycle.
As mentioned before, in this study, instead of using a constant factor in the second Masing rule, the Pyke (1980) hypothesis is implemented in the MATLAB code.

Error reduction process in nonlinear analysis
The procedure of solving a nonlinear problem by a series of linear estimations with a constant time step can lead to unacceptably inaccurate results. The two major causes of errors are utilizing tangent stiffness instead of the secant stiffness and employing a constant time step, which causes a delay in detecting reversal points in the force-deformation relationship (Chopra 1995). In order to reduce these errors, it is necessary to use iterative methods such as the Newton-Raphson. The Newton-Raphson iteration method is divided into two main categories, the regular and the modified Newton-Raphson. The difference between these two methods is in the basis point for calculating the stiffness matrix. In the regular Newton-Raphson, the stiffness is updated at every iteration. While in the modified Newton-Raphson method, the stiffness is once evaluated at the beginning of the iterations. All the consecutive iterations are calculated with the same stiffness, which reduces the computational cost and leads to faster convergence. In the present study, the modified Newton-Raphson method is implemented in the nonlinear analysis to reduce the numerical errors.

Implementing the soil variability in the ground response analysis
The variable nature of soil properties that affect the ground responses indicates the necessity of considering the stochastic aspect of input parameters. Randomizing the input parameters allows investigating the way variability of the latter is mapped onto the variability of outputs, which here are ground responses. This process will improve engineering judgment. The first step in the stochastic analysis is selecting random variables and determining their characteristics, including the mean, standard deviation, and the statistical distribution type (Johari and Amjadi 2017;Johari and Kalantari 2021).
Soil properties variation is one of the main sources of uncertainty in the ground response problem that can alter the results. Although some innovative methods, such as stochastic analysis in the scale domain, can consider new aspects of uncertainty like long-term persistence behavior in rock and soil materials (Dimitriadis et al. 2019(Dimitriadis et al. , 2021, still the common methods for evaluating the soil parameters variability are the random variable, stationary random field, non-stationary random field, and geostatistics. The random variable method is usually utilized when it aims to evaluate the variability of an output parameter, which is defined as a function of several input parameters with known variability. This method does not affect the location of the recorded data . The stationary random field method creates a variable range of parameters in a domain by considering a constant mean and standard deviation for the whole field. Although this method has been widely applied in geotechnical studies, it cannot implement the depth dependency of the mean and variance soil properties (Li et al. 2014). Therefore, the random variable and stationary random field method can not predict reasonable soil spatial variability in stochastic problems such as site response analysis. The site response analysis is significantly affected by vertical correlations and the natural trend of soil parameters (e.g., the pattern of shear wave velocity with depth in a particular site) (Rathje et al. 2010;Van Nguyen and Lee 2021). Therefore, the non-stationary random field and geostatistical method, because of their capability for considering the natural depth-dependent trend (i.e., vertical correlations) of soil parameters, are employed in this study and explained in the following.

Non-stationary random field
Site investigation data revealed the depth dependency of soil properties in natural fields (Phoon and Kulhawy 1999). The non-stationary random field describes the depth-dependent nature of the soil parameters. Li et al. (2014) proposed a non-stationary random field model for simulating shear strength variability. In the present study, the same procedure was utilized to model the soil parameters' spatial variability. In this way, nonlinear regression was conducted on the field data to find the best fit of the parameters' depth-dependent trend. Then the non-stationary random fields were generated base on the obtained trends.

Geostatistical estimation
Geostatistics, initially developed in the mining industry, is now being applied widely in geotechnical projects (Johari and Gholampour 2018;Thompson et al. 2010). The most common geostatistical estimation technique is the ordinary Kriging method (Krige 1951), which is based on the BLUE algorithm (Best Linear Unbiased Estimator). It is used as a univariate (single regionalized variable) tool to perform spatial interpolation between known borehole data. Despite some drawbacks of the Kriging method, such as uncertainty in empirical variogram fitting, as well as several degrees of freedom in the calculation of spatial predictions (Elogne et al. 2008;Marchant and Lark 2004), this method is still known as the most widely used method in estimating spatial uncertainty of soil parameters in geotechnical engineering. Furthermore, Kriging allows estimating a spatially correlated soil property at a point where no measurements are available from measurements at neighboring points. This capability of the geostatistical method can be very useful, especially when project conditions and technical limitations do not allow sampling of some desirable points. On the other hand, multivariate geostatistical analysis refers to the study of two or more co-regionalized variables. The Cokriging method can be applied to this kind of analysis in order to take the complicated cross-correlation structures into consideration (Johari and Gholampour 2018). In this paper, the Kriging method was employed to determine the variability of the dynamic soil parameters at construction location.

The procedure of stochastic analysis
In the previous sections, the procedure of deterministic nonlinear site response analysis was described. In this section, the process of applying variability of the soil parameters in the nonlinear ground response analysis is expressed.
Generally, this procedure can be divided into three main parts. The first part is the ground motion process, which includes the definition of an input acceleration, the baseline correction, and the sub-dividing of that. The second part is the random soil profile process in which, at first, the layers' soil properties are obtained by downhole and laboratory tests. Then two different methods, including non-stationary random field and geostatistical estimation, are implemented for considering the spatial variability of soil parameters. In the non-stationary random field model, estimations are based on depth-dependent trends, obtained by conducting nonlinear fitting on the average values of boreholes' data. In contrast, the geostatistical method is based on Kriging estimation to apply the boreholes' location on determining the mean and standard deviation of the dynamic soil parameters. In the next step, using the non-stationary random field and geostatistical method, "N" different random soil profiles were generated. In the third part, the nonlinear ground response analysis is carried out for each simulation. Finally, through the MCS, the statistical distribution of outputs, including peak ground displacement (PGD), peak ground acceleration (PGA), acceleration response spectra (S a ), and the amplification factor of the site, are calculated. Figure 4 shows the flowchart of the proposed stochastic procedure for considering the existing boreholes' location.

Computer program
In this study, a computer program has been developed in MATLAB for stochastic nonlinear time-domain ground response analysis based on two different methods, including nonstationary random field and geostatistical estimation. The major ability of the program is as follows: • Signal processing on strong-motion data to obtain suitable input acceleration (baseline correction and sub-dividing input acceleration).
• Calculation of the initial dynamic soil parameters of layers from imported profile data.
• Modeling the soil layers spatial variability through the non-stationary random field model and geostatistical method.  Fig. 4 Flowchart of the proposed method for the nonlinear stochastic ground response analysis • Considering the nonlinear behaviour of soil profile by the modified Darendeli model (Darendeli 2001). • Solving the equation of motion by the Newmark β method coupled with the Newton-Raphson procedure by controllable accuracy and convergence. • Calculation of stochastic peak ground motion parameters, S a , and the amplification factor for each layer by MCS.

Geotechnical and geophysical parameters of the site
To demonstrate the ability of the proposed method for considering boreholes locations effects in stochastic nonlinear site response analysis, a site with real data, which is located in Shiraz city in the Fars province of Iran, is assessed. Three boreholes with approximately 36.0 m depth from the ground surface were drilled to investigate the subsurface layers and soil properties.
In this project, due to technical limitations, it was not possible to dig a borehole at the construction location; therefore, ground response analyses were performed based on adjacent boreholes' data. The boreholes' arrangement and the hypothetical borehole at the construction location are presented in Fig. 5. The propagation velocity of shear waves was measured in each borehole every 2.0 m. In addition, unit weight and plasticity index were inferred from undisturbed samples retrieved from each borehole. The results are given in Table 1. It can be seen that the soil types are mainly of fine grain size. In these tables, PI, γ, and V s are the mean of plasticity index, unit weight, and shear wave velocity of the soil layers, respectively.
Although the 1-D ground response analysis is only appropriate for laterally homogeneous sites, even homogeneous sites always include some degree of uncertainty. Therefore, it is necessary to take advantage of a proper method for interpolating the boreholes' data at the construction location to achieve the soil spatial variability. In this study, to estimate the spatial variability of the selected stochastic soil parameters from known data of the three boreholes at the construction location (hypothetical borehole), the non-stationary random  field model and the geostatistical method were used. For this purpose, in the non-stationary random field model, after averaging between the boreholes' data, the depth-dependent trends of the mean and standard deviation of soil parameters are obtained through the nonlinear fitting. These results are summarized in Table 2.
On the other hand, in the geostatistical method, the mean and standard deviation of soil parameters are determined using Kriging estimation. In this way, an exponentially decaying correlation function was utilized. Table 3 represents the results of the Kriging interpolation at the construction location. A sample calculation process for the first layer is presented in the "Appendix".
Comparing the estimated soil parameters in Tables 2 and 3 reveal that the estimated shear wave velocities and unit weights by the non-stationary random field are lower than those by the geostatistical method. As a result, it can be stated that the soil parameters predicted by the random field method lead to softer behaviour of soil compared to the geostatistical method.
Among the statistical distributions such as uniform, normal, and log-normal distributions, researchers widely use log-normal distribution to model the variability of soil parameters because it has a zero probability of producing negative numbers (Fenton and Griffiths 2002;Kottke and Rathje 2009;Liu and Chen 2006;Van Nguyen and Lee 2021). In this study, in order to track the effect of soil properties variation on the ground responses, the log-normal non-stationary random field has been employed to generate profiles of the selected stochastic parameters (V s , γ, and PI). As well as in the geostatistical method, the selected stochastic parameters were modeled by the log-normal distribution. Figures 6, 7, 8, 9, 10 and 11 represent the variability of V s , γ, and PI at the construction location by the two methods. In each figure, the Probability Density Function (PDF) of the fourth layer's corresponding parameters is plotted. Furthermore, in these figures, the red line corresponds to the mean value of soil parameters, and the gray lines are related to 2000 realizations of the random profile. Since the boreholes' data (BH.1-3) indicates a low  level of inhomogeneity, in both non-stationary random field and geostatistical methods, the scale of fluctuation is considered 45 m and 5 m in the horizontal and vertical directions, respectively (Phoon and Kulhawy 1999). It is worth mentioning that among these two implemented models, only the geostatistical method is capable of considering existing boreholes' locations in the estimations of soil parameters at the construction location. Furthermore, regarding the issue of parsimony in the number of utilized parameters in two methods, it should be noted that the input parameters for both methods are the same, and the only different parameter is the location of the drilled boreholes, which its information is recorded in all geotechnical projects. Therefore there is no need to spend extra money or time.

Verification of the nonlinear developed code
For verifying the accuracy of the developed MATLAB code in obtaining soil layers' response, an example problem was simulated by the nonlinear analysis method of the DEEPSOIL software (Hashash et al. 2020) and the code. The mean values of soil parameters, which are provided in Table 3, were utilized in this analysis as soil properties of the layers. In order to ensure the mobilization of nonlinear soil behaviour, it was necessary to apply a relatively strong earthquake to the soil profiles. Therefore, in the following analyses, the 1994 Meymand earthquake (M w = 6.1) with the peak acceleration of 0.51 g, which occurred on June 20, at about 84 km from Shiraz city, was considered as input motion at the bottom of the soil profile. Figure 12 shows the 1994 Meymand earthquake acceleration time history.
The normalized modulus reduction and damping ratio curves corresponding to mean values of 18 layers soil properties were estimated through the Darendeli (2001) model are shown in Figs. 13 and 14.
The predicted absolute acceleration time histories, a close-up view from 3 to 7 s time window, and S a of the surface layer for the proposed method and the DEEPSOIL are compared in Figs. 15, 16 and 17, respectively. These graphs represent a good agreement between the predicted responses by the developed code and the DEEPSOIL. Figure 18 demonstrates the predicted stress-strain loops through the proposed method and DEEP-SOIL. This figure indicates that the output of the code provides a good match with DEEP-SOIL. The slight mismatch in the results is due to differences in the employed soil models in the two programs (i.e., the modified Darendeli (2001) model in developed code and the modified Kondner and Zelasko (1963) in the DEEPSOIL (Hashash et al. 2020)). As  Table 4 compares the results as obtained by means of the developed code and DEEP-SOIL, including the peak of S a , maximum and minimum values of acceleration, shear strain, and normalized shear stress (τ/σ′ 0 ). This table also emphasizes the acceptable agreement between the results of the two methods. As mentioned before, the main goal of this study is to consider the existing boreholes' location effects in stochastic nonlinear site response analysis through geostatistical methods. Although seismic ground responses are influenced by input motion, soil characteristics, and site geometry, the boreholes' location only affects the estimation of soil properties. Therefore, in the present study, loading and geometric parameters such as input motion and layers' thickness have been taken as deterministic parameters to avoid further complexity in interpreting the results. Among soil characteristics, the variations of shear modulus and damping with respect to strain are the major factors contributing to the nonlinear amplification and attenuation of soil deposits during earthquake excitations. Sensitivity analysis performed in previous studies on laboratory data of different soils illustrates that variation of soil shear modulus is more influenced by shear wave velocity, unit weight, and soil plasticity. Furthermore, the most effective parameters on soil damping are plasticity index and confining pressure (Darendeli 2001;Stokoe et al. 1999;Vucetic and Dobry 1991). Therefore, in the present study, the plastic index, shear wave velocity, and unit weight which significantly affect the nonlinear behaviour of soil were selected as stochastic soil parameters. To this end, the proposed deterministic computer program was extended to perform stochastic input parameters generation and iterative calculation via MCS. In the next step of this study, the stochastic nonlinear ground response analysis was conducted using the generated random profiles (see Sect. 7) to investigate the geostatistical method's applicability compared to the non-stationary random field.
The number of required Monte Carlo simulations depends on the desired level of confidence in the solution and the number of variables (Johari et al. 2012). Since the required number of trials based on proposed statistical equations in literature has a high computational cost, fewer simulations are typically performed in the time-consuming dynamic simulations such as the nonlinear time-domain analyses (Bazzurro and Cornell 2004). Therefore, in this study, a statistical analysis was carried out to determine the required number of stochastic simulations for MCS. The mean of PGA is considered as the target parameter for checking the stability of stochastic analysis. The employed criterion for this stability  analysis was based on a number of simulations where the average PGA variation for at least 30 consecutive numbers was less than 0.01 g. As shown in Fig. 19, when the number of simulations reaches 1675, the mean PGA of the surface layer becomes stable. Consequently, 1,800 simulations would be sufficient in stochastic analysis of the studied site. All the following stochastic results are obtained using 1,800 simulations.

Stochastic results of nonlinear ground response
The stochastic results of seismic ground response can be described in PGD, PGA, maximum and minimum shear strain, surface acceleration response spectra, and amplification factor. These stochastic results will be described and compared for the non-stationary random field and geostatistical methods in the following subsections.

Peak ground displacement
Earthquakes can cause residual and permanent deformations, which generate excess pressure and stress concentration in subsurface structures. Peak ground displacement is an important parameter in earthquake-resistant design, especially for buried structures such as pipelines and deep foundations. Figures 20 and 21 illustrate the variation of nonlinear PGD through the soil profile due to selected stochastic soil parameters by the non-stationary random field and geostatistical methods, respectively. These graphs were plotted by obtaining the PGD of the ground surface with its time and capturing displacements of sub-layers at the corresponding time. Comparing the displacement on bedrock and the ground surface indicates the high potential of the site's profile in magnifying earthquake waves. Furthermore, these figures demonstrate that the variability of input parameters, which are obtained using the non-stationary The variations of the ground surface PGD from both methods have a distribution near log-normal with the stochastic parameters summarized in Table 5. These stochastic parameters indicate that PGDs obtained using the geostatistical method have a slightly lower mean and less variation. It can be seen that the coefficient of variation for the PGD results of the geostatistical method is considerably lower than that of the non-stationary random

Variation of mobilized extrema shear strain of layers
The extrema (maxima and minima) shear strain profile determines what level of nonlinearity is mobilized at each soil layer. In this study, at each depth (the middle of each layer), the highest and lowest extrema of shear strains are extracted from the corresponding hysteresis stress-strain loops. Figures 22 and 24 show the variation of extrema shear strain, respectively, obtained from the non-stationary random field and geostatistical methods. Figures 23 and 25 depicted the calculated mean of extrema shear strain (red line) and mean plus or minus one (gray zone) and two (green zone) standard deviations for these two methods. It can be seen that the general trend of obtained shear strain is consistent with the results of previous researches (Hashash et al. 2020;Lo Presti et al. 2006). So that after reaching a peak value, it approaches zero near the ground surface.
According to these figures, it can be inferred that in the studied site, through both non-stationary random field and geostatistical methods, the highest shear strain and, as

Fig. 22
Variation of extrema shear strain of layers by nonstationary random field a result, the most level of nonlinearity mobilized in the third layer (5 m depth). As well as it can be seen that the soil in the profiles which are generated by the non-stationary random field method experiences higher shear strain values, which is due to the fact that these profiles are softer compared to the generated profiles using the geostatistical method. The obtained standard deviation of the extrema shear strain from the non-stationary random field method at all depths is greater than the corresponding value for the geostatistical method. The highest value of the extrema shear strain's standard deviation occurred at 5 m depth and equal to 0.11% and 0.08%, respectively, for the non-stationary random field and geostatistical methods.

Acceleration response spectrum
The acceleration response spectrum (S a ) is a very useful tool in designing structures under earthquake excitations. In this study, the stochastic nonlinear acceleration response spectrums have been estimated through the developed MATLAB code. Figures 26 and 28 show the variation of nonlinear S a , obtained from surface response acceleration for 5% structural damping for the non-stationary random field and geostatistical methods. Figures 27 and 29 represent the calculated mean of S a at the ground surface (red line) and mean plus or minus one (gray zone) and two (green zone) standard deviations for these two methods. Additionally, the S a of the input accelerograms at the bedrock is presented by the dashed line.
An approximation for determining the buildings' natural period is dividing the number of stories by 10. Natural periods vary from about 0.1 to 0.5 s for one to a four-story building, respectively. Other factors, such as the building's structural system, construction materials, and geometric characteristics, also affect the period, but height is the most important factor (Arnold 2007). Based on Figs. 26, 27, 28 and 29, it can be concluded that in the assessed site by the considered earthquake, buildings with short and medium height are more affected while the high-rise buildings are not influenced significantly. Figures 30 and 31 compare the variation of the mean and the COV of S a for the non-stationary random field and geostatistical methods. Figure 30 illustrates that for the geostatistical method at periods less than 0.47 s (red dash line), the mean of S a is above that of the Fig. 25 The mean of extrema shear strain of layers plus or minus one and two standard deviations by the geostatistical method 1 3 non-stationary random field. This matter is completely consistent with predicted the profiles by the two methods and the soils' nonlinear seismic behaviour. So that, by comparing Figs. 6 and 9, it can be concluded that the profiles generated by the geostatistical method are stiffer than those of non-stationary random field, so it was expected that the S a obtained from these profiles have higher values in low periods and vice versa lower values in high Fig. 26 The nonlinear stochastic S a at the ground surface by nonstationary random field Fig. 27 The nonlinear mean S a plus or minus one and two standard deviations at the ground surface by non-stationary random field 1 3 periods. In contrast, the profiles predicted by the non-stationary random field resulted in higher S a values at high periods due to being softer. Figure 31 indicates that the COV of S a obtained from the non-stationary random field method is bigger than that of the geostatistical method. By comparing Figs. 30 and 31,   Fig. 28 The nonlinear stochastic S a at the ground surface by the geostatistical method Fig. 29 The nonlinear mean S a plus or minus one and two standard deviations at the ground surface by the geostatistical method it can be inferred that at the studied site, the most critical period is about 0.08 s (green dash line). Because at this point, both the mean of S a and its COV are high. Therefore, in this site with the considered earthquake, structures with natural periods close to 0.008 s are more critical than others. Fig. 30 The mean S a of the non-stationary random field and geostatistical methods Fig. 31 The COV of S a for the non-stationary random field and geostatistical methods

Peak ground acceleration
Peak ground acceleration is defined as the maximum amplitude of absolute acceleration time history. The PGA is the most commonly used ground motion parameter in engineering applications, especially in building codes and earthquake-resistant design. Figure 32 shows the bar chart and fitted probability density function of the site PGA for the nonstationary random field and geostatistical methods with the stochastic parameters summarized in Table 5. It can be seen that the variations of the ground surface PGA from both two methods have distributions close to normal, which are in agreement with the findings of previous researchers (Bahrampouri et al. 2019;Rota et al. 2011;Tran et al. 2020). In the present study, statistical assessments revealed that the log-normal distribution could provide a good match with predicted PGA data of both methods (correlation coefficient of 0.98 and 0.97 for the non-stationary random field and geostatistical methods, respectively) as well as the Kolmogorov-Smirnov test (Lilliefors 1967) was utilized to determine whether the output PGA distributions followed the normal distribution. The results of this test show that considering the 5% significance level, only the PGA distribution obtained by the non-stationary random field method follows the normal distribution. Figure 32 also indicates that the PGA distribution from the geostatistical methods has a bigger mean value than that of the non-stationary random field. Furthermore, the dispersion of PGAs from the non-stationary random field is higher than the geostatistical method. These results are mostly due to the boreholes' location effects in the geostatistical method because in the studied site, the soil obtained from the borehole closer to the construction location was stiffer and its impact is considered in the generated stochastic profiles. Figure 33 shows the Complementary Cumulative Distribution Function (CCDF) of PGAs at the surface for the non-stationary random field and geostatistical methods. The CCDF determines how often the random variable is above a particular level. According to Table 5 and the obtained mean PGA for the non-stationary random field method (0.56 g), the probability of Fig. 32 The predicted PDF of PGA at the site surface for the non-stationary random field and geostatistical methods occurrence of PGAs can be determined using Fig. 33. This figure represents that the probability of occurring PGAs equal to or greater than 0.56 g for the non-stationary random field and geostatistical method is about 50% and 75%, respectively. In other words, the occurrence of PGAs greater than 0.56 is much more likely in the geostatistical method.
The mean, standard deviation, and coefficient of variation for both methods are presented in Table 6. The results in this table indicate that COV for the non-stationary random field is higher than the geostatistical method; therefore, the geostatistical method's results include a lower level of uncertainties.

Amplification factor
The amplification factor is determined as the ratio of displacement or acceleration response spectra (with 5% damping) at the ground surface and the input motion's corresponding response spectra at the bedrock. Figures 34 and 36 show the nonlinear amplification factor variation for the non-stationary random field and geostatistical methods. Figures 35 and 37 represented the predicted mean amplification factor (red line) and the mean plus or minus one (gray zone) and two (green zone) standard deviations for these two methods. Figures 38 and 39 compare the mean and COV variation of the predicted site's amplification factor for the non-stationary random field and geostatistical methods. Figure 38 represents that the peak of the amplification factor occurs at periods about 0.6 s and 0.4 s for the non-stationary random field and geostatistical method, respectively.

Fig. 33
The CCDF of PGA at the site surface for the non-stationary random field and geostatistical methods 1 3 It also demonstrates that the amplification factor of the geostatistical method at periods less than 0.47 s (green dash line) is above that of the non-stationary random field, and after this point, the amplification factor of the non-stationary random field is The mean amplification factor plus or minus one and two standard deviations at the ground surface by non-stationary random field 1 3 higher. This issue is resulting from the employing of boreholes' location effect (spatial effect of boreholes) in the geostatistical method. In other words, at the studied site, the BH.2 is closer to the construction location. Since BH.2 is relatively stiffer than the other two boreholes, the amplification factor of generated profiles by the geostatistical method is higher at lower periods. Furthermore, the amplification factor obtained by the Fig. 36 Stochastic amplification factor at the ground surface by the geostatistical method Fig. 37 The mean amplification factor plus or minus one and two standard deviations at the ground surface by the geostatistical method 1 3 non-stationary random field is higher at higher periods (because it is composed of softer soils). Figure 39 shows that the COV of the amplification factors obtained from the nonstationary random field method is greater than that of the geostatistical method. This Fig. 38 The mean amplification factor of the non-stationary random field and geostatistical methods Fig. 39 The COV of the amplification factor for the non-stationary random field and geostatistical methods issue demonstrates the effect of the geostatistical method in reducing the dispersion of responses.

Discussion
In this section, the obtained stochastic seismic responses from the studied site are discussed. The mean value of responses and their dispersions are examined through two different perspectives: nonlinear soil attenuation mechanism and variability of the ground responses.

Nonlinear soil attenuation mechanism
To assess the nonlinear attenuation mechanism of soil behaviour, the obtained mean values of stochastic ground responses are compared for the non-stationary random field and geostatistical methods. Table 7 represents the mean of the stochastic responses from the non-stationary random field and geostatistical methods. It can be seen that the use of the non-stationary random field method leads to larger ground displacements, shear strain, velocity response spectrum (S v ) and displacement response spectrum (S d ). In contrast, the geostatistical method provides larger surface accelerations, S a , and amplification factors, which is due to the fact that the profiles generated by the non-stationary random field are softer than those by the geostatistical method (i.e., they have less velocity and shear modulus). These results are in agreement with field observations and the findings reported by other researchers in the nonlinear attenuation behaviour of soft soils under large earthquake excitations (Idriss 1991;Suetomi and Yoshida 1998). Figure 40 was originally drawn by Idriss (1991) and later confirmed in subsequent researches (Suetomi and Yoshida 1998). This figure indicates that the occurrence of intense (i.e., PGA > 0.43 g) earthquakes at sites with soft and medium soils can activate the nonlinear attenuation mechanism (called 'deamplification'), a mechanism which leads to lower PGA and higher shear strain and surface displacement, and instigates move of the predominant period toward the long-period range by increasing the input motion intensity in soft soil sites.
In order to provide a better evaluation of the present study results, the authors added the mean of PGA through the non-stationary random field and geostatistical methods in Fig. 40. Since the peak acceleration of the considered input motion is 0.51 g, it is susceptible to be able to initiate the nonlinear attenuation mechanism. On the other hand, since the generated profiles using the non-stationary random field are softer than those of the geostatistical method, more nonlinear attenuation has occurred in the seismic responses attained from the non-stationary random field. Therefore, the mean stochastic PGA through the non-stationary random field is lower than those acquired from the geostatistical method.

Comparing the variability of the ground responses
In this section, variability in ground responses with respect to the considered stochastic soil parameters is discussed. Table 8 compares the COVs of different seismic responses at the ground surface through the non-stationary random field and geostatistical methods. This table reveals that, among these surface responses, the PGD and PGA reflect the highest and lowest dispersions due to variability of soil properties. It can also be seen that the responses attained through the geostatistical method have lower COVs than those from the non-stationary random field, which is due to considering the effect of the borehole location. It is also worth mentioning that since COV is invalid for parameters whose means are close to zero (Busing et al. 2005), the value of COV corresponding to the shear strain is not included in this comparison.

Conclusions
Nonlinear time-domain ground response analysis implementing a stepwise integration procedure provides an accurate framework for simulation of the real nonlinear behaviour of soil. On the other hand, seismic ground responses mainly depend on the soil's dynamic behaviour, which is always associated with a significant degree of heterogeneity. In this paper, a new procedure for conducting stochastic nonlinear ground response analysis considering the existing boreholes' location on the soil parameters heterogeneity through the geostatistical method, i.e., Kriging estimation, was developed and applied to a real site. Initially, a program for deterministic nonlinear site response analysis in the time-domain was developed in MATLAB and then extended to stochastic analysis. The results of the deterministic analysis demonstrate acceptable agreement between the developed program and outputs from DEEPSOIL. For the stochastic analysis, two different methods, including non-stationary random field and geostatistical estimation, were implemented to consider the spatial variability of soil parameters. The non-stationary random field model is based on depth-dependent trends obtained by nonlinear fitting of the average values of the boreholes' data. In contrast, the geostatistical method was based on Kriging estimation to apply the effect of the boreholes' location on the mean and standard deviation of the dynamic soil parameters. The obtained variability of soil parameters at the construction location was used to generate random profiles through MCS. The stochastic seismic ground responses and the amplification factor were determined through both of the aforementioned methods. Some conclusions are summarized below: • The results show that considering the effects of boreholes location can have a significant influence on seismic ground responses. In other words, the geostatistical method for prediction of the values of soil parameters at the construction location borehole may assign more weight to the boreholes which are closer to the target interpolation location. This procedure can lead to a more reasonable estimation of the variability associated with soil parameters. In the site under study, the second borehole (i.e., BH.2) was the closest to the construction location; therefore, it had the greatest impact on the estimations. Additionally, since BH.2 has a relatively stiffer soil (in terms of more velocity and shear modulus) than those attained from the other two boreholes, the generated profiles through the geostatistical method are also stiffer than those from the non-stationary random field method. As a result, in the studied site, the use of the nonstationary random field method leads to larger PGD, shear strain, S v , and S d , while the geostatistical method provides larger surface accelerations, S a , and amplification factor. • According to the obtained shear strain, it can be concluded that in the studied site the profiles generated by the non-stationary random field method may experience higher shear strain values, which can be due to the fact that these profiles are softer compared to the generated profiles using the geostatistical method. Additionally, the highest shear strain is mobilized in the third layer (i.e., 5 m depth). The obtained standard deviation of the extrema shear strain from the non-stationary random field method at all depths is greater than the corresponding value for the geostatistical method. • Based on the stochastic PGA assessment, it can be concluded that variations of the ground surface PGA from both methods have nearly normal distributions. • The results indicate that the peaks of the stochastic amplification factors occur at periods around 0.6 s and 0.4 s for the non-stationary random field method and geostatistical method, respectively. This also demonstrates that the amplification factor of the geostatistical method at periods less than 0.47 s is above that of the non-stationary random field and, after this point, the amplification factor of the non-stationary random field is higher. • Considering the borehole location in the estimations performed through the geostatistical method causes the final results to be closer to the data recorded near the construction location. Accordingly, this leads to reduction of the obtained soil parameters' variability through the geostatistical method when compared to the non-stationary random field method. The reflection of the uncertainties in the input parameters to the seismic ground responses reveals that surface responses through the geostatistical method have lower COV, and among the different responses, PGD and PGA reflect the highest and lowest dispersions due to variability in soil properties.
Placing the coordinate axis origin at the left bottom corner of the site gives the construction location coordinates x = (33, 15). Thus, the right-hand side vector M is: In the next step, the unknown weights can be obtained from the matrix equation Kβ = M (ignoring the Lagrange parameter): Since the correlation matrix is unique for PI, γ, and V s , the weights will be identical for them. Thus the best estimates at the construction location for the first layer are: The following equation can determine the estimation of variances: where x is the vector of correlation coefficients between the samples and the construction location. For the Kriging weights and given correlation structure, this yields: which gives the following individual estimation of standard deviations for the first layer: This procedure is repeated to estimate soil properties for all layers.