#### Abstract

In this research work, finite element method (FEM) and 2D Plaxis were employed to generate numerical values for bilayered soils bearing a strip footing of width B and depths, *h* and *H*, in order to predict the ultimate bearing capacity (UBC) of the strip footing underlain by layered soil profile. Several research works have tried to solve bearing capacity problems using limit equilibrium (LE) techniques. But, the LE techniques have limitations in terms of soil properties and profile arrangement. The need, however, for using constitutive models or numerical methods powered by FEM and discrete element method (DEM) has been on the rise due to the versatility and robustness of these techniques to accommodate erratic soil behaviors. Multiple numerical data were generated for the case under study and artificial intelligence (AI)-based techniques; generalized reduced gradient (GRG), genetic programming (GP), artificial neural network (ANN), and evolutionary polynomial regression (EPR) were used to predict the UBC. In order to conduct the parametric analysis to investigate the effect of different soil layers, footing width and overburden pressure at foundation level on the ultimate bearing capacity sets of finite element models were prepared. This was executed using the following soil properties: soil type of top layer (from S1 to S6), soil type of bottom layer (from S1 to S6), strip footing width (*B*) (from 1.0 to 5.0 m), thickness of the top layer (*h*) (from 0.5 B to 1.0 B), overburden pressure (from 1.0 m to 3.0 m multiplied by the *γ*′*t*), and the parameters combination of each finite model in the set is randomly selected. The results of the FEM/Plaxis parametric study produced the 2D-model, deformed shape, stress distribution, and plastic point (failure point) models. The loading produced appreciable deformation on both *x*- and *y*-axes of the soil profile with the *y*-axis showing a scattered failure configuration. The AI-based prediction produced UBC equations which performed at over 90% accuracy with ANN (99.9%; 6.3%) outperforming other techniques followed by GRG, GP, and EPR.

#### 1. Introduction

The study of soil-structure interaction has been of great interest to geotechnical engineers and soil mechanics due to its obvious connection with the ultimate bearing capacity of footings underlain by soil [1]. This is the interface between soils and superstructures. Meanwhile, ultimate bearing capacity constitutes the maximum pressure a soil can withstand prior to failure. This has always been estimated using the limit equilibrium methods in constitutive modeling exercises [2]. Constitutive models in soil mechanics have been proposed by researchers as mathematical relationships related to the stress-strain behavior of soils and the application of the geometrical and property changes with response to loading with which earthwork problems are solved. Ultimate bearing capacity, slope stability, and landslide problems have been solved by relating the stress-strain Mohr-Coulomb criterion in proposing related equations or models [1]. Limit equilibrium conditions have always governed the classical techniques that have been in use example of which include the Terzaghi method of LE, Prandtl technique, Meyerhoff, and Hansen methods of determining constitutive relationships for ultimate bearing capacity [2]. These LE techniques have been explored with the assumption that the underlying soil mass is homogeneous and infinite [2]. Obviously, however, multilayered soil conditions are experienced in practice, which is the reason for many earthwork failures and low performance of geotechnical engineering infrastructure around the world. It is important to note that the estimation of the ultimate bearing capacity of foundations underlain by multilayered soil strata is needed in the safety assessment of foundation structures [3].

There have been previous attempts made to study the bearing capacity analysis of foundations on layered soils. Mosallanezhad and Moayedi [2] studied the bearing capacity of footing on layered soils by conducting a comparative analysis of different approaches for strip footing by conducting a large number of experimental, LE, and FE analyses. It was found that the first layer thickness (*h*) to the foundation width (*B*) ratio (*h*/*B*) was the critical point and more influential than other factors. Park and Lee [1] also studied the bearing capacity for multilayered clay deposits with geosynthetic reinforcement using the discrete element method (DEM). In this stability analysis of reinforced foundation soil, a multilayered soil reinforced with a horizontal strip of geosynthetics was evaluated for its ultimate bearing capacity and the distribution of tensile forces. Several bearing capacity theories and plate loading tests were used to verify the developed model. Shoaei et al. [3] reviewed available approaches for the evaluation of the ultimate bearing capacity of two-layered soils. These works extensively presented different classical and limit equilibrium (LE) techniques that have evolved over the years. This review work by [3] focused on the importance of different failure strains and load spread angles in the analysis of foundations on multilayered soils. The computer simulation of the bearing capacity of multilayer soils in Baotou was also conducted by [4]. Various other approaches have been employed to evaluate the ultimate bearing capacity of strip foundations on multilayered soils [5–12]. However, in addition to LE and numerical approaches to determining the ultimate bearing capacity of strip foundations underlain by multilayered soils, AI-based techniques such as ANN and GP have also been employed in this case study with great results. In previous research works, ANN and some metaheuristic approaches were used to predict the ultimate bearing capacity of strip foundations over multilayered soils [13–16]. In another development, dragonfly algorithm, Harris hawks optimization, and sparse polynomial chaos expansions were employed to predict the ultimate bearing capacity of strip footing underlain by multilayered soil [17, 18] as well as other relevant research papers on the geotechnics of strip and reinforced footing on soft clays [19, 20]. Also, genetic programming (GP) was used to predict the ultimate bearing capacity of shallow foundations with good performance accuracy [21]. In another study, Dutta et al. [22] studied the bearing capacity of footing on a two-layered arrangement of sand and clay under inclined loading. It was observed that the bearing capacity of this studied arrangement compared well with lower angles of inclination and overestimated for higher ones. All the above-reviewed LE and numerical approaches to determining the ultimate bearing capacity of strip footing underlain by multilayered soils, and none was able to couple the use of LE and numerical (FE and AI) techniques to formulate numerical data and at the same time predict AI-based models for the ultimate bearing capacity (UBC). In this work, FEM/Plaxis 2D was used to generate multiple numerical data from bilayered soils selected by permutation from six (6) different soil properties and thicknesses. Subsequently, the numerical data were deployed in an AI-based predictive model exercise to predict the ultimate bearing capacity.

#### 2. Methodology

##### 2.1. FEM/PLAXIS Analysis

Geotechnical engineers often deal with layered foundation soil, which is nonhomogeneous in nature but can be simplified in representation as distinct homogeneous layers for engineering purposes. The failure mechanism of layered soil depends on the thickness and soil properties of each layer. In some cases where the top layer is relatively thick and consists of weak soil, the failure mechanism may be limited to the top layer only, and the strength of the remaining lower layers has no influence. In many other cases, however, the failure mechanism may involve two or more layers. The parametric study was carried out using PLAXIS 20 to evaluate the ultimate bearing capacity of a strip footing resting on two different soil layers.

The following Figure 1 shows the 2-dimensional plane strain model; a footing with a width B is rested on a half-space soil layer with thickness (*h*) followed by another soil layer with thickness (*H*).

Due to symmetry, only half of the problem is modeled. The model has a dimension of 10 B in both *X* and *Y* directions. The size of the finite element model is large enough to keep the boundary conditions at the bottom and the right side from restricting the soil movement due to the footing load. Furthermore, a 15-node element was used in the finite element mesh. Soil layers were modeled using Mohr-Coulomb’s constitutive law with parameters shown in Table 1.

A parametric analysis was conducted to investigate the effect of different soil layers, footing width, and overburden pressure at foundation level on the ultimate bearing capacity; therefore, sets of finite element models were arranged using the following parameter ranges: soil type of top layer (from *S*1 to *S*6), soil type of bottom layer (from *S*1 to *S*6), strip footing width (*B*) (from 1.0 to 5.0 m), thickness of the top layer (*h*) (from 0.5 B to 1.0 B), overburden pressure (from 1.0 m to 3.0 m multiplied by the *γ*′*t*), and the parameters combination of each finite model in the set are randomly selected. Samples of the generated FEM models and its output are shown in Figure 2. This shows the 2D model, deformed shape, stress distribution, and points of failure for the foundation under loading considerations.

##### 2.2. FEM/PLAXIS-Generated Database, Statistical, and Correlation Analyses

In the course of the constitutive phase of the research work, 156 PLAXIS models were developed to generate the required dataset for the regression process, and each record in the dataset contains the following fields:(i)Cohesion of top layer (Ct) kN/m^{2}(ii)Tangent of friction angle of top layer tan (*φt*)(iii)Effective density of top layer (*γ*′*t*) kN/m^{3}(iv)Thickness of top layer (*h*) m(v)Cohesion of bottom layer (Cb) kN/m^{2}(vi)Tangent of friction angle of bottom layer tan (*φb*)(vii)Effective density of bottom layer (*γ*′*b*) kN/m^{3}(viii)Width of strip footing (*B*) m(ix)Effective overburden pressure at foundation level kN/m^{2}(x)The ultimate bearing stress of strip footing kN/m^{2}

The generated records were divided into training and validation sets (110 (70%) and 46 (30%) records, respectively). Tables 2 and 3 show the statistical characteristics of the generated dataset and the correlation matrix while Figure 3 presents the statistical distribution (histograms) for input and output parameters.

**(a)**

**(b)**

##### 2.3. Research Program

Besides the well-known generalized reduced gradient (GRG) technique, another three (AI) techniques were used to predict the ultimate bearing capacity of strip footing using the generated database. The used techniques are evolutionary polynomial regression (EPR), artificial neural network (ANN), and genetic programming (GP). These models were developed to predict the ultimate bearing stress of strip footing using the soil parameters of both layers (Ct), (Cb), tan (*φt*), tan (*φb*), (*γ*′*t*), (*γ*′*b*), and geometrical parameters (*h*), (*B*), and . The accuracy of each model was evaluated using the sum of squared errors (SSE).

#### 3. Results and Discussion, and Prediction of Ultimate Bearing Capacity

##### 3.1. Model (1)—Evolutionary Polynomial Regression (EPR)

A pentagonal level EPR model was developed, for 9 inputs; there are 1287 possible combinations (792 + 330 + 120 + 36 + 8 + 1 = 1287) as follows:

The most effective 39 terms were selected by the GA technique. The developed model to predict the bearing capacity is presented in (2), while its fitness is shown in Figure 4(a). The determination factor (*R*^{2}) and average error values are 0.939 and 27.2%, respectively.

**(a)**

**(b)**

**(c)**

**(d)**

##### 3.2. Model (2)—Artificial Neural Network (ANN)

A single hidden layer (ANN) with a nonlinear activation function (hyper-tan) was trained using the backpropagation technique to predict bearing capacity. The layout of the developed ANN and their connotation weights is shown in Figure 5 and Table 4. The average error was 6.3%, and the *R*^{2} was 0.997. The relations between calculated and predicted values are illustrated in Figure 4(b).

##### 3.3. Genetic Programming (GP) Prediction of Ultimate Bearing Capacity

Five levels of complexity (GP) model (63 genes per chromosome) was developed to predict the bearing capacities values of generated database records. This model was developed using a population size of 100000 chromosomes, a survivor size of 25000 chromosomes, and 250 generations. Figure (3) presents the generated formula for the ultimate bearing capacity, while Figure 4(c) shows its fitness. The average error and (*R*^{2}) values for this model were 11.7% and 0.989, respectively.

##### 3.4. Generalized Reduced Gradient (GRG) Prediction of Ultimate Bearing Capacity

A GRG model was developed using the MS EXCEL add-in solver module. The selected structure for the bearing capacity formulas and optimized coefficients is shown in (5), and its fitness is illustrated in Figure 4(d). The corresponding (*R*^{2}) value is 0.992 while the average error is 10%. The accuracies of developed models are compared in Table 5.where

#### 4. Conclusions

This research was concerned in predicting the ultimate bearing capacity of strip footing rested on bilayered soil profile using the following four techniques: EPR, ANN, GP, and GRG. These techniques were applied to a generated dataset of 165 records, and each record contains data for soil parameters of both layers (Ct), (Cb), tan (*φt*), tan (*φb*), (*γ*′*t*), (*γ*′*b*), and geometrical parameters (*h*), (*B*), and besides the corresponding ultimate bearing capacity . The used database was generated using the well-known PLAXIS software. The output of each technique could be evaluated as follows:(i)The first technique (EPR) generated an optimized pentagonal polynomial with 39 terms selected out of 1287 possible terms. Its accuracy was lower among all developed models with an average error percent of 27.2%. Besides that, it includes only five parameters (Ct), tan (*φt*), tan (*φb*), (*h*), and (*γ*′*b*) and neglects the effect of the rest (Cb), (*γ*′*t*), (*B*), and .(ii)The second technique (ANN) showed the highest accuracy with an average error percent of 6.3%, and it utilized all the input parameters. The importance of each parameter is illustrated by the size of its input blocks in Figure 2, which indicates that all parameters have almost the same impact on the ultimate bearing capacity except (*h*), (Ct), and tan (*φt*) which have a slightly higher impact.(iii)The GP technique generated a complicated 3rd-degree exponential formula with a very good accuracy and an average error percent of 11.7%. All input parameters were used except the top layer thickness (*h*); besides that, using tan (*φt*) and tan (*φb*) as exponent assured their significant impact on the ultimate bearing capacity.(iv)Although GRG is not an advanced (AI) technique, it showed a very good accuracy with average error of 10%. The main advantage of this technique is the chosen structure of the generated formula which matched the well-known general bearing capacity equation. In addition, all input parameters were used. On a practical need, these models can be applied in predicting the ultimate bearing capacity of a similar soil layered arrangement under loading for design and foundation performance and behavior purposes.

#### 5. Recommendations

Based on the previous discussion, the following points could be recommended:(i)It is not recommended to use the EPR formula due to its low level of accuracy and the fact that ignored the impact of almost half of the input parameters.(ii)Although the ANN model showed a perfect accuracy and accounts for all input parameters, its output is a weight matrix which is very difficult to be handled manually.(iii)Both GP and GRG models showed the same level of accuracy, and both generated a closed-form formula that could be used manually. However, the absence of (*h*) parameter in the GP model is a significant shortage. On the other hand, the well-known structures of the GRG formula are an important advantage.(iv)For computerized calculations, it is recommended to use the ANN model, while the minimum of the GP and GRG models is recommended for manual calculations.(v)All generated models are valid within the considered range of input parameter values, beyond this range; the prediction accuracy should be verified.

#### Data Availability

All data, models, and code generated or used during the study appear in the submitted article.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.