New Equations to Evaluate Lateral Displacement Caused by Liquefaction Using the Response Surface Method

Few empirical and semi-empirical approaches have considered the influence of the geology, tectonic source, causative fault type, and frequency content of earthquake motion on lateral displacement caused by liquefaction (DH). This paper aims to address this gap in the literature by adding an earthquake parameter of the standardized cumulative absolute velocity (CAV5) to the original dataset for analyzing. Furthermore, the complex influence of fine content in the liquefiable layer (F15) is analyzed by deriving two different equations: the first one is for the whole range of parameters, and the second one is for a limited range of F15 values under 28% in order to the F15’s critical value presented in literature. The new response surface method (RSM) approach is applied on the basis of the artificial neural network (ANN) model to develop two new equations. Moreover, to illustrate the capability and efficiency of the developed models, the results of the RSM models are examined by comparing them with an additional three available models using data from the Chi-Chi earthquake sites that were not used for developing the models in this study. In conclusion, the RSM provides a capable tool to evaluate the liquefaction phenomenon, and the results fully justify the complex effect of different values of F15.


Introduction
When, during earthquake motion, pore water pressure rises because of applied dynamic loads, the loose saturated sand layer that is relatively close to the ground surface is liquefied.Liquefaction can be discovered through manifestations such as (1) a sharp decrease in the frequency content of a sand layer, (2) settlement, (3) flow slides, (4) sand boiling, (5) foundation failure, and (6) lateral displacement (DH).
The movements of sand blocks, which have destroyed and affected constructions and infrastructure, ranging from a few centimeters to some meters [1], have been reported.Lateral displacement can be significantly damaging for piles, piers, and pipe lines during and for a short time after earthquakes and causes more damage to structures and infrastructures than any other type of liquefaction-induced ground failure.In this phenomenon, the large blocks of soil move towards the free face or along the slope.Researchers have developed several different models and approaches to predict the DH caused by liquefaction for some decades.Some of them have proposed numerical approaches [2][3][4][5][6] such as the finite element method (FEM) and the finite difference method (FDM).Next to that, analytical approaches have been developed, for example, minimum potential energy [7] and the sliding block model [8][9][10][11][12].
Among them, due to the complicated input model parameters and difficulties in their calculations, as well as because of the complex mechanism of liquefaction, empirical and semi-empirical models are the most common models that have been performed and developed by engineers and researchers [13][14][15][16][17].However, in most cases, because of the scarcity and shortage in their database, some aspects of this phenomenon, such as geology, fault type, and the effect of near-fault sites, have been ignored, with the exception of Zhang et al. [18], who used Japanese spectral attenuation models, or Bardet et al. [19], who considered peak ground velocity (PGV) to overcome this shortage and improved the model proposed by Youd et al. [20].Nevertheless, a shortage of studies in geology and motion frequency effects still exists.
In 2006, Kramer [21] reported the result of substantial research on around 300 ground motion parameters and declared that the most efficient and sufficient intensity measure on liquefaction is one standardized form of cumulative absolute velocity (CAV), which eliminated amplitudes less than 5 cm/sec 2 and is defined as CAV5.Sufficiency defines which parameter is independent to estimate the target (increasing pore water pressure herein), and efficiency expresses which parameter is able to predict the target with lower uncertainties [22].This parameter quantifies aspects of applied frequency load, which can be affected by the near-fault region aspect and causative fault type of earthquakes.Hui et al. [23] proposed an index of PGV to peak ground acceleration (PGA) to characterize the effect of liquefaction on the piles in near-fault zones.Further, Kwang et al. [24], through performing some uniform cyclic simple shear laboratory tests, demonstrated that CAV5 provides the highest correlation with DH among ground motion parameters.While the significant correlation between CAV5 and the evaluation of liquefaction have been characterized, no attempt has been made to take it into the account when developing empirical and semi-empirical models.
Furthermore, artificial intelligence has been applied to develop models and correlations to predict DH using databases that were collected from sites [25][26][27][28].Training is organized to minimize the mean square error (MSE) function.Wang et al. [25] used a back-propagation neural network to develop a model for the prediction of lateral ground displacements caused by liquefaction.They applied the same records used by Bartlet et al. [1], along with 19 datasets of Ambraseys et al. [29].Among all datasets, 367 data points were used for the training phase, and the extra 99 datasets were used for the testing phase, while no validating phase was conducted.The model was developed using the same parameters suggested by Youd et al. [14].
Baziar et al. [26] created two subsets for training, to train a network to predict DH, and a validating phase, to prevent overtraining of the artificial neural network (ANN) model.Then, they presented an ANN model using STATISTICA software (version of Statistica 5.1, Dell Software, Round Rock, Texas, USA) to estimate DH.They inspected the performance of their model using validating subset data without considering extra available models.Furthermore, a new model was presented by Javadi et al. on the basis of genetic programming (GP) [27].They divided the dataset randomly, without paying attention to the statistical properties of the input parameters, into two subsets for the validating and training phase.Garcia et al. established a neuro-fuzzy model to use the advantages of both systems.They randomly separated their dataset into two subsets for training and testing; however, they did not take statistical aspects into account.They also compared the value predicted by their model with extra models to evaluate its performance.Baziar et al. [28] then applied ANN and GP to propose a new model.They divided their dataset randomly into two subsets for the testing and training phase; a validating process was not performed, and the statistical factors of the parameters were not considered.
Although the effects of fine content (Fc) in different values on excess pore water pressure have been investigated [30][31][32][33][34], to the best knowledge of the authors, no attempts have been made to consider the range of Fc to establish the models to predict DH.Most of the studies reveal a range of 20% to 30% for the transition of behavior of the response of sand to earthquake and liquefaction occurrence.Maurer et al. [35] investigated the Canterbury earthquakes in 2010 and 2011 through 7000 case history datasets and illustrated that a high value of Fc caused more inaccuracy in liquefaction assessments.Tao performed some laboratory tests and demonstrated that the potential of liquefaction has a significant dependency on initial relative density (Dr) when the Fc value is larger than 28% [32].
This study is based on the database of Youd et al. [20] and the addition of a new earthquake parameter of CAV5, which is CAV with a 5-cm/sec 2 threshold acceleration, through the attenuation equation presented by Kramer et al. [21].By adding CAV5, the dataset was expanded and became more capable of and efficient in considering aspects of earthquakes and geology site situations, such as earthquake motion frequency, near-fault effects and the causative fault type of an earthquake.The second dataset was created by eliminating samples with an average Fc in a liquefiable soil layer (F15) less than 28%.The response surface method (RSM) is used for the first time as a novel method to develop two equations to predict lateral displacement due to liquefaction (DH) in order to two created datasets herein.Furthermore, the meaningful and effective terms of the equations are discovered through hypothesis testing of the p-value.In this study, two ANNs with back-propagation analysis were developed to measure the coding input data of the RSM.To develop each ANN model, the main dataset is first divided into three subsets for the training, testing, and validating stages, considering statistical properties-instead of random division-to increase the capability and accuracy of the model.To achieve this goal, an attempt is made to create all three subsets with close statistical factors.Finally, the results are compared with data measured from the Chi-Chi earthquake's near fault zone of Wufeng district (Figure 1) and Nantou district (Figure 2), as well as with the predicted DH through three extra models [20,27,36] to demonstrate the accuracy and capability of the RSM models.

Review of Empirical and Semi-Empirical Models
Bartlett and Youd [1], based on factors in References [13,14,29], developed a new model to predict DH due to liquefaction; they supposed that earthquake, topographical, geological, and soil factors are the most influential parameters on DH.They studied 467 displacement vectors from the case history database.Among those vectors, 337 were from the 1964 Niigata and 1983 Nihonkai-Chubu, Japan, earthquakes; 111 were from earthquakes in the United States; and the other 19 cases were selected from Ambraseys' [29] database.In the end, they developed a new model by using multiple linear regression (MLR) for free-face and ground slope conditions [37], but they did not separate earthquakes according to their region because of a database shortage.Youd et al. revised their MLR by adding case history data from three earthquakes (1983 Borah Peak, Idaho; 1989 Loma Prieta; and 1995 Hyogoken-Nanbu (Kobe)), and they considered coarser-grained materials.They removed eight displacement sites with prevented free lateral movement and developed two equations with more accuracy given as follows: Free-face conditions: (1) where * = + (3) and = 10 ( . . ) (4) In Equations ( 1) and (2), DH is the predicted lateral ground displacement (m), Mw is the moment magnitude of the earthquake, and T15 is the cumulative thickness of saturated granular layers (m) with corrected blow counts ((N1)60) less than 15.Moreover, F15 is the average fines content of sediment within T15 (%); D5015 is the average mean grain size for granular materials within T15 (mm); S is the ground slope (%); and W is the free-face ratio (H/L), where H is the height of the free face and L is the distance from the base of the free face to the liquefied point.Finally, r is the nearest horizontal or map distance from the site to the seismic energy source (Km).
Rezania et al. [37] developed a model, based on evolutionary polynomial regression, for the assessment of liquefaction potential and lateral spreading.According to response spectral acceleration, measured through strong-motion attenuation models, Zhang et al. [38] revised the empirical model of Youd et al. [20] and demonstrated the ability of their model by comparing the predicted results with datasets from sites in Turkey and New Zealand [18].Goh et al. proposed multivariate adaptive regression splines (MARS) by using data of Youd et al. [20].They demonstrated an improvement of the original model [39].

Artificial Neural Network
An ANN is a computational process using a biological neural network structure.McCulloch et al. [40] were the first to introduce some simulations according to their neurology knowledge.Neural networks are strong approximators due to their ability to learn by samples and their independency from any algorithm or knowledge about internal features of the issue.Artificial neural networks have a number of advantages such as high accuracy in nonlinear relationships or dynamic mechanisms based on the number and effectiveness of samples.Neural networks are classically constructed in three types of layers.Layers are provided by interconnected nodes.Outlines of the network are developed via the input layer, which communicates to one or more hidden layers by performing a weighting process.Then, hidden layers link to a target (output layer).Further, learning is a supervised process that occurs with each cycle or "epoch" (i.e., each time the network is presented with a new input pattern).
Rumelhart et al. [41] introduced a back-propagation algorithm to decrease error according to the training data.The training process started with random weights to achieve minimum error.The calculation of the derivatives flows backwards through the network, which is why it is called back-propagation.The most common measure of error is the MSE: Overtraining of a neural network happens when the network trains exactly to reply to just one kind of input, which is similar to rote memorization.Therefore, learning does not occur anymore.The ANN is able to be used for problems with non-linear or dynamic correlation.

Response Surface Method
The RSM is a group of statistical and mathematical techniques to develop a capable function for a relationship of response (y) or output variable, and input variables (x), given by: = ( ) + (6) where f(x) is a vector function, consisting of powers and cross-products of input variables.This function depends on the supposed form of the response.
There are some common forms that have been used by researchers.A second-degree polynomial with cross terms is the most complicated and the strictest among them, and it is given by: The main applications of RSM are defined as follows: 1. To present an approximate relationship between input variables and an output variable or response to be able to predict the response variable.2. To discover significant factors or terms of the presented equation using RSM through hypothesis testing such as the p-value.3. To assess the optimization model to obtain a response as a maximum or minimum over a certain range of interest.

Design of Experiments
When more than one input factor is suspected of influencing an output, in order to fit physical or numerical experiments, a process by the name of design of experiments (DOE) [42] is developed.The DOE involves selecting some points according to which a response should be calculated.
In this paper, the design introduced by Box et al. [43] is used.It requires three levels to run an experiment.Furthermore, it is a special three-level design without any points at the vertices of the experiment region.This could be advantageous when the points on the corners of the cube represent level combinations that are prohibitively expensive or impossible to test because of physical process constraints.The design is applied in this study to prevent the input parameters' values from being negative.This DOE requires three levels of each input variable −1, 0, 1 (coded values) corresponding to minimum, middle, and maximum values of input parameters, respectively.

Hypothesis Test
The process in statistics science to meaningfully examine results is called a hypothesis test.During hypothesis tests, the validity of a claim, which is constructed about a population, is evaluated.This claim that is in essence on trial is called the null hypothesis.Hypothesis testing can be expressed in three steps: 1. Defining an initial assumption (null hypothesis).2. Analyzing and assessing sample data by following a formal process.3. Based on the second step, accepting or rejecting the initial assumption in the first step.
One of the main approaches to make a decision to accept or reject a null hypothesis is the p-value, defined through an α value (the probability of error is called alpha).If the p-value is smaller than α, then the null hypothesis is rejected, and contrarily, if it is larger, then the null hypothesis is accepted.
In this paper, the common value of 0.05, which researchers have used for α, is applied to analyze the meaningfulness and significance of parameters and the terms in the RSM response (equation).1. Seismic parameters-moment magnitude (MW) and horizontal distance from site to seismic energy source (r) in km. 2. Topographic parameters-free-face ratio (W) and ground slope (S), both in percent.3. Geotechnical parameters-thickness of layer with corrected blow counts (N1)60 < 15 (T15) in meters, average fines content in the T15 layer (F15) in percent, and average mean grain size in the T15 layer (D5015) in mm.

Model Proposed
Later, Youd et al. [20] eliminated eight sites' data from their dataset due to a lack of free lateral movement.Additionally, they added data from the following three earthquake sites: Borah Peak 1983, Loma Prieta 1989, and Hyogo-Ken Nanbu (Kobe) 1995.
In the present study, the main dataset was developed by adding a new parameter of CAV5, which is defined in Section 5.1.1.Kramer et al. [21] stated that CAV5 is the most efficient and sufficient earthquake intensity to evaluate liquefaction in sandy soil.Therefore, CAV5 was estimated using the attenuation equation presented by them..In this way, the causative earthquake fault types of all earthquakes in the dataset were discovered.Then, the sites with a moment magnitude range from 6.4 to 7.9 were selected in order to find the applicable magnitude range for Equation (11); the Alaska 1964 site, with a magnitude of 9.2, was thus deleted from the dataset.
Furthermore, a statistical analysis was performed to examine and estimate the coefficient of correlation (R) of all input parameters with the output (DH) one by one.The estimated values of R for r-DH and S-DH were a positive value of 0.104 and a negative value of −0.98, respectively, contrary to the supposition for them.This was possibly due to the scarcity of sites that were explored, and consequently, the shortage measured values for r and S in the main dataset.Therefore, in this study, after eliminating r and S from dataset, only the free-face condition was considered, and samples of the ground slope condition were deleted from the main dataset.Furthermore, CAV5 was added to the dataset instead of r.Figures 4 and 5 plot r versus CAV5 for range of Mw in the main dataset from 6.4 to 7 and 7 to 7.9, respectively.It should be mentioned that the bold points show more than one point coincided together.A dataset including 215 case histories with six parameters was then prepared.In addition, to investigate the complicated influence of fine content on the liquefaction, the second dataset was arranged by eliminating samples with an F15 value larger than 28%.Therefore, the second dataset included 182 samples.

Cumulative Absolute Velocity
Eed et al. utilized CAV as a criterion to evaluate the onset of structural damage for the first time.They reported it in the Electric Power Research Institute journal [44] and defined it in a mathematical framework as presented below: where a(t) is the acceleration of ground motion graph, t is the time, and tmax is the duration of the earthquake.Liyanapathirana et al. [45] studied special aspects of Australian earthquakes and found that the predominant frequency of earthquakes in Australia is much higher than in California.This is because the earthquakes in Australia are in the middle of tectonic plates, so liquefaction has not occurred.They introduced a pseudo-velocity in the same dimensions and form as CAV to quantify this difference as indicated below: By studying various Japanese codes, Orense demonstrated that seismic-induced shear stresses, calculated using empirical models containing PGA, do not have high accuracy in near-source circumstances but still have reasonable accuracy for far-source earthquakes.Next to that, he indicated threshold values of 150 gal and 20 kine for PGA and PGV respectively for the occurrence of liquefaction [46].
After inspecting approximately 300 parameters of earthquakes, Kramer et al. revealed that CAV5, as can be estimated through Equation ( 11), has better efficiency and sufficiency than other earthquake intensity parameters for liquefaction evaluation.They also utilized the strong Pacific Earthquake Engineering Research (PEER) database, consisting of 282 ground motions from 40 earthquakes to present an equation to calculate CAV5 for shallow crustal events [21].
where CAV5 is a form of CAV based on Equation (10) (m/sec), M is the moment magnitude, and r is the closest distance to the rupture (km).FN = FR = 0 for strike slip faults, FN = 1 and FR = 0 for normal faults, and FN = 0 and FR = 1 for reverse or reverse-oblique faults.The database of Equation ( 11) has a range of 4.7 to 7.4 for M, which is proposed for use for a maximum magnitude of 8, and it has a range of 1 to more than 100 km for r.In the present study, Equation (11) was applied to calculate CAV5 from the dataset by considering the causative fault type.

Artificial Neural Network Models
At first, both datasets were divided into three groups: approximately 70% for the training phase and two groups of 15% each for the testing and validating phase.The validating phase was applied to prevent the model from being overtrained.The data deviation was conducted according to statistical factors, and the three groups contained similar statistical factors, such as minimum, maximum, and mean values.Furthermore, to achieve a higher accuracy of output, the same portion of any earthquake's data was selected for the three phases, so all the sites contributed with their data in the training, testing, and validating phases with similar statistical factors.The ANN was established using a back-propagation algorithm with one hidden layer, which is the most commonly used network to drive an ANN model to predict DH.There are six inputs, including six parameters, which were considered by Youd et al. [20] and Bartlett et al. [37].In addition, there is the new measured parameter of CAV5 as well as one output as a DH.
The first dataset was divided into three groups including 151 samples for training, 32 samples for testing, and an extra 32 samples for the validating phase.Tables 1 and 2 summarize the characteristics and certificates of the first dataset and the subsequently developed ANN model.As can be seen from Table 2, the coefficient of correlation (R) given in Equation (12), which is the most common factor to assess the performance of correlations, of the first ANN model for all three groups and the main dataset was around 90%.
where n is the sample size, xi and di are the individual sample points indexed with i, and ̅ and ̅ are the mean sample sizes.To investigate the complex influence of F15, the new dataset was constructed by selecting data samples with an F15 less than critical value of 28%, which was demonstrated by Tao [32], and its characteristics are listed in Table 3.The new dataset was divided again into three groups with the same portion of each site and with similar statistical factors.Therefore, data from each earthquake contributed to the training, testing, and validating phase.In this step, around 15%, 15%, and 70% of the dataset equated to 27, 27, and 129 samples used for the testing, validating, and training processes, respectively.The characteristics of this new ANN model are presented in Table 3.Also, Table 4 presents the certificate of the second model.It can be seen in Table 4 that the R values for all groups of datasets were around 90%.

The RSM Equations for Predicting DH
First, an RSM was conducted to drive an equation on the basis of the first ANN model.Therefore, the DOE introduced by Box et al. [43] for the second-degree polynomial with cross terms was employed to provide 54 coded values to cover the full range of parameters.Through the first developed ANN model, the response value (herein referred to as DH) in coded form was calculated.Thereafter, the RSM equation with 28 terms according to the second-degree polynomial with cross terms was derived; however, through hypothesis testing considering the p-value, some terms were eliminated.Then, the RSM was applied repeatedly to achieve the final equation.The following equation was consequently developed with 22 terms to correlate the DH caused by liquefaction to the six input parameters for the whole range of parameters in this study, without any limitations on the range of the F15 value: Characteristics of the RSM equation: R 2 = 87.22%,R 2 (predicted) = 78.73%,and R 2 (adjust) = 83.99%.
The coefficient of determination (R 2 ) illustrates how well the curves fit on the data points.In addition, the R 2 (adjust) demonstrates the percentage of variation defined by the independent variables that affect the dependent variable, herein referred to as DH.Also, the R 2 (predicted) defines how well a correlation is able to predict the target for a new observation.The values of coefficients a0 to a21 are listed in Table 5.Second, the second RSM equation with the final 21 terms was derived based on the second developed ANN model in this study through the same process as that for the first RSM equation.The second-degree polynomial with cross terms was applied in conjunction with Box and Behnken's DOE on the basis of the second ANN model.Then, a hypothesis test with the same p-value was conducted to provide the second RSM equation for F15 less than 28% (Equation ( 14)).Table 6 presents the coefficients of the second RSM equation.: Characteristics of the second RSM equation: R 2 = 88.51%,R 2 (predicted) = 50.95%,and R 2 (adjust) = 78.09%It must be noted that to use the derived RSM equation, the main values of input parameters must be exchanged with coded values by the function presented in Equation (15).Those coded values must then be put into the RSM equation to achieve the results.In other words, the RSM equation declares the relationship between the coded value of parameters and responses (output); coded values have a range from −1 to 1.
The max, min, and mean values of the parameters are listed in Tables 1 and 3 for both RSM equations presented in this study.

Comparison of RSM Equations with Extra Models
Chu et al. [47] analyzed five liquefied sites during the Chi-Chi-1999 earthquake in Taiwan, all in the near-fault region from five sites in two districts of Wufeng and Natu as they are illustrated in Figures 1 and 2. Table 7 presents these parameters' values from the sites.Based on these data samples, the predicted results of the RSM equations are compared to Youd et al. [20], Javadi et Al.[27], and Rezania et al. [36].In total, 28 sites (for which the necessary parameters for the ANN model and the RSM equation are reported by Chu et.al [45]) are illustrated in Table 8.The sample numbers from 1 to 26 are from Wufeng's sites and samples of 27 and 28 are belong to Nantu's site.The capability and accuracy of the first RSM equation was demonstrated using all 28 site samples from the Chi-Chi earthquake, including samples with F15 from 13% to 48.5%.
where N is the number of samples, Xm is the measured value, and XP is the predicted value.Furthermore, all samples with an F15 greater than 28% were eliminated from the Chi-Chi earthquake cases, and 16 samples consequently remained (samples number 1 to 14 as well as numbers 27 and 28, as can be seen in Table 7).Then, the second RSM equation was validated by applying it to these samples in comparison with the extra three models.Table 9 summarizes the results of all models.

Results and Discussion
The previous sections have compared the first RSM model, which belongs to the full range of parameters, and the second RSM model, which was derived for samples whose F15 values were less than 28%, with three extra well-known models.The models were examined using new data from the Chi-Chi earthquake, which were not included in the two datasets to establish the two RSM models.As can be seen in Table 8, the RSM equation of the whole range of parameters indicated a higher R value of 0. On the other hand, as can be seen in Table 9, by considering samples with F15 less than 28%, the model of Youd et al. provided the highest R value of 0.934, closely followed by the second and the first RSM models with R values equal to 0.891 and 0.846 respectively.Table 9  The comparison between the two models developed in this study and the extra three models demonstrates that the second RSM model provided a reasonable correlation and the lowest error.The results indicated that the RSM is a highly efficient tool to perform a liquefaction hazard analysis.Furthermore, performance of the model is increased by taking into account the complex influence of Fc by eliminating an F15 larger than 28% and even by decreasing the number of samples in the dataset.
Another major advantage of the presented models is their consideration of earthquake aspects, such as the near-fault zone, the frequency of earthquake motion, and the causative fault type, by estimating and adding the CAV5 parameter to the dataset.As can be seen from Figures 6 and 7 There are some limitations for applying both first and second RSM equations as follow: 1) Both RSM models require standard penetration test SPT and laboratory tests to determine geotechnical properties parameters of T15, F15, and D5015.2) Both of the RSM models are valid for free-face conditions but not ground-slope conditions.
3) Second RSM model is valid only for F15 < 28%.4) Models are only valid for earthquakes with Mw between 6.4 and 8.0.5) Specify accuracy limits for each model.6) It is necessary to transfer all six input models' parameters measured value to a coded value using Equation ( 15) and then put the coded value in the RSM equations to predict DH.

Summary and Conclusions
The determination of lateral displacement due to liquefaction caused by an earthquake (DH) is the most important aspect of liquefaction hazard analysis.There are two main types of conditions according to the topography of the sites: free-face and sloping ground conditions.First, the parameter of corrected absolute velocity (CAV5) of sites was calculated due to it being the most efficient and sufficient parameter for the assessment of liquefaction caused by earthquakes [21], and it was added to develop the dataset to cover all aspects of earthquakes, including the frequency content of earthquake motions and the causative fault type of earthquakes.Then, a statistical parametric analysis was performed by estimating the correlation coefficient (R) between all input parameters and output as DH.To achieve a more capable and accurate model, based on the estimated values for R, the horizontal distance from a site to the seismic energy source (r) and ground slope (S) was eliminated from the original dataset due to poor correlations to the target.Therefore, the final dataset was created for free-face condition sites.
The significant aspects of earthquakes, such as the near-fault region, frequency content, and causative fault type of earthquakes, which are included in the model established by Kramer et al. [21], were considered by taking CAV5 into account.To investigate the complex effect of fine content, the main dataset was divided into two subsets.The first dataset included the whole range of parameters, and in the second one, all samples with average fine content in the liquefiable layer (F15) larger than 28% were removed from the dataset, in line with Tao [32].Furthermore, the RSM was applied to develop two equations in order to the first and the second dataset to examine its performance to assess liquefaction.In the end, the two presented models in this study were compared to three available models to demonstrate their capability and accuracy with regard to predicting DH in free-face conditions in a near fault zone case history of the Chi-Chi earthquake.
The present study highlights the importance of earthquake aspects, especially CAV5 as the most sufficient and efficient intensity to liquefaction hazard assessments.In addition, the RSM is a strong tool for the evaluation of complex non-linear phenomena such as liquefaction.
The results also confirm the complicated influence of F15 on the whole range, and they provide significant enhancements to the performance of the model by considering samples with an F15 less than 28% as a critical value defined by Tao [32].One of the most remarkable results, which shows the complex influence of fine content on evaluation of DH, is that the second model demonstrated higher accuracy and capability, even though it was developed using a database with fewer samples than the first model.

Figure 1 .
Figure 1.Location of lateral displacement and ground failure due to liquefaction during the Chi-Chi earthquake in the Wufeng district.

Figure 2 .
Figure 2. Location of lateral displacement and ground failure due to liquefaction during the Chi-Chi earthquake in the Nantou district.

Figure 3 Figure 3 .
Figure3illustrates the flow chart of the applied approach to develop the RSM equations to estimate DH in this study.Two models are presented: first one considered the whole range of the parameters and the second one was on the basis of the F15 value being less than 18%.

Figure 4 .
Figure 4. r versus CAV5 for 6.4≤Mw<7 of the main dataset applied in this study.

Figure 5 .
Figure 5. r versus CAV5 for 7≤Mw≤7.9 of the main dataset applied in this study.

Figures 6 and 7
Figures 6 and 7 visualize the comparison between both RSM equations developed in the present study and three extra models with measured data from sites of the Chi-Chi earthquake.Twenty-eight data points are evaluated in Figure6for whole range of the parameters.Meanwhile, Figure7illustrates the comparison for data points with F15 values of less than 28% at 16 data points.

Figure 6 .
Figure 6.Comparison between the first RSM equation and three extra models with 28 data points measured from sites of Chi-Chi earthquake.

Figure 7 .
Figure 7.Comparison between both RSM equations and three extra models with 16 data points including F15 of less than 28% measured from sites of the Chi-Chi earthquake.
683, in comparison with the extra models whose values were 0.433, −0.74, and 0.514.Furthermore, the RSM model comprises lower MAE and RSME values of 0.3 and 0.37, respectively, compared to 0.49 and 0.7 for Rezania et al., 1.04 and 1.19 for Javadi et al., and 3.77 and 4.37 for Youd et al.Therefore, among all of them, the RSM model provided prediction with higher accuracy.
also illustrates the MAE and RSME criteria values for all models for samples with a limited value for F15 less than 28%.The values of the MAE and RSME in the second RSM model-0.29 and 0.39, respectively-indicates the highest accuracy and performance in comparison with the others.In addition, the model of Rezania et al., with 0.42 and 0.57, illustrated lower values for the MAE and RSME, respectively.Further, Javadi et al. with 1.2 and 1.3, and Youd et al. with 4.84 and 5.34 provide less accuracy for predicting DH.
, among all models that were considered in the present study to calculate DH without any limitation on the parameters' value, the model of Youd et al. was overpredicted.Meanwhile, Youd et al.'s model provided poor and overpredicted results for samples with a limited value of F15 less than 28%.Additionally, considering samples with a limited F15 value shows the first RSM and the model of Javadi et al. present an overpredicted value for DH.Furthermore, second RSM equation and model of Rezania et al. underpredicted DH in their predictions.

Funding:
This research was funded by National Natural Sciences Foundation of China, grant number 51639002 and National Key Research & Development Plan, grant number 2018YFC1505305 and "The APC was funded by State Key Laboratory of Coastal and Offshore Engineering, Dalian University of Technology, Dalian 116024, China.Aknowledgements: The authors are grateful for the technical and financial support Provided by National Natural Sciences Foundation of China Granted No. 51639002 and National Key Research & Development Plan

Table 1 .
Characteristics of whole case histories' input parameters that were used to develop the ANN model and RSM equation.

Table 2 .
Certificates of the first ANN model for the whole dataset.

Table 3 .
Characteristics of database with F15 ≤ 28% used for the second developed ANN model and RSM equation.

Table 4 .
Certificate of the second ANN model for dataset with F15 ≤ 28%.

Table 7 .
Model parameters and measured DH at sites affected by the 1999 Chi-Chi earthquake.The results of comparison between the predicted values and measured values are summarized in terms of the root mean square error (RMSE), mean absolute error (MAE), and R in Table8.It is clear that the larger R and smaller RMSE and MAE reveal higher accuracy of predicted results.

Table 8 .
Performance certificate of first RSM equation in comparison with extra available models.

Table 9 .
Performance certificate of second RSM equation, on the basis of samples with F15 ≤ 28% in comparison with extra available models.