Seismic reliability analysis of soil slope based on Newmark sliding block model with representative slip surfaces and response surfaces

This paper presents a framework combining Monte Carlo Simulation (MCS) and the Newmark sliding block model with representative slip surfaces (RSS) (model II) and multiple response surfaces method (MRSM) to conduct seismic reliability analysis and risk assessment of soil slopes. An empirical threshold is introduced to define the limit state function to identify the failure samples in MCS and the sliding area and Newmark sliding displacement are multiplied to quantify the failure consequence. The proposed methodology is illustrated through a soil slope with multiple layers. The calculation results demonstrate that traditional Newmark sliding block model (model I) tend to underestimate the variations of yield acceleration. Both the failure probability and landslide risk exhibit decreasing trends with the increase of threshold. Significant discrepancy in failure probability and landslide risk between two models is found even for a small threshold. Therefore, the proposed methodology is highly recommended in seismic reliability analysis and risk assessment. The contributions of RSSs to the failure probability and landslide risk are insensitive to the variation of displacement thresholds.


YU et al.
Wang et al. 12 proposes a new analysis method for random dynamic response analysis and permanent displacement reliability evaluation of soil slope under mainshock-aftershock sequence.Yin et al. 13 quantitatively analyzed the influence of aftershock sequence on permanent displacement of slope.5][16] Consequently, these uncertainties are properly included in seismic reliability analysis.
7][28][29][30] However, the MCS requires a large number of executions of deterministic model to obtain the rational result at small target failure probability, and the calculation efficiency deteriorates dramatically.In order to improve the computational efficiency of MCS simulation, surrogate models have been developed to establish the relationship between the random variables and the response of geotechnical structures (e.g., the relationship between soil strength parameters and slope sliding displacement).The developed surrogate models within the literature include but not limited to response surface method (RSM), [31][32][33][34][35] artificial neural network (ANN), 36,37 support vector machine (SVM), 38,39 multivariate adaptive regression splines (MARS), 40 random forest (RF), [41][42][43] convolutional neural network (CNN).Response surface method is widely used in slope reliability analysis due to its simplicity.
In the seismic reliability of slopes, the Newmark sliding block model should be updated to consider the effect of uncertainty of soil properties on the yield acceleration k y by the consideration of representative slip surfaces.To enhance the computational efficiency of MCS, a response surface function is calibrated for each of representative slip surfaces (RSS) to predict the relationship between soil properties and k y .
This paper starts with the brief review of Newmark sliding block model with RSS, followed by the calibrations of multiple response surfaces.The seismic reliability analysis and risk assessment of slopes combining MRS are delineated.The proposed method is applied to a subgrade fill slope.The variations of k y s, failure probability and risk of the slope are fully analyzed, and finally, conclusions are drawn.

2.1
The traditional Newmark sliding block model (model I) The traditional Newmark sliding block model regards the sliding body as a rigid body and it is assumed that the sliding occurs only when the horizontal seismic acceleration a h exceeds the yield acceleration k y of the sliding body.The horizontal seismic acceleration a h , which brings the FS of the sliding body to be 1, is determined to be the yield acceleration k y .The cumulative Newmark sliding displacement, denoted by D, can be determined through double integrations over the time series when the sliding arises.

Brief review on the determination of RSS
The proposed method in Refs.33,44 is adopted in this paper to locate the RSS of slope.A brief review is present herein to make the Newmark sliding block model with RSS readable.The main idea for identifying the RSS is that among an ensemble of highly correlated potential slip surfaces, the slip surface with minimum reliability index  can be selected to represent the rest of slip surfaces.The mean value first-order second moment method 45 is used to calculate the reliability index  of slip surface.The RSS is identified using the following procedure.First, a large number of potential slip surfaces are randomly generated to cover the possible failure domain and the slip surface with minimum  are selected to be the first RSS.Second, the correlation coefficients between the current RSS and each of the residual slip surfaces.Third, the slip surfaces which are highly correlated with the current RSS are excluded from the further consideration.Finally, the next slip surface with minimum  is selected to be the next RSS and the previous three steps are repeated until all the RSS are distinguished.The readers can refer to Refs.33,44 for more details.

Newmark sliding block model with RSS (model II)
The implementation of the Newmark sliding block model with RSS is demonstrated in Figure 1, which consists of the following four steps.
1. Following the procedures described in Section 2.2, M representative sliding surfaces of the slope are selected and denoted as S 1 , S 2 , … , S M ; 2. The yield acceleration k yi , i = 1, 2, … , M, of the sliding body corresponding to S i is calculated respectively using limit equilibrium pseudo-static method; 3. The minimum value k yb is selected from k yi , i = 1, 2, … , M; 4. The sliding body corresponding to S b and k yb is used to calculate the sliding displacement of slope by traditional Newmark sliding block model.

MULTIPLE RESPONSE SURFACES METHOD (MRSM)
Since the determination of k y for a RSS is an iterative process, in order to improve the computational efficiency, a response surface function is calibrated to predict the relationship between the k y for a RSS with soil properties such as cohesion and internal friction angle.For simplicity, a second-order polynomial functions without cross terms is adopted.
Where X = (x1, x2, … , xm) represents the random variables in seismic reliability analysis of slope; K yi (X) is the yield acceleration response surface function corresponding to the ith RSS, S i , i = 1, 2, … , M; a i, is the constants,  = 0, 1, 2, 3, … , 2m.The details of the calibration of the response surface function can be referred to Refs.28,31.The calibrated K yi (X), i = 1, 2, … , M, are assembled to predict the kys corresponding to M RSS and finally to locate the k y of slope.This ensemble of K yi (X) is named multiple response surfaces in this paper.It is noted that the calibrated response surface function must have sufficiently large accuracy before it can be used in the prediction of k y .The first part aims to determine the distributions of random variables used in reliability analysis based on sufficiently large number of laboratory or in-site test data.The second part is to calculate the Newmark sliding displacement in an efficient manner.In the MCS part, N random samples are randomly generated in accordance with the characterized distributions and statistics in first part and they are denoted by X j , j = 1, 2, … , N. Consider X j for example, the model II and MRSM is adopted to locate the k y of the slope, to determine the corresponding sliding area A(X j ), and to calculate the Newmark sliding displacement D(X j ).Furtherly, If D(X j ) is greater than a displacement threshold − D, X j is a failure sample.Hence, the failure consequence related to X j , denoted by C j , equals D(X j ) × A(X j ) and the respective landslide risk R j is C j /N.otherwise, if D(X j ) is lower than − D, X j is a stable sample implying that C j and R j are zero.After N random samples have been processed in the Newmark sliding block model with RSS and MRSM, the total number of failure samples, denoted by N f , can be counted.The failure probability p f is the ratio of N f to N, and the landslide risk is the summation of R j , j = 1, 2, … , N. It must be noted that the proper choice of

Input parameters
The subgrade fill slope model is shown in Figure 3.The slope height H is 6 m and the slope angle is 18.4 • (equivalent to an inclination ratio of 1:3).The thickness of each of three soil layers is 4.5 m.The seismic stability of the slope under undrained conditions is investigated.The undrained shear strength Su and the saturated unit weight  sat are the necessary soil properties.The identical statistics in Li and Chu 33 are shown in Table 1.The first column is regarding the soil layer.The second to fifth columns are the mean value, standard deviation, and COV of Su.The last column is for distribution type of Su.EI-Centro wave is selected to be the seismic ground acceleration record.

Variation of p f
To investigate the influence of  Therefore, if a smaller threshold value is acceptable, the p f s from the two models agree fairly well with each other.Otherwise, if a greater threshold value (e.g., 0.4 m) is used, significant discrepancy (larger than 50%) will arise and the model II is highly recommended.The N f failure samples are reanalyzed.For a given failure sample, if the k y of the slope is attributed to the ith RSS(S i ), an indicator variable N i = N i + 1. N i is finally determined by checking through N f failure samples.The contribution of S i to slope p f , denoted by  i , is defined as the ratio of N i to N f .Figure 6 plots the contribution of RSS to slope p f with − D. In Figure 6, the squares are regarding the contribution of S 1 (i.e.,  1 ), the circles are regarding the contribution of S 2 (i.e.,  2 ), and the triangulars are regarding the contribution of S 3 (i.e.,  3 ).It is noticed from Figure 6 that as − D ranges between 0.1 m and 10 m,  3 varies slightly at 60%, whereas  2 fluctuates insignificantly at 40%, and  1 is negligible.This indicates that the effect of − D on the contributions of RSS to slope p f can be neglected.

Variation of landslide risk
In addition to p f , the failure consequences corresponding to N f failure samples are determined using the product of sliding area and D. The sliding area of S 1 , S 2 , and S   5B and 7B demonstrates that although the negligible discrepancy in p f between two models is observed, significant difference in R between two models has been noticed at even a small − D (e.g., m).Hence, the model II is highly recommended for seismic slope reliability analysis and risk assessment.
Like the contribution of S i to slope p f , for a given failure sample, if the k y of the slope is attributed to the ith RSS(S i ), the corresponding risk is added to an indicator variable r i .Let  i = r i /R denote the contribution of S i to risk R. Figure 8 plots the contribution of RSS to risk R with − D. In Figure 8, the squares are regarding the contribution of S 1 (i.e.,  1 ), the circles are regarding the contribution of S 2 (i.e.,  2 ), and the triangulars are regarding the contribution of S 3 (i.e.,  3 ).It is noticed from Figure 8 that as − D increases from 0.1 m and 10 m,  3 varies very slightly at 75%, whereas  2 varies very slightly at 25%, and  1 is trivial.This indicates that the effect of − D on the contributions of RSS to R can be neglected.

CONCLUSION
The necessity of conducting seismic reliability analysis of slopes using Newmark sliding displacement is warranted since the uncertainty in soil properties is inevitable.Following the methodology using MCS to conduct reliability analysis, a response surface function is calibrated to predict the yield acceleration k y of slope sliding along a given representative slip surface.The results of the case studies demonstrated that model I tends to underestimate the variations of k y .The p f exhibits a decreasing trend as the displacement thresholds increase.For a specific threshold, model I tends to underestimate the p f as opposed to model II.The discrepancy in p f between these two models shows a nonlinear growing trend with the increase of displacement threshold.If a smaller threshold (lower than 0.4 m) is adopted, the discrepancy is limited to 50%, otherwise dramatically large discrepancy is inevitable.The variation trend of R is similar to that of p f .Significant discrepancy in R is found even at a small threshold (e.g., 0.1 m).The contributions of RSSs to slope failure probability p f and landslide risk R are insensitive to the choice of displacement thresholds.

F I G U R E 1
Newmark sliding block model with RSS (model II).

Figure 2
Figure 2 demonstrates the flowchart of seismic reliability analysis and risk assessment of slope by MCS combining model II.The methodology comprises three parts, that is, (1) inherent uncertainty characterization of soil strength parameters, (2) the implementation of model II, (3) MCS.The first part aims to determine the distributions of random variables used in reliability analysis based on sufficiently large number of laboratory or in-site test data.The second part is to calculate the Newmark sliding displacement in an efficient manner.In the MCS part, N random samples are randomly generated in accordance with the characterized distributions and statistics in first part and they are denoted by X j , j = 1, 2, … , N. Consider X j for example, the model II and MRSM is adopted to locate the k y of the slope, to determine the corresponding sliding area A(X j ), and to calculate the Newmark sliding displacement D(X j ).Furtherly, If D(X j ) is greater than a displacement thresh-

−D
is crucial and also very complicated, which requires comprehensive consideration of risk acceptance degree, slope scale and other factors.A set of − D values have been selected to demonstrate the influence of − D on the results.F I G U R E 2 Flowchart of seismic reliability analysis and risk assessment of slope combining model II and MRSM.

Figure 4 F I G U R E 3 F I G U R E 4
Figure 4 plots the histogram and cumulative density curve of the Newmark sliding displacements calculated by model I and model II.In Figure 4A, the results calculated by model I are plotted in blocks with vertical lines while the results calculated by model II are plotted in blocks with inclined lines.10,000 Ds are classified into 10 sub-ranges.It can be seen from Figure 4A that the dominant range is first range in both models.The frequency of Ds in the first range is 84.95% in model II, whereas it is 94.66% in model I.The further plot in the right upper hand demonstrates that most of the Ds in the dominant range are lower than 0.7 m.The mean and standard deviation of 10,000 Ds are 1.79 m and 5.60 m (with a coefficient of variation = 5.60/1.79≈ 3.1), respectively in model II, while they are 0.54 m and 1.94 m (with a coefficient of variation =1.94/0.54≈ 3.6) in model I.The cumulative density curve of Ds is plotted in Figure 4B, where the results from model I are plotted in triangulars while those from model II are plotted in stars.Two indices, that is, the coefficient of non-uniformity C u and coefficient of curvature C c are calculated.Where, d 60 , d 30 and d 10 are the displacement values corresponding to the cumulative probability of being less than a certain displacement value of 60%, 30% and 10%, which are respectively called limited displacement, median displacement and effective displacement.It is noticed from Figure 4B

−D
on seismic reliability, a set of different values are selected form the range between 0.1 to 10 m. Figure 5A plots the variation of p f versus − D, where the results from model I are plotted in triangulars while the results from model II are plotted in stars.It is observed from Figure 5A that as − D increases from 0.1 m to 10 m, the p f calculated by model II decreases from 84.33% to 3.84%, whereas the p f by model I decreases from 71.44% to 0.51%.For a given displacement threshold, the p f obtained by the model II is greater than that from model I.

F I G U R E 6
Variation of  i with displacement threshold.The respective p f s calculated by model I and II are denoted by p 1 and p 2 .The relative discrepancy between p 2 and p 1 is quantified using  = [(p 2 −p 1 )/p 1 ] × 100%.5B plots the variation of  versus − D, where the results are plotted in triangulars.It is seen from Figure 5B that  is 18%, 182%, and 653% at − D = 0.1 m, 1 m, and 10 m, respectively showing a non-linear varying trend as − D increases from 0.1 m to 10 m.

3 are A 1 =
10.22 m 2 , A 2 = 35.49m 2 and A 3 = 86.55m 2 , respectively.Figure 7A plots the variation of R versus − D. In Figure 7A, the results from model I are plotted in triangulars while the results from model II are plotted in stars.As shown in Figure 7A, at − D = 0.1, 1.0, and 10 m, the calculated Rs by model II are 113.31,97.51, and 61.34 m 3 , while they are 45.81, 27.93, and 9.10 m 3 from model I, respectively.As − D increases from 0.1 m to 10 m, R tends to decreases in the two models.This variation trend of R with − D is similar to that of p f with − D. At a specified − D, the R from model I is significantly lower than that from model II.This observation can be attributed to the difference of 10,000 Ds between the two models as described in section 5.2.On the one hand, the model II predicts larger Ds than the other model leading to a higher p f under a given − D. On the other hand, the failure consequence corresponding to the same failure sample is magnified by the larger D.

F I G U R E 7
Variation of risk R and  with displacement threshold.(A) R. (B) .F I G U R E 8 Variation of  i with − D.Let R 1 and R 2 denote the Rs calculated by model I and II.The relative discrepancy between R 2 and R 1 is quantified by = [(R 2 −R 1 )/R 1 ] × 100%.Figure 7B plots the variation of  versus − D. It is noticed from Figure 7B that  is 147%, 249%, and 574% at − D = 0.1, 1, and 10 m, respectively exhibiting a nonlinear varying trend, which is similar to that ofin Figure 5B.The comparison between Figures