Prediction of Ultimate Bearing Capacity of Shallow Foundations on Cohesionless Soils: A Gaussian Process Regression Approach

Abstract

This study examines the potential of the soft computing technique—namely, Gaussian process regression (GPR), to predict the ultimate bearing capacity (UBC) of cohesionless soils beneath shallow foundations. The inputs of the model are width of footing (B), depth of footing (D), footing geometry (L/B), unit weight of sand (γ), and internal friction angle (ϕ). The results of the present model were compared with those obtained by two theoretical approaches reported in the literature. The statistical evaluation of results shows that the presently applied paradigm is better than the theoretical approaches and is competing well for the prediction of UBC (q_u). This study shows that the developed GPR is a robust model for the q_u prediction of shallow foundations on cohesionless soil. Sensitivity analysis was also carried out to determine the effect of each input parameter.

1. Introduction

Ultimate bearing capacity (UBC) and allowable settlement are two important criteria to consider when designing shallow foundations. The UBC is governed by the shear strength of the soil and is estimated by theories proposed by Terzaghi [1], Meyerhof [2], Hansen [3], Vesic [4], etc. However, the various bearing capacity equations reveal a wide range of variations. Furthermore, the proposed bearing capacity theories include a number of assumptions that contribute to simplifying the problems [5].

Numerous studies proposed numerical approaches for estimating bearing capacity in addition to semi-empirical solutions for determining the bearing capacity of foundations, e.g., [6,7]. Models are usually validated using the model-scale footing test. A number of researchers have investigated that how to reduce scale effects when extrapolating experimental results to full-scale footings (e.g., de Beer [8]; Yamaguchi et al. [9]). In model-scale footing experiments, Tatsuoka et al. [10] studied how particle size affects the UBC. The shearing strains vary significantly along the slip line, according to the results of large-scale footing tests on dense sand. As a result, bearing capacity formulas that use the maximum value of friction angle, (ϕ_max) tend to overstate prototype bearing capacities [11]. Therefore, the actual model-scale footing test findings differed from theoretical equations, and special consideration should be given when comparing model-scale footing test results with full-scale foundation behavior. For this reason, a new approach is needed to ensure more accurate predictions of actual bearing capacity.

Intelligent systems are typically used to model complex interactions between inputs and outputs or to explore patterns in available data. Artificial intelligence (AI)-based methods are capable of capturing the problem’s inherent nonlinearity and complex interaction between the involved variables in various domains, e.g., [12,13,14,15,16,17,18,19,20,21]. These approaches can be trained to learn the relationship between soil mechanical properties and foundation geometry with foundation bearing capacity, and no prior knowledge of the relationship’s form. Different forms of AI-based techniques have recently been used by different researchers to address the UBC problem of shallow foundations. Artificial neural networks (ANNs), fuzzy inference systems (FISs), adaptive neuro-fuzzy inference systems (ANFISs), ant colony optimization (ACO), genetic programming (GP), weighted genetic programming (WGP), soft-computing polynomials (SCP), support vector machine (SVM), random forest (RF), and relevance vector machine (RVM) have all been used to successfully estimate the UBC of shallow foundations on soil [11,22,23,24,25,26]. Soft computing methodologies are more accurate than analytical formulas, according to all of these studies. The findings revealed that the ML models mentioned above are capable of obtaining the experimental observations with acceptable accuracy. However, this field continues to be further explored.

The Gaussian process regression (GPR) approach has been successfully applied in many domains, e.g., [27,28,29], but its application in geotechnical engineering is limited based on literature surveys. Considering the improved performance of GPR, it is, however, used for the first time in this study to predict the UBC of shallow foundations on cohesionless soils. To demonstrate the efficacy of the proposed GPR-based model, the results were compared with various well-known classical formulas for calculating the UBC.

The main contributions of this paper are as follows:

• To examine the capability of the GPR model for the prediction of q_u of shallow foundation on cohesionless soil;
• To undertake a comparative study with the commonly used bearing capacity theories;
• To conduct sensitivity analyses for the determination of the effect of each input parameter on q_u.

The structure of the paper is as follows: In Section 2, the theoretical background of ultimate bearing capacity is presented. The numerical model and verification are described in Section 3. Section 4 presents the construction process of the prediction model. In Section 5, results and discussion are described. Finally, the concluding remarks are presented.

2. Theoretical Background of Ultimate Bearing Capacity

Prandtl [30] and Reissner [31] developed plastic theory-based methods for determining the ultimate bearing capacity of shallow strip footings. Over the years, other researchers significantly improved the formulation [1,2,4,22,32,33]. Terzaghi [1] improved Prandtl’s theory [30] by using the principle of superposition to account for the soil weight. Taylor [32] incorporated the effect of the overburdened soil surcharge at the foundation level into Prandtl’s formulation. Meyerhof [2] extended Terzaghi’s bearing capacity equation by incorporating different shape and depth factors. Hansen [3] later updated the Meyerhof model. Vesic [4] provided a bearing capacity prediction equation that was similar to Hansen’s equation.

Table 1. General forms of the classical prediction equations for the bearing capacity of shallow foundations.

Reference	Equation
Terzaghi [1]	$q_u=c{N}_c{s}_c+\gamma D{N}_q+\frac{1}{2}\gamma B{N}_{\gamma }{s}_{\gamma }$ $N_c=\left(N_q-1\right)\cot ɸ$ $\begin{array}{l}{N}_q=\frac{a^2}{a{\cos }^2\left(45+ɸ/2\right)}\\ a=e^{\left(0.75\pi -ɸ/2\right)\tan ɸ}\\ N_{\gamma }=\frac{\tan ɸ}{2}\left(\frac{K_{p\gamma }}{{\cos }^2ɸ}-1\right)\end{array}$
Meyerhof [2]	$q_u=c{N}_c{s}_c{d}_c+\gamma D{N}_q{s}_q{d}_q+\frac{1}{2}\gamma B{N}_{\gamma }{s}_{\gamma }{d}_{\gamma }$ $\begin{array}{l}{N}_q=e^{\pi \tan ɸ}{\tan }^2\left(45+\frac{ɸ}{2}\right)\\ N_{\gamma }=\left(N_q-1\right)\tan \left(1.4ɸ\right)\end{array}$ $N_c=\left(N_q-1\right)\cot ɸ$
Hansen [3]	$q_u=c{N}_c{s}_c{d}_c{i}_c{g}_c{b}_c+\gamma D{N}_q{s}_q{d}_q{i}_q{g}_q{b}_q+\frac{1}{2}\gamma B{N}_{\gamma }{s}_{\gamma }{d}_{\gamma }{i}_{\gamma }{g}_{\gamma }{b}_{\gamma }$ $N_q$ = same as Meyerhof above $N_c$ = same as Meyerhof above $N_{\gamma }=1.5\left(N_q-1\right)\tan ɸ$
Vesic [4]	Same as Hansen’s equation $N_q$ = same as Meyerhof above $N_c$ = same as Meyerhof above $N_{\gamma }=2\left(N_q+1\right)\tan ɸ$

Notes: q_u: ultimate bearing capacity of footing; c: cohesion; γ: average effective unit weight of the soil below and around the foundation; B: width of footing; D: depth of footing; N_c, N_q, and N_γ: non-dimensional bearing capacity factors as exponential functions of ϕ; ϕ: internal friction angle; s_c, s_q, and s_γ: non-dimensional shape factors; i_c, i_d, and i_γ: non-dimensional inclination factors; d_c, d_q, and d_γ: non-dimensional depth factors; g_c, g_q, and g_γ: non-dimensional ground factors (base on slope); b_c, b_q, and b_γ: base factors (tilted base).

summarizes the general forms of the classical equations. The form of the equations presented by the other researchers remained the same as Terzaghi’s, as can be seen in this table. Meyerhof [2], Hansen [3], Vesic [4], and other researchers proposed different shape, depth, inclination, ground, and base factors for the bearing capacity equations after thorough in situ and laboratory tests. However, one of the primary shortcomings of these traditional formulations is that they are based on some simplifying assumptions. As a result, they do not always produce accurate bearing capacity estimates [22,34].

3. Materials and Methods

3.1. Dataset

The data used in this study were adopted from Padmini et al. [11]. The five input parameters used for model development in this study were width of footing (B), depth of footing (D), footing geometry (L/B), unit weight of sand (γ), and internal friction angle (ϕ). Ultimate bearing capacity (q_u) was the single output. The data thus compiled comprised a total of 97 datasets, which consisted of results of load test data of square, rectangular, and strip footings of different sizes tested in sand beds of various densities.

The data were divided into training and testing sets, which is a method that has a substantial impact on the results when using data mining techniques [35,36]. In this case, the entire database was divided into multiple random combinations of training and testing sets until both training and testing sets had a robust representation of the entire population. A statistical analysis of the input and output parameters of the randomly selected training and testing sets was carried out to identify the most robust representation. The purpose of this analysis was to assure that the statistical properties of the data in each of the subsets were as close as possible, indicating that they represented the same statistical population. This was achieved using a trial-and-error method in this work. For the construction and testing of the GPR model, the best statistically consistent combination was selected. The data division was performed in such a way that 78 (80%) were used for training, and 19 (20%) were used for testing in all the experiments considered in this study. The parameters used in the statistical analysis included maximum, minimum, mean, and standard deviation. The results of the statistical analysis of the finally selected combinations are shown in

Table 2. Statistical aspects of the dataset.

Dataset	Statistical Parameter	B (m)	D (m)	L/B	γ (kN/m3)	ϕ (°)	qu (kPa)
Training	Minimum	0.0585	0	1	9.85	32	58.5
	Maximum	3.016	0.889	6	17.1	44.8	2847
	Mean	0.398	0.156	3.126	14.279	38.492	407.899
	Standard deviation	0.543	0.195	2.192	2.642	3.308	511.870
Testing	Minimum	0.0585	0	1	10.2	32	91.5
	Maximum	1.492	0.762	6	17.1	44.8	2033
	Mean	0.420	0.204	2.995	13.777	38.800	569.826
	Standard deviation	0.361	0.221	1.964	2.573	3.225	600.020

. It should be mentioned that due to the fact that the data contained singular, rare events that could not be replicated in all cases of the dataset, there may still have been some minor inconsistencies in the statistical parameters for the training and validation sets. The descriptive statistics (such as minimum and maximum values, mean, and standard deviations) of the selected UBC parameters with the established GPR model are provided in Table 2 (the complete database is available in

Table A1. The data used for developing the GPR model.

S. No.	B (m)	D (m)	L/B	γ (kN/m3)	ϕ (°)	qu (kPa)
1	0.6	0.3	2	9.85	34.9	270
2	0.6	0	2	10.85	44.8	860
3	0.152	0.075	5.95	17.1	42.5	335.3
4	0.0585	0.058	5.95	16.8	41.5	184.9
5	0.5	0	1	11.7	37	111
6	0.5	0	2	11.7	37	143
7 *	0.5	0.3	1	10.2	37.7	681
8	0.152	0.15	1	16.1	37	182.4
9	0.0585	0.058	5.95	15.7	34	70.91
10 *	1	0	3	11.93	40	630
11	0.152	0.15	5.95	15.7	34	122.3
12	0.152	0.15	1	16.8	41.5	361.5
13 *	0.5	0.5	4	12	40	1140
14	0.5	0.029	4	11.7	37	109
15	0.5	0	1	10.2	37.7	154
16	0.094	0.094	6	17.1	42.5	279.6
17	0.094	0.047	1	16.1	37	98.8
18	0.5	0.5	2	12.41	44	2847
19	0.5	0	4	12.41	44	797
20	0.0585	0.029	5.95	16.1	37	82.5
21 *	0.152	0.075	1	16.1	37	135.2
22	0.5	0.3	1	12.41	44	1940
23 *	0.152	0.15	1	16.5	39.5	264.5
24	0.5	0	2	10.2	37.7	195
25 *	0.6	0.3	2	10.85	44.8	1760
26	0.991	0.711	1	15.8	32	1773.7
27	0.5	0.3	1	11.7	37	446
28	0.5	0	1	10.2	37.7	165
29	0.5	0	1	11.77	37	123
30	0.5	0.49	4	12.27	42	1492
31	0.152	0.15	5.95	16.8	41.5	342.5
32	0.094	0.094	6	16.5	39.5	185.6
33	0.5	0	1	11.7	37	132
34	0.5	0	2	10.2	37.7	203
35 *	0.5	0.3	4	11.7	37	322
36	0.094	0.094	1	15.7	34	90.5
37	0.094	0.047	6	16.8	41.5	206.8
38	0.094	0.047	6	16.1	37	104.8
39	2.489	0.762	1	15.8	32	1158
40	0.094	0.094	1	16.5	39.5	191.6
41	0.152	0.075	5.95	15.7	34	98.2
42	0.094	0.047	1	16.5	39.5	147.8
43	0.5	0.127	4	11.7	37	187
44 *	0.5	0.5	4	12.41	44	2033
45 *	0.094	0.094	1	16.8	41.5	253.6
46	0.5	0.3	3	10.2	37.7	402
47	3.004	0.762	1	15.8	32	1019.4
48	0.152	0.075	1	16.8	41.5	276.3
49	0.094	0.047	6	17.1	42.5	235.6
50	0.152	0.15	5.95	16.1	37	176.4
51	0.6	0.3	2	10.2	37.7	570
52	0.5	0.3	1	11.77	37	370
53	0.5	0.3	2	10.2	37.7	530
54	0.094	0.047	1	16.8	41.5	196.8
55 *	0.094	0.094	6	16.8	41.5	244.6
56	0.094	0.094	1	17.1	42.5	295.6
57	0.152	0.15	5.95	16.5	39.5	254.5
58	0.52	0	3.85	10.2	37.7	186
59	0.152	0.15	5.95	17.1	42.5	400.6
60	1	0.2	3	11.97	39	710
61	0.0585	0.029	5.95	15.7	34	58.5
62 *	0.5	0.013	1	11.7	37	137
63	0.5	0.3	1	12.41	44	2266
64	0.0585	0.029	5.95	16.5	39.5	121.5
65	0.5	0.5	4	11.7	37	425
66	0.094	0.094	6	16.1	37	127.5
67 *	0.0585	0.029	5.95	17.1	42.5	180.5
68	0.152	0.075	1	15.7	34	91.2
69 *	0.152	0.075	5.95	16.1	37	143.3
70	0.094	0.047	1	15.7	34	67.7
71	0.152	0.075	1	17.1	42.5	325.3
72	0.094	0.047	1	17.1	42.5	228.8
73	0.5	0.3	2	10.2	37.7	542
74	0.094	0.047	6	15.7	34	74.7
75	0.094	0.047	6	16.5	39.5	155.8
76 *	1.492	0.762	1	15.8	32	1540
77	0.5	0	1	12.41	44	782
78	0.6	0	2	10.2	37.7	200
79	0.152	0.075	5.95	16.5	39.5	211.2
80	0.5	0.5	2	11.7	37	565
81	0.152	0.15	1	17.1	42.5	423.6
82 *	0.5	0.5	2	11.77	37	464
83	0.152	0.075	5.95	16.8	41.5	285.3
84 *	0.094	0.094	1	16.1	37	131.5
85	0.0585	0.058	5.95	17.1	42.5	211
86	0.152	0.075	1	16.5	39.5	201.2
87	0.0585	0.058	5.95	16.5	39.5	142.9
88	0.0585	0.058	5.95	16.1	37	98.93
89	0.152	0.15	1	15.7	34	124.4
90 *	0.5	0	3	10.2	37.7	214
91	0.0585	0.029	5.95	16.8	41.5	157.5
92 *	0.094	0.094	6	15.7	34	91.5
93 *	0.5	0	4	12	40	461
94	0.52	0.3	3.85	10.2	37.7	413
95	3.016	0.889	1	15.8	32	1161.2
96	0.5	0	2	11.77	37	134
97	0.5	0.3	1	11.7	37	406

Note: * presents the test set.

).

3.2. Correlation Analysis

The strength of the correlation between the various factors was verified using correlation coefficients (ρ) (see

Table 3. Correlation coefficients between various factors.

Parameter	B	D	L/B	γ	ϕ	qu
B	1.000
D	0.710	1.000
L/B	−0.351	−0.249	1.000
γ	−0.269	−0.125	0.340	1.000
ϕ	−0.378	−0.286	0.124	0.076	1.000
qu	0.452	0.671	−0.248	−0.238	0.258	1.000

). The following formula for ρ is given a pair of random variables (x,y):

$$\rho \left(x,y\right)=\frac{\mathrm{cov}\left(x,y\right)}{{\sigma }_x{\sigma }_y}$$

(1)

where cov is the covariance, σ_x is the standard deviation of x, and σ_y is the standard deviation of y. |ρ| > 0.8 indicates a strong correlation between x and y, values between 0.3 and 0.8 indicates a moderate connection, and |ρ| < 0.30 indicates a weak correlation [37]. A correlation is regarded as “strong” if |ρ| > 0.8, according to Song et al. [38]. Table 3 shows the correlations between B, D, L/B, γ, ϕ, and q_u, in order of moderate to weakest. As a result, no factors were removed from the UBC of the cohesionless soil estimation model. Table 3 shows that the maximum absolute value of the correlation coefficient is 0.710, and there is no “strong” correlation between the various pairs of factors.

3.3. Gaussian Processes Regression

Gaussian processes regression (GPR) is a suitable and recently proposed method that has been used in a variety of machine learning applications [39]. The GPR model’s probabilistic solution leads to the identification of generic regression problems using kernels. The applied regressor’s training process can be categorized as Bayesian framing, and the model relations are presumed to follow a Gaussian distribution to encode the previous information about the output function [40]. The Gaussian process is defined by a set of variables for which each one has a joint Gaussian allocation [41]. The overall structure of the Gaussian process is described by the following equation:

$$g\left(x\right)\sim GP\left(w\left(x\right),k\left(x,x^{\prime }\right)\right)$$

(2)

The mean function of the Gaussian process is denoted by $w\left(x\right)$, and the kernel function is denoted by $k\left(x,x^{\prime }\right)$. The mean function is typically constant, either zero or the mean of the training dataset. In this study, the mean of the training dataset was used. Consider a learning dataset with N pairs in the form of S = {(x_i, y_i)|i = 1, 2, …, N} where x is an N-dimensional input vector, and y is the corresponding target. Using x_j as a test sample, the GPR model utilizes the below formulation in determining the relationship between the provided inputs and targets [41].

$$w_j=k_j^T\left[K\left(X,X\right)+{\sigma }^2{}_{n}I\right]y$$

(3)

$${\sigma }_j^2=k\left(x_j,x_j\right)-k_j^T{\left[K\left(X,X\right)+{\sigma }_n^{2}I\right]}^{-1}{k}_j$$

(4)

where $w_j$ denotes the mean value of the most compatible predicted outputs for the test input vector (x_j). Additionally, $K\left(X,X\right)$, $k_j$, ${\sigma }_n^2$ and $y$ represent covariance matrix, the kernel distance between training and testing data, the noise variance, and the training observation, respectively. Furthermore, the produced variance $\left({\sigma }_j^2\right)$ by Equation (4) indicates a confidence measure for the acquired results. It is worth noting that this variance is inversely proportional to the confidence associated with the $w_j$ [42]. The formulas above can be combined in the form of a linear combination of the kernel function and the mean estimation $\overline{f}\left(x_j\right)$ as follows:

$$\overline{f}\left(x_j\right)={\displaystyle \sum _{i=1}^s{\left(K\left(X,X\right)+{\sigma }^2{}_{n}I\right)}^{-1}yk\left(x_j,x_j\right)}$$

(5)

GPR utilizes a number of kernel functions. GPR has a limitation in terms of selecting a suitable kernel function. Pearson VII kernel function (PUK), a widely utilized kernel function, was chosen for GPR model construction in this study, as it has shown optimized results in predicting river discharge [43].

$$PUK=\left(1/{\left[1+{\left(2{\sqrt{\Vert x_i-x_j\Vert }}^2\sqrt{2^{\left(1/\omega \right)}-1}/\sigma \right)}^2\right]}^{\omega }\right)$$

(6)

where ω and σ are kernel parameters.

3.4. Model Evaluation and Comparison

To validate and compare the developed model, four quantitative statistics were used to assess the performance of the evaluation methods: coefficient of determination (R²), the ratio of the root-mean-square error (RSR) to the standard deviation of measured data, Nash–Sutcliffe coefficient (NSE), and mean bias error (MBE). The following equations are used to express these indices:

$$R^2=1-\frac{{\displaystyle \sum _{i=1}^n{\left(y_i-{\widehat{y}}_i\right)}^2}}{{\displaystyle \sum _{i=1}^n{\left(y_i-\widehat{y}\right)}^2}}$$

(7)

$$RSR=\frac{\sqrt{{\displaystyle \sum _{i=1}^n{\left(y_i-{\widehat{y}}_i\right)}^2}}}{\sqrt{{\displaystyle \sum _{i=1}^n{\left(y_i-\overline{y}\right)}^2}}}$$

(8)

$$NSE=1-\frac{{\displaystyle {\sum }_{i=1}^n{\left(y_i-{\widehat{y}}_i\right)}^2}}{{\displaystyle {\sum }_{i=1}^n{\left(y_i-\overline{y}\right)}^2}}$$

(9)

$$MBE=\frac{1}{n}{\displaystyle {\sum }_{i=1}^n\left({\widehat{y}}_i-y_i\right)}$$

(10)

In the equations, $n$ is the number of data, $y_i$ and ${\widehat{y}}_i$ are the actual and predicted output of ith sample of the data, respectively; $\overline{y}$ is the averaged actual output of the data. The R² coefficient goes from 0 to 1, with a higher R² value indicating a more efficient model. When R² is larger than 0.8 and nearer to 1, the model is considered effective [44]. The NSE scale runs from −∞ to 1, with 1 being a perfect match. An NSE value of more than 0.65 indicates a strong correlation [45,46]. The root-mean-square error (RMSE) observations’ standard deviation ratio is calculated as the ratio of the RMSE and standard deviation of measured data. The RSR varies from an optimal value of 0 to a significant positive value. A lower RSR signifies a lower RMSE, indicating the higher predictive efficiency of the model. RSR classification ranges were described as very good, good, acceptable and unacceptable with ranges of 0.00 ≤ RSR ≤ 0.50, 0.50 ≤ RSR ≤ 0.60, 0.60 ≤ RSR ≤ 0.70 and RSR > 0.70, respectively [47]. It is self-evident that the lower the RSR criteria is, the better is the model.

4. Construction of Prediction Model

The construction process of the prediction model is shown in Figure 1. First, 80% and 20% of the UBC data were selected based on statistical consistency as training and test sets, respectively. Second, using the optimal hyperparameters configuration, the prediction model was fitted based on the training set using the trial-and-error method. We adjusted the hyperparameters so as to maximize the likelihood of the training data. We initialized the hyperparameters to random values (in a reasonable range) and then used an iterative method to search for optimal values of the hyperparameters. Numerous trials were carried out to find optimal values of primary kernel parameters—omega (ω) and sigma (σ)—and the values were 1 in the GPR model. Fourth, the test set was adopted to judge the performance of the proposed GPR model according to the four standard statistical measures: R², RSR, NSE, and MBE. Higher values of R² and NSE, and lesser values of RSR, indicate a better estimation accuracy of the proposed model. Last, the developed model was compared with the commonly used bearing capacity theories by comparing the comprehensive performance if the prediction performance of this model was acceptable, then it could be adopted for deployment and made sensitivity analysis. The entire calculation process was performed in Waikato Environment for Knowledge Analysis (WEKA) software [48]. In this study, for the GPR model, the Pearson VII function-based universal kernel [49] was used.

5. Results and Discussion

The performance of the GPR model was evaluated using root-mean-square error (RMSE), coefficient of determination (R²), the ratio of the root-mean-square error to the standard deviation of measured data (RSR), and Nash–Sutcliffe coefficient (NSE) and mean bias error (MBE). The results were compared with some popular classical methods suggested in the literature (i.e., Vesic [4]; Hansen [3]) for determining the UBC of shallow foundation on cohesionless soil. The GPR model in Section 3.3 was developed by learning from the 78 data, and the results of the training performance including R², RSR, NSE, and MBE are shown in

Table 4. Performance statistics of proposed model.

Performance Index	Training	Testing
R2	0.9552	0.9404
RSR	0.3428	0.3975
NSE	0.8825	0.8420
MBE (kPa)	−10.9661	−69.3268

. In comparison to the theoretical equations, the output of the GPR model on the testing phase has substantially higher values of coefficient of determination (R²), Nash–Sutcliffe coefficient (NSE), and lower values of MBE and RSR. In addition, the equation proposed by Vesic [50] shows the best performance in comparison to the Hansen [3] theoretical formulas. Figure 2 shows the best-fit line of estimated versus measured UBC and the associated coefficient of determination (R²). The outputs of the used theoretical formulas are more scattered than the GPR-based model, as shown in this Figure 2. From Table 4, it can be inferred that the R² and NSE are the highest (0.9404 and 0.8420), and the other criteria such as RSR and MBE are the least (0.8420 and −69.3268) for the GPR-based modeling in the testing phase. MBE (Vesic [4]: −105.2996 and Hansen [3]: −161.2911) has the largest negative score (see

Table 5. Performance statistics of bearing capacity theories.

Performance Index	Testing Set
	Vesic [4]	Hansen [3]
R2	0.9017	0.8944
RSR	0.4303	0.5227
NSE	0.8149	0.7268
MBE (kPa)	−105.2996	−161.2911

), indicating that it is the most conservative approach. It is clear that the GPR model is more accurate than analytical formulas. The model performs better when R² and NSE values are greater than 0.85. Thus, it depicts the accuracy of the finalized model.

Figure 3 shows the comparison results of the training and testing sets to characterize the difference between the actual and predicted UBC of shallow foundations on cohesionless soil. Good agreements can be also observed between the comparison results, except for the few noise points, but these results are acceptable for the proposed GPR model to predict the UBC of shallow foundation on cohesionless soil.

The proposed model was developed and validated using data inside the data range (see Table 2), but it was not tested outside of this range; hence, it is assumed that the suggested model will be valid for data inside within the data range.

The GPR model’s sensitivity results were evaluated using Yang and Zang’s [51] technique to assess the influence of input parameters on q_u. This method has been utilized in a number of studies [14,52,53,54], and it is stated as follows:

$$r_{ij}=\frac{{\displaystyle {\sum }_{m=1}^n\left(y_{im}\times y_{om}\right)}}{\sqrt{{\displaystyle {\sum }_{m=1}^n{y}_{im}{}^2{\displaystyle {\sum }_{m=1}^n{y}_{om}{}^2}}}}$$

(11)

where $n$ is the number of data values (this study used 78 data values); $y_{im}$ and $y_{om}$ are the input and output parameters. The $r_{ij}$ value ranges from zero to one for each input parameter, and the highest $r_{ij}$ values suggested the most efficient output parameter (which was q_u in this study). The $r_{ij}$ values for all input parameters are presented in Figure 4. Figure 4, shows that the depth of foundation, D $\left(r_{ij}=0.817\right)$ has the greatest effect on the UBC. The sensitivity results match well with the result of the investigations by Meyerhof [2] and Kohestani et al. [26].

6. Conclusions

In this paper, using the GPR model, values of UBC of shallow foundations on cohesionless soil were estimated. Five input variables and one output variable were used for designing the prediction model. The modeling results reveal that the GPR model has an appropriate capability for accurate estimation of the UBC of cohesionless soil. The GPR-based model also provides improved performance compared with the conventional methods considered in this study. Results of sensitivity analysis conclude that the depth of foundation is the primary essential parameter when the GPR-based model is selected for estimation of the UBC of cohesionless soil for this dataset. For this dataset, the Pearson VII kernel function-based GPR model performs better and can be used in a variety of geotechnical engineering problems with inherent uncertainties and imperfections. Furthermore, the GPR technique has the advantage of being quickly updated as new data becomes available, avoiding the requirement for expertise and time to update the old design aid or equation and/or propose a new equation.