Prediction of the Bearing Capacity of Composite Grounds Made of Geogrid-Reinforced Sand over Encased Stone Columns Floating in Soft Soil Using a White-Box Machine Learning Model

Abstract

Stone columns have been extensively advocated as a traditional approach to increase the undrained bearing capacity and reduce the settlement of footings sitting on cohesive ground. However, due to the complex interaction between the soil and the stone columns, there currently needs to be a commonly acknowledged approach that can be used to precisely predict the undrained bearing capacity of the system. For this reason, the bearing capacity of a sandy bed reinforced with geogrid and sitting above a collection of geogrid-encased stone columns floating in soft clay was studied in this research. Using a white-box machine learning (ML) technique called Multivariate Polynomial Regression (MPR), this work aims to develop a model for predicting the bearing capacity of the referred foundation system. For this purpose, two hundred and forty-five experimental results were collected from the literature. In addition, the model was compared to two other ML models, namely, a black-box model known as Random Forest (RF) and a white-box ML model called Linear Regression (LR). In terms of R² (coefficient of determination) and RMSE (Root Mean Absolute Error) values, the newly proposed model outperforms the two other referred models and demonstrates robust estimation capabilities. In addition, a parametric analysis was carried out to determine the contribution of each input variable and its relative significance on the output.

1. Introduction

Stone columns enhance the consolidation rate, decrease and regulate the soft saturated clay’s compressibility, increase the properties of the soft soil, especially the bearing capacity, and minimize differential and total settlements caused by the loading of overlying structures [1,2,3,4,5,6,7,8,9,10,11]. Stone columns get their strength and stiffness from the restriction of the soil around them, but their use in extremely soft soils has been restricted because of the low stiffness and low confinement provided by such soils [12]. The stone elements experience an increase in lateral stress proportional to the vertically applied load. As a consequence, the materials used for stone columns tend to sink or protrude excessively into the circumferential soft soil. Encasing stone columns with geotextile, in opposition to using conventional (non-encased) stone columns, affects the distribution of the lateral stress by enhancing the strain–stress behavior. This construction technique also allows for a decrease in the maximum settlement under high loading conditions because the strain–stress behavior is improved. As a result, geotextile-encased stone columns, also known as GESC for short, constitute a method for ground enhancement that can be utilized for extremely soft soils.

Recently, several experimental and numerical studies focused on the performance of GESC were performed.

Fattah and collaborators [9,13,14,15] assessed the performance of soft soil after it was reinforced with GESC and with original stone columns (OSC). They used different inter-column distances, length/diameter ratios, soil shear strengths, and area substitution ratios to explore the effect of these parameters on the improvement of the bearing capacity and on the reduction of the settlement. Kwa et al. [16] conducted an experimental research on the behavior of GESC to reinforce soft clay and concluded that the final bearing capacity of GESC is 1.6 times higher than that of OSC. Malarvizhi and Ilamparuthi performed small-scale tests to assess the enhancement of clay beds stabilized with GESC [17]. The results showed that the beds reinforced with GESC presented an ultimate bearing capacity that was three times higher when compared to the unreinforced bed.

In addition, Gu et al. [7] observed that the GESC substantially increased the final loading capacity of the soft soil and that the effective length of the encasement was 3 to 4 times higher than the diameter of the stone columns. Orekanti and Dommaraju [18] performed compression tests with restricted strains on GESC. The findings indicated that composite clays with GESC in end-bearing type and floating type allowed for a significant increase in the load-carrying capacity in comparison with clay, by 2.44 and 2.01 times, respectively. Hajiazizi and Nasiri [19] conducted an experimental study to investigate the behavior of ordinary stone columns (OSC) and GESC to strengthen sand slopes. The observed results showed that placing GESC in the middle of a slope tends to increase the bearing capacity of the sloped crown by 2.17 times more than placing OSC in the same location.

The findings of Yoo and Lee [20] and Ayadat and Hanna [21] were based on experiments that looked into the load-bearing capacity and settlement of GESC in soft and collapsible soils, respectively. These experiments showed that the encasement provided extra confinement, which resulted not only in an increased maximum carrying capacity of the stone column but also in a decrease in soft ground settlement when compared to traditional stone columns.

Almeida et al. [22] documented their findings on the behavior of a GESC-reinforced test embankment on soft soil. As the embankment height rose and the soil consolidated, they observed an increase in the differential settlement between the soil and columns. This was the conclusion they reached after observing the behavior of the test embankment on soft soil enhanced by GESC. The stone columns were typically covered with a 0.3 m or thicker layer of granular sand or gravel to serve as a drainage path and to disperse the loads from the superstructures [23]. Shahu et al. [24] proposed a clear and concise theoretical framework for investigating granular mat-topped and granular pile-reinforced soft soil. Deb [25] proposed a mechanical method to estimate the behavior of soft soils stabilized with stone columns and covered with a granular layer. Granular layer reinforcement on the top of stone column-reinforced soil has been found to minimize stress concentration around the top of the columns. Additionally, the granular layer contributes to decreasing the maximum and differential settlement and also improving the load-carrying capacity of the stone column. This allows for enhanced performance on the soft ground. However, because of the complicated geometry and the many uncertainties associated with the geotechnical parameters, it is difficult to determine the bearing capacity of composite foundation systems built with geogrid-reinforced sand beds (GRSB) and with vertically encased stone columns (VESC).

Machine learning (ML) approaches have recently attracted much attention for developing paradigms that can precisely predict complex system parameters [12,26,27,28,29,30,31,32,33,34]. Al-Obaidy and Al-Shueli [30] predicted the bearing capacity of stone columns using MATLAB’s Neural Network Toolbox. The proposed model showed reliable results, allowing the model to be used as a forecasting tool instead of relying on the results of expensive and time-consuming field or experimental tests. Similarly, Das and Dey [28] applied the Adaptive Neuro-Fuzzy Inference System (ANFIS) to predict the bearing capacity of stone columns. The authors developed three models: ANFIS-E, which used experimental results as input; ANFIS-A, which used analytical results as input; and ANFIS-EA, which used both experimental and analytical results as input. Input data were acquired from the existing literature to train the models. For each stone column arrangement and embankment condition, Chik et al. [29] employed Artificial Neural Network systems (ANNs) to predict the settlement under embankment load. Their study utilized soft soil properties together with different geometric factors as input. The ANN model was developed using Ten-Fold Cross-Validation (TFCV) and Non-Cross Validation (NVC). The results showed superior performance of the TFCV model compared with NVC. Aljanabi et al. [31] adopted another ML technique known as Support Vector Machine (SVM) to accurately and effectively predict the bearing capacity behavior of stone columns. In [12], the SVM models were also used to estimate the bearing capacity of geogrid-reinforced stone columns. Three different kernel functions—the polynomial (POLY), radial basis (RBF), and exponential radial basis (ERBF)—served as the foundation for the development and performance evaluation of three different versions of the SVM model. With the optimum generalizability and minimum statistical error of the training and testing data, ERBF-based SVR demonstrated robust performances. Additionally, SVR-ERBF was evaluated against the ANFIS model in predicting the bearing capacity. The authors concluded that the SVR-ERBF model has better performance. In addition, in the same study, the effect of the input parameters on the output predictions from the models was assessed by performing a sensitivity analysis. To forecast the bearing capacity of geogrid-reinforced stone, Ghanizadeh et al. [35] developed a prediction model based on the Escaping Bird Search (EBS) optimization method and Multivariate Adaptive Regression Splines (MARS). For the training set, test set, and full dataset, respectively, the coefficients of determination (R²) scores of their prediction model were 0.997, 0.993, and 0.995, respectively.

This study mainly focuses on providing accurate estimations for the bearing capacity of foundation systems composed of GRSB and VESC. The input variables considered in this research were the ratio of the sandy bed thickness to the foundation diameter (t/D), the ratio of the diameter of the geogrid-reinforced layer to the diameter of the footing (d/D), the undrained shear strength of unreinforced clay (q_u), the ratio of the length of the encasement to the diameter of the stone column (h/d_sc), and the ratio of the length to the diameter of the stone column (L/d_sc). For this purpose, a databank related to 245 experimental tests is gathered from the open-source literature [12]. Then, a comprehensive and reliable model is developed by using Multivariate Polynomial Regression (MPR), which is a white-box supervised intelligence system. To the best of the authors’ knowledge, this is the first study in which MPR is used to produce mathematical formulas for determining the bearing capacity of this type of composite foundation system. In contrast to the majority of the earlier uses of ML techniques in such fields, which were black-box models and required specialized software to compute the predictions, the explicitly developed model in this study is user-friendly and can be used with a basic calculator. Several statistical and graphical criteria are used to assess the accuracy and reliability of the developed model. Finally, the effect of each input variable on the bearing capacity is quantitatively evaluated using a parametric analysis.

2. Materials and Methods

2.1. Research Methodology

This study aims to develop an efficient technique for forecasting the bearing capacity of geogrid-reinforced sand over geogrid-encased stone columns floating in soft soil. The method is outlined in Figure 1 and involves comparing the performance of the MPR method with that of Linear Regression (LR) and Random Forest (RF) models. The first step was to gather experimental data from previous studies, followed by splitting the data into training and testing sets. The training set was then subjected to a fivefold cross-validation to assess the performance of the three models. All three models were then applied to k-fold cross-validation, and their performance was evaluated using standard performance metrics. In addition, a regression equation is derived from the training dataset and applied to new data to validate the model. To ensure that the model was reliable and would produce good results with different data combinations, a parametric study was conducted as a final step. The used research methodology is depicted in Figure 1.

2.2. Multivariate Polynomial Regression

A polynomial regression can be used to fit a nonlinear relationship between independent variables and a dependent variable. The equation for a kth-order polynomial regression (with k greater than 1) is represented by the coefficients ${\omega }_1$, ${\omega }_2$, … ${\omega }_k$, the intercept ${\omega }_0$, the input variable $x$, and the outcome variable $\widehat{y}$ [36,37].

$$\widehat{y}={\omega }_0+{\omega }_{1}x+{\omega }_2{x}^2+{\omega }_3{x}^3+\cdots +{\omega }_k{x}^k$$

(1)

When multiple variables are used in a polynomial regression, it is referred to as Multivariate Polynomial Regression (MPR). For a system with lease, please find attached my revision of the article, which includes proposals to improve the text, some suggestions, and also some comments for you to analyze input variables, the kth-order multivariate polynomial regression can be written as [37,38]:

$$\begin{array}{l}\widehat{y}={\omega }_0+{\displaystyle \sum _{l_1=1}^m}{\omega }_{l_1}{x}_{l_1}+{\displaystyle \sum _{l_1=1}^m}{\displaystyle \sum _{l_2=l_1}^m}{\omega }_{l_1{l}_2}{x}_{l_1}{x}_{l_2}+\cdots +\\ {\displaystyle \sum _{l_1=1}^m}{\displaystyle \sum _{l_2=l_1}^m}\cdots {\displaystyle \sum _{l_k=l_{k-1}}^m}{\omega }_{l_1{l}_2\cdots l_k}{x}_{l_1}{x}_{l_2}\cdots x_{l_k}\end{array}$$

(2)

The input variables in the polynomial regression are represented by $x_{l_1}$, $x_{l_2}$,..., $x_{l_k}$ and the coefficients by ${\omega }_{l_1}$, ${\omega }_{l_1{l}_2}$,..., ${\omega }_{l_1{l}_2\cdots l_k}$. The MPR uses a nonlinear model to fit the data, but it is linear in terms of its coefficients. Therefore, the solution of the MPR problem is similar to that of the Multiple Linear Regression (MLR) problem, which uses the least-squares method. The least-squares technique finds the optimal polynomial coefficients by reducing the residual sum of squares, which can be mathematically represented as the sum of squared errors between the predicted outcome and the actual outcome.

$$\min \mathrm{Loss} {(\mathrm{y},\overrightarrow{\mathrm{y}})}_{\mathrm{MLR}}=\sum _{i=1}^N{\left(y_i-{\omega }_0-\sum _{J=1}^p{x}_{ij}{\omega }_j\right)}^2$$

(3)

In the multivariate regression model, $x_{ij}$ represents the ith value of the jth input feature, ${\omega }_j$ is the coefficient for the jth input feature, and $p$ is the total number of input features used in the model.

2.3. Sofware Implementation

In the development of the MPR model, a second-order polynomial with interaction terms was employed. This choice allowed for capturing the relationships between the input variables more accurately while maintaining a balance between model complexity and generalization performance. The inclusion of interaction terms enables the model to consider potential synergistic effects between different input variables, leading to improved prediction accuracy. In this study, the prediction models were constructed using the WEKA workbench toolkit (version 3.9.4) [39]. Several machine learning and data mining methods are included in WEKA software version 3.9.6. It was created for academics wishing to quickly test current approaches on brand-new datasets. It provides thorough assistance for experimental data mining, including input data preparation, statistical learning scheme evaluation, and learning outcome demonstration. However, if the input attributes in the training data are strongly correlated, the performance of the MPR model can suffer. By defaulting the feature eliminate Colinear Attributes to true, WEKA can automatically identify and delete highly correlated input variables [40].

Performance can also be severely impacted by attributes that have no connection to the output variable. By configuring the feature attribute Selection Method. WEKA can perform feature selection automatically and only choose the relevant attributes [40]. It is possible to turn off this feature, which is activated by default. Finally, ridge regularization is used in the WEKA softwareto prevent overfitting in the model. Overfitting occurs when a model is too complex and fits the training data too well, leading to poor generalization performance on unseen data. Ridge regularization adds a penalty term to the loss function that discourages the model from fitting the noise in the data. This results in a model with reduced complexity, improved generalization performance, and a reduced risk of overfitting.

2.4. K-Fold Cross Validation

Cross-Validation (CV) is a well-known resampling technique used to evaluate ML models in the case of a limited data sample. This method splits the data sample into several groups determined by the k parameter. The model validation approach of CV assesses the generalizability of statistical analysis results to a different dataset. It has a diverse range of methods, for instance, leave-one-out CV, K × 2, K-fold, and repeated random sub-sampling validation. This research used K-fold cross-validation for comparisons and assessments of the model’s performance [41,42,43]. This technique randomly divides the original samples from the dataset into K-sub samples. Following that, a validation set is held using one sub-sample for model testing, and the model is trained using the remaining k-1 sub-samples. In this method, each of the k sub-samples serves as the validation data once during the Kth iteration. After that, it is possible to construct a single estimation by averaging (or combining) the k results from the folds. The fundamental process of a K-fold CV is as follows:

• Random shuffling of the dataset.
• Splitting the dataset into k number of groups.
• For every distinct k-group:
  • Use one group as a validation set;
  • Use the remaining groups as a train set.;
  • Use the train set for model training and the test set for model validation;
  • Preserve the results from each fold and discard the predictive model.
• Using a sample of model assessment scores to summarize the skills of the model.

It is common practice to summarize the results of a K-fold cross-validation run using the mean of the model skill scores. It is also a good idea to use a variance metric, such as the standard error or standard deviation, to evaluate the model’s performance.

The value of k is typically set to 5 or 10, but this is not a strict guideline. As k increases, the size difference between the resampled subsets and the training set decreases, leading to a reduction in bias [44]. In this study, k was set to 5 for the MPR model. In addition to the MPR model, two other ML algorithms were developed for evaluation: the Random Forest (RF), which is a black-box model, and the Linear Regression (LR), which is a white-box model. These models were constructed using the WEKA software version 3.9.6 and subjected to the same 5-fold cross-validation as the MPR model.

2.5. Performance Metrics

The predictive performance of the developed model was evaluated using several statistical evaluation parameters. In this research, the evaluation parameters were computed based on the values of the actual output variable (bearing capacity) from the database and on the bearing capacity predicted by the developed MPR model. Although several evaluation parameters can be used to assess the model’s performance, this research considers two key assessment parameters: the Root-Mean-Squared Error (RMSE) and the coefficient of determination ($R^2$). These parameters have been widely used in previous studies and are considered sufficient for assessing the prediction performance of the models. RMSE is a measure of the difference between the predicted values and the actual values. It is used to evaluate the performance of a model, specifically for regression problems. It is calculated by taking the square root of the mean of the squared differences between the predicted values and the actual values. The RMSE of the model is computed using Equation (4). The smaller the RMSE, the better the model is at predicting the target variable. RMSE is sensitive to outliers, so it is important to look at other metrics as well, especially when the data set has outliers.

$$\mathrm{RMSE}=\sqrt{\frac{1}{n}{\sum }_{i=1}^n{\left(y_i^{obs}-y_l^{pre}\right)}^2}$$

(4)

It is important to note that RMSE has the same unit as the target variable, which makes it easy to interpret and compare to other models.

The coefficient of determination (R²) is a statistical indicator that reflects the extent to which the independent variable(s) in a regression model explain the variance in the dependent variable. It ranges from 0 to 1, with a value of 1 indicating that the model completely explains the variance and a value of 0 meaning that the model explains none of the variance. The R² can be determined using the following equation:

$$R^2=1-\frac{{\sum }_{i=1}^n{\left(y_i^{obs}-y_l^{pre}\right)}^2}{{\sum }_{i=1}^n{\left(y_i^{obs}-y_i^{pre}\right)}^2}$$

(5)

where $y_i^{obs}$ indicates the actual outcome of the bearing capacity, $y_i^{pre}$ represents the predicted outcome, $n$ indicates the sample count employed for modeling the bearing capacity, and $y_i^{-obs}$ is the mean value of the actual outcome.

3. Database Used

The effectiveness of the model is dependent on the variables used and the amount of data. The variables used to predict the bearing capacity were gathered from the literature [12], resulting in 245 experimental tests. All the test procedures and methods were established and described in a previous paper by the same authors [45]. To construct the models, related literature was utilized to justify the following predictors: the ratio of the sandy bed thickness to the foundation diameter (t/D), the ratio of the geogrid layer diameter to the foundation diameter (d/D), the ratio of the stone column length to its diameter (L/d_sc), the ratio of the length of the encasement to the diameter of the stone column (h/d_sc), and the shear strength of unreinforced soft clay (q_u). The models aim to predict the bearing capacity of the improved soft clay soil (q_rs). These input variables have also been used by Debnath and Dey for modeling this bearing capacity [12]. The data was then split into a training set and a testing set with an 8:2 ratio, including 197 observations for training and 48 for testing. A stratified sampling approach was used to ensure that the data had a comparable outcome distribution. Figure 2 shows the overall distribution of the variables, and Table 1 presents the statistics of the input variables for the training and testing sets.

Table 1. Model-variable descriptive statistics.

Data Category	Statistics	qu (kPa)	t/D	d/D	L/dsc	h/dsc	qrs (kPa)
Training data	Standard deviation	15.05	0.08	1.42	0.96	2.02	118.33
	Mean	38.89	0.16	1.64	5.88	5.01	197.26
	Median	45.02	0.15	2.00	6.00	6.00	178.81
	Maximum	53.82	0.40	4.00	8.00	8.00	452.21
	Minimum	3.86	0.00	0.00	2.00	0.00	15.44
	Kurtosis	−0.52	1.97	−1.13	8.58	0.80	−0.87
Testing data	Standard deviation	17.35	0.09	1.38	1.26	1.93	120.11
	Mean	36.66	0.19	1.85	5.67	4.89	195.44
	Median	44.24	0.15	2.00	6.00	6.00	178.40
	Maximum	53.82	0.40	4.00	8.00	8.00	467.01
	Minimum	4.99	0.00	0.00	2.00	0.00	14.11
	Kurtosis	−1.03	0.63	−0.89	4.44	−0.17	−0.65

Table 1. Model-variable descriptive statistics.

Data Category	Statistics	qu (kPa)	t/D	d/D	L/dsc	h/dsc	qrs (kPa)
Training data	Standard deviation	15.05	0.08	1.42	0.96	2.02	118.33
	Mean	38.89	0.16	1.64	5.88	5.01	197.26
	Median	45.02	0.15	2.00	6.00	6.00	178.81
	Maximum	53.82	0.40	4.00	8.00	8.00	452.21
	Minimum	3.86	0.00	0.00	2.00	0.00	15.44
	Kurtosis	−0.52	1.97	−1.13	8.58	0.80	−0.87
Testing data	Standard deviation	17.35	0.09	1.38	1.26	1.93	120.11
	Mean	36.66	0.19	1.85	5.67	4.89	195.44
	Median	44.24	0.15	2.00	6.00	6.00	178.40
	Maximum	53.82	0.40	4.00	8.00	8.00	467.01
	Minimum	4.99	0.00	0.00	2.00	0.00	14.11
	Kurtosis	−1.03	0.63	−0.89	4.44	−0.17	−0.65

To demonstrate the variability of the data, examples include the ratio of the sandy bed thickness to the foundation diameter, which can range from 0 to 0.4; the ratio of the geogrid layer diameter to the foundation diameter, ranging from 0 to 4; the ratio of the stone column length to its diameter, ranging from 2 to 8; the ratio of the length of the encasement to the diameter of the stone column, ranging from 0 to 8; and the shear strength of unreinforced clay, varying between 3.86 and 53.82 kPa. Moreover, the measured bearing capacity of the composite soil also demonstrated a wide range of values, falling between 14.11 and 467.01 kPa. The wide range of potential values for both input and output variables in the dataset highlights the diversity of the samples included and the ability of the models constructed with this data to be applicable to a wide range of scenarios.

4. Model Results

4.1. K-Fold Cross-Validation

The training performance of each ML model is depicted in the boxplots in Figure 3 using the R² and RMSE metrics. The boxplots show the variability and fluctuations of the performance measures by using 5-fold cross-validation on 197 actual values from the training data sets. The five-number summary of a dataset, consisting of the minimum, first quartile, median, third quartile, and maximum values, is graphically displayed in a box plot, also known as a box-and-whisker plot. These values are represented in the box plot by a box that spans from the first quartile to the third quartile and a line inside the box that represents the median. The box plot also includes “whiskers” that extend from the box to the minimum and maximum values and diamonds that represent any values that fall outside of the whiskers. Specifically, outliers are any data points that are higher than 1.5 times the Interquartile Range (IQR) and are below the first quartile (Q1) or above the third quartile (Q3). Box plots are commonly used in statistics and data analysis to quickly and easily visualize the distribution of a dataset. The minimum, maximum, median, and average for RMSE and R² for three machine learning models using fivefold cross-validation are presented in

Table 2. Cross-validation measurement results.

Statistics	MPR		RF		LR
	RMSE	R2	RMSE	R2	RMSE	R2
Mean	13.9161	0.9850	19.2054	0.9820	42.1545	0.8697
Median	13.6817	0.9856	17.0350	0.9814	42.5097	0.8674
Maximum	16.0838	0.9902	31.8202	0.9900	44.7245	0.8902
Minimum	11.3142	0.9784	11.7660	0.9764	39.6920	0.8472

.

From Figure 3 and Table 2, the MPR model exhibits better performance than LR and RF in terms of the median and mean values of the two metrics. This superiority is due to the non-linear associations between the bearing capacity and the input variables that nonlinear models can effectively capture.

Further investigation is needed to ensure the stability and dependability of the model, in addition to evaluating its prediction accuracy. The random error should have a normal distribution with a mean of 0 and a symmetric shape, with most errors being close to the center. Figure 4 displays the residuals and Gaussian functions for the first three models, which all have a normal distribution pattern. The MPR model seems to have a higher level of prediction accuracy when compared to the RF and LR algorithms, as 80% of the residuals for the MPR fall within the range of [−17.267, 17.267]. Moreover, the mean value of the fitted normal distribution of residuals for the MPR model is 0.703. It is much closer to 0; that is, the MPR model has better error randomness. So, it can be stated that the MPR algorithm is more advantageous.

4.2. Evaluation of the RF Model

The constructed MPR model underwent a comparison with the two other models through 5-fold cross-validation and was based on the evaluation. For this, a method for estimating the bearing capacity $q_{rs}$ was established using the training dataset. The formula is illustrated as follows:

$$\begin{array}{cc}{q}_{rs}=17.6798+& 5.0749\times q_u-6.6152\times \left(\frac{l}{d_{sc}}\right)-0.130\times {\left(q_u\right)}^2\hfill \\ & +0.0018\times {\left(q_u\right)}^3+0.0247\times {\left(q_u\right)}^2\times \left(\frac{d}{D}\right)\hfill \\ & -0.4547\times q_u\times {\left(\frac{d}{D}\right)}^2+0.2410\times q_u\times \left(\frac{h}{d_{sc}}\right)\times \left(\frac{d}{D}\right)\hfill \\ & -247.4930\times \left(\frac{d}{D}\right)\times {\left(\frac{t}{D}\right)}^2+29.6706\times \left(\frac{t}{D}\right)\times {\left(\frac{d}{D}\right)}^2\hfill \\ & +3.1539\times \left(\frac{t}{D}\right)\times \left(\frac{h}{d_{sc}}\right)\times \left(\frac{l}{d_{sc}}\right)\hfill \end{array}$$

(6)

The MPR model was used to determine the load-bearing capacity of new data, and the results were evaluated using pair plots. Figure 5a shows a comparison of the predicted values for the load-bearing capacity to the actual values for the test data. The predictions from the MPR model align closely with the 45-degree line, which indicates that the model can provide accurate predictions for samples that were not included in the training set. The RMSE and R² values for the test data were 15.655 and 0.983, respectively. Figure 5b displays the distribution of the predicted bearing capacity values from the MPR model, the actual values, and the errors. The mean, maximum, and minimum errors were 0.92, 38.93, and 0.56 kPa, respectively. It is worth noting that the means and standard deviations of the actual and predicted bearing capacities (q_rs) in the test data were 195.44 kPa and 120.11 kPa and 194.53 kPa and 118.29 kPa, respectively. Table 3 shows the calculated values for each experiment in the test set, as well as the corresponding actual and predicted values.

Table 3. Values for the independent variables and their associated actual vs. predicted results from the MPR model.

Point	qu (kPa)	t/D	d/D	L/dsc	L/dsc	Actual qrs (kPa)	Predicted qrs (kPa)	Individual Error (kPa)
1	7.46	0	0	6	0	15.8	9.36	6.44
2	4.99	0.1	0	6	6	14.11	11.65	2.46
3	48.48	0.1	0	6	6	127.82	135.92	−8.10
4	49.95	0.1	0	6	6	135.24	143.89	−8.65
5	5.18	0.2	0	6	6	20.99	23.75	−2.76
6	30.69	0.2	0	6	6	94.57	86.28	8.29
7	42.6	0.3	0	6	6	134.5	132.15	2.35
8	48.8	0.3	0	6	6	166.92	160.31	6.61
9	5.9	0.4	0	6	6	26.07	49.19	−23.12
10	13.89	0.4	0	6	6	50.71	73.66	−22.95
11	48.26	0.4	0	6	6	168.28	168.84	−0.56
12	52.84	0.4	0	6	6	196.5	195.44	1.06
13	53.55	0.4	0	6	6	203.37	200.12	3.25
14	41.8	0.1	4	6	6	233.57	254.81	−21.24
15	47.63	0.1	4	6	6	297.8	322.93	−25.13
16	51.24	0.1	4	6	6	349.63	372.72	−23.09
17	21.62	0.15	4	6	6	143.78	125.27	18.51
18	53.22	0.15	4	6	6	436.17	419.78	16.39
19	29.28	0.3	4	6	6	201.08	189.21	11.87
20	41.07	0.3	4	6	6	279.93	285.72	−5.79
21	45.02	0.3	4	6	6	317.3	329.10	−11.80
22	53.82	0.3	4	6	6	467.01	450.68	16.33
23	5.9	0.15	1.5	6	6	42.79	30.53	12.26
24	52.34	0.15	1.5	6	6	348.15	327.30	20.85
25	53.16	0.15	1.5	6	6	369.38	336.75	32.63
26	53.82	0.15	2	6	6	420.44	381.51	38.93
27	6.36	0.15	2.5	6	6	49.12	43.82	5.30
28	14.03	0.15	2.5	6	6	92.69	82.57	10.12
29	21.62	0.15	2.5	6	6	135.77	121.82	13.95
30	21.3	0.15	3	6	6	137.76	123.92	13.84
31	36.4	0.15	2	2	2	154.23	147.24	6.99
32	48.77	0.15	2	2	2	239.71	237.26	2.45
33	51.07	0.15	2	2	2	255.44	259.59	−4.15
34	53.82	0.15	2	2	2	281.26	289.04	−7.78
35	8.15	0.15	2	4	4	37.84	43.34	−5.50
36	36.71	0.15	2	4	4	171.08	176.82	−5.74
37	32.21	0.15	2	8	8	185.72	207.77	−22.05
38	46.2	0.15	2	8	8	309.6	336.94	−27.34
39	14.68	0.15	2	6	1.5	44.38	35.71	8.67
40	23.56	0.15	2	6	1.5	70.44	61.65	8.79
41	43.46	0.15	2	6	1.5	147.23	158.34	−11.11
42	51.55	0.15	2	6	1.5	204.22	227.99	−23.77
43	52.7	0.15	2	6	1.5	215.15	239.88	−24.73
44	52.08	0.15	2	6	3	291.01	275.32	15.69
45	29.82	0.15	2	6	4.5	154.14	135.63	18.51
46	47.96	0.15	2	6	4.5	273.42	272.09	1.33
47	51.51	0.15	2	6	4.5	316.7	310.61	6.09
48	53.39	0.15	2	6	4.5	352.47	333.01	19.46

Table 3. Values for the independent variables and their associated actual vs. predicted results from the MPR model.

Point	qu (kPa)	t/D	d/D	L/dsc	L/dsc	Actual qrs (kPa)	Predicted qrs (kPa)	Individual Error (kPa)
1	7.46	0	0	6	0	15.8	9.36	6.44
2	4.99	0.1	0	6	6	14.11	11.65	2.46
3	48.48	0.1	0	6	6	127.82	135.92	−8.10
4	49.95	0.1	0	6	6	135.24	143.89	−8.65
5	5.18	0.2	0	6	6	20.99	23.75	−2.76
6	30.69	0.2	0	6	6	94.57	86.28	8.29
7	42.6	0.3	0	6	6	134.5	132.15	2.35
8	48.8	0.3	0	6	6	166.92	160.31	6.61
9	5.9	0.4	0	6	6	26.07	49.19	−23.12
10	13.89	0.4	0	6	6	50.71	73.66	−22.95
11	48.26	0.4	0	6	6	168.28	168.84	−0.56
12	52.84	0.4	0	6	6	196.5	195.44	1.06
13	53.55	0.4	0	6	6	203.37	200.12	3.25
14	41.8	0.1	4	6	6	233.57	254.81	−21.24
15	47.63	0.1	4	6	6	297.8	322.93	−25.13
16	51.24	0.1	4	6	6	349.63	372.72	−23.09
17	21.62	0.15	4	6	6	143.78	125.27	18.51
18	53.22	0.15	4	6	6	436.17	419.78	16.39
19	29.28	0.3	4	6	6	201.08	189.21	11.87
20	41.07	0.3	4	6	6	279.93	285.72	−5.79
21	45.02	0.3	4	6	6	317.3	329.10	−11.80
22	53.82	0.3	4	6	6	467.01	450.68	16.33
23	5.9	0.15	1.5	6	6	42.79	30.53	12.26
24	52.34	0.15	1.5	6	6	348.15	327.30	20.85
25	53.16	0.15	1.5	6	6	369.38	336.75	32.63
26	53.82	0.15	2	6	6	420.44	381.51	38.93
27	6.36	0.15	2.5	6	6	49.12	43.82	5.30
28	14.03	0.15	2.5	6	6	92.69	82.57	10.12
29	21.62	0.15	2.5	6	6	135.77	121.82	13.95
30	21.3	0.15	3	6	6	137.76	123.92	13.84
31	36.4	0.15	2	2	2	154.23	147.24	6.99
32	48.77	0.15	2	2	2	239.71	237.26	2.45
33	51.07	0.15	2	2	2	255.44	259.59	−4.15
34	53.82	0.15	2	2	2	281.26	289.04	−7.78
35	8.15	0.15	2	4	4	37.84	43.34	−5.50
36	36.71	0.15	2	4	4	171.08	176.82	−5.74
37	32.21	0.15	2	8	8	185.72	207.77	−22.05
38	46.2	0.15	2	8	8	309.6	336.94	−27.34
39	14.68	0.15	2	6	1.5	44.38	35.71	8.67
40	23.56	0.15	2	6	1.5	70.44	61.65	8.79
41	43.46	0.15	2	6	1.5	147.23	158.34	−11.11
42	51.55	0.15	2	6	1.5	204.22	227.99	−23.77
43	52.7	0.15	2	6	1.5	215.15	239.88	−24.73
44	52.08	0.15	2	6	3	291.01	275.32	15.69
45	29.82	0.15	2	6	4.5	154.14	135.63	18.51
46	47.96	0.15	2	6	4.5	273.42	272.09	1.33
47	51.51	0.15	2	6	4.5	316.7	310.61	6.09
48	53.39	0.15	2	6	4.5	352.47	333.01	19.46

4.3. Comparison of the MPR Model with Previously Developed Models

In the field of ML, the development and evaluation of models is an ongoing process, and researchers are constantly striving to improve upon the state-of-the-art. In this context, it is essential to compare the performance of newly developed models with that of existing ones to determine their relative strengths and weaknesses. In this section, a comparison between the performances of the constructed white-box model in this research and previously developed black-box models is presented.

To begin with, we compared the accuracy of our model with that of the black-box models developed in previous research [12]. As it can be noticed in

Table 4. Evaluating the performance of the MPR model vs. prior black-box models for predicting the bearing capacity.

Models	MPR		SVR-ERBF		SVR-RBF		SVR-POLY		ANFIS
Metrics	Train	Test	Train	Test	Train	Test	Train	Test	Train	Test
R2	0.994	0.983	0.999	0.989	0.980	0.975	0.976	0.965	0.996	0.981
RMSE	13.010	15.655	4.30	11.45	17.63	17.51	20.63	23.93	7.90	13.80

, the proposed white-box model in this study achieved an R² value of 0.994 and a RMSE value of 13.010 for training and an R² value of 0.983 and a RMSE value of 15.655 for testing, which placed it in third position among the black-box models in terms of accuracy. While this may appear to be a modest achievement, it is important to note that black-box models require specialized knowledge and expertise to implement and use them effectively, which may limit their accessibility to non-expert users. Furthermore, the interpretability of the proposed white-box model is also compared with that of the black-box model. It can be stated that the white-box model is highly interpretable, which allows the users to gain a better understanding of how the model is making predictions. In contrast, the black-box models are more difficult to interpret, making it challenging to understand the underlying mechanisms of prediction.

Using the validated SVR-ERBF, previous research [12] proposed an empirical approach to generate a design curve and correction equations for this task, but their accuracy was limited. In this study, the novel proposed approach achieves higher accuracy than the previous empirical approach, as demonstrated by its superior R² metric. The proposed model proved to outperform previous research with a remarkable R² value of 0.986 for all data points in the dataset (using Equation (6)), surpassing the empirical approach used in [12] with a R² value of 0.976.

4.4. Parametric Analysis

A parametric analysis was performed to examine the impact of each input feature on the bearing capacity (q_rs) of the composite foundation system studied in this research and to compare the predicted results obtained through the MPR algorithm with those obtained in the laboratory. Each input variable was changed between its minimum and maximum value, while other parameters remained constant with their mean values. The best-fit curves that illustrate the contribution of each input variable on q_rs are displayed in Figure 6a–d, along with a comparison with the experimental values for q_rs and also those predicted by the MPR models.

The effectiveness of models can be evaluated by comparing their predicted values with the actual behavior of the system being analyzed. As shown in Figure 6a–d, the response of the bearing capacity (q_rs) to changes in the predictor properties, such as q_u, t/D, d/D, and h/d_sc ratios, is consistent with laboratory results when using the MPR predictions. Figure 6a shows that q_rs increases as q_u increases, as can be seen in Figure 6a. In addition, an increase in both d/D and h/d_sc ratios leads to an increase in q_rs for the foundation system, as depicted in Figure 6c–d. Conversely, increasing the t/D ratio initially results in an increase in q_rs values, followed by a decrease thereafter for a t/D value of about 0.22. In general, the results indicate that the MPR models perform well, as confirmed by the experimental results and the variation in q_rs due to changes in the factors affecting the bearing capacity of the studied composite ground. Regarding the observed disparity between the experimentally best-fit curve and the scattered experimental points in Figure 6c,d, it should be recognized that this may affect the confidence in the proposed MPR model. However, it is important to note that the primary goal of the MPR model is to capture the overall trend and behavior of the system rather than perfectly fitting every individual data point. The scattered nature of the experimental points may be attributed to various factors, including inherent variability in the material properties and experimental conditions. Despite the observed disparity, the proposed MPR model provides a reasonable approximation of the overall trend and constitutes a valuable tool for understanding the general behavior of the composite foundation system. The authors believe that further research, including refining experimental procedures and incorporating additional data, may help improve the model’s performance in the future.

Figure 6. Comparison of the best-fit curves for bearing capacity (q_rs) based on input features: (a) q_u, (b) t/D, (c) h/d_sc, and (d) d/D, with experimental values and the model predictions from MPR.

Figure 6. Comparison of the best-fit curves for bearing capacity (q_rs) based on input features: (a) q_u, (b) t/D, (c) h/d_sc, and (d) d/D, with experimental values and the model predictions from MPR.

5. Conclusions

The present study focused on the development of an accurate predictive model n for the bearing capacity (q_rs) of composite ground made of geogrid-reinforced sand over geogrid-encased stone columns floating in soft soil using the advanced white-box machine learning technique known as MPR. The MPR system was trained and validated on a comprehensive dataset of 245 experimental tests, using a set of predictors including the ratio of the sandy bed thickness to the foundation diameter (t/D), the ratio of the geogrid layer diameter to the foundation diameter (d/D), the ratio of the stone column length to its diameter (L/d_sc), the ratio of the length of the encasement to the diameter of the stone column (h/d_sc), and the shear strength of unreinforced clay soil (q_u). From the results obtained in this research, the following conclusions can be drawn:

• When tested using a cross-validation approach, the MPR model outperformed both the traditional LR (white-box) and RF (black-box) models in terms of R² and RMSE. The MPR model showed average prediction performances of 0.9850 for R² and 13.9161 for RMSE, while the RF model showed performances of 0.9820 for R² and 19.2054 for RMSE. The LR model performed the worst among all the tested models;
• Further, the MPR model was successfully used to predict the bearing capacity (q_rs) for the testing dataset, with results closely aligned to the experimental values;
• Comparing new and existing models is important in ML. The MPR model performed well, achieving high accuracy and interpretability when compared to traditional white-box and black-box models, which require specialized knowledge;
• In addition, a parametric analysis was performed to assess the impact of increasing the predictors on q_rs. The findings revealed that the MPR predictions align well with laboratory results, indicating that an increase in q_u values leads to an increase in q_rs values. Similarly, an increase in both d/D and h/d_sc ratios also results in an increase in q_rs values. However, an increase in the t/D ratio shows a different effect on q_u values, which initially increase but then decrease. These findings can help in understanding how different factors impact the bearing capacity of composite foundation systems, such as the one studied in this research.

Despite the promising results, the proposed model has some limitations that should be considered in future research. Firstly, the current model does not include water level as an independent parameter, which could play a significant role in the bearing capacity and overall performance of the foundation system. Future studies should investigate the influence of water level on the bearing capacity and consider incorporating it as an input variable in the model. Secondly, the model was developed based on a dataset consisting of 245 experimental results collected from the literature. To further improve the accuracy and generalization performance of the model, it would be beneficial to include additional data points and investigate the influence of other relevant factors on the bearing capacity of the foundation system. Lastly, the performance of the model could be further enhanced by exploring other advanced machine learning techniques, such as deep learning or ensemble methods, which may provide improved predictive capabilities and a better understanding of the complex interactions between the input variables.

In conclusion, this study presents a novel approach to predicting the bearing capacity of composite ground using a white-box machine learning model, demonstrating its potential as a valuable tool for geotechnical engineers and researchers. Future work should focus on addressing the limitations of the current model and exploring alternative modeling techniques to provide a more comprehensive understanding of the factors influencing the bearing capacity of composite grounds.