Saturated Hydraulic Conductivity Estimation Using Artiﬁcial Intelligence Techniques: A Case Study for Calcareous Alluvial Soils in a Semi-Arid Region

: The direct estimation of soil hydraulic conductivity (Ks) requires expensive laboratory measurement to present adequately soil properties in an area of interest. Moreover, the estimation process is labor and time-intensive due to the difﬁculties of collecting the soil samples from the ﬁeld. Hence, innovative methods, such as machine learning techniques, can be an alternative to estimate Ks. This might facilitate agricultural water and nutrient management which has an impact on food and water security. In this spirit, the study presents neural-network-based models (artiﬁcial neural network (ANN), deep learning (DL)), tree-based (decision tree (DT), and random forest (RF)) to estimate Ks using eight combinations of soil data under calcareous alluvial soils in a semi-arid region. The combinations consisted of soil data such as clay, silt, sand, porosity, effective porosity, ﬁeld capacity, permanent wilting point, bulk density, and organic carbon contents. The results compared with the well-established model showed that all the models had satisfactory results for the estimation of Ks, where ANN7 with soil inputs of sand, silt, clay, permanent wilting point, ﬁeld capacity, and bulk density values showed the best performance with mean absolute error (MAE) of 2.401 mm h − 1 , root means square error (RMSE) of 3.096 mm h − 1 , coefﬁcient of determination (R 2 ) of 0.940, and correlation coefﬁcient (CC) of 0.970. Therefore, the ANN could be suggested among the neural-network-based models. Otherwise, RF could also be used for the estimation of Ks among the tree-based models.


Introduction
The saturated soil hydraulic conductivity (Ks) regulates hydrological activities in soils and its accurate estimation has an important value in hydrological studies, especially for simulating infiltration, soil moisture, runoff, soil erosion, and dynamics of groundwater [1,2]. Therefore, it is essential to know this specific parameter for the management of irrigation events. In addition, the hydrodynamic properties of soils provide useful information regarding the entry and storage of precipitation waters into the soil profile, especially for calcareous alluvial soils with poor structural characteristics in arid and semi-arid environments.
backgrounds. The main motivation comes from the fact that the current soil degradation scenario and the indirect impact it has on food and water security requires inputs from the workforce belonging to a spectrum of technical backgrounds. Therefore, this study is essential to find the best alternate modeling methods for each soil input combination.
The article is organized as follows. Section 2 presents the background of the study area and its parameters, machine learning methods used in this study, input selection and model development process, and performance evaluation criteria. Sections 3 and 4 present the results and discussion, respectively, and Section 5 concludes the paper.

Study Area and Data
The soil data were collected in Çumra-Konya plain (37 • 50 54 N 32 • 43 03 E-37 • 12 17 N 33 • 07 16 E, 1011 m altitude) within an area of about 280,000 ha located 30 km away from the southeast of Konya city. The map is presented for the coordinates of 291 soil samples in the study area ( Figure 1). Based on the Köppen-Geiger climate system [35], the climate of the study area is semi-arid with usually cold and snowy winters and hot and dry summers. According to the long-term climate record, the annual average temperature, relative humidity, and precipitation are 11.4 • C, 62.14%, and 296.8 mm, respectively [36]. The soils of the study area are formed by volcanic rocks, marine, and lacustrine deposits [37]. However, they have some restrictive properties such as deep clay texture formed on alluvial parental materials, low organic matter content, high pH value, low aggregate stability, and shallow soil depth [33,38]. Therefore, the region has insufficient drainage and faces soil erosion [39].
against a mature mathematical model is to eventually provide an alternative and simpler method for Ks estimation which can be used by engineers and scientists from various technical backgrounds. The main motivation comes from the fact that the current soil degradation scenario and the indirect impact it has on food and water security requires inputs from the workforce belonging to a spectrum of technical backgrounds. Therefore, this study is essential to find the best alternate modeling methods for each soil input combination.
The article is organized as follows. Section 2 presents the background of the study area and its parameters, machine learning methods used in this study, input selection and model development process, and performance evaluation criteria. Sections 3 and 4 present the results and discussion, respectively, and Section 5 concludes the paper.

Study Area and Data
The soil data were collected in Çumra-Konya plain (37°50′54″ N 32°43′03″ E-37°12′17″ N 33°07′16″ E, 1011 m altitude) within an area of about 280,000 ha located 30 km away from the southeast of Konya city. The map is presented for the coordinates of 291 soil samples in the study area ( Figure 1). Based on the Köppen-Geiger climate system [35], the climate of the study area is semi-arid with usually cold and snowy winters and hot and dry summers. According to the long-term climate record, the annual average temperature, relative humidity, and precipitation are 11.4 °C, 62.14%, and 296.8 mm, respectively [36]. The soils of the study area are formed by volcanic rocks, marine, and lacustrine deposits [37]. However, they have some restrictive properties such as deep clay texture formed on alluvial parental materials, low organic matter content, high pH value, low aggregate stability, and shallow soil depth [33,38]. Therefore, the region has insufficient drainage and faces soil erosion [39].  All the soil samples were collected from 0-20 cm soil profiles and taken from at least five points in each plot. These five samples were combined to form a representative sample. These samples were air-dried under laboratory conditions, passed through a Water 2022, 14, 3875 4 of 20 2 mm sieve, and thoroughly mixed. Finally, these processed soil samples were stored for laboratory analyses.
Physical and chemical soil properties were measured using relevant standard laboratory analytical methods (Table 1). Based on the results, porosity and effective porosity values were calculated using Equations (1) and (2).
where P is porosity (cm 3 cm −3 ), Ps is particle density (g cm −3 ), Pb is bulk density (g cm −3 ), EP is effective porosity (cm 3 cm −3 ), and FC is field capacity by volume (cm 3 cm −3 ).  [46] Experimental Ks data was absent in the study area; therefore, the estimation of Ks values was performed by the widely used methods of Saxton and Rawls [12] and Van Genuchten et al. [17]. In the end, the correlation matrix was developed to see the relation between estimated Ks values by the Saxton and Rawls [12], and the Van Genuchten et al. [17] methods, and soil input data (Table 2). Accordingly, the sum of absolute values, except insignificant values of correlation coefficients, was 2.517 for Saxton and Rawls' [12] method, while that for Van Genuchten et al.'s [17] method was 3.766. Therefore, the result of the latter was selected to evaluate the performance metrics of the machine learning methods.

Machine Learning Methods
In this study, neural-network-based machine learning methods, artificial neural networks (ANN) and deep learning (DL), and tree-based machine learning methods such as decision tree (DT) and random forest (RF) were employed to estimate the Ks parameter.

Artificial Neural Networks
The artificial neural networks (ANN) algorithm is a mathematical model which is inspired by human nervous systems. This powerful tool can handle and solve complex and difficult problems due to its structure. Artificial neurons are the main units of a neural network that are connected by weights inside the layers. The general working principle of an ANN model is based on training the model first, then validating it for performance evaluation; this process is repeated until an acceptable level of error is encountered. Specifically, each artificial neuron receives a weighted input. These inputs are the outputs of neurons in the previous layer or input variables. After this procedure, the model sums up the inputs and adds a bias term, and then passes the results using an activation function [47]. A typical ANN structure consists of an input layer, one or more hidden layers, and an output layer. The number of hidden layers and neurons increases with the complexity of data.
Several NN types can be used depending on the requirements of the application [48]. The most used are perceptron, feed-forward (FF), convolution, recurrent, Kohonen maps, and support vector machines (SVM). Perceptron is the most basic and smallest NN that performs certain computations to detect features in the input data. Having a simple structure, they are only capable of implementing linearly separable problems. FF NNs on the other hand, find applications in more complex applications such as image processing, computer vision, and speech processing. They can be further classified into single and multi-layered NNs, where the number of layers depends on the complexity. Apart from this flexibility, they can deal with data that contain significant noise and are fast and easy to implement. In contrast, convolution NNs are complex to design and slow in performance depending on the number of hidden layers. For sophisticated applications such as text auto-suggest, grammar checking, text-to-speech, and translation, recurrent NNs are used because they are capable of modeling sequential data. However, training these NNs can be a challenging task. Kohonen maps are used in specialized applications to recognize patterns in the data, for instance in medical analysis to cluster data into different categories. SVMs, which are considered very robust for prediction applications, analyze the data for classification and regression analysis.
The ANN is a well-known and widely adopted method for modeling hydrological studies [49,50]. In particular, the method is used in many studies for the estimation of Ks [6,[8][9][10][11][12][13][14]16,17]. In this study, the ANN is implemented as a reference machine learning method to compare the performance metrics of other methods.

Deep Learning
Due to satisfactory results and great potential, deep learning (DL) was first introduced by Dechter [51] and has been increasingly used for hydrological and agricultural studies over the last years [52][53][54]. DL is an extended version of ANN [55]. The main difference between DL and ANN is that it uses more hidden layers [56]. This feature gives the possibility of higher learning skills and modeling performance, especially for complex datasets. A typical DL structure consists of input, hidden, and output layers.

Decision Tree
A decision tree (DT) is one of the widely used methods among machine learning models for solving classification and regression problems. The DT model uses a tree figure to show its correlation with the output data using the observed data in the dataset to be analyzed [57]. The leaves symbolize the output of the model, while the branches symbolize the connection of the input features of the model. One of the most important features of the DT method is to convert complex decision-making problems to simpler and understandable problems by dividing them into a collection of simple decisions [58].

Random Forest
A random forest (RF) is an ensemble model which was developed by Breiman [59] and is used for solving classification and regression problems. It is an improved version Water 2022, 14, 3875 6 of 20 of the DT model which consists of multiple decision trees. This improvement is a result of the fact that more trees there are, the more robust the forest [60]. The limitations of DT compared with RF are that the former has overfitting problems since the RF model uses each DT randomly, and the output is an average of the individual DT estimation.

Selected Inputs and Model Development
The neural-network-based (ANN and DL) and tree-based machine learning methods (DT and RF) were used to simulate the Ks value. To establish the models, soil texture (clay, silt, and sand), effective porosity, bulk density, permanent wilting point, field capacity, lime, organic carbon, and porosity were used as input variables. These input variables are associated with the Ks values of soils, and they have been applied in previous studies for the estimation of Ks [6,[8][9][10][11][12][13][14]16,17].
It is well-known that the Ks value interacts with the physical and mechanical properties of the soil. Therefore, the correlation of soil data belonging to this interaction with the Ks value is shown in the correlation matrix and the combinations were developed accordingly ( Table 2). Combination 7 consisted of inputs used in the Van Genuchten et al. [17] method, while combination 8 of input was used in the Saxton and Rawls [12] method. In this way, the performance metrics of machine learning methods by using input combinations of these methods could be emulated. The data reduction technique was applied for the development of the combinations, which aimed to take advantage of the available soil feature to estimate the Ks value at a high rate. Effective porosity was added in the first 6 combinations due to the high correlation. Sand and clay contents were taken into consideration in every combination and the effects of other soil mechanical properties were also evaluated respectively. According to this phenomenon, the most relevant input combinations were created to estimate the value of Ks. The combinations of the input variables can be seen in Table 3. All the datasets were normalized between 0 and 1 for minimizing the incoherency and redundancy of the data using the following equation: where Y norm is the normalized data of Y i , Y i is the observed data, Y max and Y min are the maximum and minimum data, respectively. The observed 291 datasets were employed to estimate the accuracy of the models using k-fold cross-validation, which is more reliable than the train-test split method [56]. In this study, the k value of 5 was used, which means the dataset was separated 5-fold, and models were trained using 4 (k − 1) fold. The trained model was tested on the remaining 1-fold. This procedure was run 5 (k) times, thus repeating the experiments 5 times and the results were averaged (Figure 2). Figure 3 shows a flowchart of the implemented model. maximum and minimum data, respectively. The observed 291 datasets were employed to estimate the accuracy of the models using k-fold cross-validation, which is more reliable than the train-test split method [56]. In this study, the k value of 5 was used, which means the dataset was separated 5-fold, and models were trained using 4 (k-1) fold. The trained model was tested on the remaining 1-fold. This procedure was run 5 (k) times, thus repeating the experiments 5 times and the results were averaged (Figure 2). Figure  3 shows a flowchart of the implemented model.

Performance Evaluation
The performance metrics of the ANN, DL, DT, and RF models were used including mean absolute error (MAE), root means square error (RMSE), coefficient of determination (R 2 ), and correlation coefficient (CC). The calculations of these metrics are shown below:

Performance Evaluation
The performance metrics of the ANN, DL, DT, and RF models were used including mean absolute error (MAE), root means square error (RMSE), coefficient of determination (R 2 ), and correlation coefficient (CC). The calculations of these metrics are shown below: where Z i is the predicted value, and Y i is the observed value. Y is the mean value of observed values and Z is the mean value of the predicted values. For the evaluation of the model performance, the high model performance can be confirmed when the RMSE and MAE values are low, and R 2 and CC values are high.

Adjustment of Input Variables
The correlation matrix of the input and output soil variables is presented in Table 2. This matrix helps to develop combinations for the models, as there is no significant relation between aggregate stability and penetration resistance variables with Ks value. Therefore, these two soil input variables were not used in the combinations. However, soil organic carbon, despite having insignificant relation with Ks value, is still used in combination 8 because this combination is used as input variables for Saxton and Rawls [12] equation. The other soil input variables such as sand, silt, clay, field capacity, permanent wilting point, porosity, effective porosity, and lime contents were used for the development of combinations. The findings demonstrated that effective porosity and Ks values (0.788) have the highest correlation among the soil input variables. The second highest correlation of 0.548 was observed between Ks and sand values. The third highest correlation of 0.226 was obtained between Ks and porosity. Among the significant values, other soil input variables such as silt, clay, bulk density, field capacity, permanent wilting point, and lime had a negative correlation with Ks. The statistical values of the soil data can be seen in Table 4.
The negative values of skewness were −0.06, −0.06, −0.57, and −0.13 mm h −1 for silt, field capacity, porosity, and effective porosity, and the positive values of kurtosis were 0.12, 0.08, 0.12 and 3.37 mm h −1 for bulk density, porosity, lime, and organic carbon, respectively. The maximum, minimum, mean, and standard deviation values of Ks were 58.91, 0.88, 16.08, and 11.33 mm h −1 , respectively. The highest value of variation coefficient was observed for Ks values with 0.66 mm h −1 . The highest value of skewness was obtained by Ks value with 0.86 mm h −1 and the lowest value of kurtosis was obtained by an effective porosity value of −0.88. Table 4. Statistical values of the Ks and soil data.

Performance of Machine Learning Methods
The performance metrics of four supervised machine learning methods using eight combinations of measured soil input variables are shown in Tables 5 and 6, and Figure 4 for estimation of Ks. Scatter and residual plots of observed and simulated Ks values are shown in Figures 5-8 for the ANN, DL, DT, and RF methods, respectively, for the eight soil input combinations. In general, the highest MAE and RMSE values were observed for the ANN6 method with soil inputs of sand, clay, effective porosity, permanent wilting point, field capacity, bulk density, porosity, and lime contents, and the highest R 2 and CC values were observed for the ANN7 method with soil inputs of sand, silt, clay, permanent wilting point, field capacity, and bulk density. It can be concluded that neural-network-based models (ANN and DL) had the best performance metrics for combination 7 with soil inputs of sand, silt, clay, permanent wilting point, field capacity, and bulk density values, while they had the lowest performance metrics for combination 2 with soil inputs of clay and effective porosity. For tree-based models (DT and RF), the lowest performance was obtained for combination 8 with soil inputs of sand, clay, bulk density, and organic carbon contents, while the best performance was obtained for combination 5 with soil inputs of sand, clay, effective porosity, and bulk density.

Performances of Neural-Network-based Machine Learning Methods
The best performance was observed when the ANN 2

Performances of Neural-Network-Based Machine Learning Methods
The best performance was observed when the ANN 2(2-3-4-4-8-6-4)-4-1, which means that the algorithm consists of two neurons of combination 1, two neurons of combination 2, three neurons of combination 3, four neurons of combination 4, four neurons of combination 5, eight neurons of combination 6, six neurons of combination 7 and four neurons of combination 8 in input layers, four neurons in the hidden layer and one output layer. The rectified linear unit (ReLU) was employed as the activation function in this study since it is the most used activation function. The ANN7 model fed with soil inputs of sand, silt, clay, permanent wilting point, field capacity, and bulk density demonstrated the highest performance metrics among the other models considering MAE (2.407 mm h −1 ), RMSE  Figure 5A. The data of the scatter plots are generally close to the reference line (1:1) for all combinations. However, the combinations of 1, 2, and 8 for the ANN model showed more scattered points than combinations of 3, 4, 5, 6, and 7. The residual plots of estimated Ks values by the ANN model with eight soil input combinations are shown in Figure 5B. The residual plot demonstrated that the most errors occurred in combination 8, while the least error occurred in combination 6 for the ANN model.
The best accuracy of the Ks value was observed when two hidden layers  and ReLU were used for the DL model in this study. As can be seen in Table 5, it is seen that the DL model produces the lowest performance for combination 2 with MAE, RMSE, R 2 , and CC equal to 4.898, 6.965, 0.707, and 0.840. The combination of DL7 with soil inputs of sand, silt, clay, field capacity, permanent wilting point, and bulk density showed the best performance with MAE of 2.167 mm h −1 , RMSE of 3.423 mm h −1 , R 2 of 0.919, and CC of 0.959. The scatter plots of estimated Ks values by the DL model with eight soil input combinations are shown in Figure 6A. From the figure, it can be seen that combinations 2 and 8 had more scattered points than other combinations. The least scattered points were observed for the combination of 7. The residual plots of estimated Ks values by the DL model with eight soil input combinations are shown in Figure 6B. The least residual errors were observed for the combination of 7, while the most residual errors were observed for the combination 8 for the DL model.

Performances of Tree-Based Machine Learning Methods
The lowest performance was obtained for the DT8 fed with soil inputs of sand, clay, bulk density, and organic carbon contents considering MAE (3.179 mm h −1 ), RMSE (5.736 mm h −1 ), R 2 (0.887), CC (0.942). The performance metrics improved significantly when adding organic carbon instead of bulk density values from the soil inputs. In that case, the highest performance was observed for the fifth combination with MAE of 2.  Figure 8A. The most scattered points were obtained from combination 8, while the least scattered points were obtained from combination 7 for the RF model. The residual plots of estimated Ks values by the RF model with eight soil input combinations are shown in Figure 8B. According to the figure, similar residual errors were observed for the combination of 3, 4, and 5. However, among them, the least residual errors occurred in combination 5 for the RF model.

Discussion
The present study evaluated neural-network-based (ANN and DL) and tree-based (DT and RF) models with different combinations of soil input variables based on the van Genuchten formula in the semi-arid environment. Zhang and Schaap [2] suggested in their studies that new statistical methods should be employed using relevant and good-quality data for the estimation of Ks. In this study, the results indicated that all the machine learning methods used have a satisfactory correlation between Ks and soil input variables. Due to solving complex and nonlinear features, machine learning methods can simulate the complicated process of soil nature because these methods do not need to know the characteristic of the implemented variables [61].
Another aim of this study was to estimate the Ks values with the least input combinations ensuring high accuracy. For this purpose, eight combinations were developed using ANN, DL, DT, and RF models to evaluate the estimation of Ks value. Accordingly, the first six combinations included effective porosity to see the performances of the models, since it has the highest impact on soil water transmission and correlates to the selected soil properties, except for penetration resistance and organic carbon. The first three combinations are developed based on clay and sand from the soil texture, as it is known that sand and clay contents are one of the most important factors that directly affect the Ks in soil [62][63][64]. When the soil contains a high amount of sand, the Ks are increased. On the contrary, when the soil contains a high amount of clay, the Ks are decreased. This finding is in agreement with Table 2. It can be seen that the performance metrics of the first three combinations have similar results since the use of clay and sand contents together offer high accuracy to the models. It has been stated in many studies that the field capacity is used in the estimation of Ks [12,65,66]. However, the field capacity does not have a direct impact on the estimation of Ks and therefore does not impact the result significantly. Since bulk density is an indicator of soil compaction, it impacted positively on the estimation of Ks. The compression of the pores starts from the effective porosity and progresses towards the micropores [14]. With the decrease of effective porosity, infiltration is decreased, and runoff is increased. Finally, bulk density increases. In this respect, in five combinations, the soil textures and bulk density have a high impact on machine learning methods for the estimation of Ks. Lime content and porosity are calculated by bulk density in this study. Therefore, the combination which was developed from these two soil parameters does not impact the estimation of Ks values. The previous study conducted in this study area showed that soil lime content plays a key role in water retention up to certain levels, but this efficiency decreases due to rising lime content [67]. The decrease in the disclosure rate in combination 8 is explained by the organic carbon value not being correlated with Ks values. When many studies are examined, it is seen that the organic carbon value affects the water movement of the soils. However, in this study, it is seen in the correlation matrix values that it has no effect, since the organic carbon in the soils of the study area is very low with an average of 0.85 and the variability in the carbon value is limited.
The performance metrics demonstrated that the ANN was superior to the DL method in all soil input combinations for the neural-network-based models. This finding can be explained by the fact that the DL method requires more experimental data for boosting its modeling performance. A similar explanation was pointed out also by Kamilaris and Prenafeta-Boldú [53] who indicated that the usage of the DL method has recently increased in soil science due to its ability to solve complicated datasets. The performance metrics showed that the RF model was a better result than the DT model in all soil input combinations for the tree-based models. RF is an improved version of the DT model since it boosts the robustness of the classification feature. However, the RF is known as a black-box model and therefore it is less explicable than the DT model, but it is seen that the RF models demonstrate better performance and stability than the estimation in previous studies.
The RF model with combination 5 had the highest performance among the tree-based models, while the ANN model with combination 7 had the highest performance among the neural-network-based models. It can be seen from Table 3 that combination 5 had four input parameters but combination 7 had six input parameters. This can be explained by the fact that the ANN model had a much more complicated architecture than the RF model, which allowed for boosting its modeling performance with more input variables.
In general, the ANN model demonstrated better performance metrics compared to other models. This observation is in agreement with the results of [68] who reported that the ANN model shows satisfactory results for the estimation of Ks. Similarly, [69] indicated that the ANN demonstrates good performance for the estimation of Ks.

Conclusions
Soil hydraulic conductivity is a fundamental parameter for the estimation of water balance at regional and global levels and the determination of groundwater recharge in the vadose zone. In recent years, mathematical methods have been applied frequently by using soil features to obtain Ks. This study demonstrated that the machine learning methods can be an alternative way to overcome these challenges.
The comparison analysis of neural-network-based models (ANN and DL) and treebased models (DT and RF) models was evaluated based on the van Genuchten equation under eight combinations of soil input data. The reason for using two group of machine learning methods is to find out the best category. In general, all the machine learning methods had satisfactory performance for estimation of Ks. Among the categories, the neural-network-based models (ANN and DL) had better performance than tree-based models (DT and RF). The overall results showed that the ANN method with sand, silt, clay, field capacity, permanent wilting point, and bulk density (ANN7) had the best performance among the other methods. The RF method with sand, clay, effective porosity, and bulk density (RF5) had the best performance among the tree-based models. These findings demonstrate that the ANN method is more applicable to the study area. Likewise, the ANN method can figure out a much more complex dataset to estimate Ks for calcareous alluvial soils in a semi-arid region. In the case of the input data, it was observed that soil texture, bulk density, and effective porosity variables have a great impact on estimating the Ks value since soil properties, such as lime penetration resistance, aggregate stability, and porosity, have a low impact on the estimation of Ks. These findings are important to understand the impact of soil parameters in a study area for estimating the Ks value.
Lastly, further studies should collect more data for enhancing modeling performance. Likewise, the machine learning methods should be tested under different environmental conditions, which means seeing the performance metrics of different specific soil input combinations for the study area. This is particularly important to see the adaptation of new statistical methods such as machine learning.