Investigating hydraulic conductivity (K) is crucial for aquifer studies and groundwater flow modelling. The main objectives of the current study are to investigate the effectiveness of artificial neural network (ANN), adaptive neuro-fuzzy inference system (ANFIS), Gaussian process regression (GPR), and random forest (RF) algorithms in estimating K using data from 270 borehole soil samples, collected along the Beas riverbank in Kangra district, Himachal Pradesh, India. For the K estimation, the study utilizes the grain size parameters, i.e., d10, d50, coefficient of uniformity (Cu), and porosity (n) as input parameters. The performance evaluation of the developed models was assessed using the statistical parameters. While the performance of each model is quite satisfactory, the present study establishes the efficacy of the GPR model during validation having a determination coefficient of 0.985. The root mean square errors for ANN, ANFIS, GPR, and RF were 0.019, 0.017, 0.00853, and 0.019, respectively. The techniques used in the study offer precise K-prediction abilities that facilitate groundwater management and contaminant transport analysis. The GPR model in the study outperforms other models in estimating the K of soil samples and serves as an efficient tool for managing soil water and solute transport.

  • The primary goal of the study is to develop hydraulic conductivity models using four data-driven techniques, i.e., ANN, ANFIS, GPR, and RF.

  • These data-driven techniques enable quick estimation of hydraulic conductivity, facilitating precise determination of groundwater recharge rate.

  • Among all these models, the GPR model performs better in predicting the hydraulic conductivity of given soil samples.

Groundwater recharge is essential for sustaining water resources, especially in arid and semi-arid regions where water resources are scarce, and ecosystems face significant stress (Noori et al. 2023). Hydraulic conductivity (K), a fundamental parameter for evaluating the drainage potential of transmitting media, plays a vital role in shaping groundwater recharge rates and, consequently, the surrounding environment and ecosystems (Mahdian et al. 2024). The investigation of the K of porous media is essential for assessing soil stability, predicting landslides, and modelling groundwater flow. The K of porous media depends on the physical properties of the flowing fluid and transmission medium such as particle size, porosity, and pore connectivity (Chandel et al. 2021). Despite the importance of K in geological and geotechnical characterization, the evaluation of K remains a challenging task due to the heterogeneity of the soil particles, difficulties in soil sampling, subsequent experimental procedure, and the large extent of flow area within soil deposit (López-Acosta et al. 2019). Measurement and prediction of soil's K value can be done using direct and indirect approaches (Mozaffari et al. 2024).

Hydraulic conductivity can be measured directly via field or laboratory methods. The laboratory methods include constant head and falling head tests, while the field methods include ring infiltrometer, instant profile, and test basins (Williams & Ojuri 2021; Chandel et al. 2023), borehole permeameter and pressure infiltrometer (Alagna et al. 2023), and auger hole and tension infiltrometer (Kargas et al. 2022). Due to temporal and geographical fluctuations, the direct measurement of K becomes impractical, costly, and time-consuming. As a result, indirect methods for estimating the K from more readily and affordably measured soil parameters, i.e., specific gravity, porosity, and the proportion of clay, silt, and sand content were developed and widely used (More et al. 2022). The indirect techniques comprise empirical equations and machine learning-based data-driven techniques for the estimation of K. The main advantage of empirical equations is that they can estimate the K value quicker than the direct measurement (Williams & Ojuri 2021). Due to the domain-specific development of the empirical equations, these are restricted to specific boundary conditions.

Recent studies demonstrate the effectiveness of data-driven techniques in addressing a range of environmental and engineering problems, including water engineering (Jafari-Asl et al. 2024), flood prediction (Ahmadi et al. 2024; Donnelly et al. 2024), climate parameters (Chatrabgoun et al. 2020; Yeganeh-Bakhtiary et al. 2023), and scouring (Habib et al. 2024). These techniques are capable of precisely and accurately performing complicated tasks, i.e., learning from the data, modelling, and deriving a relation from experimental data. The data-driven techniques can be either individually operated like support vector machines (SVMs), fuzzy, and artificial neural networks (ANNs) or they can be hybrid, i.e., adaptive neuro-fuzzy inference system (ANFIS) and genetic algorithm (GA)-ANN (Williams & Ojuri 2021). Sihag (2018) estimated the K of fly-ash-mixed soils using SVM, ANN, random forest (RF), Gaussian process regression (GPR), M5P model tree, and two other traditional models and concluded that SVM with Radial Basis Function (RBF) kernel was the most accurate technique among others. Boadu (2020) investigated Multiple Linear Regression (MLR) and support vector regression (SVR) for predicting K and concluded that the SVR model performed better for estimating K. Singh et al. (2021a, b) examined ANN, MLR, RF, and M5P tree-based models to determine K and outcomes demonstrated that all models have a reasonable ability for prediction, with RF-based models outperforming the others. Abdalrahman et al. (2022) predicted the infiltration rate of treated wastewater using ANN, multilayer perceptron model (MLP), and Elman neural network (ENN), and concluded that the MLP model exhibited superior performance. Similarly, Singh et al. (2021b) developed two hybrid machine learning-based pedotransfer functions (ML-PTFs) by integrating a genetic algorithm with multilayer perceptron (MLP-GA) and support vector machine (SVM-GA). Their results indicated that the SVM-GA PTF demonstrated greater efficiency compared with the MLP-GA algorithm.

The evaluation of the research literature revealed that several data-driven methods are available for determining the K using easily measurable soil properties. ANFIS uses the Takagi-Sugeno fuzzy inference system, an effective tool for capturing nonlinear relationships between inputs and outputs (Seifi et al. 2020). GPR is a probabilistic, nonparametric model that can model complex systems and has a favourable nonlinear mapping ability (Richardson et al. 2017). RF is a meta-estimator that controls overfitting and increases predictive accuracy through the process of fitting several decision tree classifiers to various dataset subsamples and then averages the outcomes. The analysis of the literature review infers that previous studies have not integrated the applications of ANN, GPR, ANFIS, and RF techniques for estimating the hydraulic conductivity of porous media. These techniques can establish a correlation between the dependent and independent variables, even when the dataset contains inconsistent or noisy data points. Therefore, the study focuses on developing and validating the K models using four different techniques. Moreover, the study examines the performance of models through statistical indicators to identify the most suitable model for K estimation. The study presents a detailed comparison, offering valuable insights into how data-driven techniques can enhance the prediction of hydraulic conductivity (K) in porous media.

Experimental setup

In the initial phase, experiments including specific gravity, sieve analysis, and constant head permeameter tests were conducted at NIT Hamirpur's Hydraulic Laboratory. The present investigation involves the collection of 270 borehole samples with distinct diameters of fines (<0.075 mm), fine sand (0.075–0.425 mm), coarse sand (2–4.75 mm), and fine gravel (4.75–20 mm) from the banks of River Beas in the District Kangra of Himachal Pradesh, India. The geographical coordinates of the sampling locations are 31.8818° N latitude and 76.2146° E longitude. Soil samples were sieved according to the American Society for Testing and Materials standards (ASTM 2007). From the grain size curve, d50 and d10 (grain diameter at which 50 and 10% particles are finer), the coefficient of uniformity (Cu) was determined. To determine the porosity (n) of the sample, the specific gravity test was carried out using the pycnometer method (Chandel & Shankar 2022). The vertical flow type constant head permeameter, as presented in Figure 1, was used for the calculation of ‘K’ which consists of three permeameter diameters, i.e., 5.08, 10.16, and 15.24 cm, having test length and total sample length of 46.5 and 106 cm, respectively. Four pressure-taping points were positioned at a 90° angle around the circumference of the permeameter to determine the head loss between the top and bottom of the soil sample.
Figure 1

Laboratory set-up of constant head permeameter.

Figure 1

Laboratory set-up of constant head permeameter.

Close modal

Model development

In the second phase, ANN, ANFIS, GPR, and RF techniques were employed to estimate K of the porous medium utilizing observations obtained from the experimental investigations. The accuracy of these models is checked concerning the experiment data. To verify the model's performance, statistical measures like determination coefficient (R2), root mean square error (RMSE), relative absolute error (RAE), mean absolute error (MAE), Nash–Sutcliffe efficiency (NSE), root mean squared logarithmic error (RMSLE), and Kling–Gupta efficiency (KGE) were used. An overview of the modelling techniques employed in the study is given below.

Artificial neural network

The concept of ANN was developed by Frank Rosenblatt and was inspired by neural circuitry, designed to simulate a range of real-world situations by modelling the neuronal matrix of a human brain. The ANN structure comprises interconnected artificial neurons that can process large numbers of predictors and learn from them (Ewuzie et al. 2022). Each neuron in ANN receives input signals and processes them by adding a bias and computes the weighted sum of the inputs. Thereafter, an activation function receives this weighted total as an input to produce the output. Figure 2 depicts the fundamental structure of the ANN model. ANN is widely used due to its ability to model complex, and nonlinear relationships in a dataset (Zhang et al. 2018). The performance of ANNs depends on carefully setting parameters like the number of hidden layers and neurons, but they can overfit readily, particularly with smaller datasets. However, ANNs lack built-in mechanisms to estimate prediction uncertainty, making it challenging to assess the reliability of their predictions.
Figure 2

Basic architecture of the ANN model.

Figure 2

Basic architecture of the ANN model.

Close modal
ANN modelling strategy
In ANN, the input signal propagates through the network of neurons in a forward direction from one layer to another (Ewuzie et al. 2022). In the present study, a multilayer feed-forward network was adopted, where weights and biases are further updated using the backpropagation algorithm. The number of neurons in the hidden layer is determined using trial and error (Thakur et al. 2022). The neural fitting toolbox in MATLAB R2023b was used to carry out the ANN modelling in accordance with the steps illustrated in Figure 3.
Figure 3

Flowchart of ANN modelling methodology.

Figure 3

Flowchart of ANN modelling methodology.

Close modal

Adaptive neuro-fuzzy inference system

Jang established the most effective neuro-fuzzy model, ANFIS, in 1993 (Azad et al. 2018). The behaviour of the local system is based on each fuzzy rule in the FIS. Based on the input–output dataset, it generates a fuzzy inference system (FIS) (Seifi et al. 2020). FIS is implemented by the network structure, which also trains the model using hybrid learning rules. The layers and nodes in the ANFIS's structure as shown in Figure 4 are explained below (Jalalkamali 2015).
  • Layer 1: The membership grade of an input variable is generated by each node. The fuzzy set associated with each node is categorized using a variety of membership functions (MFs), including generalized bell-shaped, Gaussian, triangle-shaped, and trapezoidal-shaped functions (Seifi et al. 2020).

  • Layer 2: The incoming signal, which specifies the extent to which a rule is firing, determines each node's output, which is a fixed node with the label ∏.

  • Layer 3: Each node computes the normalized firing strengths and is a fixed node with label N.

  • Layer 4: For each adaptive node (i), compute the input of the ith rule using the node function for predicting the model's output.

  • Layer 5: Adding together the outputs of all incoming signals gives the total output of the ANFIS.

Figure 4

Basic architecture of the ANFIS model.

Figure 4

Basic architecture of the ANFIS model.

Close modal

ANFIS combines the advantages of neural networks and fuzzy logic, making it effective for modelling complex and uncertain systems. The key benefit of ANFIS lies in its interpretability, as it uses fuzzy rules to clarify the relationship between inputs and outputs. However, ANFIS can struggle with high-dimensional datasets, and its efficacy declines, when the input data points increase. ANFIS, like ANN, does not provide a built-in way to quantify uncertainty in its predictions.

ANFIS modelling strategy
In ANFIS, the neural network's learning potential is combined with fuzzy logic's reasoning capabilities. To develop the ANFIS model, ANN learns by utilizing the training datasets that are available. At the same time, a fuzzy model is mapped into the derived solution, which is essentially a set of ‘If-Then’ rules. Modelling using ANFIS involves deciding the shape and number of MFs. A fuzzy interface system (FIS) was created using three MFs for each parameter. For each input and output parameter, three shapes of MF, i.e., Gaussian, triangular, and generalized bell, were taken. The hybrid optimization algorithm was utilized to optimize the training parameters. MATLAB R2023b's neuro-fuzzy designer toolbox was used to do ANFIS modelling by the procedures shown in Figure 5.
Figure 5

Flowchart of ANFIS modelling methodology.

Figure 5

Flowchart of ANFIS modelling methodology.

Close modal

Gaussian process regression

According to Wang (2023), a Gaussian process is a set of random variables, and any finite number of Gaussian processes has a joint Gaussian distribution. A Gaussian process can be fully described by its variance, covariance, and mean functions. The GPR is a probabilistic machine learning technique that provides not only predictions but also uncertainty intervals, making it highly effective for uncertainty quantification (UQ). Compared with methods like convolutional neural network, GPR excels in UQ as it measures prediction confidence by estimating variance concerning mean (Donnelly et al. 2022). This ability is important in situations where the reliability of prediction matters. According to Donnelly et al. (2022), the following noise dataset can be shown in Equation (1):
(1)
where N is the number of data points, x is the input, and y is the output.
For the present study, d10, d50, Cu, and n are the input parameters and X is the output of the GPR. Thus, the x value can be computed using Equation (2):
(2)
The above data can be generated from Equation (3):
(3)
where ξ is the Gaussian distribution (zero mean and variance σ2) and ε is the Gaussian noise term. Equation (4) shows the joint distribution of y.
(4)
where I is the identity matrix and K(x, x) is the kernel function. Equations (5) and (6) can be used to find the GPR's mean and variance.
(5)
(6)
where σ2 is the variance and T is the transpose.
GPR modelling strategy
A GPR is a nonparametric, probabilistic model (Seyedian et al. 2023). With MATLAB R2023b's Regression Learner toolbox, the GPR models with squared exponential, rational quadratic, matern 5/2, and exponential kernel functions were developed and validated. To avoid the model from overfitting, the model was run through a 5-fold cross-validation. Figure 6 illustrates the process used for GPR modelling.
Figure 6

Flowchart of GPR modelling methodology.

Figure 6

Flowchart of GPR modelling methodology.

Close modal

Random forest

RF is a supervised machine learning technique that combines multiple machine learning algorithms to increase prediction efficiency (Azarhoosh & Koohmishi 2023). To reduce the variation of the learning strategy, manifold decision trees are created while operating concurrently, as shown schematically in Figure 7. To train independently, the training data in this case are first split into many subsets by repeatedly taking samples, known as bootstrapping. The findings of each model are then combined and averaged to obtain the final output. RF is widely recognized for its ability to handle high-dimensional data effectively and prevent overfitting through random sampling and feature selection. This technique is efficient, easy to implement, and highly versatile, making it suitable for both classification and regression tasks. However, RF is a deterministic method and does not account for uncertainty in its predictions, making it unable to provide confidence intervals for its outputs (Tyralis & Papacharalampous 2024).
Figure 7

A typical layout representing the RF model.

Figure 7

A typical layout representing the RF model.

Close modal
RF modelling strategy
The Python programming language was used to create the RF algorithm in the present study. In this technique, n-estimators describe the total number of trees in the forest, ‘min_sample_leaf’ considers the absolute minimum of samples required to be at a leaf node, and max-depth denotes the maximum depth of the trees (Sihag et al. 2019). For separate training of each subset, the training data is first divided into many subsets utilizing repeated sampling. The results of all the models are then combined and averaged to produce the final output.
(7)
where y is the final output, fb(x) is the bootstrapped training set (bth), and B is the number of training sets.

Performance evaluation of developed models

To evaluate the accuracy of developed models, six different statistical indicators, namely RMSE, R², MAE, NSE, RAE, RMSLE, and KGE (Habib et al. 2023), were computed. Lower values of RMSE, RAE, MAE, and RMSLE that are closer to 0 and R2, KGE, and NSE that are closer to 1 suggest precise agreement between the observed and projected values (Thakur et al. 2025). The equations to compute the statistical indicators are as follows:
(8)
(9)
(10)
(11)
(12)
(13)
(14)
where N is the total number of samples, Ei is the experimentally observed K values, Pi is the predicted K values by the model, E is the experimentally obtained mean value, P is the model predicted mean values, r is the linear correlation coefficient between the experimentally obtained and model predicted values. is the ratio of standard deviation of the model predicted () and experimentally obtained () values, and is the ratio of the model predicted and experimentally obtained mean values.

The number of training datasets determines how well the model estimates K in the porous medium. A total of 270 datasets were selected for the study, out of which 200 (75%) were used for training the model, while the remaining 70 (25%) samples were utilized to test the model. In the present study, four input parameters ‘d10, d50, Cu, and n’ are used to predict the K. The initial part of the study focuses on statistical analysis and the dataset employed for model development and validation. While the subsequent part focuses on K estimation through various models, along with their quantitative assessment.

Statistical analysis

Initially, grain size analysis was conducted on the collected soil samples to compute the key characteristics, including d10, d50, and Cu. Additionally, the specific gravity and constant head permeameter tests were performed to determine the n and K of the soil samples. A statistical summary of the experimental observations including minimum, maximum, mean, and standard deviation is presented in Table 1. The observed K values for the soil samples vary from 0.342 to 0.010 cm/s. Further, a relationship between the independent parameters (d10, d50, Cu, and n) and dependent parameters (K) was analysed using a correlation matrix, as presented in Table 2. The data in Table 2 reveal that d10 and d50 have a stronger influence on K with correlation coefficients of 0.92 and 0.82, respectively, while Cu and n show a moderate correlation with K, having coefficients of −0.65 and 0.70, respectively. Furthermore, a correlation heat map (Figure 8) has been developed to visually illustrate the relationship between the dependent and independent parameters. Figure 8 depicts a strong correlation of d10 and d50 with K, indicated by intense red hues in their corresponding cells. In contrast, Cu and n exhibit a moderate correlation with K, as shown by the lighter colours.
Table 1

Statistics summary of the characteristics of soil samples

CharacteristicsMinimumMaximumMeanStandard deviation
d10 (mm) 0.100 0.409 0.201 0.093 
na 0.286 0.426 0.358 0.028 
d50 (mm) 0.289 1.450 0.660 0.309 
Cua 2.150 9.567 4.691 1.559 
K (cm/s) 0.010 0.342 0.067 0.068 
CharacteristicsMinimumMaximumMeanStandard deviation
d10 (mm) 0.100 0.409 0.201 0.093 
na 0.286 0.426 0.358 0.028 
d50 (mm) 0.289 1.450 0.660 0.309 
Cua 2.150 9.567 4.691 1.559 
K (cm/s) 0.010 0.342 0.067 0.068 

aDenotes the dimensionless characteristics.

Table 2

Correlation matrix between various soil parameters and K

d10nd50CuK
d10     
n 0.13    
d50 0.91 −0.22   
Cu −0.18 −0.96 0.17  
K 0.92 0.70 0.82 −0.65 
d10nd50CuK
d10     
n 0.13    
d50 0.91 −0.22   
Cu −0.18 −0.96 0.17  
K 0.92 0.70 0.82 −0.65 
Figure 8

Correlation heatmap between the dependent and independent parameters.

Figure 8

Correlation heatmap between the dependent and independent parameters.

Close modal

Estimation of K using ANN

In the present study, ANN models were developed using 3, 4, 5, and 6 hidden neurons. The outcomes of the ANN model's performance evaluation parameters for the training and testing datasets, as shown in Table 3. ANN 3, ANN 4, ANN 5, and ANN 6 demonstrate the ANN models with 3, 4, 5, and 6 hidden neurons, respectively. The ANN model with five hidden neurons performs better than the other ANN models. For the training dataset, the best-performing ANN 5 model achieved R2, RMSE, RAE, MAE, NSE, RMSLE, and KGE values of 0.951, 0.0157, 0.1458, 0.0080, 0.9464, 0.0129, and 0.9237, respectively. While, for the testing dataset, the corresponding values were 0.902, 0.0206, 0.1982, 0.0108, 0.8928, 0.0181, and 0.8646, respectively. Figure 9 displays scatter plots of the empirically observed K values and the model predicted K values for the training and testing datasets, respectively. During training, the ANN models demonstrated satisfactory performance, but their performance was less than optimal during testing. The scatter plot of the testing dataset (right) reveals that points tend to cluster near the agreement line when K values are lower, suggesting superior prediction performance. However, when using larger K values, the predictive capability was only moderate because of the restricted data range for higher K values, necessitating larger datasets for effective training.
Table 3

ANN models' performance evaluation parameter of the training and testing datasets

Statistical indicatorsANN 3 model
ANN 4 model
ANN 5 model
ANN 6 model
TrainingTestingTrainingTestingTrainingTestingTrainingTesting
R2 0.930 0.882 0.927 0.895 0.951 0.902 0.947 0.889 
RMSE 0.0191 0.0223 0.0195 0.0217 0.0157 0.0206 0.0159 0.0218 
RAE 0.1278 0.2236 0.2048 0.2713 0.1458 0.1982 0.1567 0.2158 
MAE 0.0070 0.0121 0.0112 0.0147 0.0080 0.0108 0.0075 0.0117 
NSE 0.9231 0.8777 0.9199 0.8842 0.9464 0.8928 0.9464 0.8838 
RMSLE 0.0159 0.0197 0.0163 0.0194 0.0129 0.0181 0.0135 0.0191 
KGE 0.8762 0.8378 0.8777 0.8594 0.9237 0.8646 0.9200 0.8575 
Statistical indicatorsANN 3 model
ANN 4 model
ANN 5 model
ANN 6 model
TrainingTestingTrainingTestingTrainingTestingTrainingTesting
R2 0.930 0.882 0.927 0.895 0.951 0.902 0.947 0.889 
RMSE 0.0191 0.0223 0.0195 0.0217 0.0157 0.0206 0.0159 0.0218 
RAE 0.1278 0.2236 0.2048 0.2713 0.1458 0.1982 0.1567 0.2158 
MAE 0.0070 0.0121 0.0112 0.0147 0.0080 0.0108 0.0075 0.0117 
NSE 0.9231 0.8777 0.9199 0.8842 0.9464 0.8928 0.9464 0.8838 
RMSLE 0.0159 0.0197 0.0163 0.0194 0.0129 0.0181 0.0135 0.0191 
KGE 0.8762 0.8378 0.8777 0.8594 0.9237 0.8646 0.9200 0.8575 
Figure 9

Scatter plot between the experimentally measured and model predicted K values for the training (left) and testing (right) stages. (a) ANN 3, (b) ANN 4, (c) ANN 5, and (d) ANN 6.

Figure 9

Scatter plot between the experimentally measured and model predicted K values for the training (left) and testing (right) stages. (a) ANN 3, (b) ANN 4, (c) ANN 5, and (d) ANN 6.

Close modal
Figure 10 presents the box plots of experimentally measured and ANN model predicted K values for the training and testing datasets. In these plots, the median is positioned near the lower quartile, indicating positively skewed data distribution, whereas a larger proportion of the data points are concentrated below the median. For the training dataset, the ANN models demonstrate strong alignment with the experimental values reflecting low variability and good agreement. For the testing dataset, the ANN 5 model exhibits close alignment with the experimental values, whereas the ANN 3, ANN 4, and ANN 5 models comprise few outliers, indicating data points that fall outside the desired range.
Figure 10

Box plot of the experimentally measured and ANN model predicted K values for the training (left) and testing (right) stages.

Figure 10

Box plot of the experimentally measured and ANN model predicted K values for the training (left) and testing (right) stages.

Close modal
For more comprehensive visualization of the data distribution, violin plots were plotted for the experimentally observed and model predicted K values, as shown in Figure 11. In the training dataset, the violin plot reveals a narrow distribution around the media indicating a strong alignment between the model's prediction and the experimentally observed data. While the violin plot for the testing dataset shows wider tails reflecting a broader spread in the predicted values. This suggests that the model's predictive accuracy is more variable in the testing phase with less accurate predictions observed for particular higher K values.
Figure 11

Violin plot of the experimentally measured and ANN model predicted K values for the training (left) and testing (right) stages.

Figure 11

Violin plot of the experimentally measured and ANN model predicted K values for the training (left) and testing (right) stages.

Close modal

Estimation of K using ANFIS

Figure 12 presents the structure of the ANFIS model, in which each input parameter is linked with three MFs. The outcomes of the ANFIS model's performance evaluation parameters for the training and testing datasets are shown in Table 4. ANFIS triangular, ANFIS Gaussian, and ANFIS Gbell demonstrate the ANFIS models with triangular MF, Gaussian MF, and generalized Bell MF, respectively. From Table 4, it is evident that the ANFIS Gaussian model performs better than the ANFIS Gbell and ANFIS triangular models, with R2, RMSE, RAE, MAE, NSE, RMSLE, and KGE values of 0.973, 0.0112, 0.0840, 0.0046, 0.9735, 0.0093, and 0.9725 for the training dataset, respectively. Whereas the corresponding values for the testing dataset were 0.924, 0.0179, 0.1612, 0.0088, 0.9215, 0.0152, and 0.9373, respectively.
Table 4

Assessment metrics for ANFIS model performance on the training and testing datasets

Statistical indicatorsANFIS triangular
ANFIS Gaussian
ANFIS Gbell
TrainingTestingTrainingTestingTrainingTesting
R2 0.965 0.919 0.973 0.924 0.970 0.918 
RMSE 0.0134 0.0197 0.0112 0.0179 0.0125 0.0189 
RAE 0.0953 0.1785 0.0840 0.1612 0.0774 0.1763 
MAE 0.0052 0.0097 0.0046 0.0088 0.0048 0.0096 
NSE 0.9621 0.9060 0.9735 0.9215 0.9681 0.9121 
RMSLE 0.0110 0.0168 0.0093 0.0152 0.0087 0.0161 
KGE 0.9293 0.8561 0.9725 0.9373 0.9553 0.9274 
Statistical indicatorsANFIS triangular
ANFIS Gaussian
ANFIS Gbell
TrainingTestingTrainingTestingTrainingTesting
R2 0.965 0.919 0.973 0.924 0.970 0.918 
RMSE 0.0134 0.0197 0.0112 0.0179 0.0125 0.0189 
RAE 0.0953 0.1785 0.0840 0.1612 0.0774 0.1763 
MAE 0.0052 0.0097 0.0046 0.0088 0.0048 0.0096 
NSE 0.9621 0.9060 0.9735 0.9215 0.9681 0.9121 
RMSLE 0.0110 0.0168 0.0093 0.0152 0.0087 0.0161 
KGE 0.9293 0.8561 0.9725 0.9373 0.9553 0.9274 
Figure 12

ANFIS model with three MFs.

Figure 12

ANFIS model with three MFs.

Close modal
Figure 13 exhibits scatter plots illustrating the comparison between the model predicted K values and the empirically observed K values for both the training and testing datasets. The training and testing of ANFIS yielded good results, suggesting that the model built with ANFIS is capable of accurately predicting hydraulic conductivity. For detailed visualization of data distribution, box and violin plots were plotted for the experimentally observed and model predicted K values, as shown in Figures 14 and 15, respectively. The plots of the testing dataset reveal that as K values increase, the prediction capability tends to decrease. Conversely, for smaller K values, the points cluster more closely with the experimental data, indicating higher prediction performance. The knowledge acquired during training affects the prediction efficiency of the ANFIS model. The performance improves with higher values of input data, ANFIS models perform better in large datasets. The training dataset for the ANFIS model contained more K values in the range of 0.006–0.15, which improved model performance in the lower region during testing.
Figure 13

Scatter plot between the experimentally measured and model predicted K values for the training (left) and testing (right) stages. (a) ANFIS triangular, (b) ANFIS Gaussian, and (c) ANFIS Gbell.

Figure 13

Scatter plot between the experimentally measured and model predicted K values for the training (left) and testing (right) stages. (a) ANFIS triangular, (b) ANFIS Gaussian, and (c) ANFIS Gbell.

Close modal
Figure 14

Box plot of the experimentally measured and ANFIS model predicted K values for the training (left) and testing (right) stages.

Figure 14

Box plot of the experimentally measured and ANFIS model predicted K values for the training (left) and testing (right) stages.

Close modal
Figure 15

Violin plot of the experimentally measured and ANFIS model predicted K values for the training (left) and testing (right) stages.

Figure 15

Violin plot of the experimentally measured and ANFIS model predicted K values for the training (left) and testing (right) stages.

Close modal

Estimation of K using GPR

The performance assessment parameters for the training and testing datasets of the GPR model are determined by using different kernel functions, as provided in Table 5. For the present study, four different GPR models such as rational quadratic, squared exponential, Matern 5/2, and exponential GPR models were developed. Each developed model shows excellent performance. The evaluation of various models demonstrates that the exponential GPR model exhibited superior performance achieving R², RMSE, RAE, MAE, NSE, RMSLE, and KGE values of 0.995, 0.0049, 0.0317, 0.0017, 0.9948, 0.0041, and 0.9874, respectively, during the training stage. While the corresponding values for the testing dataset were 0.985, 0.0085, 0.0767, 0.0042, 0.9822, 0.0075, and 0.9655, respectively. Compared with other modelling techniques, GPR is more efficient when dealing with small datasets. GPR models perform effectively in small datasets because they can handle nonlinearities and capture complicated correlations within data points. Figure 16 displays scatter plots of the model predicted K values and the empirically observed K values for the training and testing datasets, respectively.
Table 5

Performance evaluation parameter of the GPR model for the training and testing datasets

Statistical indicatorsRational quadratic GPR
Squared exponential GPR
Matern 5/2 GPR
Exponential GPR
TrainingTestingTrainingTestingTrainingTestingTrainingTesting
R2 0.970 0.954 0.974 0.951 0.979 0.962 0.995 0.985 
RMSE 0.0109 0.0144 0.0111 0.0147 0.0099 0.0127 0.0049 0.0085 
RAE 0.0943 0.1540 0.0976 0.1616 0.0773 0.1240 0.0317 0.0767 
MAE 0.0052 0.0084 0.0053 0.0088 0.0042 0.0067 0.0017 0.0042 
NSE 0.9747 0.9493 0.9742 0.9466 0.9792 0.9604 0.9948 0.9822 
RMSLE 0.0091 0.0124 0.0092 0.0127 0.0082 0.0109 0.0041 0.0075 
KGE 0.9784 0.9485 0.9782 0.9465 0.9805 0.9527 0.9874 0.9655 
Statistical indicatorsRational quadratic GPR
Squared exponential GPR
Matern 5/2 GPR
Exponential GPR
TrainingTestingTrainingTestingTrainingTestingTrainingTesting
R2 0.970 0.954 0.974 0.951 0.979 0.962 0.995 0.985 
RMSE 0.0109 0.0144 0.0111 0.0147 0.0099 0.0127 0.0049 0.0085 
RAE 0.0943 0.1540 0.0976 0.1616 0.0773 0.1240 0.0317 0.0767 
MAE 0.0052 0.0084 0.0053 0.0088 0.0042 0.0067 0.0017 0.0042 
NSE 0.9747 0.9493 0.9742 0.9466 0.9792 0.9604 0.9948 0.9822 
RMSLE 0.0091 0.0124 0.0092 0.0127 0.0082 0.0109 0.0041 0.0075 
KGE 0.9784 0.9485 0.9782 0.9465 0.9805 0.9527 0.9874 0.9655 
Figure 16

Scatter plot between the experimentally measured and model predicted K values for the training (left) and testing (right) stages. (a) Rational quadratic GPR, (b) squared exponential GPR, (c) Matern 5/2 GPR, and (d) exponential GPR model.

Figure 16

Scatter plot between the experimentally measured and model predicted K values for the training (left) and testing (right) stages. (a) Rational quadratic GPR, (b) squared exponential GPR, (c) Matern 5/2 GPR, and (d) exponential GPR model.

Close modal
Figure 17 presents the box plots of experimentally measured and GPR models predicted K values for the training and testing datasets. For the training dataset, all GPR models show strong alignment with the experimental values reflecting low variability and good agreement. For the testing dataset, the exponential GPR model has the longest whiskers compared with other GPR models reflecting the model's ability to effectively capture variations in the dataset while maintaining strong prediction efficiency. Similarly, violin plots in Figure 18 show the same trend with a narrow distribution for the training phase and wider tails for the testing phase.
Figure 17

Box plot of the experimentally measured and GPR model predicted K values for the training (left) and testing (right) stages.

Figure 17

Box plot of the experimentally measured and GPR model predicted K values for the training (left) and testing (right) stages.

Close modal
Figure 18

Violin plot of the experimentally measured and GPR model predicted K values for the training (left) and testing (right) stages.

Figure 18

Violin plot of the experimentally measured and GPR model predicted K values for the training (left) and testing (right) stages.

Close modal

Estimation of K using RF regression

The RF algorithm in this study was performed using the PYTHON programming language. The variables max_depth, min_sample_leaf, and n_estimators represent the maximum depth of the trees, the least number of samples required at a leaf node, and the total number of trees in the forest, respectively. The parameters found to be most suitable were max_depth of 10, n_estimators of 100, and mean_sample_leaf of 1. The developed model recorded R2, RMSE, RAE, MAE, NSE, RMSLE, and KGE values as 0.965, 0.0137, 0.1166, 0.0065, 0.9567, 0.0116 and 0.9265, respectively, for the training dataset. Whereas, for the testing dataset, the corresponding values are 0.94, 0.0197, 0.1962, 0.0107, 0.9060, 0.0168, and 0.8405, respectively. For both the training and testing datasets, the RMSE is relatively low as compared with other developed models, indicating the better predictive potency of the model. The selection of hyperparameters, such as the number of trees, tree depth, and features examined at each split, has a significant impact on RF performance. The scatter plots for the training and testing of the RF model are shown in Figure 19.
Figure 19

Scatter plot between the experimentally measured and model predicted K values for the training (left) and testing (right) stages for RF.

Figure 19

Scatter plot between the experimentally measured and model predicted K values for the training (left) and testing (right) stages for RF.

Close modal
Figure 20 presents the box plots for the training and testing stages of the RF model. For the training dataset, the RF model closely aligns with the experimental data demonstrating the model's ability to effectively capture the distribution of K values. Whereas, for the testing dataset, the RF model exhibits longer whiskers indicating a higher spread in predicted values compared with the experimentally observed values. Violin plots in Figure 21 further highlight this trend with a broader distribution in the training stage and more pronounced tails in the testing stage, suggesting increased variability and reduced predictive accuracy of the model for higher K values.
Figure 20

Box plot of the experimentally measured and RF model predicted K values for the training (left) and testing (right) stages.

Figure 20

Box plot of the experimentally measured and RF model predicted K values for the training (left) and testing (right) stages.

Close modal
Figure 21

Violin plot of the experimentally measured and RF model predicted K values for the training (left) and testing (right) stages.

Figure 21

Violin plot of the experimentally measured and RF model predicted K values for the training (left) and testing (right) stages.

Close modal

Sensitivity analysis

The sensitivity analysis was conducted on the best-performing models based on ANN, ANFIS, GPR, and RF techniques to evaluate the influence of various input parameters. Each developed model employs four input parameters, i.e., d10, n, d50, Cu, whose significance was assessed by systematically excluding them. From Table 6, it is clear that the most sensitive variables in model development are d10 and d50, as their removal led to a substantial decrease in R2 and an increase in RMSE values. In contrast, excluding n and Cu had a minimal impact on model performance, suggesting their lower sensitivity. These results highlight that d10 and d50 parameters have a substantial influence on the model performance, underscoring the importance of their careful selection and optimization.

Table 6

Sensitivity analysis of the developed models

VariablesANN
ANFIS
GPR
RF
R2RMSER2RMSER2RMSER2RMSE
d10, n, d50 0.880 0.0279 0.919 0.0137 0.962 0.0110 0.926 0.0219 
n, d50, Cu 0.852 0.0730 0.887 0.0557 0.921 0.0279 0.901 0.0347 
d10, d50, Cu 0.891 0.0316 0.925 0.0221 0.954 0.0103 0.914 0.0279 
n, d10, Cu 0.837 0.0654 0.853 0.0706 0.917 0.0215 0.871 0.0423 
d10, n, d50, Cu 0.902 0.0206 0.924 0.0179 0.98 0.0085 0.94 0.0197 
VariablesANN
ANFIS
GPR
RF
R2RMSER2RMSER2RMSER2RMSE
d10, n, d50 0.880 0.0279 0.919 0.0137 0.962 0.0110 0.926 0.0219 
n, d50, Cu 0.852 0.0730 0.887 0.0557 0.921 0.0279 0.901 0.0347 
d10, d50, Cu 0.891 0.0316 0.925 0.0221 0.954 0.0103 0.914 0.0279 
n, d10, Cu 0.837 0.0654 0.853 0.0706 0.917 0.0215 0.871 0.0423 
d10, n, d50, Cu 0.902 0.0206 0.924 0.0179 0.98 0.0085 0.94 0.0197 

By comparing the results of the developed models, the most effective model for determining the K of the porous media was identified. From Tables 35, it is concluded that the exponential GPR model performs better in predictions than the other models. Furthermore, three radar plots were generated to visually assess the performance of the developed models based on various statistical indicators. Figure 22 depicts the performance variation of the developed models across various statistical indicators utilized in this study. This plot provides a comprehensive overview of each model's prediction capability by displaying both error metrics and goodness-of-fit measures. Figure 23 provides a comparative visualization of developed models based on two key aspects: maximum performance indicators (R², KGE, and NSE) and error metrics (RMSE, RMSLE, MAE, and RAE). The left plot highlights the prediction efficiency of each model, where higher values of R2, KGE, and NSE indicate better alignment with the observed data. While the right plot illustrates the distribution of errors, with lower values signifying better predictive accuracy.
Figure 22

Radar plot illustrating variation in performance of the developed models based on the key statistical indicators.

Figure 22

Radar plot illustrating variation in performance of the developed models based on the key statistical indicators.

Close modal
Figure 23

Radar plot highlighting maximum performance indicators (R2, KGE, and NSE) (left) and variation in error metrics (RMSE, RMSLE, MAE, and RAE) (right) for the developed models.

Figure 23

Radar plot highlighting maximum performance indicators (R2, KGE, and NSE) (left) and variation in error metrics (RMSE, RMSLE, MAE, and RAE) (right) for the developed models.

Close modal
The variation between the experimentally observed and predicted K values derived from the developed models is shown in Figure 24. Taylor diagram (Figure 25) depicts the graphical framework for comparing the outcomes of the developed models with the experimental data. Notably, the GPR model is situated closest to the actual point, showing a higher correlation coefficient. The findings of the study demonstrate the advantage of using machine learning techniques such as GPR for precise estimation of K, especially when the data are complex and limited. These findings contribute to the development of more reliable models that can be utilized to improve groundwater and aquifer management, and support environmental sustainability initiatives.
Figure 24

Comparison of the experimentally observed and model predicted K values.

Figure 24

Comparison of the experimentally observed and model predicted K values.

Close modal
Figure 25

Taylor diagram showing different models developed for hydraulic conductivity prediction.

Figure 25

Taylor diagram showing different models developed for hydraulic conductivity prediction.

Close modal

In the present study, four distinct machine learning techniques, i.e., ANN, ANFIS, GPR, and RF, are used. Based on the correlation analysis, the parameters that were most significant in predicting the K of the porous media are Cu, n, d10, and d50. The predictive performance of the developed models was compared with measured K values using various statistical indicators (R2, RMSE, MAE, RAE, NSE, RMSLE, and KGE) and scatter plots. The performance of all the developed models was satisfactory, but GPR and RF models have higher prediction efficiency as compared with the ANN and ANFIS models. The exponential GPR model outperformed other models achieving R2, RMSE, RAE, MAE, NSE, RMSLE, and KGE values of 0.94, 0.0197, 0.1962, 0.0107, 0.9060, 0.0168, and 0.9655, respectively, during validation, making it most reliable method in the current study. The outcomes of the study highlight the potential of machine learning techniques in reducing uncertainties and overcoming traditional limitations in modelling soil hydraulic properties. The results and methods from the present study can be applied to address similar environmental problems. The developed models can be utilized for predicting the hydraulic conductivity of diverse media, supporting groundwater management, soil restoration, and agricultural planning. However, the relatively small dataset used in this study for modelling K may not adequately reflect the full range of porous media conditions. Furthermore, evaluating the machine learning techniques with a larger dataset could provide insights, and future studies may assess their performance in addressing other hydrological and environmental challenges.

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.

Abdalrahman
G.
,
Lai
S. H.
,
Kumar
P.
,
Ahmed
A. N.
,
Sherif
M.
,
Sefelnasr
A.
,
Chau
K. W.
&
Elshafie
A.
(
2022
)
Modeling the infiltration rate of wastewater infiltration basins considering water quality parameters using different artificial neural network techniques
,
Engineering Applications of Computational Fluid Mechanics
,
16
(
1
),
397
421
.
Ahmadi
S. M.
,
Saeed
B.
&
Abolfathi
S.
(
2024
)
Predicting the hydraulic response of critical transport infrastructures during extreme flood events
,
Engineering Applications of Artificial Intelligence
,
133
,
108573
.
Alagna
V.
,
Autovino
D.
,
Iovino
M.
&
Toscano
A.
(
2023
) ‘
Estimating the saturated soil hydraulic conductivity in a farm constructed wetland by the borehole permeameter infiltration method
’,
2023 IEEE International Workshop on Metrology for Agriculture and Forestry (MetroAgriFor)
.
IEEE
, pp.
596
600
.
ASTM
(
2007
)
Standard D422 – Particle-Size Analysis of Soils
.
West Conshohocken, PA, USA
:
ASTM International
.
Azad
A.
,
Karami
H.
,
Farzin
S.
,
Saeedian
A.
,
Kashi
H.
&
Sayyahi
F.
(
2018
)
Prediction of water quality parameters using ANFIS optimized by intelligence algorithms (case study: Gorganrood River)
,
KSCE Journal of Civil Engineering
,
22
(
7
),
2206
2213
.
Azarhoosh
M. J.
&
Koohmishi
M.
(
2023
)
Prediction of hydraulic conductivity of porous granular media by establishment of random forest algorithm
,
Construction and Building Materials
,
366
,
130065
.
Chandel
A.
,
Shankar
V.
&
Alam
M. A.
(
2021
)
Experimental investigations for assessing the influence of fly ash on the flow through porous media in Darcy regime
,
Water Science and Technology
,
83
(
5
),
1028
1038
.
Chandel
A.
&
Shankar
V.
(
2022
)
Evaluation of empirical relationships to estimate the hydraulic conductivity of borehole soil samples
,
ISH Journal of Hydraulic Engineering
,
28
(
4
),
368
377
.
Chandel
A.
,
Fayaz
F.
&
Shankar
V.
(
2023
)
Assessment of column to particle diameter ratio on the hydraulic conductivity of porous media: Wall effect in Darcy Regime
,
ISH Journal of Hydraulic Engineering
,
29
(
4
),
437
445
.
Chatrabgoun
O.
,
Karimi
R.
,
Daneshkhah
A.
,
Abolfathi
S.
,
Nouri
H.
&
Esmaeilbeigi
M.
(
2020
)
Copula-based probabilistic assessment of intensity and duration of cold episodes: a case study of Malayer vineyard region
,
Agricultural and Forest Meteorology
,
295
,
108150
.
Donnelly
J.
,
Abolfathi
S.
,
Pearson
J.
,
Chatrabgoun
O.
&
Daneshkhah
A.
(
2022
)
Gaussian process emulation of spatio-temporal outputs of a 2D inland flood model
,
Water Research
,
225
,
119100
.
Donnelly
J.
,
Daneshkhah
A.
&
Abolfathi
S.
(
2024
)
Physics-informed neural networks as surrogate models of hydrodynamic simulators
,
Science of the Total Environment
,
912
,
168814
.
Ewuzie
U.
,
Bolade
O. P.
&
Egbedina
A. O.
(
2022
)
Application of deep learning and machine learning methods in water quality modeling and prediction: a review
. In:
Current Trends and Advances in Computer-Aided Intelligent Environmental Data Engineering
, pp.
185
218
.
Habib
M. A.
,
O'Sullivan
J. J.
,
Abolfathi
S.
&
Salauddin
M.
(
2023
)
Enhanced wave overtopping simulation at vertical breakwaters using machine learning algorithms
,
PLoS One
,
18
(
8
),
e0289318
.
Habib
M. A.
,
Abolfathi
S.
,
O'Sullivan
J. J.
&
Salauddin
M.
(
2024
)
Efficient data-driven machine learning models for scour depth predictions at sloping sea defences
,
Frontiers in Built Environment
,
10
,
1343398
.
Jafari-Asl
J.
,
Hashemi Monfared
S. A.
&
Abolfathi
S.
(
2024
)
Reducing water conveyance footprint through an advanced optimization framework
,
Water
,
16
(
6
),
874
.
López-Acosta
N. P.
,
Espinosa-Santiago
A. L.
,
Pineda-Núñez
V. M.
,
Ossa
A.
,
Mendoza
M. J.
,
Ovando-Shelley
E.
&
Botero
E.
(
2019
)
Performance of a test embankment on very soft clayey soil improved with drain-to-drain vacuum preloading technology
,
Geotextiles and Geomembranes
,
47
(
5
),
618
631
.
Mahdian
M.
,
Noori
R.
,
Salamattalab
M. M.
,
Heggy
E.
,
Bateni
S. M.
,
Nohegar
A.
,
Hosseinzadeh
M.
,
Siadatmousavi
S. M.
,
Fadaei
M. R.
&
Abolfathi
S.
(
2024
)
Anzali wetland crisis: unraveling the decline of Iran's ecological gem
,
Journal of Geophysical Research: Atmospheres
,
129
(
4
),
e2023JD039538
.
More
S. B.
,
Deka
P. C.
,
Patil
A. P.
&
Naganna
S. R.
(
2022
)
Machine learning-based modeling of saturated hydraulic conductivity in soils of tropical semi-arid zone of India
,
Sādhanā
,
47
(
1
),
26
.
Noori
R.
,
Maghrebi
M.
,
Jessen
S.
,
Bateni
S. M.
,
Heggy
E.
,
Javadi
S.
,
Noury
M.
,
Pistre
S.
,
Abolfathi
S.
&
AghaKouchak
A.
(
2023
)
Decline in Iran's groundwater recharge
,
Nature Communications
,
14
(
1
),
6674
.
Richardson
R. R.
,
Osborne
M. A.
&
Howey
D. A.
(
2017
)
Gaussian process regression for forecasting battery state of health
,
Journal of Power Sources
,
357
,
209
219
.
Seyedian
S. M.
,
Haghiabi
A.
&
Parsaie
A.
(
2023
)
Reliable prediction of the discharge coefficient of triangular labyrinth weir based on soft computing techniques
,
Flow Measurement and Instrumentation
,
92
,
102403
.
Sihag
P.
,
Mohsenzadeh Karimi
S.
&
Angelaki
A.
(
2019
)
Random forest, M5P and regression analysis to estimate the field unsaturated hydraulic conductivity
,
Applied Water Science
,
9
(
5
),
129
.
Singh
B.
,
Sihag
P.
,
Pandhiani
S. M.
,
Debnath
S.
&
Gautam
S.
(
2021a
)
Estimation of permeability of soil using easy measured soil parameters: assessing the artificial intelligence-based models
,
ISH Journal of Hydraulic Engineering
,
27
,
38
48
.
Singh
B.
,
Sihag
P.
,
Parsaie
A.
&
Angelaki
A.
(
2021b
)
Comparative analysis of artificial intelligence techniques for the prediction of infiltration process
,
Geology, Ecology, and Landscapes
,
5
(
2
),
109
118
.
Thakur
D.
,
Chandel
A.
&
Shankar
V.
(
2022
)
Estimation of hydraulic conductivity of porous media using data-driven techniques
,
Water Practice & Technology
,
17
(
12
),
2625
2638
.
Thakur
A.
,
Chandel
A.
&
Shankar
V.
(
2025
)
Prediction of groundwater levels using a long short-term memory (LSTM) technique
,
Journal of Hydroinformatics
,
27
(
1
),
51
68
.
Tyralis
H.
&
Papacharalampous
G.
(
2024
)
A review of predictive uncertainty estimation with machine learning
,
Artificial Intelligence Review
,
57
(
4
),
94
.
Wang
J.
(
2023
)
An intuitive tutorial to Gaussian process regression
,
Computing in Science & Engineering
25
(
4
),
4
11
.
Yeganeh-Bakhtiary
A.
,
EyvazOghli
H.
,
Shabakhty
N.
&
Abolfathi
S.
(
2023
)
Machine learning prediction of wave characteristics: comparison between semi-empirical approaches and DT model
,
Ocean Engineering
,
286
,
115583
.
Zhang
Q.
,
Wu
Y. N.
&
Zhu
S. C.
(
2018
) '
Interpretable convolutional neural networks
’,
Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition
, pp.
8827
8836
.
This is an Open Access article distributed under the terms of the Creative Commons Attribution Licence (CC BY-NC 4.0), which permits copying, adaptation and redistribution for non-commercial purposes, provided the original work is properly cited (http://creativecommons.org/licenses/by-nc/4.0/).