Determining Seismic Bearing Capacity of Footings Embedded in Cohesive Soil Slopes Using Multivariate Adaptive Regression Splines

Seismic bearing capacity of strip footings in cohesive soil slopes considering various embedded depths is investigated in this study. Novel solutions using pseudo-static method and finite element limit analysis (FELA) with upper bound (LB) and lower bound (LB) theorems are presented. The influences of footing depth, slope angle, slope height, undrained shear strength and pseudo-static acceleration on bearing capacity and failure mechanisms are examined using dimensionless parameters. With the comprehensive numerical results, the multivariate adaptive regression splines (MARS) model is then utilized to simulate the sensitivity of all dimensionless input parameters (i.e., the normalized depth of footing D/B, the normalized slope height H/B, the normalized distance from top slope to edge of the footing L/B, slope angle β, the strength ratio cu/γB, and the pseudo-static acceleration factor, kh). The degree of influence of each design parameter is produced, and an empirical equation for the dimensionless output parameter (i.e., bearing capacity factor Nc) is proposed. The study results are accessible in the design charts, tables, empirical equation for design practitioners.


Introduction
One of the most commonly used foundations in constructing railway tracks, retaining walls, transmission towers, and bridge piers is the strip foundation. These buried structures are typically supported by transferring their load to soils with enough bearing capacity and acceptable settling properties. Since the strip footing stability plays a significant role in practice, several researchers have considered the stability of strip footings on slopes by determining solutions from many different techniques including semi-empirical methods (e.g., [1][2][3]), limit equilibrium techniques (e.g., [4][5][6][7]), slip-line solutions (e.g., [8,9]), limit analysis (e.g., [10][11][12][13]), finite element methods (e.g., [14][15][16][17]), finite element limit analysis [18], and discontinuity layout optimization (DLO) approaches [19][20][21]. However, there are a few works to study the influence of seismic events on the bearing capacity of footings on slopes, which should be a concern in earthquake areas due to the destructive effects of the footing during seismic situations. One conventional and widely used approach for determining the stability of embedded structures is the pseudo-static approach, where the seismic forces are simply considered as horizontal and/or vertical seismic coefficients (k h and k v ) and functions of gravity acceleration.
The seismic responses of a footing on soil have been one of important issues in geotechnical engineering [22][23][24]. Several research have been conducted to calculate the bearing capacity of strip footings on slopes with pseudo-static seismic forces considerations. The theoretical approaches, including limit equilibrium methods (e.g., [25][26][27]), lower bound solutions (e.g., [28]) and upper bound solutions (e.g., [29][30][31][32][33]), and the stress characteristic method (e.g., [34]) which delivered an effective solution to evaluate the problems. In comparison to the analytical techniques discussed previously, a prior assumption regarding the failure mechanisms is not required, therefore, it can provide excellent predicted performance with a wide range of parameters considered. Shiau et al. [35] and Raj et al. [36] investigated the seismic bearing capacity of sloped footings by employing finite element limit analysis (FELA) which provided the upper and lower bounds solutions.
By utilizing the lower bound FELA, Kumar and Chakraborty [37] also computed the bearing capacity factor N g for a rough strip footing in cohesionless slopes under seismic scenarios. Subsequently, Chakraborty and Kumar [38] and Chakraborty and Mahesh [39] have studied the seismic bearing capacity of strip footings on a sloping ground surface and embankments utilizing the same methodologies. Using the Discontinuity layout optimization (DLO) technique, Zhou et al. [21] studied the ultimate seismic bearing capacity and failure mechanisms for strip footings placed close to the cohesive-frictional soil slopes. In addition, Cinicioglu and Erkli [40] investigated the seismic bearing capacity of strip footings lying on or adjacent to a slope by employing the finite element program PLAXIS in undrained conditions. Recently, the FELA technique was employed by Luo et al. [41], Beygi et al. [42], and Zhang et al. [23] to solve the seismic bearing capacity of strip footings on cohesive and cohesive-frictional soils, in spite that their solutions are limited to the cases of footings resting on the surface of slopes. Due to its popularity, the FELA technique has also been used to many other geotechnical problems [43][44][45][46][47][48][49][50][51][52][53].
In this paper, rigorous solutions of seismic bearing capacity of strip footing in cohesive soil slope are investigated by employing the finite element limit analysis (FELA) and a pseudo-static technique to evaluate the seismic loadings. The upper bound (UB) and lower bound (LB) outcomes attained by the FELA are compared with previously published results. The influences of the seismic acceleration coefficient, soil characteristics, and geometrical parameters on the seismic bearing capacity and the associated failure mechanisms of this problem are investigated. A comprehensive set of design tables and charts are also provided for the uses in design practices. The associated sensitivities are further assessed using multivariate adaptive regression splines (MARS) model, which is capable of accurately capturing the nonlinear relationships between a set of input variables and output variables in multi-dimensions. The MARS-based design equations used for forecasting the solutions of the seismic bearing capacity of strip footings embedded in cohesive slope are finally proposed using the artificial data set generated from FELA. The MARS-based design equation of the current study may be used to perform more precise and reliable evaluations of the seismic bearing capacity of this problem, while considering the coupled influences of the seismic acceleration coefficient, soil characteristics, and geometrical configurations.

Problem Statement and Modelling Technique
The problem definition of a strip footing on a slope is shown in Fig. 1. Under plane strain conditions, the slope has an inclination (β) and a height (H). The footing is assumed to be rigid material with width (B), depth (D), and the distance from top slope to edge of the footing (L). Defining the soil to Fig. 1 Schematic diagram of the model be a rigid-perfectly plastic Tresca material with a unit weight (γ) and undrained shear strength (c u ), the horizontal seismic acceleration (k h ) is applied to the footing and the slope under seismic forces. The vertical seismic acceleration is ignored in this study.
A typical model of the footing on slope problem is shown in Fig. 2. The boundary condition is determined using the standard fixity tool in OptumG2 [54], in which all boundary conditions are created as follows: both left-hand and right-hand boundaries are kept stationary in the x direction and the bottom of the boundary is fixed in both x and y directions. The movements at other boundaries are set to be freely moved in both x and y directions. The footing is model as a rigid elastic material and the footing-soil interface is considered as a perfectly rough condition. The size of the domain of this problem is chosen to be sufficiently large so that the plastic zone is contained within the domain and does not intersect with the right and bottom boundaries. An automatically adaptive mesh refinement technique in OptumG2 is used to improve the accuracy of upper and lower bound solutions [55]. Using this technique, the number of elements is automatically increased in the zone with high shear power dissipation that requires sensitivity analyses. The five iterations of the adaptive meshing are used in this study, where the number of elements is set to be automatically increased from 5000 to 10,000 elements [56][57][58][59][60][61][62][63][64][65][66][67].
The seismic bearing capacity of strip footing on a slope can be represented by dimensionless parameters [68] as follows: where N c is the undrained seismic bearing capacity factor; q u is the ultimate bearing capacity, H/B is the normalized slope height; L/B is the normalized distance from top slope to edge of the footing; D/B is the normalized depth of footing; c u /γB is the strength ratio. k h is the horizontal seismic acceleration. The seismic bearing capacity factor can be normalized as  [40][41][42].

Comparison of Results
The  case of (c u /γB = 5, D/B = 0, H/B = 4, and L/B = 0) for the four various slope angles. Also shown in Fig. 4 is for the case of L/B = 1. All other parameters are the same as in Fig. 3. In general, the bearing capacity factor N c decreases linearly with increasing k h . The numerical comparisons have shown that the present results are slightly larger than those lower bound results in Lou et al. [41], despite the fact that they are in good agreement. Seeing the results of k h = 0.3 in Fig. 4, there is a tendency that all curves merge into one point, i.e., one N c value. One possible reason may be due to the fact that the slopes become unstable as k h increases, resulting in one small value of N c . This comparison exercise has provided good confidence in producing all later parametric results in the paper.

Results and Discussion
The relationship between N c and L/B for the various values of k h and β is presented in Fig. 5 for the case of (c u /γB = 2. For β = 15° and 30°, the value of N c is almost constant as the value of c u /γB increases. However, the relationship between N c and c u /γB becomes nonlinear as β increases (see for example, β = 45° and 60°). In particular, for β = 60°, the increase in N c stops approximately at a value of c u /γB = 3.5, after which a slight decrease of N c is attained. It can, therefore, be concluded that the effect of c u /γB on N c is conspicuous at a higher value of β. Numerical results have also shown that an increase of k h or β leads to a decrease in N c .
The next study is for the relationship between N c and β. This is shown in Fig. 7   nonlinear. The results have also shown that an increase of k h results in a decrease of N c and an increase of L/B leads to an increase in N c. Using the same data in Fig. 7, the figure presented in Fig. 8 shows the effect of k h on N c . As expected, the increase of k h is to decrease the bearing capacity factor N c . The relationship is a linear reduction, and the gradients of the lines are almost the same for all values of β and L/B.
The final study is for the variations of N c with the normalized footing depth ratio D/B. This is presented in Fig. 9 using the case of (c u /γB = 2.5, L/B = 1 and H/B = 4). As D/B increases, so as the N c . The relationship is a nonlinear one. This trend is similar to a standard bearing capacity problem of a shallow foundation. In addition, the smaller the slope angle β, the larger the bearing capacity factor N c . The larger the k h , the smaller the bearing capacity factor N c .
All numerical results of the bearing capacity factor N c corresponded to the investigated dimensionless input parameters are presented in Tables 1, 2 and 3. These data will be used for MARS study in a later section. In regard to the associated failure mechanisms, selected studies are presented in Figs. 10, 11 and 12 for the effects of c u /γB, D/B, and k h , respectively.
The upper bound shear dissipation contour plots are normally used to represent failure mechanisms of geo-stability problems. The actual values of the colored contour are not important for a perfectly plasticity constitutive model, and, therefore, the contour bars for these plots are not normally shown in a technical document. Figure 10 shows the effects of c u /γB on the associated failure mechanisms. The selected Note that, as the value of c u /γB increases, the overall area of the slip zone reduces, and the failure type transforms from a toe-failure mode to a face-failure mode. The reduction in the area of failure zone indicates an increase in seismic bearing capacity as the shear strength ratio c u /γB of the soil slope increases.
The effect of D/B on the associated failure mechanisms is presented in Fig. 11. The selected case for this study is for (β = 15°, c u /γB = 2.5, L/B = 0, H/B = 4 and k h = 0.1). The seismic bearing capacity increases as D/B increases, and there appears that the local failure mechanism is similar to a single sided Prandtl type of failures. For studying the effect of k h on the associated failure mechanisms, the case of (β = 30°, c u /γB = 2.5, L/B = 0, H/B = 4 and D/B = 1) is chosen. This is shown in Fig. 12, where a slight decrease in the area of failure zone is depicted as the value of k h increases. This makes sense, as the larger seismic forces would enable a search of the shortest path of the slip line to the slope surface, and, therefore, it leads to a decrease in the seismic bearing capacity, as well as a reduction in the area of the failure zone.

Sensitivity Study Using MARS Model
Multivariate adaptive regression splines (MARS) model is a nonlinear and non-parametric regression approach that can be used to capture nonlinear relationships between the input variables and the output results using a series of piecewise linear segments (splines) with differing gradients. The MARS technique does not require any specific assumptions to build functional correlations between the input variables and the output results. The different splines are connected using a knot representing by the end of one spline and the where x is an input variable and t is a threshold value.
MARS model produces BFs by searching in a stepwise manner, of which the knot locations will be automatically determined using the adaptive regression algorithm. MARS model is presented by a two-step procedure. The first (forward) step gives BFs and finds their potential knots to   [69,70] where Δg is the increase in GCV when ith parameter is deleted.
To demonstrate the complex relationship between input variables and output results, MARS proposed the correlation equation by merging all linear basic functions (BFs), as shown in Eq. (5), where a 0 is the constant, N is the number of BFs, g n is the nth BF, a n is the coefficient of g n .
Compare to another machine learning approach (i.e., artificial neural networks (ANN), least-square support vector regression, extreme learning machine, Gaussian process regression [71][72][73][74][75][76][77]) MARS are considered as an effective approach [74,78]. Moreover, MARS model has been successfully applied to a number of geotechnical applications (see, e.g., [79][80][81][82][83][84][85][86][87][88][89][90][91][92][93]). Further details of MARS model can be found in Zhang [94]. This aforementioned MARS model was utilized to perform sensitivity analyses of each input variables (i.e., L/B, β, c u /γB, D/B, H/B, and k h ) and an empirical prediction of N c value introduced by considering the coupling effects of input variables. All FELA numerical results presented in Tables 1,  2  In engineering practice, a careful design requires sensitivity analysis of each input variable [95,96]. This is best presented through the relative importance index (RII) for the design output, i.e., the N c values. As mentioned above, the value of RII shows the degree of influence, i.e., a RII of 100% indicates that the corresponding input variable has the most significant influence on the output N c . Figure 13 shows the RII of each dimensionless parameter from MARS analysis. Numerical results have shown that the normalized embedded depth D/B has the greatest effect on seismic bearing capacity of footings placed on slope with a RII of 100%. This is followed by k h , H/B, β, L/B, and c u /γB with RII of 39.38%, 38.83%, 35.16%, 29.15%, and 23.24%, respectively. These results indicates that, for a general shallow foundation, the width and depth of the foundation is most important to determine bearing capacity. Although the present study considers the influence of other parameters, they cannot replace the most importance of width and depth of the shallow foundation. Moreover, this RII study has improved our understanding on the level of importance of each design parameter for the problem considered. The confidence level in practical design can, therefore, be enhanced greatly with these RII values. Table 4 presents an empirical prediction equation provided by MARS model where the 30 BFs are listed. They can be written as in the following equation:

Conclusions
This study has examined the seismic bearing capacity performance of a strip footing resting on undrained cohesive slopes using the robust finite element limit analysis with upper and lower bound theorems. The following conclusions are drawn based on the study.  3. An empirical equation with good accuracy (R 2 = 92.21%) based on MARS model was proposed to determine the seismic bearing capacity factor N c . This study has paved the road for future geo-stability research to include sensitivity analysis of multi-variable problems with the useful relative importance index (RII) and design equation. This has many practical implications in the seismic design of soil structures in geotechnical engineering. However, it still has some limitations and should be investigated further in the future. For example, the proposed equation for the seismic bearing capacity