Geotech Geol Eng
https://doi.org/10.1007/s10706-024-02767-8

ORIGINAL PAPER

Prediction of Stress‑Dependent Soil Water Retention Using
Machine Learning
Seyed Farid Fazel Mojtahedi · Ali Akbarpour
Ali Golaghaei Darzi · Hamed Sadeghi ·Martinus Theodorus van Genuchten·

Received: 24 August 2023 / Accepted: 13 February 2024
© The Author(s), under exclusive licence to Springer Nature Switzerland AG 2024

Abstract The soil water retention curve (SWRC)
provides information for a wide range of geoenvironmental problems, such as analyses of transient twophase flow, the bearing capacity and shear strengthof unsaturated soils. Many past studies have shown
experimentally the effects of stress on the SWRC.
Unfortunately, direct stress-dependent water retention
measurements are relatively time-consuming and generally require special equipment and a certain level ofexpertise. This study primarily aimed to develop anovel predictive framework within the context of soft
computing to capture the dependency of the SWRC
on several variables, with an emphasis on stress andsoil type. To achieve this, the three shape parametersof van Genuchten’s water retention model were estimated using a comprehensive database of 102 SWRC
tests retrieved from the literature. In this study, 60%
S. F. Fazel Mojtahedi
Department of Infrastructure Engineering, Universityof Melbourne, Parkville, Australia
A. Akbarpour
School of Civil Engineering and Environmental Science,
University of Oklahoma, Norman, USA
A. G. Darzi · H. Sadeghi (*)
Department of Civil Engineering, Sharif Universityof Technology, Tehran, Irane-mail: hsadeghi@sharif.edu
M. T. van Genuchten
Department of Nuclear Engineering, Federal Universityof Rio de Janeiro, Rio de Janeiro, Brazilof the datasets were employed for model training,with an additional 20% being designated for validation, while the remaining 20% were set aside for testing the model’s performance. The data were analyzed
using two machine learning techniques: the groupmethod of data handling and multi-layer perceptron
approaches. Results showed excellent performanceof the two methods. A sensitivity analysis was conducted to explore the relative significance of the different variables. Interestingly, net stress was found tobe almost as significant as soil type. The introduced
artificial intelligence based predictive framework provided a very effective method of integrating theory
and practice.
Keywords Soil water retention · Net stress ·Machine learning · Group method of data handling ·Multilayer perceptron
1 Introduction
The soil–water retention curve (SWRC) provides
key information about the hydraulic and mechanicalresponse of unsaturated soils to varying field conditions. Therefore, a sound understanding of retention
characteristics can result in a more reliable analysisof stability under transient two-phase flow conditions,among many other geotechnical problems (Garakani
et al. 2021). The SWRC expresses the relationshipbetween the volumetric (θ) or gravimetric (w) water
Vol.: (0123456789)

13

Geotech Geol Eng

content, or alternatively the degree of saturation (S),as a function of the capillary pressure or suction (ψ)(Genuchten 1980; Bazargan et al. 2015; Mahmoodabadi and Bryson 2021; Roy and Rajesh 2020; Vahdani
et al. 2022; Onyelowe et al. 2022). Previous research
has shown that several factors such as soil type, void
ratio, and net stress significantly affect the SWRC
(Zhang et al. 2019; Jiang et al. 2020; Sadeghi andDarzi 2021a; Dias et al. 2023). Any change in thesevariables, including net stress, results in variationsin the soil hydraulic properties. To characterize thestress-dependent SWRC, special experiments suchas unsaturated oedometer, stress-controlled pressure
plate, and unsaturated triaxial cell are required (Zhou
and Ng 2014; Mendes and Marinho 2020; Raheem
2023; Sadeghi and Darzi 2021b; 2023). Since theexperimental methodology is costly and time-consuming, a need exists to introduce alternative predictive frameworks incorporating basic soil properties.
Accordingly, many predictive and regressionbased models have been proposed to follow the mathematically sigmoidal form of SWRCs (Assouline andOr 2013; Wang et al. 2021; Zhao et al. 2023). The
first-order continuous three-parameter function of van
Genuchten (1980) has, for this reason, become popular for modeling the SWRC over the past four decades. The three model parameters (α, N, and M) arecharacteristics of the air entry suction, the pore size
distribution, and the overall symmetry of a SWRC,
respectively. Of particular interest is how the influence of net stress leads to a higher air entry suction
and a steeper desorption rate, resulting in a higher
retention capacity (Eldin and Senouci 1994; Mojtahedi et al. 2019b; Rezaee et al. 2021; Sadeghi et al.
2024). A critical issue is that even if the water retention parameters of a given soil under a specific net
stress are known, retention behavior at other stressstates cannot be predicted unless additional measurements are taken (Zhou and Ng 2014; Satyanaga et al.
2023). Hence, a need exists for a more flexible predictive framework that considers explicitly the effects
of net stress on the SWRC of different soils.
In recent years, soft computing has found ways tooptimize and solve various problems in soil, civil, andgeotechnical engineering (Behzadafshar et al. 2019;
Mojtahedi et al. 2019a; 2019b; Rezaee et al. 2021).
At the beginning of this era, machine learning techniques such as neural networks (NN) were preferredsince they did not require an a-priory predetermined
Vol:. (1234567890)

13

mathematical model while also tolerating specificmeasurement errors (Eldin and Senouci 1994; dosSantos Pereira et al. 2023). Other methods, such as
adaptive probabilistic neural networks (APNNs),
probabilistic neural networks (PNNs), and fuzzy
polynomial neural networks (FPNNs), were developed later to improve the performance, precision, andspeed of NNs. Nonetheless, applications of artificial
neural networks (ANNs) to a wide range of problems were found to be very beneficial. For instance,
ANNs have been used successfully to model properties related to the soil-atmosphere continuum such as
the soil hydraulic conductivity (Schaap et al. 1998),
soil water retention (Haghverdi et al. 2012; Khlosi
et al. 2016), various soil chemical soil properties(Kovačević et al. 2010), soil temperature (Feng et al.
2019), soil hydraulic parameters (Twarakavi et al.
2009), and the soil water holding capacity (Kaingo
et al. 2018).
ANNs have been used successfully to describe
SWRCs (Adeyemi et al. 2018; Bayat et al. 2013; Janaet al. 2008; Karandish and Šimůnek 2016; dos Santos Pereira et al. 2023). However, few, if any, previous models based on soft computing considered thedependency of SWRC on net stress, even though thesignificant role of net stress has been known for morethan 20 years. In addition, recent advances in ANNs,
such as the group method of data handling (GMDH),
have enhanced the accuracy of existing methods andhelped optimize network configurations. They canalso be used effectively to identify data interrelations
and improve sensitivity analyses (Mulashani et al.
2022). In GMDH algorithms, the collection of neurons in each layer is described by means of a quadratic polynomial, which produces new neurons in a
subsequent layer. GMDH techniques help to resolve
the modeling difficulty of single outputs of multiple
data (Ivakhnenko et al. 1999). Moreover, GMDH, asan algorithm modeling method, can be used to represent nonlinear relationships between the input andoutput variables (Menad et al. 2019).
The purpose of our study is to formulate a new
predictive framework for the SWRC within the context of machine learning. For the predictions, weuse an extensive database of water retention data asreported by several research groups worldwide. To
our knowledge, no published studies exist in which
the capabilities of input–output mapping of MLP and
GMDH based NN techniques have been explored forGeotech Geol Eng

modeling the stress-dependent SWRC of different
soils.
2 Methods
2.1 Neural Networks
In order to learn the dominant data classification
patterns and process the information, neural networks (NN) have been used extensively as one of themachine learning techniques. NNs can be utilized to
predict many tasks, such as statistical modeling andresolving mathematical problems using various datatypes. The inspiration for a NN comes from humanbiological neural systems. One of the main components of NN is neurons, which are the processor
members that make the connection between components feasible by the execution of specified mathematical functions. Similar to a human neural network,neural networks conduct training through the use of
initial instances. Therefore, the governing pattern canbe determined through the training phase (Behzadafshar et al. 2019; Mojtahedi et al. 2019a).
2.1.1 Multi‑layer Perceptron Neural Network(MLP‑NN)The multi-layer perceptron (MLP) presented firstby Hornik et al. (1989) is one of the most reputableneural networks. An MLP is a structure made up of
several hidden layers of neurons in which the output
from one layer of neurons becomes the input of thenext layer’s neurons. MLP is comprised of three layers covering output, hidden, and input. Each layer
contains several neurons linked to each other using
layers by weights. Because each unit (or neuron) of alayer is linked to all units of the next layer but not to
the neurons of the same layer, a multi-layer perceptron (MLP) has a basic architecture. The number ofneurons in each layer and the number of hidden layers are the parameters of architecture. For example,
to approximate several functions in high-dimensional
spaces, MLP results have been found to perform
effectively in terms of accuracy and error (Geva andSitte 1992). MLPs propagate input data in a forward
direction. Randomly selected weight values are thenused between the input and hidden neurons and also
between the hidden and output neurons. Multiplying

each of the precedent layer’s input signals (xi) by aconstant adaptive weight (wij) produces the weighted
signals, which will feed the hidden neurons. The
weighted signals then go into a summation function.
To produce the final output of the system, this process
is repeated for each layer. According to MLP definitions, the output value is as follows:

nethj =

n∑(5)

(wij × xi + bj )i=1

where nethj is the output of the network, wij is theweight between each input neuron and output neuron, xi is the input signal, and bi is the bias for eachnode. Various node transfer functions can be used
depending on the application. Moreover, based on
the complexity of the problem, an activation function is selected. A sigmoid function such as a tangent sigmoid or log-sigmoid can be used for mostnonlinear problems (Hosseini et al. 2020; Mojtahedi
et al. 2023). The same training procedure as network training occurs, which can be executed using
different algorithms. The Levenberg–Marquardt(LM) approach is one of the most common optimization algorithms known for reducing the squareerror between the target variable and the output. The
convergence speed of this variable is relatively high(Jahangir and Rezazadeh Eidgahee 2021).
2.1.2 Group Method of Data Handling Neural
Network (GMDH‑NN)
The group method of data handling (GMDH) is basedon a self-organizing inductive approach to estimating black-box models with unknown relationships
between the variables (Vissikirsky et al. 2005). Ivakhnenko et al. (1999) developed this method as one of
the first GMDH approaches to design nonlinear relationships systematically. The optimal structure of themodel within the GMDH algorithm is obtained fromsuccessive generations of partial descriptions of thedata.
In the GMDH algorithms, simple polynomials are
generated at the beginning; these polynomials allow
a rough estimation of the aimed systems to be generated, but as they become more complex, more accurate models can be generated. Polynomials typicallytake the Kolmogorov–Gabor form, which is a discrete
version of the Volterra functional (Ivakhnenko 1971):
Vol.: (0123456789)

13

Geotech Geol Eng

𝜙 = a0 +m∑
i=1

ai xi +mm ∑∑i=1 j=1

aij xi xj +

m ∑
m
m ∑
∑
i=1 j=1 k=1

aijk xi xj xk + …

(2)

To generate a comprehensive system, multiple
layers of partial descriptions are built into a network
based on the multi-layer perceptron (MLP), which
has the characteristics of a higher-order function(Ivakhnenko 1971):)(
𝜙 = f x1 , x2 , x3 , … , xn
(3)

The regular style of a trained GMDH network anda fully connected MLP are manifested in Fig. 1. In
both models, input values are represented by squared
nodes, while the circular nodes of the output and hidden layers are neurons or processing elements (PE).
The most prevalent MLP type is characterized bynodes in each layer being connected to all nodes ofthe previous layer. On the other hand, in the mostcommon type of GMDH, each node contains one
partial description, with the outputs of two different nodes of the previous layer forming its input.
Since all neurons in the GMDH network are trained

(a)

(b)

to estimate the same function, which is the desired
output of the network, they should have one single
output; however, the MLP-NN could have more than
one output. If the same inputs are required for several
outputs, separate GMDH-NNs must be generated foreach output.
Using a training algorithm based on the number of
inputs, the number of partial descriptions that form
the neurons of the first layer is determined. The number of partial descriptions is denoted using the number of the n inputs of the network taken two by two,as represented in Eq. 4. Two inputs are taken by eachpartial description, which delivers one output (Ivakhnenko 1971).

C2n =n2 − nn!
=
2!(n − 2)!
2

Partial descriptions (neurons) of the first layer will
be m = C2n, where n is the number of inputs of the
network:()
y11 = f11 x1 , x2()
y12 = f12 x1 , x3
(5)
⋮⋮
(
)
y1m = f1m xn−1 , x1

In Eq. 5, f1k is the partial description of neuron k inlayer 1, and y1k represents the output of that neuron,
where k = [1, 2, …, m].
After training the first layer, some neurons are
removed according to a selection criterion so that
only those neurons that meet that criterion are
retained. As manifested by Eqs. (6) and (7), the collection of chosen m̂ partial descriptions is a subset of
the total number of partial descriptions:
]
] [
[
�
y11 , �y12 , �y13 , … , �y1�m ⊂ y11 , y12 , y13 , … , y1m
(6)

[
] [
]
�f11 , �f12 , �f13 , … , �f1�m ⊂ f11 , f12 , f13 , … , f1m

Fig. 1  Comparison between a a typical GMDH neural network and b a fully connected multi-layer perceptron (MLP)
Vol:. (1234567890)

13

(4)

(7)

where m̂ < m.
A second layer is added after selecting the mosteffective nodes in the first layer. This new layer has
p = C2m̂ neurons, where m̂ is the number of nodes
in the first layer following the selection procedureaccording to Eq. (8):Geotech Geol Eng

Fig. 2  Graphical examples of the training process for the group method of data handling (GMDH)()
y21 = f21 ŷ 11 , ŷ 12()
y22 = f22 ŷ 11 , ŷ 13⋮⋮(

y2p = f2p ŷ 1(m−1), ŷ 1m̂
̂)(8)

The method of selecting neurons and adding layers
continues until one of the stopping rules is fulfilled, atwhich point the procedure is terminated.
The formation and training processes of the GMDH
Polynomial Neural Networks are depicted in Fig. 2 andsummarized in the following stepwise procedure:
1. Start by creating the first layer of the network.
Every node serves as a partial description. Figure 2a indicates a layer created with four inputs;
2. Compute the coefficients of the nodes;
3. Evaluate the performance and eliminate thoseneurons that provide the poorest results. The
lighter color nodes in Fig. 2a show the removed
neurons in Fig. 2a;4. Add one more layer that uses the output of theprevious layer as input (Fig. 2b);5. Select and train the nodes of the new layer. Add
another layer after selection (Fig. 2c);
6. Stop the training process if any layer has only one
node after selection. (Fig. 2d);
7. Train until the best performing value is lower
than that of the previous layer, except for layerone. When this occurs, the algorithm chooses oneof the best nodes from the remaining nodes in thelast layer and discards the others;
8. Remove all nodes that are not contributing tothe output of the network in all previous layers(Fig. 2f).
3 Case Study
3.1 Model Performance
The following statistical indicators were used toevaluate the performance of models:
1. Root mean squared error (RMSE):

Vol.: (0123456789)

13

Geotech Geol Eng

√
√
Mt
√ 1 ∑
(h − t )2
RMSE = √Mt i=1 i i(9)

2. Mean of absolute percent error (MAPE):

MAPE =∑Mt

�hi − ti �
× 100
∑Mt
h
i=1 i

i=1

(10)

3. Coefficient of determination (R2):
2
��
⎤
⎡
�
∑Mt �
⎥
⎢
hi − hi ti − ti
i=1⎥
R2 = ⎢ �
�2 ∑ �
⎥
⎢ ∑ �
�
2
Mt
Mt
⎢
hi − hiti − t i ⎥
i=1
i=1
⎦
⎣

(11)

where Mt is the total number of data, ti and hi are predicted, and actual values of the output, respectively,
and t̅i and h̅i are the averages of the predicted andactual output, respectively. We note that lower values of MAPE and RMSE obtained for a model indicate better performance of the model. The RMSE isa measure representing the difference between modeled and measured data, but because it is squaredin the formulation, it significantly penalizes outliers. The MAPE indicates the average or mean of theabsolute percentage errors of the forecast. Error is thedifference between the observed or actual value andthe predicted value. This measure is easy to comprehend because the error is expressed in percentages.
The problem associated with canceling negative andpositive errors is also avoided by using absolute percentage errors. In fact, for a precise prediction model,
the expected values of RMSE and MAPE are zero.
On the other hand, R2 indicates the total variance inthe observed data that the model explains; hence, theideal value of R2 is unity.
3.2 Database
The selected database contains 102 test results
obtained from 27 publications (Albadri et al. 2021;
Elkady et al. 2015; Garakani et al. 2015; Haeri et al.
2017; 2014a, b; Karube and Kawai 2001; Jafarzadeh et al. 2021; Khosravi et al. 2016a; 2016b; Kimet al. 2010; Kit 2007; Lee et al. 2005; Leung and Ng
2016; Li et al. 2007; Ng et al. 2016; 2020; Ng andPang 2000; Oh and Lu 2014; Sadeghi et al. 2016;
Vol:. (1234567890)

13Song et al. 2020; Sun et al. 2007; Tavakoli Dastjerdi
et al. 2014; Uchaipichat and Khalili 2009; Vanapalli et al. 2015; Sadeghi and Nasiri 2021; Yao et al.
2021; Ahmadinezhad et al. 2019; Zhang et al. 2020).
Details of the retrieved water retention data are summarized in Table 1.
Input parameters comprised the optimum watercontent (wopt), the maximum dry density (γdmax),
plastic limit (PL), liquid limit (LL), specific gravity (Gs), net stress (σ), and soil type. The output values are the three shape parameters (α, N, and M) in
van Genuchten’s soil water retention model given by(Genuchten 1980):[]M
Sr = 1 + (𝛼𝜓)N(12)where Sr is relative saturation and ψ the soil suction(or capillary pressure, here given in kPa). The parameters α, N, and M represent the air entry suction (ψa,approximately the inverse of α), the pore size distribution, and the overall symmetry of the water retention curve, respectively. We note here that the threeparameter form of van Genuchten’s SWRC was usedin this study, although the two-parameter form (whereM = 1–1/N) could be an alternative for fitting theexperimental data. An advantage of the three-parameter form is that the downward and upward curvaturesof the retention curve are independent of each other.
It’s important to mention that this study utilizes fitting analysis on the collected datasets to determine
the actual values of parameters a, n, and m in the VG
model. As mentioned earlier, the dataset includes
102 experimental data for the SWRC drying curve.
The general form of the SWRC is shown in Fig. 3.
In accordance with the collected database evidencing the dependency of SWRC on net stress, this figure
illustrates the expected change in the SWRC due toan increase in net stress. An increase in the net stressleads to a lower average pore volume and a more
uniform pore size distribution (Sadeghi et al. 2015).
Consequently, both the air entry suction (ψa) and theslope of the curve will increase.
Statistical properties of inputs and outputs are
summarized in Tables 2 and 3. Furthermore, Fig. 4
shows the variety and number of soil types weused as input in this study. It should be mentionedthat the net stress is a tensor, but for the data sets
used in this paper, this refers to the normal component only.

Geotech Geol Eng
Table 1  Details of theselected dataset

Data Source

No. of data

Type of soil

Range ofnet stress(kPa)

Ng and Yung (2008)
Haeri et al. (2014b)
Kit (2007)
Ng and Pang (2000)
Ng et al. (2020)
Oh and Lu (2014)
Khosravi et al. (2016b)
Li et al. (2007)
Tavakoli et al. (2014)
Zhang et al. (2020)
Ng et al. (2016)
Karube and Kawai (2001)
Garakani et al. (2015)
Yao et al. (2021)
Uchaipichat and Khalili (2009)
Albadri et al. (2021)
Elkady et al. (2015)
Tavakoli and Habibagahi (2013)
Khosravi et al. (2016a)
Khosravi et al. (2016b)
Lee et al. (2005)
Kim et al. (2010)
Haeri et al. (2017)
Vanapalli et al. (1999)
Haeri et al. (2014a)
Leung and Ng (2016)
Sun et al. (2007)
Total Number of Data:

2
4
3
3
3
4
3
2
3
3
3
2
3
16
4
3
3
3
3
3
4
5
4
3
4
3
5
102

CL
ML
SP
ML
CL
SM
ML
CH
CL—ML
ML
CL
ML
ML
SC
CL—ML
SP
CH
CL – ML
ML
SP
SM
SP-SM
ML
CL
ML
CL
CL

110—300
50—200
0—80
0—80
0—80
0—200
125- 225
10—4000 – 250
100—400
0—150
19.6—98
50 – 200
0—90
50—200
0—200
7—600
0—250
40—90050—2000—300
0 – 100
50 -200
50—200
150—400
0—80
49—588Degree of saturation (Sr)

3.3 Data Pre‑processing

Zero net stress SWRC
Increasing net stress SWRC

Soil Suction ( )Fig. 3  Schematic representation of the effect of stress on the
SWRC​

One of the problems commonly associated with large
datasets is the existence of missing values. Missing
values in a dataset can create biased estimates anddecrease the statistical power of a model, which leads
to inapplicable conclusions. Therefore, there is a needto impute missing values before training the model.
Unfortunately, using inaccurate methods for missing data imputation may engender several-fold larger
errors in the model predictions (Song et al. 2020).
Because the datasets for this study were derived fromprevious publications, missing values were apparent for several input parameters. Missing values forGs (specific gravity), PL (plastic limit), LL (liquid
Vol.: (0123456789)

13

Geotech Geol Eng
Table 2  Statistical analysisof data used (inputs)

Variable

Abbreviation

Min

Max

Mean

Sd

Sample size

Optimum water content (%)
Maximum dry density (g/m3)
Liquid limit (%)
Plastic limit (%)
Specific gravity
Net stress (kPa)

wopt
γdmax
LL
PL
Gs
-

3.6
12
14.0
4.0
2.60
0

26
18.94
137.0
60.0
2.75
900

13.07
17.33
41.3
23.0
2.68
129

5.10
1.89
22.2
10.0
0.03
152

90
91
82
82
84
102

Table 3  Statistical analysis of data used (outputs)
VG water retention
shape factors

Min

Max

Mean

Sd

α
N
M

0.15
0.2
0.006

3.29
104
11.24

0.26
7.5
0.51

0.52
16
1.44

Fig. 4  Statistical report of the different soil types employed inthis study based on USCS

limit), γdmax (Maximum dry density), and wopt (optimum water content) were 14.2, 20.9, 20.9, 10.4, and14.2%, respectively. In order to impute missing values, the following three approaches were exploredand evaluated:(1) Mean imputation Mean imputation is the usual
method of imputing missing data. In this method,
the missing value of a specific input parameteris replaced by the mean value of all available
values of that parameter. Mean substitution presumes that the mean is an acceptable estimate foran irregularly picked observation drawn from a
Vol:. (1234567890)

13

normal distribution. While popular, this imputation approach may result in conflicting bias when
missing data are not entirely random, especiallywhen there is a significant disparity in the number of missing values for the different variables.
Mean imputation also provides that no new information is obtained but only increases the sample
size, resulting in an underestimation of the errors
involved (Kang 2013).
(2) K-nearest neighbor imputation (KNN) KNN isa part of the data mining classification processthat finds patterns depending on how far away
the parameters are from their closest neighbors.
Based on the value of Euclidean distance, KNN
calculates the distance to the nearest neighbor(Sallaby and Azlan 2021). Recent research suggests that KNN performs better than hot-nearest
neighbor imputation and Naïve Bayes imputation
algorithms (Fadillah and Muchlisoh 2019).
(3) DataWig This method tries to impute missing
data in tables by building machine learning models. DataWig supports the imputation of both
numerical and categorical columns. To determine
the best fit for missing data, DataWig combines
automatic hyperparameter optimization withdeep learning feature extractors (Biessmann et al.
2018).
Figure 5 shows root mean squared error (RMSE)
values calculated for the three imputation algorithms(i.e., mean imputation, k-nearest neighbor imputation, and DataWig). Results indicate that the error
rates of mean imputation are higher than those of theother two algorithms. Mean imputation distorts the
distribution and underestimates the standard deviation, mostly by ignoring correlations between theattributes of the data set. As shown in Fig. 5, KNN
imputation produced lower RMSE values than mean
imputation, but they were still relatively high. Clearly,Geotech Geol Eng

Fig. 5  Root mean squared errors (RMSEs) of Mean Imputation, K-Nearest Neighbor Imputation, and DataWig for different parameters. a Gs, b γdmax, c PL, d LL, and e wopt

DataWig yielded the best results for the examineddata compared with the other two algorithms. The
performance of all approaches became worse withan increase in the missing rate. This is due to the factthat more significant missing rates are associated witha smaller amount of data available for imputation.
Vol.: (0123456789)

13

Geotech Geol Eng

4 Results and Discussion
4.1 GMDH ModelingAs noted earlier, the GMDH methodology has several layers. In the beginning, input data go throughthe first layer. After processing, output from thisfirst layer becomes input for the second layer, etc.
This process continues until the error of the resultsof layer n + 1 is lower than that of layer n, at whichtime the calculations terminate. The iterative natureof this process ensures progressive model refinement,balancing complexity and accuracy to prevent overfitting, a notable consideration in machine learningmodels.
The soil water retention data in this paper wereanalyzed using van Genuchten’s three-parametermodel with α, N, and M as the unknowns. Two steps
play a pivotal role during the optimization process.
First, a suitable performance index should be selected
to evaluate the performance of the artificial intelligence algorithm as applied to the selected model (inthis case, Eq. 12). In our research, the performanceof the algorithm was evaluated in terms of the mean
absolute percent error (MAPE). MAPE was chosen due to its precise representation of prediction
accuracy, which is essential for handling diverse
geotechnical datasets. Second, to obtain an optimal
form of the algorithm, its control parameters should
be defined in the most effective way. The maximum number of layers (MNL), the maximum number of neurons in a layer (MNNL), and the selection pressure (SP) are the control parameters thatplay a crucial role in finding the most suitable solution, including proper convergence of the algorithm.
The selection pressure, a dimensionless parameter
between 0 and 1, plays a key role in the sensitivity ofthe modeling errors. Values of SP close to zero willlead to more output. If this parameter is set to zero,
all generated answers are transferred to the next layer,
which is unacceptable. The strategic choice of control
parameters, particularly the SP value of 0.6, reflects
a balance between model complexity and information
redundancy. In general, no single unique equation
exists to calculate control parameters; in most cases,
they are determined based on expert opinion andtrial and error. Based on the previous studies (Guido
et al. 2020; Morosini et al. 2020), a value of 0.6 was
Vol:. (1234567890)

13

selected for SP, which implies that only 0.4 of the
produced answers are transferred to the next layer.
Values of MNNL and MNL were determined using
trial and error. For this purpose, we considered values of 5, 10, 20, 30, and 50 for MNNL and 5, 10, 20,
and 30 for MNL. Finally, 60% of the data were considered training data, 20% were used as validation
data, and the remaining 20% as test data. According
to the control parameters, 20 models generated withthe calculated results are reported in Tables 4, 5, 6 forthe inverse of the VG shape parameters α, N, and M.
Since model 18 produced the lowest MAPE, this formulation was selected as the best-developed model.
The selection of Model 18, with its low MAPE, highlights the effectiveness of the optimization process inachieving accurate predictions for practical geotechnical applications. A flowchart of the GMDH-NN
model is presented in Fig. 6. This visual representation aids in comprehending the complex structure anditerative development process of the model. The flowchart of the GMDH-NN model is presented in Fig. 6.
4.2 Comparison of Results
Comparisons between the actual and predicted valuesof the VG parameters α, N, and M are shown in Fig. 7
for the training data set, in Fig. 8 for the validation
data, and in Fig. 9 for the testing data. Results indicate that the developed models are very accurate. The
close alignment between actual and predicted values
in these figures underscores the models’ robustness
in capturing the complex relationships in soil–water
retention data. During the training process, the
GMDH model automatically learns the factors that
dominate the structure variables. The optimum neuron structure is created automatically by removing
redundant neurons from the network and minimizing
values of the prediction error criteria. This optimization is crucial for enhancing the model’s predictive
accuracy and efficiency, as seen in the progressively
lower error rates in the iterative training process.
Table 7 shows the results of the performance indices for all developed models. The RMSE and MAPE
results for the training, validation, and testing data
sets show that GMDH performed slightly better than
MLP. This performance superiority of GMDH is
attributed to its ability to capture and model the nonlinear characteristics inherent in soil–water retention
data more effectively. This is reflected in Figs. 10, 11,

Geotech Geol Eng
Table 4  The effect ofdifferent values of twocontrol parameters on theaccuracy of predicting α
in training, validation, andtesting datasets

Table 5  The effect ofdifferent values of twocontrol parameters on theaccuracy of predicting n
in training, validation, andtesting datasetsmodel no

SP

MNL

MNNL

MAPE oftraining

MAPE ofvalidation

MAPE of testing

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

0.6
0.6
0.6
0.6
0.6
0.6
0.6
0.6
0.6
0.6
0.6
0.6
0.6
0.6
0.6
0.6
0.6
0.6
0.6
0.6

5
5
5
5
5
10
10
10
10
10
20
20
20
20
20
30
30
30
30
30

5
10
20
30
50
5
10
20
30
50
5
10
20
30
50
5
10
20
30
50

9.24
7.74
5.94
11.34
4.44
13.44
4.64
8.64
5.84
12.14
10.84
11.54
12.64
10.04
3.94
4.44
13.04
3.44
7.58
7.14

10.73
9.35
6.98
12.26
5.12
14.09
5.69
9.36
7.42
13.21
12.42
12.54
13.54
10.92
5.12
5.51
14.59
4.09
9.01
8.52

14.12
12.82
10.19
15.29
8.51
17.46
8.84
12.39
10.94
16.35
15.84
15.61
16.65
13.98
8.89
8.63
17.99
7.05
12.35
11.83

Model no

SP

MNL

MNNL

MAPE oftraining

MAPE ofvalidation

MAPE of testing

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

0.6
0.6
0.6
0.6
0.6
0.6
0.6
0.6
0.6
0.6
0.6
0.6
0.6
0.6
0.6
0.6
0.6
0.6
0.6
0.6

5
5
5
5
5
10
10
10
10
10
20
20
20
20
20
30
30
30
30
30

5
10
20
30
50
5
10
20
30
50
5
10
20
30
50
5
10
20
30
50

10.27
8.32
6.23
12.32
4.65
14.55
5.26
9.62
6.78
13.11
11.65
11.76
13.42
10.49
4.48
5.02
13.68
3.65
7.81
7.98

15.20
12.29
9.16
18.29
6.79
21.63
7.70
14.24
9.97
19.47
17.28
17.45
19.94
15.54
6.52
7.34
20.33
5.28
11.52
11.77

17.28
16.06
13.62
19.81
12.27
21.03
12.89
16.84
13.88
20.57
19.23
19.04
20.64
17.85
11.93
12.33
20.83
10.93
15.57
15.31
Vol.: (0123456789)

13

Geotech Geol Eng
Table 6  The effect ofdifferent values of twocontrol parameters on theaccuracy of predicting m
in training, validation, andtesting datasets

Model no

SP

MNL

MNNL

MAPE oftraining

MAPE ofvalidation

MAPE of testing

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

0.6
0.6
0.6
0.6
0.6
0.6
0.6
0.6
0.6
0.6
0.6
0.6
0.6
0.6
0.6
0.6
0.6
0.6
0.6
0.6

5
5
5
5
5
10
10
10
10
10
20
20
20
20
20
30
30
30
30
30

5
10
20
30
50
5
10
20
30
50
5
10
20
30
50
5
10
20
30
50

8.16
6.80
5.16
10.09
3.59
11.69
4.03
7.60
5.04
10.79
9.60
10.27
11.20
8.89
3.07
3.87
11.59
2.94
6.67
6.28

15.44
13.79
11.78
17.81
9.86
19.77
10.40
14.76
11.64
18.66
17.21
18.03
19.17
16.34
9.22
10.20
19.64
9.07
13.63
13.16

15.70
13.89
11.71
18.29
9.62
20.42
10.20
14.96
11.55
19.22
17.63
18.52
19.77
16.68
8.92
9.98
20.28
8.75
13.72
13.21

and 12 by the higher R2 values obtained with GMDH
as compared to MLP. For instance, R2 values of thetest data set using the GMDH model for α, N, and M
were 0.9961, 0.972, and 0.9899, respectively, slightly
higher than those for the MLP model (0,9950, 0.9725,
and 0.950, respectively). These values, especially inthe testing dataset, are indicative of the model’s generalizability and predictive strength. Hence, we conclude that the proposed GMDH model provides moreaccurate results than MLP; however, both approaches
performed satisfactorily in terms of predicting stressdependent SWRCs. The satisfactory performanceof both models highlights the potential of machinelearning techniques in enhancing predictive capabilities in geotechnical engineering, opening avenues formore sophisticated and nuanced soil analysis.
4.3 Further Evaluation of the Models
To evaluate the effectiveness of the proposed machinelearning models, a blind dataset consisting of 18 data
points drawn from four separate studies was utilized(Cai et al. 2023; Jadar et al. 2023; Ng and Pang 2000;
Seth et al. 2023). Details of the used blind data set
Vol:. (1234567890)

13

can be found in Table 8. Figure 13 illustrates theresults obtained from this analysis. The MAPE values for parameters α (a), m (b), and n (c) consistently
demonstrated lower errors when the GMDH approach
was applied, as opposed to the MLP. Specifically, the
GMDH method outperformed the MLP across most
test instances, indicating its superior accuracy in predicting the given outcomes. This trend was particularly pronounced in the case of parameter α, where
the GMDH model exhibited a significantly reduced
MAPE. Additionally, it is noteworthy that the MAPE
values for α at test instances 4, 5, and 6 exhibit a significant discrepancy under the MLP method. This
divergence is attributed to the α values for these
instances falling outside the scope of the range
employed during the training phase.
In the presented research, an approach has been
introduced in the field of geotechnical engineering by employing machine learning techniques,
specifically the group method of data handling(GMDH) and multi-layer perceptron (MLP), formodeling the stress-dependent soil water retention
curve (SWRC). Based on the results in this section, this approach is recognized as innovative dueGeotech Geol Eng

Fig. 6  Flowchart of proposed GMDH-NN

to its provision of a subtle understanding of how
different stress levels impact soil–water retention, a
critical aspect in geotechnical analysis and design.
For example, the model can predict how the SWRC
of clay soil changes under increased compaction
stress, which is vital for understanding its behavior in foundation engineering. The adaptability ofthe model to various soil types is a distinguishingfeature. A wide range of soil physical properties
has been incorporated into the machine learningmodels, resulting in a framework that is both versatile and highly applicable in real-world scenarios
where soil types vary significantly. However, it has
been noted that the performance of the model isconfined to the specific range of input and output
parameters represented in the dataset. As shownin Fig. 13, while robust predictions are demonstrated within this range, accuracy may diminish
when applied to soil types or stress conditions thatdeviate significantly from those included in the
study. This limitation underscores the necessity for

expanding the range of input and output parametersin future research, thereby enhancing the model’s
generalizability and applicability to a broader spectrum of geotechnical engineering challenges.
4.4 Discussion About the Wetting SWRC​
The selected database for wetting path contains 49
test results obtained from 13 publications (Ng andPang 2000; Tavakoli and Habibagahi 2013; Karube
and Kawai 2001; Tse 2007; Oh and Lu 2014; Ng
et al. 2016; 2020; Leung and Ng 2016; Khosravi et al.
2016b; Zhang et al. 2020; Albadri et al. 2021; Sethet al. 2023; Cai et al. 2023), as shown in Table 9. Dueto the limited data for the wetting path, the models
were impeded in their learning and validation capabilities, which necessitated the prioritization of thedrying path where the dataset size was in alignment
with the recognized standards for a ’decent dataset’(Vabalas et al. 2019). This decision was corroborated
Vol.: (0123456789)

13

Geotech Geol Eng

Fig. 7  Comparison between measured and predicted values of
the training data set for the VG model parameters of a α, b N,
and c M

Fig. 8  Comparison between measured and predicted values ofthe validation data set for the VG model parameters of a α, bN, and c M

by attempts to model the wetting path, which did notmeet the established threshold for reliability.
Furthermore, Table 10 highlights that both MLP
and GMDH models’ error metrics—RMSE and
MAPE—are quite high, especially in the validationand test phases. This could be due to the relatively
small dataset size, which may not be sufficiently representative to allow the models to generalize well. It’s
a common challenge in machine learning that limited
data can lead to models that perform well on training

Vol:. (1234567890)

13

Geotech Geol Eng
Table 7  Obtained results of RMSE and MAPE for the constructed MLP and GMDH models
Metrics (model) VG parameter Training Validation Test
RMSE (MLP)
MAPE (MLP)
RMSE(GMDH)
MAPE(GMDH)RMSE (MLP)
MAPE (MLP)
RMSE(GMDH)
MAPE(GMDH)RMSE (MLP)
MAPE (MLP)
RMSE(GMDH)
MAPE(GMDH)

Fig. 9  Comparison between measured and predicted values of
the testing data set for the VG model parameters of a α, b N,
and c M

data but fail to predict new, unseen data accurately.
The noticeable high value in error emphasizes theneed for a larger dataset to improve model reliability
and performance.

α

n

m

0.073
8.75
0.009

0.038
9.13
0.023

0.027
11.30
0.023

3.44

4.09

7.05

1.630
7.94
1.220

2.453
10.76
0.399

1.581
13.60
1.385

3.65

5.28

10.93

0.120
8.53
0.068

0.162
13.46
0.086

0.048
16.19
0.021

2.94

9.07

8.75

In addition, Figs. 14, 15 and 16, depict the relationship between actual and estimated values forthe parameters α, n, and m, across training, validation, and testing datasets for the wetting path using
MLP and GMDH models. These scatter plots withfitted lines show a consistent trend: the R-squared
values are relatively low across all datasets andparameters, indicating a moderate to weak predictive relationship. Notably, the R-squared values
decrease from training to testing datasets, which
can be a sign of overfitting during the training
phase due to the small dataset size. This is further
evidenced by the scatter of points around the fitted line being broader in the test dataset figures,
suggesting higher variability and lower prediction
accuracy when the models encounter new data.
The consistency of this pattern for both MLP and
GMDH models across all parameters highlights thechallenge of limited data in developing robust predictive models for the wetting path.
The small size of the dataset posed significant
challenges for effectively training and validatingmachine learning models for the wetting path. It is
well-documented that the correlation between dataset size and a model’s ability to generalize is strong,Vol.: (0123456789)

13

Geotech Geol Eng

Fig. 10  Predicted versus measured values for the training data
set. a α-MLP, b α-GMDH, c N-MLP, d N-GMDH, e M-MLP,
f M-GMDH. The black solid line depicts the linear relationshipbetween predicted and observed values. The dashed red line isa one-to-one line and provides a benchmark for comparison

meaning that models trained on insufficient data
tend to perform poorly when used on new examples(Uddin and Lu 2024).

Therefore, it is suggested that further research be
conducted to explore the potential benefits of thewetting path. Efforts to increase the dataset for the

Vol:. (1234567890)

13

Geotech Geol Eng

Fig. 11  Actual versus estimated values for the validation data
set. a α-MLP, b α-GMDH, c N-MLP, d N-GMDH, e M-MLP,
f M-GMDH. The black solid line depicts the linear relationshipbetween predicted and observed values. The dashed red line isa one-to-one line and provides a benchmark for comparison

wetting path through additional experiments are currently underway. As more data is collected, the plan isto refine the existing models and expand the scope of

the predictive framework. This approach is believed
to maintain the integrity of the current findings whileproviding a clear direction for future investigation.
Vol.: (0123456789)

13

Geotech Geol Eng

Fig. 12  Actual versus estimated values for the testing data set.
a α-MLP, b α-GMDH, c N-MLP, d N-GMDH, e M-MLP, f
M-GMDH. The black solid line depicts the linear relationship

Vol:. (1234567890)

13

between predicted and observed values. The dashed red line isa one-to-one line and provides a benchmark for comparison

Geotech Geol Eng
Table 8  Details of theblind data set

Data Source

No. of data Type of soil Range of net Gs
stress (kPa)

PL

LL

γdmax

wopt

Ng et al. (2000)
Jadar et al. (2023)
Cai et al. (2023)(Seth et al. 2023)
Seth et al. (2023)

3
3
3
4
5

33.4
28.45
14.3
29
24

55.4
84.3
27.9
50
44

16.03
15.99
18
15.48
16.96

22
18.25
14.6
25
20

SM/SC
Ben-Sand
Loess
CH
CL

4.5 Sensitivity Analysis
We next performed a sensitivity analysis to better
understand how uncertainty in the model prediction
output can be attributed to various sources of uncertainty in the input parameters. Values of the sensitivities (SA) of each input parameter to the output were
calculated using:

1 ∑
SA =
Nt i=1Nt(

𝛿Output%
𝛿Input%

)

× 100

(13)

i

where δinput% is the percentage change in the input,δoutput% is the percentage change in the output, Nt is
the number of model inputs, and i is the sample point.
This indicates that the values of parameters are varied
from a maximum to a minimum value. When a specific input value highly affects the predicted output,
the corresponding SA value of that input is high. On
the other hand, when the input value has less effect on
the generated output, the corresponding SA value isrelatively low.
The sensitivities of the input parameters on
the three VG output parameters are displayed inFig. 17. Results show that soil type and net stresshave by far the highest effects on the output. On theother hand, the specific gravity (Gs) had a very marginal effect on the output parameters. Furthermore,
based on the sensitivity analysis, it was found thatthe maximum dry density (γdmax) and the optimum
water content (­ wopt) are not among the most significant parameters affecting the output (Fig. 17). It issomehow surprising that optimum compaction conditions were revealed not to affect the predictions
much. The reason behind is that majority of database for water retention tests do not comprise the

0—80
10—20050—250
0—120
0–200

2.6
2.6
2.68
2.75
2.67

optimum conditions during the retention tests. In
other words, many of the available test results wereobtained for samples with a reduced dry density
compared to that of optimum conditions simply due
to the fact that the corresponding natural soil strata
had a relatively low dry density such as naturalloess, lateritic soils or some residual soils. Hence,the soft computing-based predictive framework successfully distinguished between irrelevant factors
such as specific gravity, and influential parameterssuch as net stress on the generated results. This
gives much credit to the effectiveness and robustness of our study’s newly formulated AI-basedsoil–water retention model. It should be mentionedthat other influencing factors, such as the pore fluid
chemistry, may also lead to a change in the retention
properties of fine-grained soils (Sadeghi and Nasiri
2021; Hedayati-Azar and Sadeghi 2022; Darzi et al.
2023). As suggested by these references, the potential impact of pore fluid chemistry indicates thecomplexity of soil–water interactions and points toan area for further exploration and refinement of themodel.
5 Conclusions
A new model based on soft computing was derived
to predict stress-dependent soil water retention characteristics. The widely used van Genuchten (VG) soil
water retention model was implemented within themodeling framework to illustrate variations in the
SWRC caused by different net normal stresses andsoil physical properties. Accordingly, the three shapeparameters of van Genuchten’s model (i.e., a, N, andM) were determined using a relatively comprehensive

Vol.: (0123456789)

13

Geotech Geol Eng

Vol:. (1234567890)

13

Geotech Geol Eng◂Fig. 13  Comparing the obtained MAPE value for MLP and

GMDH methods, a α, b m, c ndatabase of 102 experiments taken from the literature.
Two different machine learning models (GMDH and
MLP) were employed to predict the α, N, and M asa function of soil type and net stress. Specifically, net
stress (σ), optimum water content (wopt), maximum
dry density (γdmax), plastic limit (PL), liquid limit(LL), specific gravity (Gs), and soil type were used as
input soil parameters to predict the output in terms ofthe three VG shape parameters. Results showed that
GMDH and MLP both performed excellently, with
the GMDH approach delivering slightly better results.
The sensitivity analysis indicated that net stress andsoil type were the two most significant variables governing the soil water retention characteristics of soil
and geomaterials.

Table 9  Details of theselected dataset for wetting
path

Table 10  Obtained results of RMSE and MAPE for the constructed MLP and GMDH models for wetting path
Metrics(model)VG parameter Training Validation Test

RMSE (MLP)
MAPE (MLP)
RMSE(GMDH)
MAPE(GMDH)RMSE (MLP)
MAPE (MLP)
RMSE(GMDH)
MAPE(GMDH)RMSE (MLP)
MAPE (MLP)
RMSE(GMDH)
MAPE(GMDH)

α

n

m

9.51
72.02
7.10

11.75
98.44
14.87

14.92
101.96
15.01

53.35

72.32

100.47

18.31
88.72
25.34

20.42
114.6616.34

37.81
121.80
16.91

70.21

88.36

91.47

0.18
89.48
0.13

0.27
109.270.11

0.25
126.62
0.13

88.44

84.09

84.31

Data Source

No. of data

Type of soil

Range ofnet stress(kPa)

Ng and Pang (2000)
Tavakoli and Habibagahi (2013)
Karube and Kawai (2001)
Tse (2007)
Oh and Lu (2014)
Ng et al. (2016)
Leung and Ng (2016)
Khosravi et al. (2016b)
Ng et al. (2020)
Zhang et al. (2020)
Albadri et al. (2021)
Seth et al. (2023)
Cai et al. (2023)
Total Number of Data:

3
3
2
3
4
3
3
3
4
6
3
9
3
49

SM – SCCL – ML
CL
SP
SM
CL
CL
ML
CL
ML
SP
CH / CL
Loess

0—80
0—250
19.6—98
0—80
0—200
0—150
0—80
125- 225
0—120
0—400
0—200
0—200
50—250Vol.: (0123456789)

13

Geotech Geol Eng

Fig. 14  Predicted versus measured values for the training data set of wetting path. a α-MLP, b α-GMDH, c N-MLP, d N-GMDH, e
M-MLP, f M-GMDH

Vol:. (1234567890)

13

Geotech Geol Eng

Fig. 15  Actual versus estimated values for the validation data set for wetting path. a) α-MLP, b α-GMDH, c N-MLP, d N-GMDH, e
M-MLP, f M-GMDH

Vol.: (0123456789)

13

Geotech Geol Eng

Fig. 16  Actual versus estimated values for the testing data set for wetting path. a α-MLP, b α-GMDH, c N-MLP, d N-GMDH, e
M-MLP, f M-GMDH

Vol:. (1234567890)

13

Geotech Geol Eng

Fig. 17  Sensitivity analysis of different input parameters: soil
type, net stress, Plastic limit (PL), Liquid limit (LL), Maximum dry density (γdmax), Optimum water content (­wopt), specific gravity (Gs)
Acknowledgements The authors are grateful for the financial
support provided by the Research Grant Office at Sharif University Technology. The corresponding author is also gratefulto the Iran’s National Elites Foundation for the financial support provided to him by way of “Dr Kazemi-Ashtiani Award”.
Data Availability The authors declare that they have noknown competing financial interests or personal relationshipsthat could have appeared to influence the work reported in thispaper.
Declarations
Conflict of interest The authors declare that they have noknown competing financial interests or personal relationshipsthat could have appeared to influence the work reported in thispaper.

References
Adeyemi O, Grove I, Peets S, Domun Y, Norton T (2018)
Dynamic neural network modelling of soil moisture content for predictive irrigation scheduling. Sensors (switzerland) 18:3408. https://​doi.​org/​10.​3390/​s1810​3408
Ahmadinezhad A, Jafarzadeh F, Sadeghi H (2019) Combination of water head control and axis translation techniques in new unsaturated cyclic simple shear tests. Soil
Dyn Earthq Eng 126:105818. https://​doi.​org/​10.​1016/j.​
soild​yn.​2019.​105818
Albadri WM, Noor MJM, Alhani IJ (2021) The relationship between the shear strength and water retention
curve of unsaturated sand at different hydraulic phases.
Acta Geotech 16:2821–2835. https://​doi.​org/​10.​1007/​
s11440-​021-​01189-7
Assouline S, Or D (2013) Conceptual and parametric representation of soil hydraulic properties: a review. Vadose
Zo J. https://​doi.​org/​10.​2136/​vzj20​13.​07.​0121

Bayat H, Neyshaburi MR, Mohammadi K, Nariman-Zadeh N,
Irannejad M, Gregory AS (2013) Combination of artificial neural networks and fractal theory to predict soil
water retention curve. Comput Electron Agric 92:92–103.
https://​doi.​org/​10.​1016/j.​compag.​2013.​01.​005
Bazargan A, Sadeghi H, Garcia-Mayoral R, McKay G (2015)
An unsteady state retention model for fluid desorption
from sorbents. J Colloid Interface Sci 450:127–134.
https://​doi.​org/​10.​1016/j.​jcis.​2015.​02.​036
Behzadafshar K, Sarafraz ME, Hasanipanah M, Mojtahedi SFF,
Tahir MM (2019) Proposing a new model to approximate
the elasticity modulus of granite rock samples based onlaboratory tests results. Bull Eng Geol Environ 78:1527–
1536. https://​doi.​org/​10.​1007/​s10064-​017-​1210-5
Biessmann F, Salinas D, Schelter S, Schmidt P, Lange D.
(2018) Deep learning for missing value imputation intables with non-numerical data. In: Int Conf Inf Knowl
Manag Proc 2017–2026. https://​doi.​org/​10.​1145/​32692​
06.​32720​05
Cai G, Liu Q, Li K, Zhang J, Liu Y, Zhou A (2023) Volume
change behavior of compacted loess under drying/wetting and freezing/thawing cycles. Eng Geol. https://​doi.​
org/​10.​1016/j.​enggeo.​2023.​107329
Darzi AG, Sadeghi H, Zhou C (2023) A soil-brine retention
model for wetting processes considering the hysteresis
effects. Transp Geotech 41:101032. https://​doi.​org/​10.​1016/j.​trgeo.​2023.​101032
Dias AS, Hughes PN, Toll DG (2023) Soil-water retention
curve prediction for compacted london clay subjected tomoisture cycles. Geotech Geol Eng. https://​doi.​org/​10.​
1007/​s10706-​023-​02706-z
dos Santos Pereira SA, Silva Junior AC, Mendes TA, Gitirana
Junior GDFN, Alves RD (2023) prediction of soil-water
characteristic curves in bimodal tropical soils using artificial neural networks. Geotech Geol Eng. https://​doi.​
org/​10.​1007/​s10706-​023-​02716-x
Eldin NN, Senouci AB (1994) Measurement and predictionof the strength of rubberized concrete. Cem Concr Compos 16:287–298. https://​doi.​org/​10.​1016/​0958-​9465(94)​90041-8
Elkady TY, Al-Mahbashi AM, Al-Refeai TO (2015) Stressdependent soil-water characteristic curves of limetreated expansive clay. J Mater Civ Eng 27:04014127.
https://​doi.​org/​10.​1061/​(asce)​mt.​1943-​5533.​00009​95
Fadillah IJ, Muchlisoh S (2019) Perbandingan metode hotdeck imputation dan metode knni dalam mengatasi
missing values. Semin Nas off Stat 2019(1):275–285.
https://​doi.​org/​10.​34123/​semna​soffs​tat.​v2019​i1.​101
Feng Y, Cui N, Hao W, Gao L, Gong D (2019) Estimation ofsoil temperature from meteorological data using different
machine learning models. Geoderma 338:67–77. https://​doi.​org/​10.​1016/J.​GEODE​RMA.​2018.​11.​044
Garakani AA, Haeri SM, Khosravi A, Habibagahi G (2015)
Hydro-mechanical behavior of undisturbed collapsible
loessial soils under different stress state conditions. Eng
Geol 195:28–41. https://​doi.​org/​10.​1016/j.​enggeo.​2015.​
05.​026
Garakani AA, Birgani MM, Sadeghi H (2021) An effective
stress-based parametric study on the seismic stabilityof unsaturated slopes with implications for preliminary

Vol.: (0123456789)

13

Geotech Geol Eng
microzonation. Bull Eng Geol Environ 80:7525–7549.
https://​doi.​org/​10.​1007/​s10064-​021-​02440-x
Geva S, Sitte J (1992) A constructive method for multivariate
function approximation by multilayer perceptrons. IEEE
Trans Neural Netw 3:621–624. https://​doi.​org/​10.​1109/​72.​
143376
Guido G, Haghshenas SS, Haghshenas SS, Vitale A, Gallelli
V, Astarita V (2020) Development of a binary classification model to assess safety in transportation systems
using GMDH-type neural network algorithm. Sustain
12(17):6735. https://​doi.​org/​10.​3390/​SU121​76735
Haeri SM, Khosravi A, Garakani AA, Ghazizadeh S (2017)
Effect of soil structure and disturbance on hydromechanical behavior of collapsible loessial soils. Int J Geomech
17:04016021.
https://​doi.​org/​10.​1061/​(asce)​gm.​1943-​
5622.​00006​56
Haeri SM, Garakani AA, Khosravi A, Meehan CL (2014a)
Assessing the hydro-mechanical behavior of collapsible
soils using a modified triaxial test device. Geotech Test J.
https://​doi.​org/​10.​1520/​GTJ20​130034
Haeri SM, Khosravi A, Ghaizadeh S, Garakani AA, Meehan
CL (2014b) Characterization of the effect of disturbanceon the hydro-mechanical behavior of a highly collapsible
loessial soil. In: Unsaturated soils: research and applications - proceedings of the 6th international conference on
unsaturated soils, UNSAT 2014. pp 261–266. https://​doi.​
org/​10.​1201/​b17034-​35
Haghverdi A, Cornelis WM, Ghahraman B (2012) A pseudocontinuous neural network approach for developing water
retention pedotransfer functions with limited data. J
Hydrol 442–443:46–54. https://​doi.​org/​10.​1016/j.​jhydr​ol.​
2012.​03.​036
Hedayati-Azar A, Sadeghi H (2022) Semi-empirical modelling of hydraulic conductivity of clayey soils exposed todeionized and saline environments. J Contam Hydrol
249:104042.
https://​doi.​org/​10.​1016/j.​jconh​yd.​2022.​
104042
Hornik K, Stinchcombe M, White H (1989) Multilayer feedforward networks are universal approximators. Neural
Netw 2:359–366. https://​doi.​org/​10.​1016/​0893-​6080(89)​90020-8
Hosseini SAA, Mojtahedi SFF, Sadeghi H (2020) Optimisationof deep mixing technique by artificial neural network basedon laboratory and field experiments. Georisk 14:142–157.
https://​doi.​org/​10.​1080/​17499​518.​2019.​16125​26
Ivakhnenko AG (1971) Polynomial theory of complex systems.
IEEE Trans Syst Man Cybern 1:364–378. https://​doi.​org/​
10.​1109/​TSMC.​1971.​43083​20
Ivakhnenko AG, Wunsch D, Ivakhnenko GA (1999) Inductive
sorting-out GMDH algorithms with polynomial complexity for active neurons of neural network. In: Proceedingsof the international joint conference on neural networks.
IEEE, pp 1169–1173. https://​doi.​org/​10.​1109/​ijcnn.​1999.​
831124
Jadar CM, Rajesh S, Roy S (2023) The effect of stress-dependent SWRC on the load carrying capacity of the slope subjected to the drying-wetting path. In: E3S web of conferences. https://​doi.​org/​10.​1051/​e3sco​nf/​20233​82120​04
Jafarzadeh F, Ahmadinezhad A, Sadeghi H (2021) Effects ofinitial suction and degree of saturation on dynamic properties of sand at large strain. Sci Iran 28(1):156–174
Vol:. (1234567890)

13

Jahangir H, Rezazadeh Eidgahee D (2021) A new and robust
hybrid artificial bee colony algorithm—ANN model for
FRP-concrete bond strength evaluation. Compos Struct257:113160. https://​doi.​org/​10.​1016/j.​comps​truct.​2020.​
113160Jana RB, Mohanty BP, Springer EP (2008) Multiscale Bayesian neural networks for soil water content estimation.
Water Resour Res. https://​doi.​org/​10.​1029/​2008W​R0068​
79
Jiang X, Wu L, Wei Y (2020) Influence of fine content on thesoil-water characteristic curve of unsaturated soils. Geotech Geol Eng 38:1371–1378. https://​doi.​org/​10.​1007/​
s10706-​019-​01096-5
Kaingo J, Tumbo SD, Kihupi NI, Mbilinyi BP (2018) Prediction of soil moisture-holding capacity with support vector
machines in dry subhumid tropics. Environ Soil Sci Appl.
https://​doi.​org/​10.​1155/​2018/​92632​96
Kang H (2013) The prevention and handling of the missing
data. Korean J Anesthesiol 64:402–406. https://​doi.​org/​10.​4097/​kjae.​2013.​64.5.​402
Karandish F, Šimůnek J (2016) A comparison of numerical andmachine-learning modeling of soil water content with limited input data. J Hydrol 543:892–909. https://​doi.​org/​10.​1016/j.​jhydr​ol.​2016.​11.​007
Karube D, Kawai K (2001) The role of pore water in themechanical behavior of unsaturated soils. Geotech Geol
Eng 19:211–241. https://​doi.​org/​10.​1023/A:​10131​88200​
053
Khlosi M, Alhamdoosh M, Douaik A, Gabriels D, Cornelis
WM (2016) Enhanced pedotransfer functions with support vector machines to predict water retention of calcareous soil. Eur J Soil Sci 67:276–284. https://​doi.​org/​10.​1111/​ejss.​12345
Khosravi A, Salam S, McCartney JS, Dadashi A (2016b) Suction-induced hardening effects on the shear modulus of
unsaturated silt. Int J Geomech 16:1–10. https://​doi.​org/​
10.​1061/​(asce)​gm.​1943-​5622.​00006​14
Khosravi A, Rahimi M, Shahbazan P, Pak A, Gheibi A (2016a)
Characterizing the variation of small strain shear modulus for silt and sand during hydraulic hysteresis. E3S Web
Conf. https://​doi.​org/​10.​1051/​e3sco​nf/​20160​914018
Kim BS, Shibuya S, Park SW, Kato S (2010) Application ofsuction stress for estimating unsaturated shear strengthof soils using direct shear testing under low confiningpressure. Can Geotech J 47:955–970. https://​doi.​org/​10.​1139/​T10-​007
Kit T (2007) Influence of stress states on soil-water characteristics, conjunctive surface-subsurface flow modelling andstability analysis (Doctoral dissertation)
Kovačević M, Bajat B, Gajić B (2010) Soil type classification and estimation of soil properties using support vector machines. Geoderma 154:340–347. https://​doi.​org/​10.​1016/j.​geode​rma.​2009.​11.​005
Lee IM, Sung SG, Cho GC (2005) Effect of stress state on theunsaturated shear strength of a weathered granite. Can
Geotech J 42:624–631. https://​doi.​org/​10.​1139/​t04-​091
Leung AK, Ng CWW (2016) Field investigation of deformation characteristics and stress mobilisation of a soil
slope. Landslides 13:229–240. https://​doi.​org/​10.​1007/​
s10346-​015-​0561-x

Geotech Geol Eng
Li J, Sun D, Sheng D, Sloan SW, Fredlund DG (2007) Preliminary study on soil-water characteristics of Maryland clay.
In: Proc. 3rd Asian conf. unsaturated soils c, pp 569–574
Mahmoodabadi M, Bryson LS (2021) Direct application of thesoil-water characteristic curve to estimate the shear modulus of unsaturated soils. Int J Geomech 21:04020243.
https://​doi.​org/​10.​1061/​(asce)​gm.​1943-​5622.​00018​93
Menad NA, Noureddine Z, Hemmati-Sarapardeh A, Shamshirband S, Mosavi A, Chau K (2019) Modeling temperature dependency of oil - water relative permeability
in thermal enhanced oil recovery processes using groupmethod of data handling and gene expression programming. Eng Appl Comput Fluid Mech 13:724–743. https://​doi.​org/​10.​1080/​19942​060.​2019.​16395​49
Mendes RM, Marinho FAM (2020) Soil water retention curves
for residual soils using traditional methods and MIP. Geotech Geol Eng 38:5167–5177. https://​doi.​org/​10.​1007/​
s10706-​020-​01354-x
Mojtahedi SFF, Ebtehaj I, Hasanipanah M, Bonakdari H,
Amnieh HB (2019a) Proposing a novel hybrid intelligent model for the simulation of particle size distribution
resulting from blasting. Eng Comput 35:47–56. https://​doi.​org/​10.​1007/​s00366-​018-​0582-x
Mojtahedi SFF, Tabatabaee S, Ghoroqi M, Soltani Tehrani
M, Gordan B, Ghoroqi M (2019b) A novel probabilistic
simulation approach for forecasting the safety factor ofslopes: a case study. Eng Comput 35:637–646. https://​doi.​
org/​10.​1007/​s00366-​018-​0623-5
Mojtahedi SFF, Hosseini SAA, Sadeghi H (2023) An artificial intelligence based data-driven method for forecasting
unconfined compressive strength of cement stabilized soil
by deep mixing technique. Geotech Geol Eng 41(1):491–
514. https://​doi.​org/​10.​1007/​s10706-​022-​02297-1
Morosini AF, Haghshenas SS, Haghshenas SS, Geem ZW(2020) Development of a binary model for evaluatingwater distribution systems by a pressure driven analysis(PDA) approach. Appl Sci. https://​doi.​org/​10.​3390/​app10​
093029Mulashani AK, Shen C, Nkurlu BM, Mkono CN, Kawamala M(2022) Enhanced group method of data handling (GMDH)
for permeability prediction based on the modified Levenberg Marquardt technique from well log data. Energy
239:121915. https://​doi.​org/​10.​1016/j.​energy.​2021.​121915
Ng CWW, Pang YW (2000) Influence of stress state on soilwater characteristics and slope stability. J Geotech Geoenviron Eng 126:157–166. https://​doi.​org/​10.​1061/​(asce)​
1090-​0241(2000)​126:​2(157)
Ng CWW, Yung SY (2008) Determination of the anisotropic
shear stiffness of an unsaturated decomposed soil. Geotechnique 58:23–35. https://​doi.​org/​10.​1680/​geot.​2008.​
58.1.​23
Ng CWW, Pang WY, Ng CWW, Pang YW (2000) Influenceof stress state on soil-water characteristics and slope stability. J Geotech Geoenviron Eng 126(2):157–166. https://​doi.​org/​10.​1061/​(ASCE)​1090-​0241(2000)​126:​2(157)
Ng CWW, Owusu ST, Zhou C, Chiu ACF (2020) Effects ofsesquioxide content on stress-dependent water retention
behaviour of weathered soils. Eng Geol 266:105455.
https://​doi.​org/​10.​1016/j.​enggeo.​2019.​105455

Ng CWW, Sadeghi H, Hossen SB, Chiu CF, Alonso EE, Baghbanrezvan S (2016) Water retention and volumetric characteristics of intact and re-compacted loess. Can GeotechnJ 53:1258–1269. https://​doi.​org/​10.​1139/​cgj-​2015-​0364
Oh S, Lu N (2014) Uniqueness of the suction stress characteristic curve under different confining stress conditions.
Vadose Zo J 13:1–10. https://​doi.​org/​10.​2136/​vzj20​13.​04.​
0077
Onyelowe KC, Mojtahedi FF, Azizi S, Mahdi HA, Sujatha ER,
Ebid AM, Darzi AG, Aneke FI (2022) Innovative overview of SWRC application in modeling geotechnical engineering problems. Designs 6:69. https://​doi.​org/​10.​3390/​
desig​ns605​0069
Raheem AM (2023) Developing a new soil-water characteristic
model using nondestructive electrical properties for unsaturated soil. Geotech Geol Eng 41:1589–1601. https://​doi.​
org/​10.​1007/​s10706-​022-​02355-8
Rezaee M, Mojtahedi SFF, Taherabadi E, Soleymani K, Pejman M (2021) Prediction of shear strength parameters ofhydrocarbon contaminated sand based on machine learning methods. Georisk 15:317–335. https://​doi.​org/​10.​
1080/​17499​518.​2020.​18616​33
Roy S, Rajesh S (2020) Simplified model to predict features ofsoil-water retention curve accounting for stress state conditions. Int J Geomech 20:04019191. https://​doi.​org/​10.​1061/​(asce)​gm.​1943-​5622.​00015​91
Sadeghi H, Darzi AG (2021a) Modelling of soil-water retention curve considering the effects of existing salt solution
in the pore fluid. MATEC Web Conf 337:02001. https://​doi.​org/​10.​1051/​matec​conf/​20213​37020​01
Sadeghi H, Darzi AG (2021b) A review of different approaches
to analytical‌modeling of soil-water retention curve. SJCE
37(3.1):111–123.
https://​doi.​org/​10.​24200/​J30.​2021.​
56932.​2870
Sadeghi H, Darzi AG (2023) A soil-water retention model
incorporating pore-fluid osmotic potential. Int J Geomech
23(11):04023209.
https://​doi.​org/​10.​1061/​IJGNAI.​
GMENG-​8648
Sadeghi H, Nasiri H (2021) Hysteresis of soil water retention
and shrinkage behaviour for various salt concentrations.
Géotechnique Lett 11(1):21–29. https://​doi.​org/​10.​1680/​
jgele.​20.​00047
Sadeghi H, Hossen SB, Chiu AC, Cheng Q, CWW N (2016)
Water retention curves of intact and recompactedloessat different net stresses. Jpn Geotech Soc Spec Publ
2(4):221–225.https://​doi.​org/​10.​3208/​jgssp.​HKG-​04
Sadeghi H, Heydari A, Zhou C (2024) Insights into the fluid
retention and shrinkage of a lean clay under the combined influence of microstructure, solute concentrationand salt species. Acta Geotech. https://​doi.​org/​10.​1007/​
s11440-​023-​02206-7
Sallaby AF, Azlan A (2021) Analysis of missing value imputation application with K-nearest neighbor (K-NN) algorithm in dataset. IJICS (int J Inf Comput Sci) 5(2):141–
144. https://​doi.​org/​10.​30865/​ijics.​v5i2.​3185
Satyanaga A, Rahardjo H, Zhai Q, Moon SW, Kim J (2023)
Modelling particle-size distribution and estimation ofsoil–water characteristic curve utilizing modified lognormal distribution function. Geotech Geol Eng. https://​doi.​org/​10.​1007/​s10706-​023-​02638-8

Vol.: (0123456789)

13

Geotech Geol Eng
Schaap MG, Leij FJ, van Genuchten MT (1998) Neural network analysis for hierarchical prediction of soil hydraulic properties. Soil Sci Soc Am J 62:847–855. https://​doi.​org/​10.​2136/​sssaj​1998.​03615​99500​62000​40001x
Seth OT, Charles WWN, Chao Z (2023) An experimental
study of stress-dependent water retention behaviour of
two lateritic clays with different minerals. Eng Geol.
https://​doi.​org/​10.​1016/j.​enggeo.​2023.​107398
Song I, Yang Y, Im J, Tong T, Ceylan H, Cho IH (2020)
Impacts of fractional hot-deck imputation on learning andprediction of engineering data. IEEE Trans Knowl Data
Eng 32:2363–2373. https://​doi.​org/​10.​1109/​TKDE.​2019.​
29226​38
Sun D, Sheng D, Xu Y (2007) Collapse behaviour of unsaturated compacted soil with different initial densities. Can
Geotech J 44:673–686. https://​doi.​org/​10.​1139/​T07-​023
Tavakoli MH, Habibagahi G. (2013) Influence of confining
stress on the effective stress parameter. In: Poromechanics V: Proceedings of the fifth biot conference on poromechanics. https://​doi.​org/​10.​1061/​97807​84412​992.​165
Tavakoli Dastjerdi MH, Habibagahi G, Nikooee E (2014)
Effect of confining stress on soil water retention curveand its impact on the shear strength of unsaturated soils.
Vadose Zo J. https://​doi.​org/​10.​2136/​vzj20​13.​05.​0094
Tse MK (2007) Influence of stress states on soil-water characteristics, conjunctive surface-subsurface flow modellingand stability analysis (Doctoral dissertation)
Twarakavi NKC, Šimůnek J, Schaap MG (2009) Developmentof pedotransfer functions for estimation of soil hydraulicparameters using support vector machines. Soil Sci Soc
Am J 73:1443–1452. https://​doi.​org/​10.​2136/​sssaj​2008.​
0021
Uchaipichat A, Khalili N (2009) Experimental investigationof thermo-hydro-mechanical behaviour of an unsaturated
silt. Geotechnique 59:339–353. https://​doi.​org/​10.​1680/​
geot.​2009.​59.4.​339
Uddin S, Lu H (2024) Dataset meta-level and statistical features affect machine learning performance. Sci Rep
14(1):1670. https://​doi.​org/​10.​1038/​s41598-​024-​51825-x
Vabalas A, Gowen E, Poliakoff E, Casson AJ (2019) Machinelearning algorithm validation with a limited sample size.
PLoS ONE 14(11):e0224365. https://​doi.​org/​10.​1371/​
journ​al.​pone.​02243​65
Van Genuchten MT (1980) A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil
Sci Soc Am J 44:892–898. https://​doi.​org/​10.​2136/​sssaj​
1980.​03615​99500​44000​50002x
Vanapalli SK, Fredlund DG, Pufahl DE (1999) The influence of soil structure and stress history on the soil-water

Vol:. (1234567890)

13

characteristics of a compacted till. Geotechnique 49:143–
159. https://​doi.​org/​10.​1680/​geot.​1999.​49.2.​143
Vanapalli SK, Fredlund DG, Pufahl DE (2015) The influenceof soil structure and stress history on the soil–water characteristics of a compacted till. Géotechnique 49:143–159.
https://​doi.​org/​10.​1680/​geot.​1999.​49.2.​143
Vissikirsky VA, Stepashko VS, Kalavrouziotis IK, Drakatos
PA (2005) Growth dynamics of trees irrigated with wastewater: GMDH modeling, assessment, and control issues.
Instrum Sci Technol 33:229–249. https://​doi.​org/​10.​1081/​
CI-​20004​8085
Wang C, Li SY, He XJ, Chen Q, Zhang H, Liu XY (2021)
Improved prediction of water retention characteristic
based on soil gradation and clay fraction. Geoderma
404:115293.
https://​doi.​org/​10.​1016/j.​geode​rma.​2021.​
115293
Yao Y, Ni J, Li J (2021) Stress-dependent water retention of
granite residual soil and its implications for ground settlement. Comput Geotech 129:103835. https://​doi.​org/​10.​1016/j.​compg​eo.​2020.​103835
Zhang T, Zhang J, Jiang T, Wang X, Jia H, Wang L (2019)
SWCCs of silt in Yudong zone and its prediction
under different drying-wetting cycle conditions. Geotech Geol Eng 37:1977–1986. https://​doi.​org/​10.​1007/​
s10706-​018-​0738-x
Zhang D, Wang J, Chen C, Wang S (2020) The compression and collapse behaviour of intact loess in suctionmonitored triaxial apparatus. Acta Geotech 15:529–548.
https://​doi.​org/​10.​1007/​s11440-​019-​00829-3
Zhao Y, Rahardjo H, Satyanaga A, Zhai Q, He J (2023) A general best-fitting equation for the multimodal soil-water
characteristic curve. Geotech Geol Eng 41:3239–3252.
https://​doi.​org/​10.​1007/​s10706-​023-​02447-z
Zhou C, Ng CWW (2014) A new and simple stress-dependent
water retention model for unsaturated soil. Comput Geotech 62:216–222. https://​doi.​org/​10.​1016/j.​compg​eo.​2014.​
07.​012Publisher’s Note Springer Nature remains neutral with regard
to jurisdictional claims in published maps and institutional
affiliations.
Springer Nature or its licensor (e.g. a society or other partner)
holds exclusive rights to this article under a publishingagreement with the author(s) or other rightsholder(s); author
self-archiving of the accepted manuscript version of this article
is solely governed by the terms of such publishing agreement
and applicable law.