Reliability-Based Stability Analysis of a Baltic Cliff by the Combined Response Surface Method

The study presents a probabilistic stability analysis of a Baltic cliff in Jastrzębia Góra, Poland. Progressive slope erosion is a threat to adjacent buildings, so safety assessment of the slope is essential. The cliff shows a compound, multi-layered geological structure, which makes the analysis of its reliability a complex multivariate problem. A simple, straightforward computational procedure was proposed, incorporating the Response Surface Method (RSM) linked with the standard Monte Carlo (MC) simulation method and the Point Estimate Method (PEM). PEM samples make it possible to analyse the sensitivity of the cliff’s stability to variation in subsoil parameters and to reduce the number of random variables of the problem. The proposed methods were tested in two cases: a high failure probability (undrained state) and a moderate failure probability (drained case). The so-called combined Response Surface Method (CRSM) proposed here may be successfully applied in geotechnical computations characterized by dispersion and uncertainty of soil data as well as a relatively high damage probability.

A milestone in probabilistic analysis of slopes was the application of the Finite Element Method (FEM) in geotechnical problems. Various FEM-based probabilistic methods have been proposed. Developments in both hardware and software have made the MC simulation dominant. The FEM-based MC application is usually called the Random Finite Element Method (RFEM), which fully incorporate spatial correlation and averaging (Griffiths and Fenton 2004). This method has been applied to slope reliability (Griffiths et al. 2009;Cho 2007;Huang et al. 2010Huang et al. , 2013Wang et al. 2011;Li et al. 2014;Farah et al. 2015;Liu et al. 2015). In the case of complicated slopes, the application of RFEM involves extensive computational effort. In order to improve the computational efficiency of slope reliability assessment, research was conducted to simplify the slope stability performance function and to reduce the sample for Monte Carlo simulations.
A method addressing slope stability problems is the Response Surface Method (RSM) (see Wong 1985;Cho 2010;Huang et al. 2010;Li et al. 2011Li et al. , 2013Tan et al. 2013Tan et al. , 2016Ji and Low 2012;Zhang et al. 2011Zhang et al. , 2013Jiang et al. 2014;Li et al. 2016;Zhou and Huang 2018). RSM is usually applied to relatively simple cases; it is difficult to find a literature example of its application to multi-layered slopes with more than a few uncertain soil parameters.
The Point Estimate Method (PEM) proposed by Rosenblueth (1975) and its modified procedures have also attracted interest in many engineering fields. Various authors apply this method in probabilistic slope design, e.g. Gibson (2011), Wang and Huang (2012), Ahmadabadi and Poisel (2015).
Solving the same analytical problem by different probabilistic approaches may often produce diverse results. The differences are even clearer in non-linear cases. On the other hand, confirming the results by at least two methods is required to obtain a better insight into the problem. An MC simulation is usually a reference here. Moreover, random variable and random field generation by the MC method is relatively easy and widely used in great many engineering applications. Thus it seems reasonable to combine traditional MC sampling with other algorithms, such as PEM or RSM, to accelerate computations and verification of computational results.
The paper aims at a directed approach to RSM, the so-called combined Response Surface Method (CRSM), which is based on simple random methods (MC, PEM and RSM) and combines their major advantages. Note that CRSM is not universal; it is dedicated to cases of high or moderate failure probability. Thus, it is particularly suitable for analyses of landslides and soils with a complex structure.
The applicability of CRSM is discussed here with regard to a stability study of a Baltic cliff in Jastrzębia Góra, Poland. The cliff has a complex multi-layered geological structure with over a dozen uncertain soil parameters. Complex investigations of the critical section of the cliff (approximately 750 m long) conducted over the years have been addressed in numerous works and reports, e.g. Tejchman et al. (1995a), Zabuski and Korzec (2017). Some works, e.g. Tejchman et al. (1995b), have been presented at conferences, and some include a probabilistic stability analysis of the cliff, e.g. Tejchman et al. (1996). All these publications provide a background for the indepth analysis presented here. It is also worth noting that a damaged building is situated on the crest of the cliff (in the landslide-prone zone), and numerous attempts have recently been made to repair it and restore it to use. The paper complements a feasibility study of this project, and therefore the computations regard both the present situation (the case of high failure probability) and the hypothetical post-drainage situation (the case of moderate failure probability).

Combined Response Surface Method (CRSM)
The proposed computational CRSM algorithm draws on three popular reliability assessment methods: Monte Carlo simulation (MC), Point Estimate Method (PEM) and Response Surface Method (RSM).
The direct MC approach estimates the probability of failure P f in the entire domain of x by means of the expression applying the so-called indicator function I 0=1 where LSðxÞ is the limit state function, and NS denotes the number of samples used in the approximation. Cornell's reliability index b (Cornell 1971) is estimated in the form in which UðÁÞ denotes the standard Gaussian cumulative distribution function.
The main idea behind the PEM is to properly discretize the continuous random variable, see Rosenblueth (1975). In the case of n-dimensional random vector x s ¼ x i f g; i ¼ 1; 2; 3; . . .; n, the structural response function is given as where values of the response function yðxÞ are computed at relevant combinations of PEM discretization points. For example, two combinations related to the variable x i defined by the mean value l i and standard deviation r i take the form In turn, the mean value l y and the standard deviation r y of the response function yðxÞ are defined as Hence, the reliability index b reads b ¼ l y r y ð8Þ and the probability of failure P f is obtained by inversion of Eq. (3). The aim of RSM is to determine a function to approximate the actual structural response induced by random variation of input parameters. The response surfaceŷðxÞ maps a set of discrete values of actual response yðxÞ corresponding to a finite number of discrete realizations of the random vector x: where x 1 ; x 2 ; x 3 ; . . . ; x n are realizations of basic random variables, and e is the assessment error of the actual structural response. If the response is approximated in a small subspace only or its detected skewness is low, the first-order response surface (RS) model may be applied. On the other hand, cases of spacewise response mapping, large curvature of the actual response or an inaccurate first-order RS model approximation require higherorder RS models, e.g. the second-order model, defined aŝ In order to assess surface slope factors b ðÁÞ given in Eq. (10), standard approximation techniques are commonly used, e.g. the least square method or the ANOVA tabular variance reduction technique, see Montgomery (1997). On the basis of the surface slope factors determined in RSM-Win (Winkelmann and Górski 2014), the reliability index b is estimated as the shortest distance from the system origin to the socalled design point x Ã on the failure surface in the reduced (standardized) space (Hasofer and Lind 1974). In general, MC results are usually used as reference values, especially in non-linear and multidimensional random variable cases. However, this method requires multi-sample computations and may thus become extremely inefficient. Computational effort, even in multidimensional cases, is significantly reduced by PEM. Unfortunately, due to the small number of samples used in PEM, its accuracy is frequently unsatisfactory, therefore it is recommended to verify its results by other methods. On the other hand, PEM samples may be successfully introduced in the sensitivity assessment of structural response to variations in input parameters.
The idea behind CRSM is to combine the advantages of the three abovementioned methods: MC, PEM and RSM. Samples generated by the MC method and PEM may be directly used in surface approximation by RSM. CRSM enables parallel comparison of the MC, PEM and RSM results. Their convergence is checked at every approximation step, so that CRSM computations can be terminated as soon as the required convergence level has been achieved. Thus, a CRSM solution is obtained in a relatively short time.
MC computations in the course of CRSM are not independent, since this would involve a high computational expenditure due to a large sample space. MC samples are generated primarily to be used in the approximation of further, higher-accuracy RSM functions. However, MC-generated samples may be distant from the limit state boundary, so surface approximation by RSM is not precise in this region. Hence CRSM is dedicated to cases of moderate failure probability. The introduction of MC samples imposes limitations on the methodology, but this approach is highly accessible, requiring neither a vast theoretical background on probability, nor specialized, advanced software. Such requirements are inherent to the target sampling methods that iteratively search for the limit state boundary. CRSM seems simple and straightforward to use in engineering applications, but its simplicity does not exclude non-linear, multi-variable unknown limit state functions.
CRSM is applied here to analyse the stability of a Baltic cliff in a region where progressive erosion threatens a building situated on the crest of the cliff.

The Cliff in Jastrzębia Góra
The Jastrzebia Góra cliff ( Fig. 1) is a several-kilometre-long section of the Polish Baltic cliff coast, which extends over 100 km. The cross-section of the cliff, descending steeply into the sea, rises 31 m above sea level (Fig. 2). The complex structure of this cliff stretch is dominated by fluvioglacial sediments prone to active landslides. In this case, the initial layers of the cliff, composed mainly of clays, induce the sliding of the other beds.
Investigations of coastline movement started at the end of the 19th century. They were initially carried out using historical topographic maps, then with the use of aerial photography and finally, in recent years, by terrestrial laser scanning (TLS), capable of the most accurate monitoring of movement. According to longterm observations, the average landward movement of the coastline between 1875 and 1976 was 0.3 m/year, causing the cliff coastline to move 80 m landwards. Later, in the years 1977-1990, the movement accelerated substantially to 0.94 m/year. Since 1991, the average velocity of displacements has decreased to 0.5 m/year. This slowdown has been achieved by erecting a gabion barrier at the base of the cliff, at the sea level (Fig. 2). In spite of that, several small landslides have developed. In 2006, one of the landslides caused the collapse of a seaward wall of the building (Fig. 1). Therefore, remedial works were conducted to replace the colluvium with geogridreinforced gravel and erect a shelf-like structure between the building and the cliff's edge (Fig. 2).
An extensive research programme was conducted (Tejchman et al. 1995b) in the most active section of the cliff. Based on deep borings in the crest of the cliff, up to a depth of 30 m, three cross-sections were identified in a 750 m range. All boreholes were also used to install piezometers and inclinometers. Physical and strength parameters of particular soils composing taken from these sections. The majority of laboratory tests were carried out with a simple shear apparatus, although some additional complex tri-axial tests were also performed. The primary objective was to estimate the strength parameters of soils lying in the vicinity of potential slip surfaces. The in situ investigations and numerical computations revealed the dominant role of clays in the stability analysis. The results are collected in Fig. 2 and Table 1, which are complemented by other data presented in this paper. Two formations have been distinguished on this basis. The first one, up to an altitude of 15 m a.s.l., is composed of relatively strong impermeable soils, i.e. clays, loams and silts. The second, upper formation is a complex of sandloam soils (Zabuski and Korzec 2017).
An urban project has recently been initiated to restore the damaged building ( Fig. 1). Thus, it was necessary to verify the cliff's stability in the vicinity of the building. Due to the lack of data on the considered cross-section, the arrangement of soil layers, water conditions and geomechanical parameters of particular soils were assumed on the basis of the three nearest cross-sections available, included the closest located 70 m from the building (Tejchman et al. 1995b).
Landslide processes in the Jastrzębia Góra cliff are also triggered by unfavourable hydrogeological conditions, mostly by temporary increases in the groundwater level. Measurements of the ground water table (GWT) by all piezometers proved that the GWT was located between the permeable granular soil layer and the impermeable clay layer (Tejchman et al. 1996). Thus the same GWT location, i.e. at the Table 1 Soil layers of the slope, their geotechnical parameters, the adopted random variables and their statistical characteristics (see Fig. 2 Loamy sand 2.18 x 11 c L6 ½ l 11 ¼ 21:8 r 11 ¼ 6:54 x 12 / L6 ½ l 12 ¼ 27:0 A cross-section of the Jastrzębia Góra cliff. A detailed description of the highlighted soil layers is presented in Table 1 contact between layers L4 and L3, was assumed in the section analysed, see Fig. 2. The continuous GWT is relatively high and appears on the cliff face in the form of seepage springs. Summing up, it is reasonable to extrapolate the data available for the nearest cross-sections to the section analysed in the paper. Furthermore, a number of redundant inquiries may be avoided this way.

Deterministic Stability Study of the Slope
Slope stability was investigated with the Itasca FLAC (FLAC 4.0, 2000) commercial software in a twodimensional explicit model based on the Finite Difference Method integrated with the strength reduction technique. The safety of the cliff is assessed on the basis of the safety factor F, which is calculated using the strength reduction technique fundamentals (Dawson et al. 1999). It is typically used to calibrate the factor during a stepwise reduction in the material shear strength in order to reach a state of limit equilibrium of the slope. The method is applied with the Mohr-Coulomb failure criterion. A simulation sequence is performed using trial values of the F factor to modify cohesion c, friction angle / and tensile strength r t . An incremental reduction (or increase) in strength is performed until a failure state is reached, and on this basis the selected technique finds strength values to match the final safety factor F. It is then normalized appropriately, and the safe situation means that F [ 1.
Simulations performed to assess the safety factor F of the cliff (Fig. 2) in its natural (undrained) state, using strength parameters presented in Table 1, give the safety factor F undr ¼ 1:06055. The slip, expressed by the maximum shear strain increments (SSI), encompasses the entire slope (Zabuski and Korzec 2017). Hence, the safety margin is extremely low. The potential landslide zone is shown in Fig. 3. The building at the crest of the cliff is threatened by the slip, as it is situated within the critical zone.
After drainage, which is one of the possible stabilization methods, the safety factor of the cliff increases significantly to F dr ¼ 1:29492.

Probabilistic Analysis of the Slope
To properly take into account the variation in strength parameters, the deterministic cliff analysis is extended by a probabilistic approach. It is aimed at determining both the safety margin F and slope reliability.
The basic random variables of the problem are related to the parameters of the assumed soil layers forming the slope, as shown in Fig. 2. Eight different soil layers can be distinguished in the slope: six cohesive and two non-cohesive layers. For cohesive soils, both soil cohesion c and the internal friction angle / are assumed random, while for non-cohesive soils the internal friction angle / is the only variable. Thus, a total of 14 random parameters are adopted here. The volumetric weight of each layer is assumed deterministic.
On the basis of several hundred tests classified into several groups (Tejchman et al. 1995a), the coefficients of variation of cohesion m c and coefficients of Fig. 3 Location of the potential slip zone of the cliff in its natural (undrained) state (Itasca FLAC 4.0, 2000) variation of the internal friction angle m / were defined for the soil of the highest impact on stability, i.e. clay: m c was assumed to range from 0.20 to 0.27, and m / was assumed to fall within the interval from 0.05 to 0.11. In the case of silts, these coefficients were m c ¼ 0:30 and m / ¼ 0:075 (Tejchman et al. 1996). Eventually, the values of m c ¼ 0:30 and m / ¼ 0:10 were adopted in the subsequent probabilistic analysis. Identical coefficients of variation were applied to assess the strength parameters of all geotechnical layers of the cliff, consistently with remarks found in the literature, see e.g. (Christian et al. 1994). All random parameters are assumed Gaussian. Other random variable types are applicable here as well, e.g. log-normal ones or more advanced (Brejda and Moorman 2000;Fenton and Griffiths 2003;Tobutt 1982), but the insufficient database makes the Gaussian choice reasonable. The statistical characteristics of the random variables related to particular soil layers of the slope are presented in Table 1. The high diversity in soil cohesion c i and the internal friction angle / i is intended to produce a high dispersion in cliff response results. Moreover, a uniform variation of an entire layer was assumed here, whereas a real slope may display a random field-like character, see e.g. Griffiths et al. (2009), Metya andBhattacharya (2011), Liu et al. (2014), Zhu et al. (2019). Neglecting the two-dimensional parameter variation leads to an overestimation of safety factor variance. Thus the simplifications concerning strata variability are justified from the engineering viewpoint.

Analysis of the Slope in Its Present State
(Under Undrained Conditions)

Parameter Sensitivity Analysis
First of all, test computations were performed to determine the sensitivity of the slope safety factor to random variation in individual input parameters. For each random variable, three computational cases were considered: one in which the parameter is set to its mean value (x i ¼ l i ) and two cases in which the parameter is changed by a single standard deviation (x i ¼ l i AE r i ). The remaining random parameters were all set to their respective expected values. This approach is coincident with PEM (Ahmadabadi and Poisel 2015). The results corresponding to these samples were also adopted later in estimating the slope reliability indices by CRSM. Figure 4 presents safety factor values obtained in the preliminary sensitivity analysis of the undrained slope. An overall safety factor variation d (given in [%]), induced by random variation of a given parameter, is estimated as follows where F is the mean value of the safety factor F. The impact of each random parameter on the safety factor value expressed by Eq. (11) is given in Table 2 in the column ''Present state''.
In the cases considered, the extreme safety factors were estimated at F max;undr ¼ 1:09180 ½À due to the increment in the cohesion c of layer L3 and at F min;undr ¼ 0:99414 ½À due to the decrement in the cohesion c of layer L7 (see Table 1 and Fig. 2). While the mean slope safety factor is F undr ¼ 1:06055 ½À, it can also be concluded that the dispersion of the results is not symmetric, and hence their distribution is probably not Gaussian, which is a key issue in subsequent computations.
Moreover, this variable yields a non-linear response, which also needs to be considered in the later choice of the approximation model. Thus, the character of the structural response is shown to be complex, and a higher-order approximation given in Eq. (10) is expected to minimize the fit error of response surface mapping.
The preliminary analysis presented above shows that 7 out of 14 variables are decisive to the slope response. The significance threshold was set to an impact level of d i ¼ 1:0%, and a large-gap criterion was assumed (the difference between the least significant variable and the most irrelevant variable should be greater than 25%). For example, the difference between the impacts of the variable x 14 (the least significant one) and the variable x 6 (the first irrelevant one) is 1:11 À 0:74 ð Þ =1:11 ½ Á 100% % 33%. Since the above-mentioned criteria are subjective, it is easy to choose a different number of variables (6, 8 or 11) if a different acceptance level of the sensitivity significance criterion is adopted. Moreover, it should be emphasized that variables assumed insignificant in the preliminary sensitivity analysis should not be treated as deterministic in the subsequent probabilistic analysis. This validates a full-scale analysis involving 14 random variables, rather than a smaller number of selected variables. In addition, this approach minimizes extra errors triggered by reducing the range of the random vector. On the other hand, in many cases, an excessive number of variables substantially increases the time of computations or even makes them impracticable.

PEM and MC Approach
The MC approach employs 100 samples (NS ¼ 100) to optimize the numerical quality of the analysis while maintaining an acceptable computational time. All 14 parameters of the slope were assumed random variable, and their parameters are presented in Table 1. All 100 MC samples, i.e. the generated sets of cliff soil parameters, are then introduced into distinct numerical models computed by the FLAC software. The software provides a safety factor value for each sample along with a graphical presentation of the response of the cliff.
First of all, estimators of the mean value F undr ¼ 1:00131 and the standard deviation r F undr ¼ 0:08338 of the safety factor F were computed. Convergence of these parameters is presented in Fig. 5.
With statistical distribution adopted in the form of F À F ð Þ ffiffiffiffiffiffi NS p =r F , the confidence interval was assumed according to the safety factor where t a=2; NSÀ1 is the Student's t variable argument to be exceeded with a probability of a=2. Assuming a 99% confidence interval, corresponding to a ¼ 0:01 and t 0:005; 99 ¼ 2:5674, the result was Numerical output proves that the randomness of layer parameters and layer-wise relations affect both the stability of the building and the shape of the slipping body, jointly expressed by the safety factor F. In most cases, the cliff-top building is located within the landslip zone (Fig. 6). Moreover, specific safety factor F values are attributed to various slipping surface modes.
The MC-based probability of failure P f MC;undr ¼ 0:44 was estimated by Eq. (1), and its corresponding reliability index b MC;undr ¼ 0:15 was obtained from Eq. (3). The estimated reliability indices correspond to the anticipated high failure probability cases. In the MC case, the reliability indices do not converge sufficiently well (Fig. 7). However, the direct MC variant is not a separate background procedure, but only intended to support the major probabilistic CRSM routine. It is necessary to point out that the non-linear problem type involves the computation of 100 samples, which is time-consuming beyond the limitations of a standard geotechnical design process. CRSM is introduced to estimate the safety factor F accurately enough using a smaller sample population.
Applying PEM to all 14 variables, 28 samples were taken for computations, see Eq. (5). The estimated reliability index of the slope is b PEM;undr ¼ 0:25 ]. Although MC results are more reliable, PEM will not be further referenced for this random problem. The PEM samples were already used in the sensitivity analysis of parameters essential for cliff stability. For this reason, they were also employed in CRSM computations.

Analysis by Combined Response Surface Method (CRSM)
The samples from the preliminary analysis were reused to create an approximate response surface of the phenomena by means of CRSM with a combined sample database. The values of the slope's probability of failure and reliability index were estimated using all 129 samples from direct calculations (one central sample as the starting point of the analysis, 2 9 14 PEM samples and 100 MC samples). Since a highly non-linear response surface was expected here, the second-order model was applied to assess the key values, see Eq. (10). Figure 7 shows the convergence of reliability indices checked by the direct MC method and variants of RSM approximations based on the PEM The convergence of the results obtained by CRSM and the MC method is an advantageous, decisive feature. The CRSM procedure is halted if two criteria are simultaneously met: • The CRSM solution is assumed to converge if in 10 subsequent computational runs a difference of less than 5% is detected between the computed b or P f value and its corresponding mean value based exclusively on these 10 samples. The convergence criterion in the undrained case is met after 40 software runs (69 samples, the first run uses 29 samples at once), see Fig. 7. Regarding the interval between 30 and 40 runs, the traced differences range from 2.1 to 4.9%; • The MC results are assumed to exhibit satisfactory proximity to convergent CRSM results when the difference between the mean RSM-based value of b or P f in an interval of 10 subsequent samples and 10 direct MC solutions is not greater than 15%. As shown in Fig. 7, the first 10 MC samples that fulfil the proximity criterion after CRSM convergence are samples 66-75 (in their case, the differences range from 6.7 to 14.2%).
Thus, Fig. 7 suggests that an approximate number of software runs to achieve a sufficient convergence of the CRSM solution is 75. The obtained reliability index is b CRSM;undr ¼ 0:15, so the probability of stability loss is P f CRSM;undr ¼ 0:44. These results coincide with the case of 100 direct MC calculations (b MC;undr ¼ 0:15). The relative error of the reliability index b assessed by the MC method with respect to the result obtained by CRSM was 0.64%.
Although CRSM re-uses MC samples, the failure probability P f is actually assessed by two different methodologies. This comparison is possible only in cases of a high P f , when the direct MC method does not require a large sample domain. In the case of a small P f , even several hundred MC samples may not suffice for a proper reliability assessment, so the proximity criterion may strongly limit the efficiency of CRSM. In such cases, other methods are recommended, e.g. the Targeted Random Sampling (TRS) technique, see Shields and Sundar (2015), or the subset simulation, see Jiang and Huang (2016).
The application of PEM samples along with MC samples led to expected results. The initial outline of the response surface curvature provided by PEM points proves sufficient to map the problem domain Àr i ; þr i h iaccurately, so that the introduction of only a few strictly random computational points makes the surface equation converge quickly. It is interesting to compare the results from the two RSM variants (RSM with MC samples and CRSM) (Fig. 7). It seems that CRSM accurately identifies the initial curvature of the approximated surface thanks to the application of the PEM-based sample set, next, this curvature is further refined using the MC sample set.

Analysis of the Slope Under Drained Conditions
One possible method of stabilizing the cliff is to permanently lower the groundwater table to the sea level. It is worth noting that this solution is not merely theoretical, since in dry weather the water table depth was proven to coincide with the sea level. Nevertheless, a random re-computation of a more reliable drained slope is a perfect opportunity to test the proposed CRSM computational procedure and determine its limitations.

Sensitivity Analysis
Similarly to the previous analysis, test computations were first conducted to determine the sensitivity of the slope to random variation of individual input parameters. Extreme safety factor values were F max;dr ¼ 1:35742½À, due to the increment in the cohesion c of layer L3, and F min;dr ¼ 1:14648½À, due to the decrement in the cohesion c of layer L3 (see Fig. 2). This suggests a shift of the analytical focus to marginal values of the key variable parameters of the soil strata. The impact of each random parameter on the safety factor, defined by Eq. (11), is shown in Table 2. The decisive variable contribution is 8.14%, which becomes crucial to the structural response variability. In this case, only 5 out of 14 variables affect the slope response, fewer than in the undrained slope.
The layers with the biggest contribution to the failure scenario considered here are generally different from those affecting the stability of the undrained slope, see Table 2. As the L3 clays constitute the main slip surface, the impact of the related variables on the final cliff response is almost three times as large as that of the rest. Thus, it is essential to ensure a high cohesion of the third layer of the slope during and after drainage procedures.
Referring to prior sensitivity analysis, conclusions were drawn on the further conduct of investigations. Due to the apparent diversity of the combined impact of significant and insignificant variables, displayed in Table 2, the number of variables may eventually be reduced to 5. Moreover, in view of a possible reduction in the number of random variables, the efficiency of the proposed CRSM algorithm in a relatively small random problem may be proven.
While the impact of most variables on the safety factor in the range AEr i is linear, the key layer parameters with the highest impact on the output result are clearly non-linear, so it seems reasonable to assume the second-order response surface model (Eq. 10).

CRSM Analysis
The CRSM analysis is conducted in parallel by both the direct MC approach and RSM approximations based on 11 PEM samples (one central sample and 2 9 5 starting points) and MC samples. The MC samples are subsequently added, and the CRSM convergence check is performed with respect to the cliff reliability index. The computations are terminated once a satisfactory CRSM convergence has been achieved. The course of convergence is shown in Fig. 8.
The CRSM approach clearly stabilizes the approximation process. According to Fig. 8, convergence is achieved after about 40 runs (51 samples in total: 1 central sample, 10 PEM samples and 40 MC samples). The CRSM convergence is reached after the 30 th run (differences traced for 10 calculated b values versus the mean value range from 0.2 to 0.4%), whereas the MC result proximity criterion is fulfilled for the following 10 MC samples (samples 31-40; differences traced for 10 MC-calculated b values versus the RSMbased value on this integral range from 0.3 to 4.8%). The CRSM variant results in b CRSM;dr ¼ 1:567 (P f CRSM;dr ¼ 0:059), while the 40-sample MC variant yields b MC;dr ¼ 1:645 (P f MC;dr ¼ 0:05). The direct MC results cannot serve as a reference level for CRSM because of the small population, since the probability of failure is relatively low here (about 6%). However, the low P f value obtained is evidence that drainage significantly affects slope stability. The drained case reveals a moderate probability of failure (P f CRSM;dr ¼ 0:059), so CRSM converges relatively fast, as some samples are still located in the vicinity of the limit state function. If the anticipated probability of failure decreases, convergence is expected to decrease too, resulting in reduced efficiency of the CRSM approach.

Conclusions and Remarks
This paper proposes a combined response surface approximation approach (CRSM). It is based on the sum of specific PEM samples and random MC samples. CRSM may be successfully employed in cases of a moderate or high failure probability, i.e. P f [ 0:05, representing various tasks in geotechnics.
The main feature of CRSM is multiple use of the same samples: • The CRSM combines sensitivity analysis with the probabilistic approach. Incorporation of PEM-like sampling routines made it possible to perform a basic sensitivity analysis in the course of identifying the linear or non-linear curvature of the structural response surface, which is essential for assuming a proper order of the RS approximation model. The sensitivity analysis highlights variables with a small impact on structural response in order to consider them as deterministic in computations. Moreover, PEM-like samples contribute to further reliability estimation, as they constitute the initial run in the CRSM approximation of the mechanical response of the slope. • The CRSM-postulated joint use of MC simulations with the widespread RSM software makes it possible to perform independent, parallel computations re-using the same samples. After addition of each MC sample, the solution convergence is checked, and the computations are terminated when CRSM both converges in its own right and fulfils the criterion of MC result proximity. MC sampling is merely auxiliary, so no computational convergence is required due to this method alone. • CRSM makes it possible to reduce the sample population from several hundred to several dozen, which greatly facilitates the use of probabilistic analysis in standard real-life geotechnical calculations when the value of P f is relatively large, e.g. P f [ 0:05).
In both cases presented here (i.e. an undrained cliff and a drained cliff), the parallel MC-RSM algorithm works well, each time resulting in a reasonable number of software runs. The time savings due to the possibility of sample re-usage and parallel computing is the greatest advantage of this novel approach. The estimated reliability indices b and probabilities of failure P f result in randomization of a deterministically computed slope stability factor F. In the undrained cliff, the factor was estimated at close to 1.0, but the CRSM results indicate its potential high variability, and hence a high instability of the cliff. Thus, the probabilistic computations confirm the preliminary deterministic analysis and the subsequent conclusion that the repair and operation of the building at the cliff top will be possible only after reinforcement of the cliff.
The CRSM computations prove that the proposed drainage-based cliff stabilization is effective. The probability of stability loss decreases from 44 to 6%. However, the idea of overall drainage is merely theoretical. On the other hand, the preliminary sensitivity analysis performed by CRSM identifies strata with the highest contribution to slope stability, so an alternative solution, such as a partial drainage of the key strata, may be proposed. This case needs to be reanalysed, so the variability parameters of particular soil layers should be verified to evaluate the repair works performed.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.