*Article* **Developing Two Hybrid Algorithms for Predicting the Elastic Modulus of Intact Rocks**

**Yuzhen Wang 1,2, Mohammad Rezaei 3,\*, Rini Asnida Abdullah <sup>4</sup> and Mahdi Hasanipanah <sup>4</sup>**


**Abstract:** In the primary and final designs of projects related to rock mechanics and engineering geology, one of the key parameters that needs to be taken into account is the intact rock elastic modulus (E). To measure this parameter in a laboratory setting, core samples with high-quality and costly tools are required, which also makes for a time-consuming process. The aim of this study is to assess the effectiveness of two meta-heuristic-driven approaches to predicting E. The models proposed in this paper, which are based on integrated expert systems, hybridize the adaptive neurofuzzy inference system (ANFIS) with two optimization algorithms, i.e., the differential evolution (DE) and the firefly algorithm (FA). The performance quality of both ANFIS-DE and ANFIS-FA models was then evaluated by comparing them with ANFIS and neural network (NN) models. The ANFIS-DE and ANFIS-FA models were formed on the basis of the data collected from the Azad and Bakhtiari dam sites in Iran. After applying several statistical criteria, such as root mean square error (RMSE), the ANFIS-FA model was found superior to the ANFIS-DE, ANFIS, and NN models in terms of predicting the E value. Additionally, the sensitivity analysis results showed that the P-wave velocity further influenced E compared with the other independent variables.

**Keywords:** elastic modulus; ANFIS; differential evolution; firefly algorithm

### **1. Introduction**

When planning most projects pertinent to geotechnical issues and rock engineering, it is of high importance to properly analyze how the intact rock behaves and carefully estimate its associated mechanical properties. The intact rock elastic modulus (E) has substantial effects on both the initial and final steps of designing geoscience-related projects, which include planning tunnels; designing blasting operations in rock materials; analyzing the constancy of rock slopes; and designing rock pillars, roads, dams, bridges, etc. Moreover, E is the most significant parameter applied to analyzing the stress-strain chart of rock specimens in a laboratory. E also plays an important role in analyzing the deformations and breakage of rocks surrounding underground excavation projects. As a result, inaccurate predictions of E can result in serious damages, leading to economic issues and severe safety problems due to the breakage probability during construction processes [1,2]. Thus, it is necessary to determine the E value quickly and accurately in order to correctly plan geoengineering structures, accurately design mining- and civil engineering-related projects, and enhance the general safety level and effectiveness of operations at hand.

In general, the E value can be obtained using direct or indirect methods. The former are typically carried out within rock mechanics laboratories, where rock core specimens are subjected to experiments in a variety of conditions [1–3]. In contrast, indirect methods make

**Citation:** Wang, Y.; Rezaei, M.; Abdullah, R.A.; Hasanipanah, M. Developing Two Hybrid Algorithms for Predicting the Elastic Modulus of Intact Rocks. *Sustainability* **2023**, *15*, 4230. https://doi.org/10.3390/ su15054230

Academic Editor: Jianjun Ma

Received: 30 January 2023 Revised: 22 February 2023 Accepted: 23 February 2023 Published: 26 February 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

use of predictive equations or models to estimate E. The direct methods have accuracy, but at the same time, they suffer from some drawbacks. First, it is not easy to provide the required specimens during the coring process with a high level of accuracy, especially in jointed, layered, and weakened rock structures. Second, it is both difficult and time consuming to prepare the core specimens with the appropriate geometry for the purpose of carrying out laboratory E tests. Such issues hinder the use of direct methods unless there is a high necessity [3–5].

Due to the above-mentioned challenges, various indirect methods have been introduced in the literature on the basis of predictive models/algorithms and equations to determine the E value of intact rocks. These methods have been normally configured based on arithmetical and smart intelligent models. Statistical models generally make use of simple or multiple regression models aiming to develop a number of empirical equations between the E value and effective mechanical and physical rock properties. The literature consists of numerous empirical equations formed based on analyzing petrology and mineralogy in addition to values estimated using the rock physical and mechanical properties, such as Schmidt hammer numbers [6], porosity of rock [7], slake durability of rock [8], and compressional/primary wave velocity [9].

In recent years, scholars have made several efforts to develop artificial intelligence (AI) models applicable to mining and rock engineering problems. Such efforts have resulted in a number of novel models proposed to estimate E and some other rock mechanical properties. These are mostly based on probable and intelligent methods, such as particle swarm optimization (PSO), fuzzy inference systems (FISs), genetic algorithms (GAs), Bayesian methods, adaptive neuro-fuzzy inference systems (ANFIS), tree models, extreme gradient-boosting (XGB), and artificial neural networks (ANNs), as well as their hybridized forms [10–16]. Sarkhani Benemaran et al. [17] employed an XGB model in combination with several optimization algorithms to predict the resilient modulus of flexible pavement foundations. They concluded the effectiveness of PSO-XGB models in this field. In another study, conducted by Shahani et al. [18], different AI models such as XGB, gradient-boosted tree regressors (GBTRs), Catboost, and light gradient-boosting machines were used to predict E. According to their results, the performance of GBTR was better than that of the other developed models. Recently, Tsang et al. [19] predicted the E values through some other models, i.e., extreme gradient-boosting trees, ANNs, random forests, and classification and regression trees. The results showed that the extreme gradient-boosting trees predicted the E value with the highest accuracy.

Such applications show that intelligent algorithms typically outperform the traditionally used statistical methods regarding E value prediction.

The present study is carried out to assess the potential of applying two hybrid evolutionary models to predict E. The proposed models are based on the integrated expert systems comprising ANFIS with two optimization algorithms, i.e., the firefly algorithm (FA) and differential evolution (DE). To check the effectiveness of FA and DE, the results of ANFIS-FA and ANFIS-DE are then compared with the ANFIS and NN results. The rest of this study is organized as follows. In Section 2, the source of the database is described. Then, the methodologies used in this paper and their implementations are explained in detail in Section 3. Finally, Sections 4 and 5 present the results/discussions and conclusions of this study, respectively.

#### **2. Source of Database**

An inclusive database is needed to be formed for E modeling by means of indirect intelligent methods. Such data were obtained through performing laboratory experiments on the core specimens provided from the excavation drilling processes carried out in two under-construction dam sites, namely the Bakhtiari and Azad dams located in Iran. The precise location of the Azad dam site is in the west of Iran, 40 km away from the western city of Sanandaj in the Kurdistan state. It is situated on the Sanandaj–Marivan cities road inside Kurdistan, with the 46◦3205700 to 35◦1905900 geographical coordinates of eastern and

northern latitudes, respectively. The construction of this dam is currently in progress. It is mainly aimed at supplying electrical energy and producing power plant storage. The dimensions (length, height, and width) and water storage capacity of this dam are 595 m, 115 m, 11 m, and 260,000,000 m<sup>3</sup> , respectively. It is mainly aimed at supplying electrical energy and producing power plant storage. The dimensions (length, height, and width) and water storage capacity of this dam are 595 m, 115 m, 11 m, and 260,000,000 m<sup>3</sup> , respectively. The Bakhtiari dam is located in the Zagros Highlands, 65 km southwest of the town

precise location of the Azad dam site is in the west of Iran, 40 km away from the western city of Sanandaj in the Kurdistan state. It is situated on the Sanandaj–Marivan cities road inside Kurdistan, with the 46° 32′ 57″ to 35° 19′ 59″ geographical coordinates of eastern and northern latitudes, respectively. The construction of this dam is currently in progress.

*Sustainability* **2023**, *15*, 4230 3 of 24

The Bakhtiari dam is located in the Zagros Highlands, 65 km southwest of the town of Dorud in the Lorestan state, and 70 km northeast of the town of Andimeshk in Khuzestan, Iran. The position of the dam site is at the 48◦45034.8700 to 32◦57023.5800 geographical coordinates of eastern and northern latitudes, respectively (Figure 1). The dam was built upon the Bakhtiari River, aiming to provide adequate water for many purposes, such as drinking, electrical power generation, flood control, and agricultural activities. The dam's body is at an elevation of 840 m. In addition, in the case of this dam, the peak elevation, crown width, crown length, and foundation width are 325 m, 10 m, 434 m, and 30 m, respectively [20]. The situations of both case studies (the Azad and Bakhtiari dams) on the Iran map are illustrated in Figure 1. of Dorud in the Lorestan state, and 70 km northeast of the town of Andimeshk in Khuzestan, Iran. The position of the dam site is at the 48° 45′ 34.87″ to 32° 57′ 23.58″ geographical coordinates of eastern and northern latitudes, respectively (Figure 1). The dam was built upon the Bakhtiari River, aiming to provide adequate water for many purposes, such as drinking, electrical power generation, flood control, and agricultural activities. The dam's body is at an elevation of 840 m. In addition, in the case of this dam, the peak elevation, crown width, crown length, and foundation width are 325 m, 10 m, 434 m, and 30 m, respectively [20]. The situations of both case studies (the Azad and Bakhtiari dams) on the Iran map are illustrated in Figure 1.

**Figure 1.** Situations of both case studies (the Azad and Bakhtiari Dams). **Figure 1.** Situations of both case studies (the Azad and Bakhtiari Dams).

The Azad dam comprises a series of common structures, such as tailraces, higher deposits, different caverns, gurgitation storages, and access tunnels. Geologically, this dam is situated in Iran's famous formation, Sanandaj–Sirjan, with an alternation of schist, sandstone, limestone, and phylite rocks. The bedrock of the dam mainly comprises sandstone with a low degree of metamorphic, phyllite, and schist. Additionally, within the highland areas, lenses of limestone are also observable. From the stratigraphy point of view, rock outliers from the higher Cretaceous period to the present can be observed The Azad dam comprises a series of common structures, such as tailraces, higher deposits, different caverns, gurgitation storages, and access tunnels. Geologically, this dam is situated in Iran's famous formation, Sanandaj–Sirjan, with an alternation of schist, sandstone, limestone, and phylite rocks. The bedrock of the dam mainly comprises sandstone with a low degree of metamorphic, phyllite, and schist. Additionally, within the highland areas, lenses of limestone are also observable. From the stratigraphy point of view, rock outliers from the higher Cretaceous period to the present can be observed within the investigated region. Such rocks consist of four types from the past to the current session:

within the investigated region. Such rocks consist of four types from the past to the current

(1) metamorphic rock related to the Cretaceous period that includes a combination of clay and shale; (2) phyllite formation rocks related to the participation of the Cretaceous Paleocene periods, containing limestone and shale with sand; (3) rocks related to the participation of the Paleocene and Eocene periods, comprising sandstone, shale, limestone lenses, and volcanic rocks; and (4) formations related to the Quaternary period, consisting of shallow terraces and debris. From a tectonic viewpoint, the Sarvabad, Kargineh, and Satileh faults are situated 23, 4.5, and 32 km to the south, east, and northeast of the Azad dam, respectively [21]. The geological conditions and faults of the Azad dam site are shown in Figure 2. session: (1) metamorphic rock related to the Cretaceous period that includes a combination of clay and shale; (2) phyllite formation rocks related to the participation of the Cretaceous Paleocene periods, containing limestone and shale with sand; (3) rocks related to the participation of the Paleocene and Eocene periods, comprising sandstone, shale, limestone lenses, and volcanic rocks; and (4) formations related to the Quaternary period, consisting of shallow terraces and debris. From a tectonic viewpoint, the Sarvabad, Kargineh, and Satileh faults are situated 23, 4.5, and 32 km to the south, east, and northeast of the Azad dam, respectively [21]. The geological conditions and faults of the Azad dam site are shown in Figure 2.

*Sustainability* **2023**, *15*, 4230 4 of 24

**Figure 2.** Geological conditions of the Azad dam site [15].

The Bakhtiari dam's bedrock is made of separate limestone and limestone combined with marl, which incorporates chert nodes. The limestone sections might be synthesized by a combined dolomite substance. From the perspective of geological structure, the dam area is positioned within the pleated Zagros, a portion of the tectono-sediment region of the Zagros. In the lowest northern portion, the area is restricted by pushed Zagros, while in the southwest, it is confined by the Khuzestan plain. With regard to the age of the compressed reservoirs of the area, they date back to between the Triassic and Pliocene eras, and then would have been wrinkled from the Plio-Pleistocene via the latest Alpine organic phase. A number of syncline and anticline sets have been created through such tectonism procedures. Primarily, the above arrangements have been identified by perpendicular axial levels related to the lots of pushed faults in the Zagros area. Additionally, key bed-rocks of the investigated area are made of limestone siliceous related to the current famous formation, Sarvak. This formation (Sarvak) belongs to the Bangestan collection of the middle period of the Cretaceous [20]. For a better review, the geological cross-section of the Bakhtiary dam and plant site is shown in Figure 3. by a combined dolomite substance. From the perspective of geological structure, the dam area is positioned within the pleated Zagros, a portion of the tectono-sediment region of the Zagros. In the lowest northern portion, the area is restricted by pushed Zagros, while in the southwest, it is confined by the Khuzestan plain. With regard to the age of the compressed reservoirs of the area, they date back to between the Triassic and Pliocene eras, and then would have been wrinkled from the Plio-Pleistocene via the latest Alpine organic phase. A number of syncline and anticline sets have been created through such tectonism procedures. Primarily, the above arrangements have been identified by perpendicular axial levels related to the lots of pushed faults in the Zagros area. Additionally, key bedrocks of the investigated area are made of limestone siliceous related to the current famous formation, Sarvak. This formation (Sarvak) belongs to the Bangestan collection of the middle period of the Cretaceous [20]. For a better review, the geological cross-section of the Bakhtiary dam and plant site is shown in Figure 3.

The Bakhtiari dam's bedrock is made of separate limestone and limestone combined with marl, which incorporates chert nodes. The limestone sections might be synthesized

*Sustainability* **2023**, *15*, 4230 5 of 24

**Figure 2.** Geological conditions of the Azad dam site [15].

**Figure 3.** Geological cross-section of the Bakhtiary dam and plant site [15]. **Figure 3.** Geological cross-section of the Bakhtiary dam and plant site [15].

#### *Database Description Database Description*

To create inclusive datasets, adequate core samples with NX sizes, i.e., 54 mm in diameter, perpendicular cylindrical shapes, and ratios of height to diameter in the middle To create inclusive datasets, adequate core samples with NX sizes, i.e., 54 mm in diameter, perpendicular cylindrical shapes, and ratios of height to diameter in the middle of

of 2:1 to 2.5:1 were used based on the process recommended by the ISRM [22,23]. The samples were arranged from the two dam sites introduced before. When the core speci-

2:1 to 2.5:1 were used based on the process recommended by the ISRM [22,23]. The samples were arranged from the two dam sites introduced before. When the core specimens were prepared, their various characteristics, such as E value, porosity, density, and durability index, were measured in a laboratory. Furthermore, in the course of the coring operation, each sample's coring depth was recorded for the purpose of evaluating its impact on the rocks' geomechanical properties. The laboratory experiments in this study were carried out totally based on the ISRM and ASTM standard methods [22,23]. In this regard, a total of 50 test series were done successfully, and the outputs were documented in the cases of all variables noted above. As a result, 50 datasets were provided, aiming to construct the ANFIS-FA, ANFIS-DE, ANFIS, and NN models. Then, a sorting approach was adopted to divide the available database into training (constructing) and testing datasets. Roughly 20% of the database was determined as the testing dataset in order to be used later in the process of evaluating the models built in this paper.

Note that the ratio of 80 to 20 for training and testing groups has been widely suggested by many scholars, such as Ye et al. [24], Fang et al. [25], Nguyen et al. [26], and Zhou et al. [27]. Aside from that, we also tested the ratio of 70:30. Nevertheless, the 80:20 ratio had better performance; thus, this ratio was used in this study.

The statistical characteristics of all variables used in this study are shown in Table 1. For a better view, the frequency histogram of all input and output variables are depicted in Figures 4 and 5. For example, in Figure 4, regarding the depth of coring variables, 11, 21, 4, and 4 data were varied in the range of 0–50 m, 50–100 m, 100–150 m, and 150–250 m, respectively. In addition, Figure 6 illustrates the Pearson correlation plots for all variables.


**Table 1.** Modelling variables and the statistical characteristics datasets.

DC: depth of coring, ρ: density, n: porosity, DI: durability, ν: Poisson ratio, Vp: P-wave velocity, E: elastic modulus.

**Figure 4. Figure 4.**Frequency histogram of all Frequency histogram of all variables for the training group. variables for the training group.

**Figure 5. Figure 5.** Frequency histogram of all variables for the testing group. Frequency histogram of all variables for the testing group.

**Figure 6.** Pearson correlation plots for all input variables vs. the output (E). **Figure 6.** Pearson correlation plots for all input variables vs. the output (E).

#### **3. Methodology 3. Methodology**

This section explains how ANFIS combined with the FA and DE algorithms is implemented. Additionally, the modeling process of the NN model is explained in this section. In the aforementioned models, of the total 50 datasets, 40 were used for the training phase and 10 were used for the testing phase. For a better overview, a schematic flowchart of the ANFIS-FA and ANFIS-DE proposed in this study is shown in Figure 7. It is worth mentioning that MATLAB@2018 was used to encode the proposed hybrid models. This section explains how ANFIS combined with the FA and DE algorithms is implemented. Additionally, the modeling process of the NN model is explained in this section. In the aforementioned models, of the total 50 datasets, 40 were used for the training phase and 10 were used for the testing phase. For a better overview, a schematic flowchart of the ANFIS-FA and ANFIS-DE proposed in this study is shown in Figure 7. It is worth mentioning that MATLAB@2018 was used to encode the proposed hybrid models.

**Figure 7.** A schematic flowchart of the ANFIS-FA and ANFIS-DE proposed in this study. **Figure 7.** A schematic flowchart of the ANFIS-FA and ANFIS-DE proposed in this study.

#### *3.1. ANFIS Combined with DE 3.1. ANFIS Combined with DE*

Differential evolution (DE), which was originally proposed by Storn and Price [28], is an effective evolutionary algorithm that works on the basis of a global optimization approach. In general, DE offers three benefits: (1) a simple structure, (2) high-quality solutions achieved, and (3) easy implementation [29]. As a result, it is applied to a variety of conditions. In the present study, DE is used for the aim of minimizing the function of fitness using the amounts of optimized variable. By definition, the function of fitness refers to the root mean square error (RMSE) between the estimated and target datasets. DE, as an innovative algorithm, was implemented in order to adjust the functions of membership amounts of the ANFIS model and, consequently, enhance its overall prediction capability. Figure 4 illustrates the schematic presentation of the hybrid ANFIS-based DE algo-Differential evolution (DE), which was originally proposed by Storn and Price [28], is an effective evolutionary algorithm that works on the basis of a global optimization approach. In general, DE offers three benefits: (1) a simple structure, (2) high-quality solutions achieved, and (3) easy implementation [29]. As a result, it is applied to a variety of conditions. In the present study, DE is used for the aim of minimizing the function of fitness using the amounts of optimized variable. By definition, the function of fitness refers to the root mean square error (RMSE) between the estimated and target datasets. DE, as an innovative algorithm, was implemented in order to adjust the functions of membership amounts of the ANFIS model and, consequently, enhance its overall prediction capability. Figure 4 illustrates the schematic presentation of the hybrid ANFIS-based DE algorithm.

rithm. To model ANFIS-DE, four parameters must be specified, namely the number of the iteration, crossover probability, mutation probability, and population size. For the selection of the optimal mutation probability, various values were examined, which can be seen in Table 2. The table also shows that by setting the mutation probability to 0.3, the optimum performance with the highest rank related to the testing phase (i.e., the maximum *R<sup>2</sup>* values) was attained. To obtain the best crossover probability, different values were examined (see Table 3). The table shows that the maximum *R<sup>2</sup>* values were attained when the crossover probability was fixed at 0.75. Different population sizes were also tested, as can be observed in Table 4. The table clearly demonstrates that when the population size was set to 250, the best result (the maximum *R<sup>2</sup>* values) was achieved. In these To model ANFIS-DE, four parameters must be specified, namely the number of the iteration, crossover probability, mutation probability, and population size. For the selection of the optimal mutation probability, various values were examined, which can be seen in Table 2. The table also shows that by setting the mutation probability to 0.3, the optimum performance with the highest rank related to the testing phase (i.e., the maximum *R <sup>2</sup>* values) was attained. To obtain the best crossover probability, different values were examined (see Table 3). The table shows that the maximum *R <sup>2</sup>* values were attained when the crossover probability was fixed at 0.75. Different population sizes were also tested, as can be observed in Table 4. The table clearly demonstrates that when the population size was set to 250, the best result (the maximum *R <sup>2</sup>* values) was achieved. In these tests, the smallest amount of error was fixed at <sup>1</sup> <sup>×</sup> <sup>10</sup>−<sup>5</sup> , and the peak repetition was set to 500. Accordingly, the crossover probability, mutation probability, and scope of population were fixed to 0.75, 0.3, and 250, respectively. It is worth mentioning that the bolded amounts in Tables 2–4 are related to the best results (highest rank) obtained from the developed models.


**Table 2.** Selection of the most optimum mutation rate value in implementing the ANFIS-DE model.

**Table 3.** Selection of the most optimum crossover value in implementing the ANFIS-DE model.


**Table 4.** Selection of the most optimum population size value in implementing the ANFIS-DE model.


Figure 8 also shows the ANFIS-DE flowchart used in this study. 500 0.951 0.950 7

**Figure 8. Figure 8.**ANFIS ANFIS-DE Flowchart. -DE Flowchart.

#### *3.2. ANFIS Combined with FA*

This section introduces a hybridized model combining ANFIS and FA, called ANFIS-FA, with the objective of optimizing the premise parameters of ANFIS. To initiate the modeling process, there is a need to first determine the input-target variables/parameters. When the assembly of the input-target variables/parameters of the model is determined, it is time to determine the training samples. The reason for this is that a model should be capable of predicting the target parameter for those samples that have no effect on the model training process. As a result, the entire dataset was divided into two categories: training/construction and testing/examination samples. As previously mentioned, from among a total of 50 datasets, 20% (*n* = 10) were chosen in a random way and assigned to the testing group, while the remaining 80% (*n* = 40) were assigned to the training samples. In addition, it was required to initialize the ANFIS and FA parameters prior to ANFIS modeling in order to predict the target variables. With the use of a trial-and-errorbased approach, the optimal amount of membership functions (MFs) was archived as 6. According to the literature [30–32], the *α*, *β*0, *γ*, number of variables, and the population are the FA parameters. To improve the ANFIS performance, it was necessary to select the most appropriate values for the aforementioned FA parameters. Table 2 clearly shows that the number of variables is equal to six. As stated in the literature [31,32], in some cases, the value of 1 is suitable for the *β*<sup>0</sup> parameters. Therefore, in the modelling of ANFIS-FA, the value of *β*<sup>0</sup> was set to 1.

To select the most appropriate values for the *α* and *γ* parameters, various amounts of these parameters were examined, as given in Tables 5 and 6. Considering these tables, the most appropriate values (the highest *R 2* ) for the *α* and *γ* parameters were obtained with *α* = 0.6 and *γ* = 1.5. As a result, the values of 0.6 and 1.5 were used for the *α* and *γ* parameters in ANFIS-FA modelling. By setting the number of iterations to 1000, different values were also tested to select the most appropriate value for the population size, as shown in Table 7. The table shows that the best performance was attained with population = 200. Based on the above descriptions, the values of 6, 1, 0.6, 1.5 and 200 were set as the number of variables, *β*0, *α*, *γ*, and population size, respectively. It is worth mentioning that the bolded amounts in Tables 5–7 are related to the best results (highest rank) obtained from the developed models. In this step, the most appropriate value of the number of iterations needed to be determined. According to the results, after the 15th iteration, no significant change was observed in the ANFIS-FA performance. In other words, after 15 iterations, the performance for different populations was constant. Accordingly, the number of iterations in ANFIS-FA modelling used in this study was set to 15. When the user-based defined parameters in the investigated models (FA and ANFIS) were determined, then the ANFIS training process was begun with the use of the training samples. To this end, the FA algorithm was used to optimize the primary part of the fuzzy If-Then rules, and the least-square method was applied to the optimization of the linear consequent fuzzy rules.

Additionally, the preliminary light strength corresponding to the primary generation was computed, and then each firefly's attractiveness level was measured. With the use of the movement equation, those fireflies that had a lower level of attraction were pushed toward the brighter firefly. Afterward, the light strength and individual firefly's position were updated, and the function of fitness was computed again. All steps involved in ANFIS-FA are displayed in Figure 9.


**Table 5.** Selection of the optimum *α* value in implementing the ANFIS-FA.

**Table 6.** Selection of the optimum *γ* value in implementing the ANFIS-FA.


**Table 7.** Selection of the optimum population value in implementing the ANFIS-FA.


**Figure 9.** Implementing the ANFIS-FA used in this study. **Figure 9.** Implementing the ANFIS-FA used in this study.

#### *3.3. Neural Network (NN)*

Neural Networks (NNs), especially the Multi-Layer Perceptron (MLP), are widely used in prediction models applied to different engineering problems [31,32]. MLP, which is employed in the present study, contains three layers: input, hidden, and output layers. Therefore, as can be seen in Figure 6, the nodes that exist within the input layer correspond to DC, *ρ*, n, DI, ν, and Vp, while those in the output layer correspond to E. Based on the trial-and-error approach, we considered the number of nodes within the hidden layer. The evaluation results showed that the existence of seven nodes within the hidden layer can result in a higher reliability. As can be observed in Figure 10, the hidden layer with seven nodes resulted in the optimal performance of NN (with the maximum *R 2* ). It is worth mentioning that, to select the suitable number of nodes inside the hidden layer, different numbers were tested. As a result, the NN structure in the present research was built on the basis of six nodes within the first/input layer, seven nodes within the second/hidden layer, and one node within the last/output layer. *Sustainability* **2023**, *15*, 4230 16 of 24 numbers were tested. As a result, the NN structure in the present research was built on the basis of six nodes within the first/input layer, seven nodes within the second/hidden layer, and one node within the last/output layer.

**Figure 10.** Use of various numbers of nodes in the second/hidden layer with their *R<sup>2</sup>* values. **Figure 10.** Use of various numbers of nodes in the second/hidden layer with their *R <sup>2</sup>* values.

#### **4. Results and Discussion**

tions:

**4. Results and Discussion** The present study was aimed at examining the effectiveness of the FA and DE algorithms in optimizing ANFIS for the prediction of E. The results obtained from the proposed ANFIS-FA and ANFIS-DE models were compared to those of ANFIS and NN models. Here, the models' prediction capabilities were assessed regarding RMSE, mean of av-The present study was aimed at examining the effectiveness of the FA and DE algorithms in optimizing ANFIS for the prediction of E. The results obtained from the proposed ANFIS-FA and ANFIS-DE models were compared to those of ANFIS and NN models. Here, the models' prediction capabilities were assessed regarding RMSE, mean of average percentage error (MAPE), mean of absolute error (MAE), variance account for (VAF), *A10-index*, and performance index (PI) [33–37], as presented in the following equations:

$$\text{MAE} = \frac{1}{n} \sum\_{i=1}^{n} |A\_i - P\_i| \tag{1}$$

$$\text{RMSE} = \sqrt{\frac{\sum\_{i=1}^{n} \left(A\_i - P\_i\right)^2}{n}} \tag{2}$$

(1)

(2)

(5)

(6)

̅ signify the actual,

$$\text{MAPE} = \left[ \frac{1}{n} \sum\_{i=1}^{n} \frac{|A\_i - P\_i|}{A\_m} \right] \times 100\tag{3}$$

$$\text{VAF} = \left[1 - \frac{var \left(A\_i - P\_i\right)}{var \left(A\_i\right)}\right] \times 100\tag{4}$$

(

10 − =

1 ̅

=

)

 + 1

estimated, and average of actual E values, respectively. Additionally, m10 is the number of data with values of rate actual/predicted values (ranging from 0.9 to 1.1), and R in Equation (6) is the correlation coefficient. Table 8 presents the MAPE (%), MAE, RMSE,

As can be observed in Table 8, the lowest MAPE (%), MAE, RMSE, and PI values were determined for the ANFIS-FA model as 6.254%, 1.1, 1.152, and 0.031, respectively. In addition, the highest VAF(%) and *A*10*-index* values were determined for the ANFIS-FA model as 98.778% and 0.9, respectively. These values were calculated for the ANFIS-DE model as 1.827, 1.662, 9.452%, 0.049, 96.957%, and 0.6, respectively; for the ANFIS model as 2.557, 2.456, 13.965%, 0.070, 92.514%, and 0.4, respectively; and for the NN model as

10 

$$A10 - index = \frac{m10}{n} \tag{5}$$

$$PI = \frac{1}{\overline{A\_i}} \frac{RMSE}{R+1} \tag{6}$$

,

, and

where *n* stands for the number of data (*n* = 50), and

VAF(%), and *A*10*-Index* values attained by the developed models.

where *n* stands for the number of data (*n* = 50), and *A<sup>i</sup>* , *P<sup>i</sup>* , and *A<sup>i</sup>* signify the actual, estimated, and average of actual E values, respectively. Additionally, m10 is the number of data with values of rate actual/predicted values (ranging from 0.9 to 1.1), and R in Equation (6) is the correlation coefficient. Table 8 presents the MAPE (%), MAE, RMSE, VAF(%), and *A10-Index* values attained by the developed models.


**Table 8.** Performance of the models to predict E by using six statistical criteria.

As can be observed in Table 8, the lowest MAPE (%), MAE, RMSE, and PI values were determined for the ANFIS-FA model as 6.254%, 1.1, 1.152, and 0.031, respectively. In addition, the highest VAF (%) and *A10-index* values were determined for the ANFIS-FA model as 98.778% and 0.9, respectively. These values were calculated for the ANFIS-DE model as 1.827, 1.662, 9.452%, 0.049, 96.957%, and 0.6, respectively; for the ANFIS model as 2.557, 2.456, 13.965%, 0.070, 92.514%, and 0.4, respectively; and for the NN model as 2.781, 2.639, 15.006%, 0.076, 90.985%, and 0.4, respectively. According to Table 8, the highest total rank values for both the training and testing groups were obtained by the ANFIS-FA model. It is worth mentioning that the bolded amounts in Table 8 are related to the best results (highest rank) obtained from the ANFIS-FA model. For a better overview, the predicted E values provided by all models in the testing phase are depicted in Figure 11. Additionally, Figure 12 shows the amount of error for each model related to the testing phase. According to these two figures, the prediction of E by the ANFIS-FA model is very accurate and closer to measured E values. In addition, Figures 13 and 14 demonstrate the scatter plots of actual versus estimated E values with the use of all predictive models. The figures show that the ANFIS-FA model obtained a greater value for the coefficients of determination (*R* 2 ). The *R* 2 values of 0.988, 0.970, 0.928, and 0.913 were obtained by the ANFIS-FA, ANFIS-DE, ANFIS, and NN models, respectively. Accordingly, FA was more effective in comparison with DE in regard to the ANFIS improvement. Furthermore, the absolute error of ANFIS-FA, ANFIS-DE, ANFIS, and NN models in predicting E for testing datasets (ten datasets) is depicted in Figure 15. According to this Figure, the orange-coloured line, which was obtained by the ANFIS-FA model, yields the lowest absolute error for all ten datasets. Moreover, the Taylor diagrams for both training and testing groups are shown in Figure 16. The results show that the ANFIS-FA has a stronger potential to predict E than the others. In this study, a sensitivity analysis was also performed. For this work, the effect of removing each input variable on E for the ANFIS-FA was calculated. In this regard, six new models based on the combination of input variables were constructed, as follows:

Model 1: inputs: all variables given in Table 1.

Model 2: inputs: all variables given in Table 1 except the depth of coring.

Model 3: inputs: all variables given in Table 1 except density.

Model 4: inputs: all variables given in Table 1 except porosity.

Model 5: inputs: all variables given in Table 1 except durability.

Model 6: inputs: all variables given in Table 1 except Poisson ratio.

Model 7: inputs: all variables given in Table 1 except P-wave velocity.

*Sustainability* **2023**, *15*, 4230 18 of 24

**Figure 11.** Comparing the values of E predicted by the four models. **Figure 11.** Comparing the values of E predicted by the four models. **Figure 11.** Comparing the values of E predicted by the four models.

**Figure 12.** Amount of error for each model, as related to the testing phase. **Figure 12.** Amount of error for each model, as related to the testing phase.

**Figure 13.** Comparison of the actual E value with those predicted, obtained by (**a**) NN, (**b**) ANFIS, (**c**) ANFIS-DE, and (**d**) ANFIS-FA for the training group. **Figure 13.** Comparison of the actual E value with those predicted, obtained by (**a**) NN, (**b**) ANFIS, (**c**) ANFIS-DE, and (**d**) ANFIS-FA for the training group. **Figure 13.** Comparison of the actual E value with those predicted, obtained by (**a**) NN, (**b**) ANFIS, (**c**) ANFIS-DE, and (**d**) ANFIS-FA for the training group.

**Figure 14.** Comparison of the actual E value with those predicted, obtained by (**a**) NN, (**b**) ANFIS, (**c**) ANFIS-DE, and (**d**) ANFIS-FA for the testing group. **Figure 14.** Comparison of the actual E value with those predicted, obtained by (**a**) NN, (**b**) ANFIS, (**c**) ANFIS-DE, and (**d**) ANFIS-FA for the testing group. **Figure 14.** Comparison of the actual E value with those predicted, obtained by (**a**) NN, (**b**) ANFIS, (**c**) ANFIS-DE, and (**d**) ANFIS-FA for the testing group.

**Model**

**Model**

**Train; Rank**

**Train; Rank**

*Sustainability* **2023**, *15*, 4230 20 of 24

**Figure 15.** Absolute error values obtained by the models using the testing phase. **Figure 15.** Absolute error values obtained by the models using the testing phase. **Figure 15.** Absolute error values obtained by the models using the testing phase.

**Figure 16.** Taylor diagrams obtained by the predictive models. **Figure 16.** Taylor diagrams obtained by the predictive models. **Figure 16.** Taylor diagrams obtained by the predictive models.

**Table 9.** Performance of all seven ANFIS-FA models. **Statistical Criteria Total Rank RMSE MAE MAPE (%) VAF (%) Test; Rank Train; Rank Test; Rank Train; Rank Test; Rank Train; Rank Test; Rank Train Test Table 9.** Performance of all seven ANFIS-FA models. **Statistical Criteria Total Rank RMSE MAE MAPE (%) VAF (%) Test; Rank Train; Rank Test; Rank Train; Rank Test; Rank Train; Rank Test; Rank Train Test** Model 1 **0.909; 7 1.152; 7 0.865; 7 1.100; 7 3.899; 7 6.254; 7 98.962; 7 98.778; 7 28 28** The results of the above models are presented in Table 9, which shows that models 1 and 7 had the highest total rank (i.e., the best performance) and lowest total rank (i.e., the worst performance), respectively (Figure 17). Note that, the results of model 1 is bolded in Table 9. The results of presented in Table 9 indicated that once the P-wave velocity was removed from the modeling, the worst performance was obtained; thus, P-wave velocity can be determined as the most effective variable in the modeling.

Model 1 **0.909; 7 1.152; 7 0.865; 7 1.100; 7 3.899; 7 6.254; 7 98.962; 7 98.778; 7 28 28** Model 2 1.776; 3 2.306; 3 1.683; 4 2.204; 3 7.583; 4 12.532; 3 95.979; 3 95.848; 3 14 12

Model 2 1.776; 3 2.306; 3 1.683; 4 2.204; 3 7.583; 4 12.532; 3 95.979; 3 95.848; 3 14 12


**Table 9.** Performance of all seven ANFIS-FA models.

**Figure 17.** Total ranks of all seven constructed ANFIS-FA models to evaluate the sensitivity analysis.

#### **5. Conclusions**

The elastic modulus (E) is considered one of the most significant factors in the primary and ultimate plans of projects related to the geo-engineering field. As a result, it is highly necessary to predict E with a high accuracy level. This paper examined the use of two hybrid evolutionary models, namely ANFIS-FA and ANFIS-DE, to predict E. Additionally, the traditional ANFIS and NN models were developed for comparison aims. In total, 50 datasets were collected during the drilling process in the Azad and Bakhtiari underconstruction dams in Iran. Out of the 50 datasets, 40 were used to construct the models, and the remaining datasets were used to test them. The input parameters considered in the construction of the models were porosity, density, depth of coring, Poisson's ratio, compressional/primary wave velocity, and durability, which were assigned as the input variables, whereas E was the output/target variable. Finally, some statistical indices were designed in order to demonstrate the capacity of the models in the prediction of E. According to the findings, the following results and remarks can be briefly listed:


these two algorithms can be effectively used to address other predicting problems in rock engineering fields.


**Author Contributions:** Conceptualization, M.R. and M.H.; methodology, Y.W. and M.H.; validation, Y.W., M.R., and M.H.; investigation, M.H. and M.R.; data curation, M.R.; writing—original draft preparation, Y.W., M.R., M.H., and R.A.A.; writing—review and editing, Y.W., M.R., M.H., and R.A.A.; supervision, M.R. and M.H. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Project of Tackling Key Problems of Science and Technology in Henan Province (222102320164).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data used in this study may be available on request from the corresponding author.

**Conflicts of Interest:** The authors have no conflicts of interest to declare that are relevant to the content of this article.

### **Nomenclature**


### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

**Zhvan Baqi Qader 1,\*, Zuheir Karabash <sup>2</sup> and Ali Firat Cabalar <sup>1</sup>**


**Abstract:** The aim of this study is to analyze and model the geotechnical characteristics of soils in Erbil city using Geographic Information Systems (GIS) and Artificial Neural Networks (ANNs). The study used GIS to analyze the geotechnical properties of soils by collecting data from 102 boreholes in three different depth levels (1.5 m–3.5 m, 3.5 m–6.5 m and 6.5 m–9.5 m) to visualize and analyze soil characteristics such as fines content, moisture content, soil plasticity, shear strength parameters, compressibility, Standard penetration test (SPT), and bearing capacity. The paper also establishes the prediction of SPT-N value and bearing capacity based on geotechnical properties of soils using ANN methods and made correlations between SPT values and shear strength parameters with the bearing capacity of the soil. The results analyzed via GIS indicated that the soil classification was silty clay with a small amount of sandy gravel (CL) in most of the study area. According to the SPT–N values, most of the soils in Erbil City ranged between 33 and 50; a higher SPT value generally indicates denser and stronger soil. The value of the shear strength parameter for the maximum friction angle of the soil layers was found to be 36◦ , and the predominant cohesion was approximately 100 kPa. The compression index of soils ranged between 0.11 to 0.31. The results showed that the ANN models were able to accurately predict the geotechnical parameters of the soil types in the study area. In addition, the use of GIS and ANN techniques allowed for a comprehensive analysis of the geotechnical characteristics of the soils in Erbil, providing valuable information for future construction and development projects.

**Keywords:** Erbil; geotechnical characterization; GIS; ANN

### **1. Introduction**

One of the most important steps before constructing infrastructure is the geotechnical site investigation. It provides information on the site suitability for design criteria and possible construction problems such as time and resources. There are many methods for site investigations and in-situ tests, including pressure meter test, dilatometer test, SPT, cone penetration test (CPT), and plate load test [1]. For the construction of multistory buildings, highways, bridges, and industrial facilities, a soil survey is required to determine the type of soil, consistency, index properties, relative density, groundwater level, shear strength parameters, (SPT) value and bearing capacity [2]. It is necessary to know the bearing capacity of the soil layers for design, the choice of the foundation type, and the foundation depth for any superstructure [2–4]. A geotechnical investigation provides valuable information on the physical and mechanical properties of the soil and rock at a site, which is necessary for safe and durable engineering structures. The information collected from a geotechnical investigation is used to make informed decisions about the design and construction of the foundation, and to identify potential hazards, such as soil liquefaction or instability of slopes, that could compromise the safety of the structure. Therefore, to obtain the geotechnical parameters required for the calculation of the soil bearing capacity

**Citation:** Qader, Z.B.; Karabash, Z.; Cabalar, A.F. Analyzing Geotechnical Characteristics of Soils in Erbil via GIS and ANNs. *Sustainability* **2023**, *15*, 4030. https://doi.org/10.3390/ su15054030

Academic Editors: Jian Zhou, Mahdi Hasanipanah and Danial Jahed Armaghani

Received: 12 January 2023 Revised: 14 February 2023 Accepted: 18 February 2023 Published: 22 February 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

and settlement, in situ testing is required in addition to the collecting of disturbed and undisturbed specimens at different depths. Thus, several geotechnical experiments are conducted on these specimens to determine various parameters that are typically used to design the foundations [5]. Researchers have studied the reliability of SPT to determine the bearing capacity of soil [6–8]. Currently, the SPT test is used to evaluate the bearing capacity to design foundations [9,10].

GISs and ANNs could be used together in analyzing the geotechnical characteristics, and to predict the shear strength, settlement, and bearing capacity of the soil from the index properties of soils. The GISs provide an analytical function that is time–consuming for developing model entry data at different spatial scales [10,11]. A GIS is an organization of data that people interact with to integrate, analyze, and visualize data, to identify relationships, patterns and trends, and to resolve complicated issues; GIS has been used by many researchers to analyze various data [12–15]. ArcGIS was designed to capture data, store, update, process and present data, and to conduct analyses [16]. GIS can help to recognize possible challenges to the completion of the project early in the design process, which can help to avoid time losses for a construction project. Therefore, a GIS is a modular instrument that can be used to support geotechnical site assessments. It has been used to guide land preparation and to integrate field data with existing data [17].

ANN is one of the prevalent algorithms among researchers nowadays, specifically in geotechnical issues. ANN holds three significant advantages: first, the counting speed is high. Second, it has a strong fault-tolerant capability. Third, it is proficient in dealing with problems with complicated problem-solving rules [18]. The technique of utilizing ANN could be a suggestion for predictions, especially in cases where theoretical modeling does not give foreseen outcomes [19]. ANN aims to model the behavior of the nervous system in the human brain. ANN is an adequate solution for solving complex and nonlinear data modeling. Ref. [20] presents the estimation of standard penetration test values on cohesive soil using an artificial neural network without data normalization. Some previous studies investigated the assessment of geotechnical properties and determination of shear strength parameters by unitized ANN [21–23]. In the geotechnical domain, the development situations generally have multiple variables, making them challenging to model employing conventional mathematics [24].

In this study, the test results of 102 boreholes were gathered, categorized, and analyzed and modeled using ArcGIS10.7 software. These data were used to construct models using ANNs to predict the SPT values and bearing capacity of soils. The data represent the area of Erbil City in Iraq, covering a depth of 9.5 m below the ground surface. Data were collected from the Andrea Engineering Test Laboratory and the construction laboratory in Erbil. Therefore, for a geotechnical engineer, this information can be used to classify areas into zones according to GIS results. The purpose of this study is to analyze and model the geotechnical parameters such as the fines content, moisture content, soil plasticity, shear strength parameters, compressibility parameters, SPT–N values and soil bearing capacity via GIS to create a group of maps in different layers. In addition, the prediction of SPT-N value and bearing capacity based on geotechnical parameters were modeled by ANN methods. Correlations between shear strength parameters, SPT values and bearing capacities of soils were made by Minitab 17 programming.

By combining the capabilities of GIS and ANNs, researchers can develop models to analyze geotechnical issues at different spatial scales, producing results that are more accurate and efficient compared to manual methods In summary, GIS and ANNs complement each other by providing an analytical function that is efficient for developing model data entry at different spatial scales for geotechnical issues, allowing for more accurate predictions and a better understanding of the relationships between soil properties and other factors. To the best of our knowledge, this is the first study to cover mapping and modelling all the geotechnical characteristics of soil in Erbil city.

#### **2. Study Area** belt of northern Iraq in the structural trough with a NW–SE axial trend, and within the

*Sustainability* **2023**, *15*, x FOR PEER REVIEW 3 of 39

Erbil is located in the northwestern region of Iraq. Geologically, it is in the low–folded belt of northern Iraq in the structural trough with a NW–SE axial trend, and within the foothill zone, which is part of the stable shelf tectonic unit of Iraq. Erbil City has an area of approximately 250 km<sup>2</sup> and GPS coordinates of 36◦11027.4<sup>00</sup> N 44◦00033.7<sup>00</sup> [25]. A location map of Erbil City is shown in Figure 1. From a geomorphological perspective, the area is flat with uncommon low–lying hills. In addition, Erbil City is stratigraphically covered by quaternary and Pleistocene deposits, which are dominated by clay, silt, and sand [26]. Erbil City is mainly covered by soils such as gravel and conglomerates with sand, clay, and silt. Conglomerates cover more than 80% of the study area [27]. foothill zone, which is part of the stable shelf tectonic unit of Iraq. Erbil City has an area of approximately 250 km2 and GPS coordinates of 36°11′27.4″ N 44°00′33.7″ [25]. A location map of Erbil City is shown in Figure 1. From a geomorphological perspective, the area is flat with uncommon low–lying hills. In addition, Erbil City is stratigraphically covered by quaternary and Pleistocene deposits, which are dominated by clay, silt, and sand [26]. Erbil City is mainly covered by soils such as gravel and conglomerates with sand, clay, and silt. Conglomerates cover more than 80% of the study area [27].

Erbil is located in the northwestern region of Iraq. Geologically, it is in the low–folded

**Figure 1.** Location of Erbil City, Iraq. **Figure 1.** Location of Erbil City, Iraq.

In the last decade, Erbil city had extensive development in the construction of railways, internal and ring roads. Therefore, collecting data, building a database, producing GIS maps for soil properties, and developing a model for soil behaviors and the bearing capacity of the foundations would be very useful for site engineers to make immediate decisions regarding the selection of project positions. Finally, the geotechnical properties at different depths were collected from the study area (Erbil city), analyzed and tabulated In the last decade, Erbil city had extensive development in the construction of railways, internal and ring roads. Therefore, collecting data, building a database, producing GIS maps for soil properties, and developing a model for soil behaviors and the bearing capacity of the foundations would be very useful for site engineers to make immediate decisions regarding the selection of project positions. Finally, the geotechnical properties at different depths were collected from the study area (Erbil city), analyzed and tabulated for the 102 boreholes. The locations of the boreholes were selected to ensure a uniform distribution throughout the study area. The borehole locations are shown in Figure 2.

for the 102 boreholes. The locations of the boreholes were selected to ensure a uniform distribution throughout the study area. The borehole locations are shown in Figure 2.

*Sustainability* **2023**, *15*, x FOR PEER REVIEW 4 of 39

**Figure 2.** Borehole locations in the study area. **Figure 2.** Borehole locations in the study area. **3. Methodology**  The methodology of the study involved collecting data from the field in Erbil city

#### **3. Methodology 3. Methodology** center and analyzing it using both GIS and ANNs. The data was processed and analyzed

The methodology of the study involved collecting data from the field in Erbil city center and analyzing it using both GIS and ANNs. The data was processed and analyzed using these tools to gain insights and conclusions about the study area. The flow chart represents in Figure 3 the methodology of this study as a tool to help readers understand the process used. A flow chart can show the different steps involved in integrating GIS, ANN, and lab analyses, making it easier for readers to follow and comprehend the study's The methodology of the study involved collecting data from the field in Erbil city center and analyzing it using both GIS and ANNs. The data was processed and analyzed using these tools to gain insights and conclusions about the study area. The flow chart represents in Figure 3 the methodology of this study as a tool to help readers understand the process used. A flow chart can show the different steps involved in integrating GIS, ANN, and lab analyses, making it easier for readers to follow and comprehend the study's methodology. using these tools to gain insights and conclusions about the study area. The flow chart represents in Figure 3 the methodology of this study as a tool to help readers understand the process used. A flow chart can show the different steps involved in integrating GIS, ANN, and lab analyses, making it easier for readers to follow and comprehend the study's methodology.

**Figure 3.** Flow chart showing in a simple way the methodology adopted in the present study. **Figure 3. Figure 3.** Flow chart showing in a simple way the me Flow chart showing in a simple way the methodology adopted in the present study. thodology adopted in the present study.

#### *3.1. Data Collection*

In this study, the results of site investigations and numerous series of soil laboratory tests were collected from 102 borehole locations that cover the main part of the region. Soil investigation included drilling boreholes, taking both disturbed and undisturbed soil samples at 1–meter intervals from 1.0 m to 9.5 m depth, and various field tests. Laboratory work included a series of geotechnical tests to determine the soil index properties, sieve analysis, compressibility, settlement, and shear strength.

A sample of the data collected from the laboratory and field tests is presented in Table 1. This table includes some statistical information (e.g., min. values, max. values, average values, and standard deviation of the input and output data). All data of the study could find in Appendix A as a Table A1.


**Table 1.** Inputs and output of the present study.

SD \* = Standard deviation.

#### *3.2. Geographical Information Systems (GISs)*

The results of the soil investigation and field and laboratory tests were employed to create a digital database for the study region. A database of geotechnical properties was used to provide the input values for the mapping software. In this study, the data were analyzed and presented as maps using ArcGIS (10.7) software. Deterministic methods (inverse distance weighting) were used to create maps and models of spatial data, which rely on probability and uncertainty. Deterministic methods use a fixed set of rules or algorithms to create maps and models, unlike geostatistical methods used for analyzing and modeling spatial data [11]. These methods are commonly used to analyze patterns, relationships, and trends in large, complex datasets.

Some of the most commonly used deterministic methods in ArcGIS include:


The resulting digital maps illustrate the soil formation patterns, distribution, and geotechnical properties of soils at different depths. These maps simplify and help designers and site engineers make the right decisions in the construction of projects. The aim of drawing digital maps using the GIS method was to illustrate the bearing capacity of foundations at three different depths. In this study, the bearing capacity was estimated using two methods. The first Meyerhof method (1963) used shear strength parameters for (10 × 10 m). In the second method, standard penetration numbers were used for bearing capacity estimation [28].

#### *3.3. Statistical Analysis*

To make a correlation between the geotechnical properties, the data obtained from the soil investigation of all boreholes in the study area were correlated by MINITAB 17 software. For instance, correlations were made between the SPT values of soil strata with the shear strength parameters of soils and the ultimate bearing capacity. This method is used to find a correlation between the response (Y) and predictor (X) using regression analysis, which is an extensively used method for analyzing multifactor data.

### *3.4. Neural Network Model*

The process for creating an artificial neural network is assumed by using the Matlab application. This study aims to make models by an artificial neural network. The network model developed was formed from data collection of geotechnical properties of soils in the study area. The ANN analysis result aims to predict SPT N-value and bearing capacity using the identical algorithm, the Back-propagation algorithm, and the same activation function. The network architecture was chosen using hidden layers and varying the number of neurons in the hidden layer. The relation number of neurons in the hidden layer is between 15 and 18 according to previous researches [29,30]. The network performance that has the smallest error and the correlation coefficient value that is proximate to 1 is most suitable for data predictions. Root Mean Squared Error (RMSE) is a commonly used evaluation metric in ANN models. RMSE measures the difference between the predicted and actual values, and it is expressed in the same units as the target variable. A low RMSE value indicates that the predictions are close to the actual values, while a high RMSE value displays that the predictions are far from the actual values. In ANN models, RMSE is used to evaluate the performance of the model and determine the quality of the predictions. A lower RMSE value shows a better fit between the predictions and the actual values and a more accurate model. In this research, the R<sup>2</sup> , RMSE, and MAE values of the estimated and actual target parameters are computed in the implementation evaluation of regression models. The R<sup>2</sup> , RMSE, and MAE represented mathematically as Equations (1)–(3):

$$\mathbf{R}^2 = 1 - \frac{\sum\_{\mathbf{I}=1}^{\mathbf{N}} \left( \mathbf{y}\_{\text{mea}} - \mathbf{y}\_{\text{pre}} \right)^2}{\sum\_{\mathbf{I}=1}^{\mathbf{N}} \left( \mathbf{y}\_{\text{mea}} - \mathbf{y}\_{\text{m}} \right)^2} \tag{1}$$

$$\text{RMSE} = \sqrt{\frac{\sum\_{I=1}^{N} \left(\mathbf{y}\_{\text{mean}} - \mathbf{y}\_{\text{pre}}\right)^{2}}{N}} \tag{2}$$

$$\text{MAE} = \frac{1}{\text{N}} \sum\_{\text{I}=1}^{\text{N}} \left| \mathbf{y}\_{\text{mean}} - \mathbf{y}\_{\text{pre}} \right| \tag{3}$$

*Sustainability* **2023**, *15*, x FOR PEER REVIEW 7 of 39

MAE =

1

୍ୀଵ

<sup>N</sup>หy୫ୣୟ <sup>−</sup> y୮୰ୣห

where ymea, ypre, and ym represent the average of existing output, predicted output, and

with R2 immediacy to 1. RMSE and MAE are utilized to assess the model's prediction

where ymea, ypre, and y<sup>m</sup> represent the average of existing output, predicted output, and actual output, respectively. N denotes for all number of data. The degree of fitting is raised with R<sup>2</sup> immediacy to 1. RMSE and MAE are utilized to assess the model's prediction capability. For the RMSE and MAE, the prediction model will be more exact and its accuracy will be higher with a smallish value. capability. For the RMSE and MAE, the prediction model will be more exact and its accuracy will be higher with a smallish value. Figure 4 illustrates the structure of neural network models to predict SPT N value as output with two models: (a) using input as (LL%, PL%, PI%, WC, cohesion, ϕ, Fine con-

(3)

Figure 4 illustrates the structure of neural network models to predict SPT N value as output with two models: (a) using input as (LL%, PL%, PI%, WC, cohesion, φ, Fine content) and (b) using inputs as (LL%, PL%, PI%, WC, φ, Fine content). The structures of neural network models to predict ultimate bearing capacity are presented in Figure 5. The parameters were used as input in two models: (c) (LL%, PL%, PI%, WC, cohesion, φ, Fine content) and (d) (LL%, PL%, PI%, WC, φ, Fine content). tent) and (b) using inputs as (LL%, PL%, PI%, WC, ϕ, Fine content). The structures of neural network models to predict ultimate bearing capacity are presented in Figure 5. The parameters were used as input in two models: (c) (LL%, PL%, PI%, WC, cohesion, ϕ, Fine content) and (d) (LL%, PL%, PI%, WC, ϕ, Fine content).

**Figure 4.** Structure of neural network models (a) and (b) to predict SPT-N value. **Figure 4.** Structure of neural network models (a) and (b) to predict SPT-N value.

**Figure 5.** Structure of neural network models (c) and (d) to predict Q-ultimate. **Figure 5.** Structure of neural network models (c) and (d) to predict Q-ultimate.

#### **4. Results and Discussions**

### *4.1. Modeling of Soil Properties Using GIS Maps*

**4. Results and Discussions**  *4.1. Modeling of Soil Properties Using GIS Maps*  GIS maps of the soil properties of the study area were produced. These maps included the distribution of the soil fines content, natural water content, liquid limit, plastic limit, cohesion, angle of internal friction, compression index, rebound index, SPT–N values, and bearing capacity of foundations at three different levels (1.5 m–3.5 m, 3.5 m–6.5 GIS maps of the soil properties of the study area were produced. These maps included the distribution of the soil fines content, natural water content, liquid limit, plastic limit, cohesion, angle of internal friction, compression index, rebound index, SPT–N values, and bearing capacity of foundations at three different levels (1.5 m–3.5 m, 3.5 m–6.5 m, and 6.5 m–9.5 m). The soil characteristics at different depths were interpolated for the survey area to show the distribution of these properties in a clear way. In general, the soil characteristics in all the maps were divided into six major legends, each of which was represented by a unique color.

### m, and 6.5 m–9.5 m). The soil characteristics at different depths were interpolated for the 4.1.1. Fines Content Model

4.1.1. Fines Content Model

survey area to show the distribution of these properties in a clear way. In general, the soil characteristics in all the maps were divided into six major legends, each of which was represented by a unique color. Fines content in soils is one of the most significant parameters that affect soil behaviors such as shear strength, compressibility, plasticity, and indirectly, the bearing capacity of foundations [31]. The results presented in Figure 6 indicate that the majority of the study area in Erbil city center has a high fines content. According to the figure, the fines content

of foundations [31]. The results presented in Figure 6 indicate that the majority of the study area in Erbil city center has a high fines content. According to the figure, the fines content in most of the area is greater than 81%, meaning that the proportion of soil particles that are smaller than 0.075 mm in diameter is high. This high fine content will likely have significant effects on the soil's engineering properties, such as its shear strength, compressibility, and plasticity. The impact of soil properties on its ability to support loads, resist deformation, and transmit loads must be carefully considered in future construction and development projects in the area. The results from Figure 6 provide crucial information for engineers and planners in Erbil city center, emphasizing the need to consider

in most of the area is greater than 81%, meaning that the proportion of soil particles that are smaller than 0.075 mm in diameter is high. This high fine content will likely have significant effects on the soil's engineering properties, such as its shear strength, compressibility, and plasticity. The impact of soil properties on its ability to support loads, resist deformation, and transmit loads must be carefully considered in future construction and development projects in the area. The results from Figure 6 provide crucial information for engineers and planners in Erbil city center, emphasizing the need to consider the fines content of soil in decision-making for development projects. The compression level of soil becomes crucial when large particles are replaced by fine particles, and the impact of fines content is more pronounced when the soil is near saturation. This has been noted by other researchers in the field [31–33]. The relationship between natural water content and fines content is intertwined and cannot be evaluated separately. Sometimes, improving natural water content weakens the effect of fines content on soil shear strength, due to the sensitivity of fines to changes in natural water content. In dry conditions, fines do not significantly affect soil behavior due to the influence of suction [34,35]. *Sustainability* **2023**, *15*, x FOR PEER REVIEW 9 of 39 the fines content of soil in decision-making for development projects. The compression level of soil becomes crucial when large particles are replaced by fine particles, and the impact of fines content is more pronounced when the soil is near saturation. This has been noted by other researchers in the field [31–33]. The relationship between natural water content and fines content is intertwined and cannot be evaluated separately. Sometimes, improving natural water content weakens the effect of fines content on soil shear strength, due to the sensitivity of fines to changes in natural water content. In dry conditions, fines do not significantly affect soil behavior due to the influence of suction [34,35].

**Figure 6.** Fines content at depths (**a**) 1.5–3.5 m. (**b**) 3.5–6.5 m. (**c**) 6.5–9.5 m. **Figure 6.** Fines content at depths (**a**) 1.5–3.5 m. (**b**) 3.5–6.5 m. (**c**) 6.5–9.5 m.

Atterberg limits can be used to characterize soil behavior and classification, including

liquid limit and plastic limit [36]. Soil water content affects its consistency. A high water content in clay makes it behave like a liquid due to the reduction of attraction between

4.1.2. Atterberg Limits

#### 4.1.2. Atterberg Limits

Atterberg limits can be used to characterize soil behavior and classification, including the swelling potential of expansive soils, the consistency and plasticity of the soil. Two essential index properties can be obtained from the values of the natural water content: liquid limit and plastic limit [36]. Soil water content affects its consistency. A high water content in clay makes it behave like a liquid due to the reduction of attraction between clay particles caused by the excess water between the particles [37]. The liquid limit was used to determine the consistency of the fine-grained soil. This measure of soil consistency is useful for estimating soil consolidation properties and calculating the acceptable bearing capacity and settlement of foundations [38].

The variations in the soil liquid limit throughout the study area at the three different depths are illustrated in Figure 7. The results of the analysis of liquid limit values in the study area show that a significant portion of the soils have a liquid limit range between 40% and 52%. This range indicates the presence of low plastic clay, which is a type of soil that has low resistance to deformation when subjected to stress and is prone to collapsing. The presence of low plastic clay in high percentages in the study area is an important factor to consider in construction and development projects, as it may impact the stability and integrity of structures built on these soils. Additionally, the observation of soils with high liquid limit values (greater than 53%) in relatively central regions is also important. These central zones are considered critical points in the study area, as soils with high liquid limit values are prone to deformation and instability. These critical points need to be carefully evaluated and addressed in any future construction and development projects in the area to ensure the stability and integrity of the structures built on these soils. Numerous researchers have investigated the relationship between the liquid limit of soils and swelling potential. Some types of clay minerals with a high cation exchange capacity (CEC) suffer from expansion and an increase in the volume of the available water [38–41].

Plasticity is one of the most important features of clay, and the crystallinity of clay minerals is the primary source of this plasticity [39–41]. Soil is plastic when the water content is below the liquid and plastic limits. The plastic range, which is the difference between the two limit values, is called the plasticity index [42]. The plastic limit provides geotechnical engineers with indirect information about the activity, toughness index, and optimum moisture content of soils. Figure 8 shows the plastic limit variation in the study area at different depths. The analysis of the figure reveals that half of the study area has a plastic limit of soils that ranges between 19% and 24%, while the other half has a plastic limit ranging between 25% and 30%. This indicates that the soils in the study area have different levels of plasticity, which is a measure of the soil's ability to change shape without breaking under stress. The presence of soils with low plastic limits (Figure 8a) in small zones in the southeast direction is also an important observation. Soils with low plastic limits are less plastic and more brittle, making them more prone to cracking and failure under stress. These areas need to be carefully assessed and addressed in future construction projects in the area to ensure the stability and reliability of the structures built on these soils.

### 4.1.3. Natural Water Content Model

The natural water content can be considered a parameter that profoundly affects the geotechnical properties [43,44]. Figure 9 shows the variations in the natural water content in Erbil City. The figure presents the water content at depths of 1.5 m–3.5 m, 3.5 m–6.5 m, and 6.5 m–9.5 m. The results show that the natural water content in the study area is mostly found to range between 16% and 20% at the three levels of investigation depths. The natural water content of soils is an important factor to consider in construction and development projects, as it affects the soil's stability, bearing capacity, and compressibility. Soils with high natural water content are more susceptible to instability, while soils with low natural water content are more prone to drying and cracking.

**Figure 7.** Liquid limit of soils at depths (**a**) 1.5–3.5 m. (**b**) 3.5–6.5 m. (**c**) 6.5–9.5 m. **Figure 7.** Liquid limit of soils at depths (**a**) 1.5–3.5 m. (**b**) 3.5–6.5 m. (**c**) 6.5–9.5 m.

This is consistent with the findings of many researchers working in this region [43,44]. The water table, during the time of exploration, was in very high depths below the natural ground level (NGL), meaning the groundwater was relatively shallow. The water table fluctuates seasonally, with an increase during spring. The soil above the water table affects its strength and compressibility, as more moisture results in decreased strength and increased compressibility. Saturated soil below the water table creates settling issues as the consolidation process reduces the natural water content under stable load.

**Figure 8.** Plastic limit of soils at depths (**a**) 1.5–3.5 m. (**b**) 3.5–6.5 m. (**c**) 6.5–9.5 m. **Figure 8.** Plastic limit of soils at depths (**a**) 1.5–3.5 m. (**b**) 3.5–6.5 m. (**c**) 6.5–9.5 m.

4.1.4. Shear Strength Parameters

4.1.3. Natural Water Content Model The natural water content can be considered a parameter that profoundly affects the geotechnical properties [43,44]. Figure 9 shows the variations in the natural water content in Erbil City. The figure presents the water content at depths of 1.5 m–3.5 m, 3.5 m–6.5 m, and 6.5 m–9.5 m. The results show that the natural water content in the study area is mostly found to range between 16% and 20% at the three levels of investigation depths. The natural water content of soils is an important factor to consider in construction and development projects, as it affects the soil's stability, bearing capacity, and compressibility. Soils with high natural water content are more susceptible to instability, while soils with low natural water content are more prone to drying and cracking. This is consistent with the findings of many researchers working in this region [43,44]. The water table, during the time of exploration, was in very high depths below The shear strength of soils is a crucial parameter in many foundation engineering designs. It refers to the ability of soil to resist forces that cause slipping or sliding along a plane within the soil. The shear strength of soils is important in determining the bearing capacity of shallow and deep foundations, slope stability, tunnels, and lateral pressure on structures [28]. The soil's shear strength comes from its cohesive strength (c) and frictional strength, represented by the angle of internal friction (φ). Cohesive strength is due to the bonding force between soil grains and the binding material, while frictional strength arises from the friction, interlocking, and rolling of soil grains [44–47]. The strength of any soil decreases as the shear strain and expansion or contraction increase or decrease, respectively, with respect to the soil density due to applied loads. Shear strength parameters are widely utilized by different standard equations in the design of foundations, particularly in empirical equations. [48].

the natural ground level (NGL), meaning the groundwater was relatively shallow. The water table fluctuates seasonally, with an increase during spring. The soil above the water table affects its strength and compressibility, as more moisture results in decreased load.

**Figure 9.** Natural water content of soils at depths (**a**) 1.5–3.5 m. (**b**) 3.5–6.5 m. (**c**) 6.5–9.5 m. **Figure 9.** Natural water content of soils at depths (**a**) 1.5–3.5 m. (**b**) 3.5–6.5 m. (**c**) 6.5–9.5 m.

4.1.4. Shear Strength Parameters The shear strength of soils is a crucial parameter in many foundation engineering designs. It refers to the ability of soil to resist forces that cause slipping or sliding along a plane within the soil. The shear strength of soils is important in determining the bearing capacity of shallow and deep foundations, slope stability, tunnels, and lateral pressure on structures [28]. The soil's shear strength comes from its cohesive strength (c) and frictional strength, represented by the angle of internal friction (ϕ). Cohesive strength is due to the bonding force between soil grains and the binding material, while frictional strength arises from the friction, interlocking, and rolling of soil grains [44–47]. The strength of any Figures 10 and 11 show the variations in the cohesion and angle of internal friction, respectively. The combination of these two parameters produces the shear strength of the soil, its variation throughout the study area, and the soil depth. Most of the soils in the study area at shallow depths (1.5 m–3.5 m) had cohesion values between 76 kPa and 100 kPa. However, there were relatively small regions with cohesion values greater than 100 kPa and less than 50 kPa. The area covered with soil had cohesion values greater than 100 kPa, which increased with the depth. The results of the figures are in agreement with the distribution of the fines content in the study area. The fines content of soils in the majority of the study area was found to be greater than 80%, and this high fine content is likely contributing to the high cohesion values in the area.

soil decreases as the shear strain and expansion or contraction increase or decrease, respectively, with respect to the soil density due to applied loads. Shear strength parameters

strength and increased compressibility. Saturated soil below the water table creates settling issues as the consolidation process reduces the natural water content under stable

contributing to the high cohesion values in the area.

are widely utilized by different standard equations in the design of foundations, particu-

Figures 10 and 11 show the variations in the cohesion and angle of internal friction, respectively. The combination of these two parameters produces the shear strength of the soil, its variation throughout the study area, and the soil depth. Most of the soils in the study area at shallow depths (1.5 m–3.5 m) had cohesion values between 76 kPa and 100 kPa. However, there were relatively small regions with cohesion values greater than 100 kPa and less than 50 kPa. The area covered with soil had cohesion values greater than 100 kPa, which increased with the depth. The results of the figures are in agreement with the distribution of the fines content in the study area. The fines content of soils in the majority of the study area was found to be greater than 80%, and this high fine content is likely

larly in empirical equations. [48].

**Figure 10.** Cohesion of soils (kPa) at depths (**a**) 1.5–3.5 m. (**b**) 3.5–6.5 m. (**c**) 6.5–9.5 m. **Figure 10.** Cohesion of soils (kPa) at depths (**a**) 1.5–3.5 m. (**b**) 3.5–6.5 m. (**c**) 6.5–9.5 m.

Soils with high fines content produce higher cohesion values [49]. The angle of internal friction is a parameter of the soil shear strength and is employed in bearing capacity estimation, slope stability analyses, and estimation of soil lateral earth pressures. The soils in the study area were found to have an angle of internal friction between 2◦ and 6◦ , which was found to be similar at all three depths in this investigation. The east-south part of the study area had soils with an angle of internal friction ranging from 7◦ to 12◦ . These results indicate that the angle of friction values found in the study area are consistent with those found in similar studies in the region. The angle of internal friction is an important parameter in determining the shear strength of the soil, and these results suggest that the soils in the study area have moderate to low shear strength values. The results of the figures can be used to identify zones of high and low shear strength, which are important in determining the suitability of the soil for different types of structures. The information can also be used to design and construct structures that are appropriate for the soil conditions in the area [47,50].

**Figure 11.** Angle of internal friction of soils at depths (**a**) 1.5–3.5 m. (**b**) 3.5–6.5 m. (**c**) 6.5–9.5 m. **Figure 11.** Angle of internal friction of soils at depths (**a**) 1.5–3.5 m. (**b**) 3.5–6.5 m. (**c**) 6.5–9.5 m.

4.1.5. Consolidation Parameter Model

in the area [47,50].

Soils with high fines content produce higher cohesion values [49]. The angle of internal friction is a parameter of the soil shear strength and is employed in bearing capacity estimation, slope stability analyses, and estimation of soil lateral earth pressures. The soils in the study area were found to have an angle of internal friction between 2° and 6°, which was found to be similar at all three depths in this investigation. The east-south part of the study area had soils with an angle of internal friction ranging from 7° to 12°. These results indicate that the angle of friction values found in the study area are consistent with those found in similar studies in the region. The angle of internal friction is an important parameter in determining the shear strength of the soil, and these results suggest that the soils in the study area have moderate to low shear strength values. The results of the figures can be used to identify zones of high and low shear strength, which are important in The consolidation parameters (compression index and swelling index) of saturated clayey soil should be checked during the analysis and design of the foundations [51]. In this study, the distribution of the compression index and swelling index throughout the study area and their variation with the soil depth were investigated. Figure 12 shows that the compression index of soils in the study area, Erbil city center, has a range of values between 0.17 and 0.22, with a lower value found at greater depths. The compression index is an important factor in determining soil compression and consolidation. The results suggest that the shallow soil strata in the study area have high consolidation and settlement potential, while the settlement potential is expected to decrease with depth due to the variation in compression indices.

determining the suitability of the soil for different types of structures. The information can also be used to design and construct structures that are appropriate for the soil conditions

4.1.5. Consolidation Parameter Model

variation in compression indices.

The consolidation parameters (compression index and swelling index) of saturated clayey soil should be checked during the analysis and design of the foundations [51]. In this study, the distribution of the compression index and swelling index throughout the study area and their variation with the soil depth were investigated. Figure 12 shows that the compression index of soils in the study area, Erbil city center, has a range of values between 0.17 and 0.22, with a lower value found at greater depths. The compression index is an important factor in determining soil compression and consolidation. The results suggest that the shallow soil strata in the study area have high consolidation and settlement potential, while the settlement potential is expected to decrease with depth due to the

**Figure 12.** Compression index, Cc, of soils at depths (**a**) 1.5–3.5 m. (**b**) 3.5–6.5 m. (**c**) 6.5–9.5 m. **Figure 12.** Compression index, Cc, of soils at depths (**a**) 1.5–3.5 m. (**b**) 3.5–6.5 m. (**c**) 6.5–9.5 m.

The rebound index is an important parameter in geotechnical engineering as it represents the soil's tendency to expand or contract under changes in moisture content. A high rebound index value designates that the soil is more susceptible to expansion and contraction, whereas a low value indicates that the soil is more stable. The results shown The rebound index is an important parameter in geotechnical engineering as it represents the soil's tendency to expand or contract under changes in moisture content. A high rebound index value designates that the soil is more susceptible to expansion and contraction, whereas a low value indicates that the soil is more stable. The results shown in Figure 13 suggest that the soil in the study area has moderate to low rebound index values in the range 0.015–0.078, indicating that the soil is relatively stable and less likely to expand or contract under changes in moisture content. This indicated that there was no soil swelling or shrinkage potential in the study area. The minimum values showed a consistent trend across all levels in the study area (Figure 13a), with the highest swelling index values located north of Erbil at the first level and in the south at the second and third levels. Most of the study area had moderate parameter values of 0.015–0.078. The consolidation parameters were consistent with previous studies [47,50,52].

in Figure 13 suggest that the soil in the study area has moderate to low rebound index values in the range 0.015–0.078, indicating that the soil is relatively stable and less likely to expand or contract under changes in moisture content. This indicated that there was no soil swelling or shrinkage potential in the study area. The minimum values showed a consistent trend across all levels in the study area (Figure 13a), with the highest swelling index values located north of Erbil at the first level and in the south at the second and third levels. Most of the study area had moderate parameter values of 0.015–0.078. The consol-

idation parameters were consistent with previous studies [47,50,52].

**Figure 13.** Rebound index, Cr, of soils at depths (**a**) 1.5–3.5 m. (**b**) 3.5–6.5 m. (**c**) 6.5–9.5 m. **Figure 13.** Rebound index, Cr, of soils at depths (**a**) 1.5–3.5 m. (**b**) 3.5–6.5 m. (**c**) 6.5–9.5 m.

4.1.6. Standard Penetration Test Model

4.1.6. Standard Penetration Test Model The SPT is one of the many standard in situ tests used to identify soil type, stratigraphy, and relative strength measures during site investigations [53,54]. The eastern part of Erbil had higher SPT values, which is attributed to the higher unit weight of soils and the presence of stiffer and stronger soil layers as shown in Figure 14. The SPT values of soils at a depth of 6.5 m–9.5 m was mostly between 35 and 60, and the range increased from The SPT is one of the many standard in situ tests used to identify soil type, stratigraphy, and relative strength measures during site investigations [53,54]. The eastern part of Erbil had higher SPT values, which is attributed to the higher unit weight of soils and the presence of stiffer and stronger soil layers as shown in Figure 14. The SPT values of soils at a depth of 6.5 m–9.5 m was mostly between 35 and 60, and the range increased from the west to the east of the study area. The increase in the SPT values with increasing depth is due to the influence of several factors. Soil type affects the SPT values as different soil types have different characteristics such as density, porosity, and strength, which all impact the SPT results. Unit weight, or the weight per unit volume of soil, is also a factor, as a higher unit weight results in higher SPT values. Confining pressure, the pressure applied to the soil from the surrounding material, also increases with depth, leading to higher SPT values. Overall, the results of SPT tests provide a good indication of the soil strength and its variation throughout the study area. To confirm the relationship between the SPT values and shear strength parameters (cohesion and angle of internal friction), data obtained from the boreholes were correlated. Figures 15 and 16 illustrate the correlation between the

SPT values and shear strength parameters for the soils within the study area. The figures evidently show that there is a good correlation between the SPT values, cohesion, and angle of internal friction. Generally, most values of SPT–N distribution from the center to the west of Erbil City were found to be between 17 and 48, while in the northeast and southeast, the values were higher than 50. As mentioned in various studies [55–57], the results of SPT–N showed that the soil in Erbil city is medium to very dense. tion between the SPT values and shear strength parameters for the soils within the study area. The figures evidently show that there is a good correlation between the SPT values, cohesion, and angle of internal friction. Generally, most values of SPT–N distribution from the center to the west of Erbil City were found to be between 17 and 48, while in the northeast and southeast, the values were higher than 50. As mentioned in various studies [55– 57], the results of SPT–N showed that the soil in Erbil city is medium to very dense.

the west to the east of the study area. The increase in the SPT values with increasing depth is due to the influence of several factors. Soil type affects the SPT values as different soil types have different characteristics such as density, porosity, and strength, which all impact the SPT results. Unit weight, or the weight per unit volume of soil, is also a factor, as a higher unit weight results in higher SPT values. Confining pressure, the pressure applied to the soil from the surrounding material, also increases with depth, leading to higher SPT values. Overall, the results of SPT tests provide a good indication of the soil strength and its variation throughout the study area. To confirm the relationship between the SPT values and shear strength parameters (cohesion and angle of internal friction), data obtained from the boreholes were correlated. Figures 15 and 16 illustrate the correla-

*Sustainability* **2023**, *15*, x FOR PEER REVIEW 18 of 39

**Figure 14. Figure 14.** SPT-N values in site at depths ( SPT-N values in site at depths (**aa** ) 1.5–3.5 m. ( ) 1.5–3.5 m. ( **b b** ) 3.5–6.5 m. ( ) 3.5–6.5 m. (**cc**) 6.5–9.5 m. ) 6.5–9.5 m.

### 4.1.7. Bearing Capacity

The design of a footing depends on its soil's bearing capacity. Many methods for estimating soil bearing capacity exist, relying on factors such as soil shear strength, footing type, and SPT value [58]. In this study, the soil strata's ultimate bearing capacity was estimated using two methods, one based on shear strength parameters using Meyerhof's equation, and the other based on SPT-N values. The variation in the soil strata's ultimate bearing capacity based on Meyerhof's equation is shown in Figure 17.

*Sustainability* **2023**, *15*, x FOR PEER REVIEW 19 of 39

**Figure 15.** Correlation of SPT-N values with cohesion of soils in the study area. **Figure 15.** Correlation of SPT-N values with cohesion of soils in the study area.

2.5 3.0 3.5 4.0 4.5 5.0 5.5 **ϕ° Figure 16.** Correlation of SPT-N values with angle of internal friction (ϕ) of soils in the study. **Figure 16.** Correlation of SPT-N values with angle of internal friction (φ) of soils in the study.

The design of a footing depends on its soil's bearing capacity. Many methods for

The design of a footing depends on its soil's bearing capacity. Many methods for

Figure 18 shows the ultimate bearing capacity of the soil strata in the study area at

Figure 18 shows the ultimate bearing capacity of the soil strata in the study area at

three different depths that were estimated from the SPT values. Changes in the ultimate bearing capacity were observed throughout the study area. The ultimate bearing capacity of the soil strata in the study area is a significant factor in the design of shallow and deep foundations. The results from Figure 18 indicate that the majority of the study area at

type, and SPT value [58]. In this study, the soil strata's ultimate bearing capacity was estimated using two methods, one based on shear strength parameters using Meyerhof's equation, and the other based on SPT-N values. The variation in the soil strata's ultimate

equation, and the other based on SPT-N values. The variation in the soil strata's ultimate

three different depths that were estimated from the SPT values. Changes in the ultimate bearing capacity were observed throughout the study area. The ultimate bearing capacity of the soil strata in the study area is a significant factor in the design of shallow and deep foundations. The results from Figure 18 indicate that the majority of the study area at

bearing capacity based on Meyerhof's equation is shown in Figure 17.

bearing capacity based on Meyerhof's equation is shown in Figure 17.

4.1.7. Bearing Capacity

4.1.7. Bearing Capacity

**Figure 16.** Correlation of SPT-N values with angle of internal friction (ϕ) of soils in the study.

estimating soil bearing capacity exist, relying on factors such as soil shear strength, footing type, and SPT value [58]. In this study, the soil strata's ultimate bearing capacity was es-

sign and load-bearing capacity of foundations in the study area.

**Figure 17.** Ultimate bearing capacity according to Meyerhof's equation at depths (**a**) 1.5–3.5 m. (**b**) 3.5–6.5 m. (**c**) 6.5–9.5 m. **Figure 17.** Ultimate bearing capacity according to Meyerhof's equation at depths (**a**) 1.5–3.5 m. (**b**) 3.5–6.5 m. (**c**) 6.5–9.5 m.

shallow depths had an ultimate bearing capacity be-tween 170 kPa and 940 kPa. However, there were some small areas with lower ultimate bearing capacities. The ultimate bearing capacity increased with increasing depth, which can be attributed to the increase in soil confining pressure and soil unit weight. These findings are useful in determining the de-

Figure 19 presents the correlation between the ultimate bearing capacity from Meyerhof's equation and the SPT values of the soils within the study area. This figure indicates a good correlation between the SPT values and the ultimate bearing capacity of the foundations. The correlation between the ultimate bearing capacity from the SPT values of the soil strata within the study area and soil cohesion (shear strength parameter). The ultimate bearing capacity from Meyerhof's equation and the ultimate bearing capacity from the SPT values of the soil strata within the study area were correlated, and the results are presented in Figure 20. The estimation of ultimate bearing capacity from the Standard Figure 18 shows the ultimate bearing capacity of the soil strata in the study area at three different depths that were estimated from the SPT values. Changes in the ultimate bearing capacity were observed throughout the study area. The ultimate bearing capacity of the soil strata in the study area is a significant factor in the design of shallow and deep foundations. The results from Figure 18 indicate that the majority of the study area at shallow depths had an ultimate bearing capacity be-tween 170 kPa and 940 kPa. However, there were some small areas with lower ultimate bearing capacities. The ultimate bearing capacity increased with increasing depth, which can be attributed to the increase in soil confining pressure and soil unit weight. These findings are useful in determining the design and load-bearing capacity of foundations in the study area.

tigraphy, and loading conditions [59].

**Figure 18.** Ultimate bearing capacity from SPT-N value at depths (**a**) 1.5–3.5 m. (**b**) 3.5–6.5 m. (**c**) 6.5– 9.5 m. **Figure 18.** Ultimate bearing capacity from SPT-N value at depths (**a**) 1.5–3.5 m. (**b**) 3.5–6.5 m. (**c**) 6.5–9.5 m.

Penetration Test (SPT) values is a relatively simple method compared to other methods that require more experimental tests and complex equations established on the shear strength parameters of the soil. The estimation of ultimate bearing capacity from SPT values is based on empirical correlations and has been widely used in the field of geotechnical engineering. The advantage of this method is its simplicity; however, it may not accurately reflect the actual soil conditions, which can be affected by factors such as soil type, stra-

Figure 19 presents the correlation between the ultimate bearing capacity from Meyerhof's equation and the SPT values of the soils within the study area. This figure indicates a good correlation between the SPT values and the ultimate bearing capacity of the foundations. The correlation between the ultimate bearing capacity from the SPT values of the soil strata within the study area and soil cohesion (shear strength parameter). The ultimate bearing capacity from Meyerhof's equation and the ultimate bearing capacity from the SPT values of the soil strata within the study area were correlated, and the results are presented in Figure 20. The estimation of ultimate bearing capacity from the Standard Penetration Test (SPT) values is a relatively simple method compared to other methods that require more experimental tests and complex equations established on the shear strength parameters of the soil. The estimation of ultimate bearing capacity from SPT values is based on empirical correlations and has been widely used in the field of geotechnical engineering. The advantage of this method is its simplicity; however, it may not accurately reflect the actual soil conditions, which can be affected by factors such as soil type, stratigraphy, and loading conditions [59].

area.

*Sustainability* **2023**, *15*, x FOR PEER REVIEW 22 of 39

**Figure 19.** Correlation of ultimate bearing capacity (Meyerhof) with SPT-N values of soils in study **Figure 19.** Correlation of ultimate bearing capacity (Meyerhof) with SPT-N values of soils in study area.

**Figure 20.** Correlation of ultimate bearing capacity (SPT-N) with soil cohesion of soils in study area. **Figure 20.** Correlation of ultimate bearing capacity (SPT-N) with soil cohesion of soils in study area.

#### **Figure 20.** Correlation of ultimate bearing capacity (SPT-N) with soil cohesion of soils in study area. *4.2. Artificial Neural Network Models*

*4.2. Artificial Neural Network Models*  ANN models are commonly used for regression and classification tasks, including *4.2. Artificial Neural Network Models*  ANN models are commonly used for regression and classification tasks, including prediction problems. ANN models consist of interconnected nodes, or artificial neurons, ANN models are commonly used for regression and classification tasks, including prediction problems. ANN models consist of interconnected nodes, or artificial neurons, that process and transmit information. The nodes are organized into layers, and the

ious fields due to their ability to handle complex relationships and make accurate predictions. The relationship between the predicted original SPT values and bearing capacity in the training data and test data was demonstrated in four models that have an acceptable

tions. The relationship between the predicted original SPT values and bearing capacity in the training data and test data was demonstrated in four models that have an acceptable

prediction problems. ANN models consist of interconnected nodes, or artificial neurons,

that process and transmit information. The nodes are organized into layers, and the con-

connections between them can be adjusted during the training process to minimize the error between the predicted and actual outcomes. ANN models have been widely used in various fields due to their ability to handle complex relationships and make accurate predictions. The relationship between the predicted original SPT values and bearing capacity in the training data and test data was demonstrated in four models that have an acceptable R<sup>2</sup> and the smallest error. After input and output data are gathered and structured, training and test sets were established. 70% of the data were used for training and 15% used for testing and 15% for validation of the total data of boreholes. Predictions for SPT-N value from two models were developed.

#### 4.2.1. Validation of Interpolations Based on Semivariograms

The transformed data's spatial autocorrelation is modeled using semivariograms/ covariance modeling for SPT-N values at depth 1.5–3.5 m, depicting the similarity decrease between data points as their distance increases. Binned values (red dots) are generated by grouping semivariograms/covariance points using square cells, while average points (blue crosses) are generated by binning empirical points in angular sectors. Binned points show local variation, while average values show smooth variation. A stable type model (dark blue line) is fitted to the empirical variogram for the measured data points.

At h = 0, the semivariogram should be 0. However, at an infinitesimally small separation distance, the semivariogram often exhibits a nugget effect, which is some value greater than zero. In this case, the nugget effect exists which is zero. The range is the distance at which the model levels out (5750.634). Locations closer than the range are spatially auto correlated, while farther locations are not. The partial sill is the sill minus the nugget (462.1716). The lag size is the distance class size (674.7548) with 12 lags. Semivariogram values are shown in Figure 21 with higher values in orange/red and lower values in blue/green.

Kriging estimates unknown spatial values. The search neighborhood step involves selecting nearby points with significant influence on the prediction location, determined by spatial auto-correlation. The method eliminates irrelevant points and weights nearby points using a search neighborhood of adjacent points, radius, and number of sectors to estimate values at the unknown location. Accurate neighborhood identification and selection of nearby points are crucial for successful kriging. As shown in Figure 22, five neighboring points are selected and a circle with four sectors is selected. The points highlighted in the data give an indicator of the weights associated with each point, and these weights are used to estimate the value at the unknown location, which is at the center of the crosshair.

#### 4.2.2. Prediction for SPT-N Value

In Figures 23 and 24 the results of the ANN model (a) show a good agreement between the predicted and measured SPT-N values, with an R<sup>2</sup> of 0.92 for the training data and 0.81 for the testing data. The model uses Atterberg limit values, water content, cohesion, and internal friction as input variables and predicts the SPT value as its output. On the other hand, to predict SPT N value as output with (LL%, PL%, PI%, WC, φ, Fine content) as input, the result showed R<sup>2</sup> values of 0.90 and 0.8 for training and testing, respectively. As mentioned in previous studies [60], the predictions of the SPT values in model (a) and model (b) were conducted with more superficial R<sup>2</sup> values that give a significant agreement for using ANN modeling in geotechnical engineering that helps the engineering to be utilized in the design of infrastructures.

are used to estimate the value at the unknown location, which is at the center of the cross-

R2 and the smallest error. After input and output data are gathered and structured, training and test sets were established. 70% of the data were used for training and 15% used for testing and 15% for validation of the total data of boreholes. Predictions for SPT-N

The transformed data's spatial autocorrelation is modeled using semivariograms/covariance modeling for SPT-N values at depth 1.5–3.5 m, depicting the similarity decrease between data points as their distance increases. Binned values (red dots) are generated by grouping semivariograms/covariance points using square cells, while average points (blue crosses) are generated by binning empirical points in angular sectors. Binned points show local variation, while average values show smooth variation. A stable type model (dark blue line) is fitted to the empirical variogram for the measured data points.

At h = 0, the semivariogram should be 0. However, at an infinitesimally small separation distance, the semivariogram often exhibits a nugget effect, which is some value greater than zero. In this case, the nugget effect exists which is zero. The range is the distance at which the model levels out (5750.634). Locations closer than the range are spatially auto correlated, while farther locations are not. The partial sill is the sill minus the nugget (462.1716). The lag size is the distance class size (674.7548) with 12 lags. Semivariogram values are shown in Figure 21 with higher values in orange/red and lower values in

value from two models were developed.

blue/green.

4.2.1. Validation of Interpolations Based on Semivariograms

**Figure 21.** Semivariograms/Covariance Modeling. **Figure 21.** Semivariograms/Covariance Modeling. hair.

**Figure 22.** Searching Neighborhood Dialog Box of Ordinary Kriging. **Figure 22.** Searching Neighborhood Dialog Box of Ordinary Kriging.

In Figures 23 and 24 the results of the ANN model (a) show a good agreement between the predicted and measured SPT-N values, with an R2 of 0.92 for the training data and 0.81 for the testing data. The model uses Atterberg limit values, water content, cohesion, and internal friction as input variables and predicts the SPT value as its output. On the other hand, to predict SPT N value as output with (LL%, PL%, PI%, WC, ϕ, Fine con-

tively. As mentioned in previous studies [60], the predictions of the SPT values in model (a) and model (b) were conducted with more superficial R2 values that give a significant agreement for using ANN modeling in geotechnical engineering that helps the engineer-

4.2.2. Prediction for SPT-N Value

ing to be utilized in the design of infrastructures.

*Sustainability* **2023**, *15*, x FOR PEER REVIEW 25 of 39

**Figure 23.** Model (a): relations between (LL%, PL%, PI%, WC, c, φ, Fine content) and SPT-N.

**Figure 24.** Model (b): relations between (LL%, PL%, PI%, WC, ϕ, Fine content) and SPT-N. **Figure 24.** Model (b): relations between (LL%, PL%, PI%, WC, <sup>φ</sup>, Fine content) and SPT-N.

**Figure 24.** Model (b): relations between (LL%, PL%, PI%, WC, ϕ, Fine content) and SPT-N.

#### 4.2.3. Prediction of Ultimate Bearing Capacity

The prediction of the ultimate bearing capacity of soil using ANNs is a machine learning approach that involves training a neural network model with a dataset of soil properties and corresponding ultimate bearing capacity values. The trained model can then be used to predict the ultimate bearing capacity of new soil samples based on their properties. Advantages of ANN prediction include the handling of non-linear relationships between soil properties and ultimate bearing capacity, the ability to incorporate complex relationships between soil properties and ultimate bearing capacity, and the ability to handle large datasets with many input variables.

To determine the best prediction results for the ultimate bearing capacity, two models were used. The inputs and outputs of the models are listed in Table 2 along with the value of R<sup>2</sup> , which is a measure of the goodness of fit of the model to the data. R<sup>2</sup> ranges from 0 to 1, with a higher value indicating a better fit between the model and the data. The R <sup>2</sup> value provides an indication of how well the models are able to predict the ultimate bearing capacity based on the inputs.

**Table 2.** ANN models for Q-Ultimate prediction.


The results of bearing capacity prediction using ANN modeling show similarities with the results of previous studies [60]. In Figures 25 and 26 the relationship between the predicted ultimate bearing capacity and the original ultimate bearing capacity on the training and testing data are illustrated. *Sustainability* **2023**, *15*, x FOR PEER REVIEW 27 of 39

**Figure 26.** Model (d): relations between (LL%, PL%, PI%, WC, ϕ) and Q ultimate. **Figure 25.** Model (c): relations between (LL%, PL%, PI%, WC, c, <sup>φ</sup>, Fine content) and Q ultimate.

**Figure 25.** Model (c): relations between (LL%, PL%, PI%, WC, c, ϕ, Fine content) and Q ultimate.

**Figure 25.** Model (c): relations between (LL%, PL%, PI%, WC, c, ϕ, Fine content) and Q ultimate.

**Figure 26.** Model (d): relations between (LL%, PL%, PI%, WC, φ) and Q ultimate.

**Figure 26.** Model (d): relations between (LL%, PL%, PI%, WC, ϕ) and Q ultimate. Model (c) has a high R<sup>2</sup> value, indicating a strong correlation between the predicted and actual values, and a low variance in the residuals. The high R<sup>2</sup> value for the training data (0.91) indicates that the model is able to fit the training data well, while the R<sup>2</sup> value of 0.82 for the testing data suggests that the model has good generalization ability and can predict unseen data with a certain degree of accuracy. In Model (d), the decrease in R<sup>2</sup> value is likely due to the absence of the cohesion value as an input. Cohesion is an important factor that affects the bearing capacity of soil, so its absence in the input could result in a decrease in the accuracy of the model. This highlights the importance of considering all relevant factors in the input variables of the model to improve its accuracy and predictability.

Soil cohesion is a measure of the shear strength of soil, which determines its ability to support loads. The results of research in geotechnical engineering have shown that soil cohesion is the most important factor in estimating soil bearing capacity. This is because it determines the resistance of soil to sliding or deformation under load, and is essential for ensuring the stability of structures built on the soil. Therefore, accurate determination of soil cohesion is crucial for safe and effective design of geotechnical structures.

#### 4.2.4. Percentage Error of ANN Models

Figures 27 and 28 demonstrate the error percentage lines for model (a) and (b) outcomes referring to the difference between the predicted values and the actual values. It is a measure of the accuracy of the model's predictions. A lower error percentage indicates a more accurate model, while a higher error percentage means that the model is less accurate. The acceptable error percentage depends on the specific application and the acceptable level of error for that particular field. In some cases, a low error percentage, such as 1–2%, may be acceptable, while in others, a higher error percentage may be acceptable if the cost of a lower error is too high [61]. It is important to note that no model can be 100% accurate and some level of error is always present. The goal is to reduce the error to the lowest possible level while still making predictions that are useful and relevant.

Figures 27 and 28 demonstrate the error percentage lines for model (a) and (b) outcomes referring to the difference between the predicted values and the actual values. It is a measure of the accuracy of the model's predictions. A lower error percentage indicates a more accurate model, while a higher error percentage means that the model is less accurate. The acceptable error percentage depends on the specific application and the acceptable level of error for that particular field. In some cases, a low error percentage, such as 1–2%, may be acceptable, while in others, a higher error percentage may be acceptable if the cost of a lower error is too high [61]. It is important to note that no model can be 100% accurate and some level of error is always present. The goal is to reduce the error to

Figures 27 and 28 demonstrate the error percentage lines for model (a) and (b) outcomes referring to the difference between the predicted values and the actual values. It is a measure of the accuracy of the model's predictions. A lower error percentage indicates a more accurate model, while a higher error percentage means that the model is less accurate. The acceptable error percentage depends on the specific application and the acceptable level of error for that particular field. In some cases, a low error percentage, such as 1–2%, may be acceptable, while in others, a higher error percentage may be acceptable if the cost of a lower error is too high [61]. It is important to note that no model can be 100% accurate and some level of error is always present. The goal is to reduce the error to

the lowest possible level while still making predictions that are useful and relevant.

the lowest possible level while still making predictions that are useful and relevant.

**Figure 27.** Error percentage lines of model (a) and (b) for prediction SPT-N value. **Figure 27.** Error percentage lines of model (a) and (b) for prediction SPT-N value. **Figure 27.** Error percentage lines of model (a) and (b) for prediction SPT-N value.

4.2.4. Percentage Error of ANN Models

4.2.4. Percentage Error of ANN Models

*Sustainability* **2023**, *15*, x FOR PEER REVIEW 28 of 39

**Figure 28.** Error percentage lines of model (c) and (d) for prediction Q ultimate. **Figure 28.** Error percentage lines of model (c) and (d) for prediction Q ultimate.

#### **Figure 28.** Error percentage lines of model (c) and (d) for prediction Q ultimate. 4.2.5. Analysis of Models

4.2.5. Analysis of Models Sensitivity analysis is an important aspect of model evaluation, as it helps to assess the impact of individual model parameters on the final results. The absence of a sensitivity analysis in the study you are referring to could limit the understanding of the model's behavior and the confidence in the results obtained. Analysis of models is particularly important when using ANOVA, as the parameters of these models can have a significant 4.2.5. Analysis of Models Sensitivity analysis is an important aspect of model evaluation, as it helps to assess the impact of individual model parameters on the final results. The absence of a sensitivity analysis in the study you are referring to could limit the understanding of the model's behavior and the confidence in the results obtained. Analysis of models is particularly important when using ANOVA, as the parameters of these models can have a significant Sensitivity analysis is an important aspect of model evaluation, as it helps to assess the impact of individual model parameters on the final results. The absence of a sensitivity analysis in the study you are referring to could limit the understanding of the model's behavior and the confidence in the results obtained. Analysis of models is particularly important when using ANOVA, as the parameters of these models can have a significant impact on the spatial structure of the data. By conducting a sensitivity analysis, the researchers could determine the robustness of the model results. The sensitivity analysis to evaluate the effect of each model parameter based on the semivariograms shown in Table 3.


**Table 3.** Analysis of models and effect each parameter on the output of prediction.

A *p*-value less than 0.05 is commonly used as a threshold for determining statistical significance in hypothesis testing. In this context, a *p*-value less than 0.05 for a particular model parameter means that the results suggest a statistically significant relationship between the parameter and the outcome variable being modeled.

In this study, a *p*-value less than 0.05 for a particular model parameter would indicate that the parameter has a significant impact on the model results. This information can be useful in understanding the underlying relationships in the data and in guiding further analysis or interpretation of the results. However, it is important to consider other factors such as the sample size, the quality of the data, and the overall fit of the model in evaluating the reliability and robustness of the results.

### **5. Conclusions**

This study focuses on developing maps for soil geotechnical properties that are widely utilized by geotechnical engineers in foundation design capacity. The main conclusions are as follows. An artificial neural network (NN) model was established for estimating SPT-N value and ultimate bearing capacity:


structures and in designing and constructing structures that are appropriate for the soil conditions in the area.


**Author Contributions:** Conceptualization, K.Z., Z.B.Q. and A.F.C.; methodology, Z.B.Q.; software, Z.B.Q.; validation, Z.B.Q., A.F.C. and Z.K.; formal analysis, Z.B.Q.; investigation, Z.B.Q.; resources, Z.B.Q.; data curation, Z.B.Q.; writing—original draft preparation, Z.B.Q.; writing—review and editing, A.F.C., Z.B.Q. and Z.K.; visualization, Z.B.Q.; supervision, A.F.C. and Z.K.; project administration, A.F.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data are available from the conforming author upon request.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **Appendix A**

The data of all boreholes are provided demonstrated in Table A1.

**Table A1.** Full data of all boreholes of study area.









### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
