Environmental Modelling and Software 150 (2022) 105341

Contents lists available at ScienceDirect

Environmental Modelling and Software
journal homepage: www.elsevier.com/locate/envsoft

Balancing exploitation and exploration: A novel hybrid global-local
optimization strategy for hydrological model calibration
Giuseppe Brunetti a, *, Christine Stumpp a, Jiří Šimůnek bab

University of Natural Resources and Life Sciences, Vienna (BOKU), Institute for Soil Physics and Rural Water Management, Muthgasse 18, 1180, Vienna, Austria
University of California, Riverside, Department of Environmental Sciences, CA, 92521, USA

A R T I C L E I N F O

A B S T R A C T

Keywords:
Optimization
Model calibration
HYDRUS
Environmental modeling

Optimization problems in hydrological modeling are frequently solved using local or global search strategies,
which either maximize exploitation or exploration. Thus, the elevated performance of one strategy for one classof problems is often offset by poor performance for another class. To overcome this issue, we propose a hybrid
strategy, G-CLPSO, that combines the global search characteristics of the Comprehensive Learning Particle
Swarm Optimization (CLPSO) with the exploitation capability of the Marquardt-Levenberg (ML) method andimplement it into the hydrological model, HYDRUS. Benchmarks involving optimizing non-separable unimodaland multimodal functions demonstrate that G-CLPSO outperforms CLPSO in terms of accuracy and convergence.
Synthetic modeling scenarios involving the inverse estimation of soil hydraulic properties are used to compare
the G-CLPSO against the original HYDRUS ML solver, the gradient-based algorithm PEST, and the stochastic SCEUA strategy. Results demonstrate the superior performance of the G-CLPSO, suggesting a potential use in otherenvironmental problems.

1. Introduction
Due to the poor transferability of laboratory-estimated parameters(e.g., Ritter et al., 2003), hydrological models are generally calibrated
against transient field measurements for predictive purposes. The main
goal is to inversely estimate model parameters by minimizing the de­viations between model predictions and observations (i.e., the objectivefunction). The same strategy can also be used for designing purposes [e.
g., irrigation system optimization (Roy et al., 2019)]. These optimization
problems are frequently dealt with using local or global search strate­gies, which either maximize exploitation or exploration. Thus, theelevated performance of one strategy for one class of problems is often
offset by poor performance for another class.
A classic example is the hydrological model HYDRUS (Šimůnek et al.,2016), one of the most widespread mechanistic models for analyzingwater flow and transport processes in the vadose zone. HYDRUS in­cludes the Marquardt–Levenberg (ML) algorithm (Levenberg, 1944;
Marquardt, 1963) to solve a nonlinear least-squares minimization
problem arising from the model calibration (Šimůnek and Hopmans,
2002). The ML is a gradient-based local search method, which has beenused extensively in the past in soil physics and vadose zone hydrology (e.

g., Kool and Parker, 1988; Šimůnek et al., 1998; Šimůnek and Van
Genuchten, 1996; van Dam et al., 1994) due to its efficiency forwell-posed low-dimensional unimodal inverse problems. However,calibration problems are often ill-posed and non-separable due tomeasurement errors and model inadequacy, thus leading to response
surfaces characterized by long banana-shaped narrow valleys or multi­
ple local optima (e.g., Abbaspour et al., 2004; Arora et al., 2012; Bru­
netti et al., 2016; Šimůnek and Van Genuchten, 1996). In these
conditions, the ML algorithm becomes extremely sensitive to the initialstarting point and can easily converge to a local minimum (Transtrum
and Sethna, 2012). To deal with irregularities in the response surface inhigh dimensions, the PEST software (White et al., 2020) couples an ML
solver with a subspace parameter dimensionality reduction based on the
singular value decomposition. The resulting optimization strategy
generated significant computational gains and was successfully appliedin combination with HYDRUS (Watson et al., 2013). Nevertheless, its
local search behavior can still lead to a modest performance in complex
optimization tasks (e.g., Sedaghatdoost et al., 2019).
To overcome this limitation, multiple studies have externallycoupled HYDRUS with global optimization strategies such as theannealing-simplex method (Pan and Wu, 1998), genetic algorithms (Ines* Corresponding author.
E-mail address: giuseppe.brunetti@boku.ac.at (G. Brunetti).
https://doi.org/10.1016/j.envsoft.2022.105341
Received 29 June 2021; Received in revised form 28 January 2022; Accepted 28 January 2022
Available online 1 February 2022
1364-8152/© 2022 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).

G. Brunetti et al.

Environmental Modelling and Software 150 (2022) 105341

and Droogers, 2002), ant colony optimization (Abbaspour et al., 2001),
shuffled complex methods (Duan et al., 1993; Vrugt et al., 2003), andparticle swarm optimization (Brunetti et al., 2016, 2018), among manyothers. Global optimization algorithms outperform local search methods
for high-dimensional complex problems due to their high exploration
capability. However, according to the so-called “no free lunch” theorem(Wolpert and Macready, 1997): “For any algorithm, any elevated perfor­
mance over one class of problems is offset by poor performance over anotherclass.” Thus, the high performance of global search strategies on com­plex multimodal problems is counterbalanced by their slow convergenceon unimodal problems. Nevertheless, since the shape of the response
surface (i.e., objective function) is not known a priori when solving
real-world problems, it is preferable to apply global optimization
methods.
Particle Swarm Optimizers (PSOs) are a promising tool due to theirease of implementation and good performance on moderately highdimensional inverse problems. To improve their efficiency, multiple
studies have coupled them with local search methods (e.g., Cao et al.,2019; Han and Liu, 2015; Noel, 2012; Zhao et al., 2008). For instance, inNoel (2012), the steepest descent scheme is used at every iteration or
periodically to survey the vicinity around the swarm’s best solution(exploitation) to complement the PSO’s global exploration (Shi andEberhart, 1998). However, if the local search is performed too early, theswarm can lose diversity and converge to a local optimum. Recently, Caoet al. (2019) coupled the Comprehensive Learning Particle Swarm
Optimization (CLPSO) (Liang et al., 2006) with both the Nelder-Meadand Broyden–Fletch–Goldfarb–Shannon algorithms and proposed tostart the local search only when all particles in the CLPSO enter the
global optimum basin. By doing so, both the convergence time and ac­
curacy of the CLPSO improve. However, the limited interaction between
the CLPSO and the local search during the optimization process suggeststhat there is still room for improvement, especially for unimodal
problems.
Starting from this background, the main goal of this study is todevelop an optimization strategy that combines the exploitation capa­
bility of the ML algorithm with the exploration features of the CLPSO,
which was shown to outperform other algorithms for multimodal
problems (Liang et al., 2006). We aim to obtain an algorithm that is
flexible enough to be applied to multiple optimization tasks without
deteriorating its performance. The developed approach is implemented
in the HYDRUS source code as an inverse solver to extend its applica­
bility to moderately high-dimensional multimodal optimization prob­
lems. The internal coupling overcomes the usability limitations ofexternal coupling routines by making the approach readily available toall HYDRUS users. The problem is addressed in the following way. First,
we describe the ML and CLPSO algorithms. Next, we explain their
combination in the Gradient-Based Comprehensive Learning Particle
Swarm Optimization (G-CLPSO) and its implementation into the HYD­
RUS source code. Non-separable and ill-conditioned unimodal andmultimodal test functions are then used to compare the original CLPSO
and newly developed strategies and assess achieved improvements.
Finally, synthetic modeling scenarios involving the inverse estimation ofa layered lysimeter’s soil hydraulic properties are used to test the
G-CLPSO performance and assess its benefit against the original HYD­
RUS ML solver, the gradient-based PEST strategy (White et al., 2020),
and the stochastic search algorithm SCE-UA (Duan et al., 1993).

minimize ​ f (x1 , x2 , …, xD )subjectto ​ xli < xi < xiu ​ with ​ i = 1, 2, …, D

(1)where xil and xiu are the lower and upper bounds of the ith parameter,
respectively. Since the optimization strategy includes the ML algorithm,
the fitness f is the Sum of Squared Residuals (SSQ).
2.1.1. Comprehensive Learning Particle Swarm Optimization
The PSO is a metaheuristic and gradient-free global search strategy
based on a social-psychological metaphor involving individuals of aswarm that interact with each other in a social world to achieve an
optimum state (Kennedy and Eberhart, 1995; Shi and Eberhart, 1998).
By making no assumptions on the problem being optimized, it is espe­
cially suited for black-box optimization problems often encountered inenvironmental sciences (e.g., Brunetti et al., 2018, 2016; Chu andChang, 2009; Xi et al., 2017).
The swarm consists of S individuals or particles (i.e., potential so­lutions). Each particle j is characterized by three D-dimensional vectors→
and a scalar value. These are its current position X j , its best position
best→→
X j , its velocity V j , and its best objective function fjbest . The current
position is a set of coordinates in the parameters’ space, which is
updated at each iteration using a specific learning strategy that involves
an information exchange with some neighboring individuals or thewhole swarm. This communication allows the swarm to progressively
locate better solutions and eventually converge to the global optimum.
Intuitively, the core of the PSO is social interaction. In the original
PSO, the particle’s position is updated by random averaging of the best
position encountered by the particle so far and the best position found by
the neighboring particles. However, this trade-off between individual
and social interests can reduce the swarm’s diversity, thus leading topremature convergence to a local optimum in complex multimodal
problems. The CLPSO (Liang et al., 2006) overcomes this limitation byusing a different learning strategy that preserves diversity and reducesswarm behavior. The jth particle’s movement in the ith dimension is
described as follows:)({
Vj,i ←ωVj,i + Ui (0, 1)⋅c⋅ Xfbest− Xj,i ​j (i),i(2)
Xj,i ←Xj,i + Vj,iwhere Vj,i and Xj,i are the velocity and position of the ith parameter of thejth particle, respectively; ω is an inertia weight that is reduced as thenumber of iterations grows to favor exploitation, Ui is a random numberin the range [0, 1], c is a learning parameter typically set to 1.4995, and→best
Xfbestdefines which particles’ X j the particle j should follow. Xfbest
canj (i),ij (i),i→best
be the corresponding dimension of any particle’s X j including its own,

and the decision depends on probability Pc, referred to as the learning
probability, which takes different values for different particles:()1)
exp 10(j−− 1
S− 1
Pcj = 0.05 + 0.45
(3)exp(10) − 1
For each dimension of the particle j, a random number is generated.
If this number is lower than Pcj, the corresponding dimension learns
from other particles based on a tournament selection procedure;→best
otherwise, it will learn from its own X j . The learning process is
continuously monitored. If a particle ceases to improve for a certain
number of iterations (i.e., refreshing gap), m, the exemplar from which
the particle is learning, is reassigned. A thorough description of the se­lection procedure is reported in Liang et al. (2006).
As problem dimensionality increases, the particles tend to leave thesearch space and exhibit unwanted roaming behavior. To address thisissue, both the particle’s velocity and position should be corrected. In

2. Materials and methods
2.1. Optimization strategy
In this section, we first provide a general description of the CLPSO
and ML algorithms and then explain their combination in the G-CLPSO,
as well as its implementation into the HYDRUS source code. Consideredis the following optimization problem in a bounded parameter space ofdimension D:
2

G. Brunetti et al.

Environmental Modelling and Software 150 (2022) 105341

the CLPSO, the particle’s velocity in each dimension is bounded by a
⃒ ⃒maximum magnitude, Vmax. If ⃒Vj,i ⃒ exceeds a positive constant value,⃒ ⃒Vmax,j , then the velocity of that dimension is set to sign(⃒Vj,i ⃒)Vmax,j .

by Stanovov et al. (2019), who concluded that rank and tournament
selection could significantly enhance the performance of differential
evolution strategies by increasing their exploitation capabilities while
not affecting exploration significantly. However, one inherent danger
with high selective pressure is the rapid loss of diversity within the
population. Thus, it is reasonable to assume that the coupling with alocal search algorithm will further increase selective pressure, withdetrimental effects for multimodal fitness landscapes.
In the Gradient-based Comprehensive Learning Particle Swarm
Optimization (G-CLPSO) developed here, we propose a loose couplingstrategy between the CLPSO and ML. The CLPSO is used first for NL it­
erations. Then, one random individual rand is selected from the swarm
and Xrand is used as the starting point for the ML local search. If thecalculated fitness valuef(Xrand ), is lower than the corresponding per­
best
best
sonal best frand, then Xrandis replaced by the optimum found by the local→best
Simultaneously, the fitness value of a particle is calculated and its X j

is updated only if its position is within the limits of the search range. Thisguarantees that the particle will return in the search range in thefollowing iterations.
The main advantage of the CLPSO learning strategy is that all par­ticles can be potentially used to guide the search direction of other in­dividuals, and this learning process can be different for each dimensionof the particle. Such reduced emphasis on the global best position pre­serves the swarm’s diversity and enhances the algorithm’s performanceon complex multimodal problems (Liang et al., 2006). However, thishigh exploration capability delays the convergence for unimodal
problems.

newbest. By doing so, the new Xrand
can enter the CLPSO tournament
search, Xrand
selection procedure and improve the swarm without significantly
reducing the diversity of the swarm. The main advantage is the possi­
bility to start the ML algorithm after every NL iterations at different
points, with practical implications for different types of problems:

2.1.2. Marquardt-Levenberg algorithm
The Marquardt-Levenberg algorithm is a gradient-based techniquethat combines the steepest descent and Gauss-Newton methods (Lev­
enberg, 1944; Marquardt, 1963) to solve nonlinear least-squares mini­
mization problems. In the classic steepest descent scheme, the algorithm
is started at a selected position in the search space and small steps are
made against the direction of the gradient according to:) →( )( (→( ))→(
X k + 1 = X k − β∇ C X k(4)- Multimodal problems: This allows to progressively identify different
local minima, which will then enter and improve the tournament
selection procedure without hampering the swarm’s diversity.
- Unimodal problems: The ML will calculate similar fitness values fordifferent starting points, thus driving the entire swarm towards theglobal optimum in a few iterations.

→where X (k +1) is the approximation of the local minimum at the iter­→
ation k+1, β is the step size, and ∇(C( X (k))) is the gradient of the→function evaluated at X (k). The gradient consists of first-order de­rivatives, which are calculated using the forward difference approxi­
mation. This step can be relatively computationally intensive when the
inverse problem’s dimensionality is high and can lead to inaccuracies
when the response surface is topologically complex. Large β will result ina faster convergence towards the local minimum, yet oscillations will→appear once the algorithm reaches the sub-optimal regions as X (k)overshoots the minimum. On the other hand, a small step size will
improve stability but reduce the algorithm’s convergence rate.
The ML algorithm solves these issues by using the steepest descentmethod when the objective function is far from its minimum, and theGauss-Newton method as the minimum is approached. This switch isoperated by the parameter λ, which is set to a modest value at the
beginning (e.g., 0.02) and then reduced as the solution approaches theminimum. The algorithm stops when a certain number of iterations,
Nmax, is reached, or no further reductions of the fitness value are ob­tained. Results are then used to calculate the correlation matrix between
estimated parameters, and their confidence intervals, thus providing astatistical basis for uncertainty assessment.
The ML algorithm outperforms stochastic search algorithms forunimodal response surfaces characterized by a single well-defined op­
timum, or when the initial point lies near the global optimum. In suchcircumstances, the high exploration ability of stochastic search methods
is redundant. However, when the problem’s dimensionality is high, andthe response surface is noisy and topologically complex, the ML is likely
to converge early to a local optimum (Transtrum and Sethna, 2012).

The algorithm stops when the total number of function evaluations
exceeds a user-defined threshold value, NFES
max , or the difference between
the best and worst objective function values is below a user-definedtolerance value, tol. The latter criterion has proven to be reliable formultiple problems (Zielinski and Laur, 2007). The pseudocode of the
G-CLPSO is shown in Fig. 1.
2.2. Algorithm testingThis section describes the test functions and synthetic case studies
used to assess the performance of the G-CLPSO for both widely usedbenchmark functions and classic vadose zone inverse estimation prob­
lems. In particular, test functions are used only to highlight numericalimprovements of the G-CLPSO against the original CLPSO and investi­
gate the influence of G-CLPSO input parameters, S and NL, on the al­
gorithm’s performance. Conversely, the synthetic scenarios serve tocompare the G-CLPSO, the original HYDRUS inverse solver, thegradient-based PEST optimization strategy (White et al., 2020), and thestochastic search algorithm SCE-UA (Duan et al., 1993) for specific hy­
drological problems. Algorithms are evaluated in terms of accuracy,
which is the capability to approximate better the global optimum, andconvergence speed, which is the number of function calls required to
identify the optimum or meet the termination criteria.
The use of test functions and synthetic scenarios is an often-usedpractice in hydrological modeling for testing newly developed algo­rithms. For instance, Duan et al. (1993) tested the shuffled complex
evolution approach against six benchmark functions and one synthetic
dataset obtained using the SIXPAR hydrological model. In anotherstudy, Duan et al. (1994) used a synthetic sequence of streamflow to test
the SCE-UA optimization strategy. Multiple numerical experiments wereused by Vrugt and Robinson (2007) and Ter Braak and Vrugt (2008) toassess the performance of the AMALGAM and Differential Evolution
Monte Carlo algorithms, respectively. The a priori knowledge of the
global optimum and topological features of the response surface gua­rantees a fair assessment of the proposed optimization strategy.

2.1.3. Gradient-based Comprehensive Learning Particle Swarm
Optimization
The critical point in hybrid global-local search strategies is to find a
trade-off between their conflicting explorative and exploitative in­terests. One of the main advantages of the CLPSO over other PSOs is theuse of selective pressure in the learning strategy. Sutton et al. (2007)
noted how the application of selective pressure to the optimization of anon-separable function leads to significant improvements in differential
evolution strategies by focusing the search. This was further confirmed
3

G. Brunetti et al.

Environmental Modelling and Software 150 (2022) 105341

Fig. 1. The pseudocode of the G-CLPSO.

2.2.1. Test functions
As the G-CLPSO is supposed to outperform the original CLPSO in
unimodal problems while preserving its performance for multimodal
functions, we selected one unimodal and two multimodal functions fromthe list of benchmark problems used at the IEEE Congress on Evolu­
tionary Computation (Liang et al., 2014) to test the proposed strategy.
The choice of the test functions is not casual but intentionally directed
towards the inclusion of mathematical functions whose landscapesclosely resemble the response surfaces of real hydrological optimization
problems. All functions are non-separable since model parameters often
exhibit high interactions in hydrological inverse problems. To this aim,
different rotation matrices, M, are assigned to each function, which isalso shifted by using a vector o. The search range is [− 100,100]D for all
problems, which are tested in ten and thirty dimensions. The maximum
FES
number of function evaluations Nmax
and the convergence tolerance, tol,are set to 10,000D and 0.001, respectively. The properties and formulasof selected functions are listed in Table 1, while their fitness landscapes
are shown in Fig. 2.
Twenty-five combinations of S and NL are generated. For each ofthem, the average final fitness and the total number of function evalu­
ations for each test function are calculated by averaging the results of ten
independent G-CLPSO optimizations. Finally, a color map is used to
explain the sensitivity of the G-CLPSO performance to changes in the
population size and local search frequency. The better combination is
identified, and then it is used in the following comparison between the
G-CLPSO and the CLPSO.
2.2.2. Implementation of G-CLPSO into the HYDRUS code
The main reason to implement the G-CLPSO directly into the HYD­
RUS code is to extend the capabilities of the existing HYDRUS inverse
4

G. Brunetti et al.

Environmental Modelling and Software 150 (2022) 105341

2.2.3. Synthetic modeling scenarios and model settings
We use three synthetic modeling scenarios focused on the inverse
estimation of a layered lysimeter’s soil hydraulic properties to test theperformance of the G-CLPSO in vadose zone-related problems and assess
its benefit against the original HYDRUS ML solver, the gradient-based
PEST optimization strategy (White et al., 2020), and the stochastic
search algorithm SCE-UA (Duan et al., 1993). The synthetic lysimeter isa widely used numerical experiment for the inverse estimation of soil
hydraulic properties (e.g., Durner et al., 2008; Schelle et al., 2012).
A 150 cm deep monolithic lysimeter installed in Gumpenstein(Austria) has been used as a reference (Herndl et al., 2011). The
lysimeter consists of three soil layers, classified from top to bottom asloam, silt loam, and silt. The vegetation consists of grass with homo­
geneous roots distributed to a depth of 20 cm. Climatic data collected bya weather station are used to calculate daily potential evapotranspira­tion (Hargreaves, 1994). The evapotranspiration demand is then parti­
tioned into potential evaporation and transpiration fluxes using the Leaf
Area Index (Sutanto et al., 2012), which is constant during the simula­tion period and equal to 2.0. Precipitation is measured using a tipping
bucket rain collector placed nearby. Data from January to September
2012 are used in the numerical simulations.
HYDRUS-1D simulates variably-saturated water flow and root water
uptake in the lysimeter by numerically solving the Richards equation(Richards, 1931) and the Feddes macroscopic model (Feddes et al.,1978), respectively. The van Genuchten-Mualem functions (VGM) (vanGenuchten, 1980) describe the unsaturated soil hydraulic properties of
the three layers, and their parameters are initially set according to Carsel
and Parrish (1988) to generate synthetic data (Table 2). Field volumetricwater content measurements at a pressure head of − 15,000 cm are used
to adjust the residual water content, θr. The domain is discretized into
100 finite elements. An atmospheric boundary condition is applied at
the soil surface, while a seepage face boundary condition is set at thebottom.
Warm-up simulations covering a period of 3 months (1st Jan – 31st
Mar 2012) are performed to reduce the influence of the initial condition,
which is set constant and equal to a pressure head of − 100 cm. Water
flow from 1st Apr to 30th Sept 2012 is calculated to generate the syn­thetic data used in the following inverse modeling scenarios. In partic­
ular, simulated daily seepage outflow and volumetric water contents atthree different depths (i.e., z = [− 10, − 45, − 130] cm) are retrieved andcombined in different modeling scenarios.
Forward simulation results are perturbed with a normally distributed
measurement error with zero mean and standard deviation σ. We usedσ θ = 1.1 × 10− 4 cm3cm− 3 and σq = 0.01 cm for the water content andcumulative outflow, respectively (Schelle et al., 2012). Synthetic per­turbed data are then combined in three different modeling scenarios

Table 1
Formulas and properties of the test functions.
ID

Function

F1

Rotated High
Conditioned
Elliptic
Function

F2

F3

Shifted andRotatedRosenbrock’s
Function

Shifted andRotated
Rastrigin’s
Function

Formulaf1 (x) =Properties

D
∑
i=1

i− 1
(106 )D − 1 x2i

F1 (x) = f1 (M(x − o1 ))f2 (x) =D− 1∑
2[100(x2i − xi+1 ) + (xi − 1)2 ]
i=1

( (
)
)
2.048(x − o2 )F2 (x) = f2 M
+1100

f3 (x) =D∑[y2i − 10 cos(2πyi ) + 10]
i=1

( (
))
5.12(x − o3 )
F3 (x) = f3 M
100

- Unimodal- NonSeparable- Quadratic illconditioned- Multimodal- NonSeparable- Having anarrow valley
from a local
optimum to a
globaloptimum- Multimodal- NonSeparable- The numberof local
optima ishugesolver and overcome the usability limitations of the external couplingroutines. Furthermore, the addition of a few new parameters will not
require extensive changes to the HYDRUS graphical user interface.
The new G-CLPSO algorithm has been implemented in a separate
FORTRAN file (GCLPSO.for) to preserve the original HYDRUS code’s
structure as much as possible. The file contains only four subroutines:
GCLPSO, GRADIENT, TOURNAMENT, and FITNESS. The first two
names are self-explanatory, and the subroutines contain the instructions
to run the GPSO and ML algorithms as described above in Sections 2.1.1
and 2.1.2. TOURNAMENT includes the CLPSO learning strategy, while
FITNESS is a subroutine used to call the HYDRUS model and calculate
SSQ. The main file is modified to load the G-CLPSO parameters from anexternal input file (i.e., FIT.in), which are then passed to the GCLPSO
subroutine. At this stage of development, an integer, iOptimization,
loaded from an external input file, is used to choose between the original
HYDRUS inverse solver (i.e., iOptimization = 0) or the new G-CLPSO
algorithm (i.e., iOptimization = 1). Modifications of the Graphical User
Interface are currently ongoing to facilitate the input. At each G-CLPSO
iteration k, the number of HYDRUS calls, the lowest SSQ calculated thus
far, and its corresponding best position in the parameter space are savedin an external output file (i.e., GCLPSO.out).

Fig. 2. The fitness landscapes of the three test functions described in Table 1.
5

G. Brunetti et al.

Environmental Modelling and Software 150 (2022) 105341

Table 2
Soil hydraulic parameters and their search bounds in the inverse estimation.
Parameter

Parameter description

Soil Layer 1 (0 < z < -25 cm; Loam)θr1
Residual water content[cm3cm¡3]
θs1
Saturated water content[cm3cm¡3]
α1 [cm¡1]
VGM shape parametern1 [-]VGM shape parameter
Ks1 [cm
Saturated hydraulicconductivityday¡1]
Soil Layer 2 (-25 < z < -55 cm; Silt loam)θr2
Residual water content[cm3cm¡3]
θs2
Saturated water content[cm3cm¡3]¡1
α2 [cm ]VGM shape parameter
VGM shape parametern2 [-]
Ks2 [cm
Saturated hydraulicday¡1]
conductivity
Soil Layer 3 (-55 < z < -150 cm; Silt)θr3
Residual water content[cm3cm¡3]θs3
Saturated water content[cm3cm¡3]
α3 [cm¡1]
VGM shape parametern3 [-]VGM shape parameterKs3 [cm
Saturated hydraulicday¡1]
conductivity

Original
parameters

Table 3
Termination criteria for G-CLPSO, ML, PEST, and SCE-UA in scenarios M1, M2,
and M3.

Parameter
bounds

0.07

0.07

0.44

(0.4, 0.6)

0.036
1.56
24.96

(0.001, 0.1)
(1.1, 2.5)
(1, 100)

0.15

0.15

0.45

(0.4, 0.6)

0.02
1.41
10.8

(0.001, 0.1)
(1.1, 2.5)
(1, 100)

0.07

0.07

0.46

(0.4, 0.6)

0.016
1.37
6.00

(0.001, 0.1)
(1.1, 2.5)
(1, 100)

G-CLPSOaML

PEST

SCE-UA

Maximum
number of
HYDRUS
executions =10,000
Minimum
relativedifference
between best
and worst SSQ
value in theswarm = 0.01
Minimum
relative SSQ
reduction
betweenconsecutive
iterations =0.005

Maximum
number of
iterations = 100

Maximum
number of
iterations = 100

Minimum relative
SSQ reduction
betweenconsecutive
iterations = 0.005

Minimum relative
SSQ reduction
betweenconsecutive
iterations = 0.005

Maximum
number of
HYDRUS
executions =10,000
Minimum relative
SSQ reduction
betweenconsecutive
iterations = 0.005

Maximum
number ofconsecutive
iterations forwhich theminimum
reduction isachieved = 5
Maximum
number ofconsecutive
iterations forwhich there is no
SSQ reduction =5

Maximum
number ofconsecutive
iterations forwhich theminimum
reduction isachieved = 5
Maximum
number ofconsecutive
iterations forwhich there is no
SSQ reduction =5

Maximum
number ofconsecutive
iterations forwhich theminimum
reduction isachieved = 5
Maximum
number ofconsecutive
iterations forwhich there is
no SSQ
reduction = 5

focused on the inverse estimation of the saturated water content andconductivity, θs and Ks, respectively, and the VGM shape parameters, αand n, for the three soil layers (i.e., a total of 12 soil hydraulic param­
eters). Scenarios M1 and M2 use separately synthetic seepage outflow
and volumetric water contents at three different depths, respectively,
while scenario M3 combines the two time series.

Maximum
number ofconsecutive
iterations forwhich theminimum
reduction isachieved = 5a
For G-CLPSO termination, only iterations encompassing a gradient-based
search are considered.

2.2.4. Algorithms’ settings
The PEST software (White et al., 2020) is based on a
Gauss-Levenberg-Marquardt algorithm similar to that implemented in
HYDRUS, which calculates the gradient through numerical differentia­
tion. The main differences are the treatment of the parameter λ and the
PEST use of regularization to restrict the parameter search to identifi­
able parameters either by adding additional constraints to the parame­ters (Tikhonov Regularization) or separating identifiable parametersfrom non-identifiable parameters (Subspace Regularization). In thepresent study, PEST with subspace regularization based on singularvalue decomposition is used. Optimization control parameters are set
according to the PEST documentation (PEST, 2021). The PEST software
suite also implements the global optimizers SCE-UA, which uses thesame instructions and control files. Therefore, this implementation isapplied in the present study. Termination criteria for all algorithms arereported in Table 3.
PEST is connected to a Python (Python Software Foundation, 2013)batch file that externally executes HYDRUS and checks its convergence.
If the simulation is not convergent, HYDRUS output files read by PEST
are replaced with existing similar files having all outputs of interest set
to a very large value. This will lead to a very high calculated SSQ value
for this simulation, which can slow the calculations but avoids theirinterruption. To conduct a fair comparison between multiple optimiza­tion strategies, we compare them by averaging the results of ten inde­
pendent optimizations (e.g., Duan et al., 1994; Raj Shrestha and Rode,
2008; Sorooshian et al., 1993). In particular, the local search methods (i.
e., ML and PEST) are started at ten different random points in the search
space to resemble the initialization of the particles in the G-CLPSO,
which are randomly spread in the parameters’ bounds. Similarly, thestochastic SCE-UA is run with ten different seeds.

3. Results and discussion
3.1. G-CLPSO: sensitivity analysis
Fig. 3 shows the mutual influence of the swarm size, S, and local
search frequency, NL, on the normalized final SSQ and FES (Function
Evaluations) for each 10D test function. Results indicate that NL has asmall effect on the final fitness value (SSQ) for both the Elliptic andRosenbrock functions, while it influences more the outcomes for the
Rastrigin problem. The combined effect of S and NL on the latter isnonlinear, suggesting that specific parameter combinations can yield
better results for multimodal problems. However, in general, a more
frequent local search (1 < NL < 10) leads to better results for both
unimodal and multimodal problems independently from the population
size.
This also holds for the total number of function evaluations (FES),
though the behavior is more variegated among different functions. A
frequent local (ML) search significantly accelerates the convergence of
the algorithm for the Rosenbrock problem, while its effect is morenuanced for the Elliptic and Rastrigin functions. This is intuitive sincefrequently initializing the ML solver at different points in the curved
narrow valley of the Rosenbrock function improves the swarm’s
convergence towards the global optimum. This behavior is less criticalbut still appreciable for the unimodal elliptic function characterized by amore regular fitness landscape.
These findings are exacerbated when the dimensionality (30D) of the
problem increases (Fig. 4). On the one hand, a more frequent local
search is again enormously beneficial for both the Elliptic and Rose­
nbrock functions, which are only minimally affected by the swarm size.
6

G. Brunetti et al.

Environmental Modelling and Software 150 (2022) 105341

Fig. 3. Interpolated color maps showing the mutual influence of the swarm size, S, and iterations’ number, NL, on the normalized final SSQ (top) and FES (bottom)
for three analyzed functions (Elliptic (left), Rosenbrock (middle), and Rastrigin (right)) for each 10D test function. SSQ and FES are calculated by averaging theresults of five independent G-CLPSO optimizations. The ticks correspond to the values of S and NL used in the G-CLPSO.

On the other hand, the interaction between S and NL is nonlinear for themultimodal Rastrigin problem.
The main findings of the sensitivity analysis are:not observe this problem in the G-CLPSO. This may be mainly related
to the effects of the ML solver on the search behavior of the particles
and the corresponding reduction in their roaming behavior. How­ever, no general conclusions can be drawn since the maximum
investigated swarm size was 40.
- The analysis suggests that performing the local search every 1 to 10
G-CLPSO iterations is beneficial in terms of convergence andcomputational cost for low- and moderately high-dimensional
unimodal and multimodal problems. We argue that this is again
related to the local search’s magnifying effect on selective pressure in
the tournament selection procedure of the CLPSO.

- It is not possible to identify a single combination S-NL that guarantees
the best performance for all problems. Conversely, the algorithm
consistently exhibits good performance for multiple parameters’
combinations for all benchmark functions, thus suggesting theoverall robustness of the optimization approach. The use of only two
parameters confirms the flexibility of the proposed search strategy.
- The algorithm shows limited sensitivity to the swarm size for lowand moderately high-dimensional unimodal and multimodal prob­
lems. Piotrowski et al. (2020) compared the CLPSO with other PSOs
and concluded that a swarm of 20 individuals is optimal for the
CLPSO. With a larger population, particles’ movement in the CLPSO
becomes too chaotic, reducing the convergence to any optima. We do

In the following comparison, S and NL in the G-CLPSO are set to 10
and 1, respectively. This choice is a compromise based on the sensitivity
analysis results and is intended to highlight the effect of a very frequentlocal search on the algorithm performance. The population size in theFig. 4. Interpolated color map showing the mutual influence of the swarm size, S, and iterations’ number, NL, on the normalized final SSQ (top) and FES (bottom) for
three analyzed functions (Elliptic (left), Rosenbrock (middle), and Rastrigin (right)) for each 30D test function. SSQ and FES are calculated by averaging the results offive independent G-CLPSO optimizations. The ticks correspond to the values of S and NL used in the G-CLPSO.
7

G. Brunetti et al.

Environmental Modelling and Software 150 (2022) 105341

CLPSO is set to 20, according to Piotrowski et al. (2020), who carried out
an extensive analysis on the effect of the population size in multiple
particle swarm optimization strategies.

Table 4Means and standard deviations of the final SSQ for each strategy (G-CLPSO,
CLPSIO, and ML) for three analyzed functions (F1, F2, and F1) and three hypo­thetical scenarios (M1, M2, and M3) calculated from 10 independent
optimizations.

3.2. Test functions: comparison between G-CLPSO and CLPSO
The convergence characteristics of the G-CLPSO (black) and CLPSO(grey) on the selected test functions in both 10 and 30 dimensions areshown in Fig. 5, which reports the fitness value vs. the number offunction evaluations for ten independent optimizations. The mean andthe standard deviation of the final fitness value for both competing al­
gorithms are reported in Table 4. Results indicate that the G-CLPSO
significantly outperforms the CLPSO in terms of both accuracy andconvergence speed for low- and moderately high-dimensional unimodal
problems (i.e., F1). The local search quickly focuses the search of theswarm on the global optimum basin and reduces the roaming behaviorof the CLPSO. This also happens, although to a lower extent, for themultimodal Rosenbrock problems, where more function evaluations areneeded to identify the global minimum. Nevertheless, the accuracy gainof the G-CLPSO against the CLPSO remains significant also for thisfunction (Table 4).
Similar behavior is encountered for the multimodal Rastrigin prob­lem, for which the G-CLPSO identifies a better minimum. Interestingly,
the use of the local search does not hinder the exploration of the in­dividuals in the early stage of the optimization, while it focuses theirsearch towards the end when the global optimum basin is likely to be
identified. This is a desirable property to avoid premature convergencein local minima. Finally, Table 4 confirms the better performance of the
G-CLPSO for all problems, thus indicating good flexibility of the devel­
oped strategy over a wide range of potential applications.

Function/
Scenario

G-CLPSO

F1 (10D)

0 (0)aF2 (10D)

0 (0)F3 (10D)

3.14 (1.27)

F1 (30D)

0 (0)F2 (30D)

0 (0)F3 (30D)
M1

48 (10.15)
4.5e-3(0.9e-3)
3.2e-5(2.02e-5)
3.1e-3(1.3e-3)

M2
M3
a

CLPSO

ML

PEST

SCE-UA

–

–

–

–

–

–

–

–

–

–

–

–

–

–

–

–
0.45
(1.04)
0.31(0.54)
0.65(0.62)–
0.54(1.12)
0.02(0.06)
0.75(1.31)–
0.017(3.5e-3)
5.9e-3(0.01)
0.016(0.01)

Mean SSQ (standard deviation)
0.05(0.02)
3.28(2.21)
21.02(2.89)
2.5e-4(2.8e-4)
101(14.85)112 (7.58)–––

Values lower than 1e-10 are reported as zero.

value vs. the number of HYDRUS evaluations for ten independent op­
timizations. The mean and standard deviation of the final fitness value
for all competing algorithms are reported in Table 4. On the first in­spection, it is evident that the G-CLPSO outperforms other solvers interms of accuracy and convergence speed in all scenarios. The local
search strategies, ML and PEST, have similar behaviors and exhibit early
convergence in local minima due to their sensitivity to the initialization
point. Their accuracy is comparable (Table 4), with PEST outperforming
ML in scenario M2. Conversely, ML is generally less computationallyexpensive than PEST, and this can be related to the scaling of theparameter λ during the optimization process.
The global stochastic optimization strategy SCE-UA ranks second interms of accuracy (Table 4) but exhibits slow convergence (Fig. 6).

3.3. Synthetic modeling scenarios: Comparison between G-CLPSO, ML,
PEST, and SCE-UA
The convergence characteristics of the ML (grey), PEST (red), SCEUA (blue), and G-CLPSO (black) for the three inverse modeling sce­narios (M1, M2, and M3) are shown in Fig. 6, which reports the fitness

Fig. 5. Convergence analysis of the G-CLPSO (black) and CLPSO (grey) for three analyzed functions (Elliptic (left), Rosenbrock (middle), and Rastrigin (right)) for10D (top) and 30D (bottom) test functions.
8G. Brunetti et al.

Environmental Modelling and Software 150 (2022) 105341

Fig. 6. Convergence analysis of the ML (grey), PEST (red), SCE-UA (blue) and G-CLPSO (black) for the three inverse modeling scenarios (M1, M2, and M3).

Multiple optimization runs were terminated because the algorithm
exceeded the maximum number of allowed model executions (i.e.,
10,000). The algorithm can consistently identify the global optimum
basin, but then a significant number of model executions is required toachieve small accuracy gains. This behavior was already observed byWöhling et al. (2008), and it is common for global optimization
methods, which generally lack exploitation capabilities. Compared to
ML and PEST, SCE-UA has a much higher computational cost, hindering
its applicability when the model execution time increases (e.g., inreactive solute transport models).
On the other hand, the G-CLPSO always reaches the global optimum
region in a reasonable amount of function evaluations. In particular, thenumber of HYDRUS model executions is, on average, only slightly
higher than for ML and PEST and significantly lower than for SCE-UA.
This is extremely important since the computational cost is generally abottleneck in the calibration of environmental models, especially whenusing global optimization strategies. Table 4 confirms the superior
performance of the G-CLPSO, the convergence characteristics of which
exhibit a recurrent pattern characterized by early exploration followedby an exploitation phase induced by the local search. After this, minor
improvements are achieved, suggesting that the optimal region has beenreached and that this region is likely to be flat and similar to the Rose­
nbrock fitness landscape (Fig. 2). This shape of the objective function isfrequent in inverse hydrological problems (Abbaspour et al., 2004;
Šimůnek and Van Genuchten, 1996). As Fig. 5 shows, more iterations areneeded to identify precisely the global minimum in these cases. Morestringent convergence criteria would have further reduced the spread of
the final optimized value among different runs. However, this is notparticularly important for practical applications since the modeler ismore interested in knowing the fitness landscape’s shape near the op­
timum than a single value characterizing the parameter’s uncertainty.
This is automatically done in the G-CLPSO, which calculates the pa­
rameters’ confidence interval around the final global optimum.
3.4. Strengths and limitations
We consider this study to be innovative in three main aspects:- Performance: The analysis demonstrates that the newly developed
strategy outperforms both its “parent” algorithms, the CLPSO andML, and the gradient-based and stochastic search strategies PEST and
SCE-UA, respectively, in terms of accuracy and convergence speed.
9G. Brunetti et al.

Environmental Modelling and Software 150 (2022) 105341

In particular, it preserves the ML solver’s exploitation and robustness
without sacrificing the global search capabilities of the CLPSO. This
allows the algorithm to be successfully used in moderately highdimensional unimodal and multimodal problems frequently
encountered in vadose zone hydrology and environmental sciences.
- Computational cost: The G-CLPSO does not replace the well-tested androbust ML solver but extends its capabilities without dramaticallyincreasing the run time. This is crucial for several applications, but
especially for multi-dimensional problems and problems involving
reactive solute transport, for which high computational costs makeglobal search strategies impractical.
- Ease of implementation and coupling: The G-CLPSO needs only a few
additional parameters (i.e., S and NL) compared to the original
HYDRUS local search solver. This makes it easy to use and simplifies
its implementation into the HYDRUS graphical user interface, whichmakes the ML algorithm to be favored by HYDRUS users due to its
internal coupling into the HYDRUS source code.

the work reported in this paper.
Acknowledgments
Partial financial support from the Austrian Academy of Sciences(ÖAW) is acknowledged for funding the project RechAUT (Variability of
Groundwater Recharge and its Implication for Sustainable Land Use inAustria) within the Earth System Sciences research programme. A Py­
thon implementation of the developed algorithm is available uponrequest from the corresponding author.
References
Abbaspour, K.C., Johnson, C.A., van Genuchten, M.T., 2004. Estimating uncertain flow
and transport parameters using a sequential uncertainty fitting procedure. Vadose
Zone J. 3, 1340–1352. https://doi.org/10.2113/3.4.1340.
Abbaspour, K.C., Schulin, R., van Genuchten, M.T., 2001. Estimating unsaturated soil
hydraulic parameters using ant colony optimization. Adv. Water Resour. 24,
827–841. https://doi.org/10.1016/S0309-1708(01)00018-5.
Arora, B., Mohanty, B.P., McGuire, J.T., 2012. Uncertainty in dual permeability model
parameters for structured soils. Water Resour. Res. 48 https://doi.org/10.1029/
2011WR010500.
Brunetti, G., Porti, M., Piro, P., 2018. Multi-level numerical and statistical analysis of thehygrothermal behavior of a non-vegetated green roof in a mediterranean climate.
Appl. Energy 221, 204–219. https://doi.org/10.1016/J.APENERGY.2018.03.190.
Brunetti, G., Šimůnek, J., Piro, P., 2016. A comprehensive numerical analysis of thehydraulic behavior of a permeable pavement. J. Hydrol. 540, 1146–1161. https://doi.org/10.1016/j.jhydrol.2016.07.030.
Cao, Y., Zhang, H., Li, W., Zhou, M., Zhang, Y., Chaovalitwongse, W.A., 2019.
Comprehensive learning particle swarm optimization algorithm with local search formultimodal functions. IEEE Trans. Evol. Comput. https://doi.org/10.1109/
TEVC.2018.2885075.
Carsel, R.F., Parrish, R.S., 1988. Developing joint probability distributions of soil water
retention characteristics. Water Resour. Res. 24, 755–769. https://doi.org/10.1029/
WR024i005p00755.
Chu, H.-J., Chang, L.-C., 2009. Applying particle swarm optimization to parameter
estimation of the nonlinear Muskingum model. J. Hydrol. Eng. 14, 1024–1027.
https://doi.org/10.1061/(ASCE)HE.1943-5584.0000070.
Duan, Q., Sorooshian, S., Gupta, V.K., 1994. Optimal use of the SCE-UA global
optimization method for calibrating watershed models. J. Hydrol. 158, 265–284.
https://doi.org/10.1016/0022-1694(94)90057-4.
Duan, Q.Y., Gupta, V.K., Sorooshian, S., 1993. Shuffled complex evolution approach foreffective and efficient global minimization. J. Optim. Theor. Appl. 76, 501–521.
https://doi.org/10.1007/BF00939380.
Durner, W., Jansen, U., Iden, S.C., 2008. Effective hydraulic properties of layered soils at
the lysimeter scale determined by inverse modelling. Eur. J. Soil Sci. https://doi.org/
10.1111/j.1365-2389.2007.00972.x.
Feddes, R.A., Kowalik, P.J., Zaradny, H., 1978. Simulation of Field Water Use and Crop
Yield. PUDOC, Wageningen.
Han, F., Liu, Q., 2015. An improved hybrid PSO based on ARPSO and the quasi-Newton
method. In: Lecture Notes in Computer Science (Including Subseries Lecture Notes in
Artificial Intelligence and Lecture Notes in Bioinformatics. https://doi.org/10.1007/
978-3-319-20466-6_48.
Hargreaves, G.H., 1994. Defining and using reference evapotranspiration. J. Irrigat.
Drain. Eng. 120, 6. https://doi.org/10.1061/(asce)0733-9437, 1132, 1994.
Herndl, M., Pötsch, E.M., Bohner, A., Kandolf, M., 2011. Lysimeter als Bestandteil eines
technischen Versuchskonzeptes zur Simulation der Erderwärmung im Grünland. In:
14. Gumpensteiner Lysimetertagung “Lysimeter in Der Klimafolgenforschung Und
Wasserwirtschaft. ” LFZ Raumberg-Gumpenstein, pp. 119–126.
Ines, A.V.M., Droogers, P., 2002. Inverse modelling in estimating soil hydraulicfunctions: a Genetic Algorithm approach. Hydrol. Earth Syst. Sci. Discuss. 6, 49–66.
Kennedy, J., Eberhart, R., 1995. Particle swarm optimization. In: Proceedings of ICNN’95- International Conference on Neural Networks. IEEE. https://doi.org/10.1109/
ICNN.1995.488968, 1942–1948.
Kool, J.B., Parker, J.C., 1988. Analysis of the inverse problem for transient unsaturated
flow. Water Resour. Res. https://doi.org/10.1029/WR024i006p00817.
Levenberg, K., 1944. A method for the solution of certain non-linear problems in leastsquares. Q. Appl. Math. https://doi.org/10.1090/qam/10666.
Liang, J.J., Qin, A.K., Suganthan, P.N., Baskar, S., 2006. Comprehensive learning particle
swarm optimizer for global optimization of multimodal functions. IEEE Trans. Evol.
Comput. 10, 281–295. https://doi.org/10.1109/TEVC.2005.857610.
Liang, J.J., Qu, B.Y., Suganthan, P.N., 2014. Problem definitions and evaluation criteria
for the CEC 2014 special session and competition on single objective real-parameter
numerical optimization, 2014 IEEE Congr. Evolu. Comput.
Marquardt, D.W., 1963. An algorithm for least-squares estimation of nonlinear
parameters. J. Soc. Ind. Appl. Math. 11, 431–441. https://doi.org/10.1137/
0111030.
Noel, M.M., 2012. A new gradient based particle swarm optimization algorithm foraccurate computation of global minimum. Appl. Soft Comput. J. https://doi.org/
10.1016/j.asoc.2011.08.037.

Despite these major benefits, there are still some aspects that must be
further investigated:- The synthetic lysimeter experiment is limited and does not represent
other inverse modeling scenarios frequently encountered in vadosezone hydrology. Therefore, the algorithm should be further tested onsynthetic and real case studies to identify potential issues and thenresolve them.
- The effect of the local search on selective pressure in the tournament
procedure is empirically hypothesized but not mathematicallyexplained. More numerical tests are needed to fully understand theinteraction between the local and global search, particularly how therefreshing gap m of the CLPSO might interact with NL.
4. Conclusions
This study’s main aim was to extend the capabilities of the HYDRUS
inverse solver to moderately high-dimensional multimodal calibration
problems by combining its exploitation capability with the global searchcharacteristics of the PSO. To this aim, we developed a new loosely
coupled strategy, the G-CLPSO, which couples the CLPSO and ML al­
gorithms. The G-CLPSO has been internally coupled into the HYDRUS
source code, facilitating its use and implementation in the HYDRUS
graphical user interface. Multiple numerical and synthetic tests
confirmed that the hybrid strategy outperforms both the CLPSO and ML
on unimodal and multimodal problems by appropriately balancing
exploration and exploitation. Furthermore, synthetic modeling sce­narios on the inverse estimation of soil hydraulic properties showed GCLPSO outperforming the original HYDRUS ML solver, the gradientbased strategy PEST, and the stochastic global optimization method
SCE-UA in terms of accuracy and convergence speed.
Despite the promising results presented in this study, more researchis needed to test the developed optimization strategy further, better
understand its structure, and improve its performance. While this is
generally accomplished by using a probabilistic framework, the inclu­
sion of an approximate appraisal of the parameters’ uncertainty in the
numerical framework is desirable to better characterize the model pre­dictive performance. Nevertheless, this represents an important step
toward efficiently embodying global search strategies in vadose zone
model calibration. Finally, we would like to emphasize that while wehave developed the G-CLPSO strategy for the HYDRUS model, the al­
gorithm can be easily adapted for other hydrological or environmental
models.
Declaration of competing interest
The authors declare that they have no known competing financial
interests or personal relationships that could have appeared to influence
10

G. Brunetti et al.

Environmental Modelling and Software 150 (2022) 105341

Pan, L., Wu, L., 1998. A hybrid global optimization method for inverse estimation of
hydraulic parameters: annealing-simplex method. Water Resour. Res. 34,
2261–2269. https://doi.org/10.1029/98WR01672.
PEST, 2021. PEST - Model-independent Parameter Estimation User Manual.
Piotrowski, A.P., Napiorkowski, J.J., Piotrowska, A.E., 2020. Population size in particle
swarm optimization. Swarm Evol. Comput. 58, 100718. https://doi.org/10.1016/j.
swevo.2020.100718.
Python Software Foundation, 2013. Python Language Reference, Version 2.7. Python
Softw. Found.
Raj Shrestha, R., Rode, M., 2008. Multi-objective calibration and fuzzy preference
selection of a distributed hydrological model. Environ. Model. Software 23,1384–1395. https://doi.org/10.1016/j.envsoft.2008.04.001.
Richards, L.A., 1931. Capillary conductions of liquids through porous mediums. Physics
1, 318. https://doi.org/10.1063/1.1745010.
Ritter, A., Hupet, F., Muñoz-Carpena, R., Lambot, S., Vanclooster, M., 2003. Using
inverse methods for estimating soil hydraulic properties from field data as analternative to direct methods. Agric. Water Manag. https://doi.org/10.1016/S03783774(02)00160-9.
Roy, P.C., Guber, A., Abouali, M., Nejadhashemi, A.P., Deb, K., Smucker, A.J.M., 2019.
Simulation optimization of water usage and crop yield using precision irrigation. In:
Lecture Notes in Computer Science (Including Subseries Lecture Notes in Artificial
Intelligence and Lecture Notes in Bioinformatics). Springer Verlag, pp. 695–706.
https://doi.org/10.1007/978-3-030-12598-1_55.
Schelle, H., Iden, S.C., Fank, J., Durner, W., 2012. Inverse Estimation of Soil Hydraulicand Root Distribution Parameters from Lysimeter Data. Vadose Zo. J. https://doi.
org/10.2136/vzj2011.0169.
Sedaghatdoost, A., Ebrahimian, H., Liaghat, A., 2019. An inverse modeling approach tocalibrate parameters for a drainage model with two optimization algorithms on
homogeneous/Heterogeneous soil. Water Resour. Manag. 33, 1383–1395. https://doi.org/10.1007/s11269-019-2191-x.
Shi, Y., Eberhart, R., 1998. A modified particle swarm optimizer. In: 1998 IEEE Int. Conf.
Evol. Comput. Proceedings. IEEE World Congr. Comput. Intell. (Cat. No.98TH8360)
69–73. https://doi.org/10.1109/ICEC.1998.699146.
Šimůnek, J., Hopmans, J.W., 2002. 1.7 parameter optimization and nonlinear fitting. In:
Methods of Soil Analysis, Part 1, Physical Methods, pp. 139–157. https://doi.org/
10.2136/sssabookser5.4.c7.
Šimůnek, J., Van Genuchten, M.T., 1996. Estimating unsaturated soil hydraulicproperties from tension disc infiltrometer data by numerical inversion. Water
Resour. Res. https://doi.org/10.1029/96WR01525.
Šimůnek, J., van Genuchten, M.T., Šejna, M., 2016. Recent developments andapplications of the HYDRUS computer software packages. Vadose Zone J. 15, 25.
https://doi.org/10.2136/vzj2016.04.0033.
Šimůnek, J., van Genuchten, M.T., Wendroth, O., 1998. Parameter estimation analysis of
the evaporation method for determining soil hydraulic properties. Soil Sci. Soc. Am.
J. 62, 894. https://doi.org/10.2136/sssaj1998.03615995006200040007x.
Sorooshian, S., Duan, Q., Gupta, V.K., 1993. Calibration of rainfall-runoff models:application of global optimization to the sacramento soil moisture accounting model.
Water Resour. Res. 29, 1185–1194. https://doi.org/10.1029/92WR02617.

Stanovov, V., Akhmedova, S., Semenkin, E., 2019. Selective Pressure Strategy indifferential evolution: exploitation improvement in solving global optimization
problems. Swarm Evol. Comput. https://doi.org/10.1016/j.swevo.2018.10.014.
Sutanto, S.J., Wenninger, J., Coenders-Gerrits, A.M.J., Uhlenbrook, S., 2012. Partitioningof evaporation into transpiration, soil evaporation and interception: a comparison
between isotope measurements and a HYDRUS-1D model. Hydrol. Earth Syst. Sci.
16, 2605–2616. https://doi.org/10.5194/hess-16-2605-2012.
Sutton, A.M., Lunacek, M., Whitley, L.D., 2007. Differential evolution and nonseparability: using selective pressure to focus search. In: Proceedings of GECCO
2007: Genetic and Evolutionary Computation Conference. https://doi.org/10.1145/
1276958.1277221.
Ter Braak, C.J.F., Vrugt, J.A., 2008. Differential Evolution Markov chain with snooker
updater and fewer chains. Stat. Comput. 18, 435–446. https://doi.org/10.1007/
s11222-008-9104-9.
Transtrum, M.K., Sethna, J.P., 2012. Improvements to the Levenberg-Marquardt
Algorithm for Nonlinear Least-Squares Minimization.
van Dam, J.C., Stricker, J.N.M., Droogers, P., 1994. Inverse method to determine soil
hydraulic functions from multistep outflow experiments. Soil Sci. Soc. Am. J.
https://doi.org/10.2136/sssaj1994.03615995005800030002x.
van Genuchten, M.T., 1980. A closed-form equation for predicting the hydraulicconductivity of unsaturated Soils1. Soil Sci. Soc. Am. J. 44, 892–898. https://doi.
org/10.2136/sssaj1980.03615995004400050002x.
Vrugt, J.A., Gupta, H.V., Bouten, W., Sorooshian, S., 2003. A Shuffled Complex Evolution
Metropolis algorithm for optimization and uncertainty assessment of hydrologic
model parameters. Water Resour. Res. 39 https://doi.org/10.1029/2002WR001642
n/a-n/a.
Vrugt, J.A., Robinson, B.A., 2007. Improved evolutionary optimization from genetically
adaptive multimethod search. Proc. Natl. Acad. Sci. Unit. States Am. 104, 708–711.
https://doi.org/10.1073/pnas.0610471104.
Watson, T.A., Doherty, J.E., Christensen, S., 2013. Parameter and predictive outcomes ofmodel simplification. Water Resour. Res. 49, 3952–3977. https://doi.org/10.1002/
wrcr.20145.
White, J.T., Hunt, R.J., Fienen, M.N., Doherty, J.E., 2020. Approaches to Highly
Parameterized Inversion: PEST++ Version 5, a Software Suite for Parameter
Estimation, Uncertainty Analysis, Management Optimization and Sensitivity
Analysis, Techniques and Methods. https://doi.org/10.3133/tm7c26.
Wöhling, T., Vrugt, J.A., Barkle, G.F., 2008. Comparison of three multiobjective
optimization algorithms for inverse modeling of vadose zone hydraulic properties.
Soil Sci. Soc. Am. J. 72, 305–319. https://doi.org/10.2136/sssaj2007.0176.
Wolpert, D.H., Macready, W.G., 1997. No free lunch theorems for optimization. IEEE
Trans. Evol. Comput. https://doi.org/10.1109/4235.585893.
Xi, M., Lu, D., Gui, D., Qi, Z., Zhang, G., 2017. Calibration of an agricultural-hydrologicalmodel (RZWQM2) using surrogate global optimization. J. Hydrol. 544, 456–466.
https://doi.org/10.1016/J.JHYDROL.2016.11.051.
Zhao, S.Z., Liang, J.J., Suganthan, P.N., Tasgetiren, M.F., 2008. Dynamic multi-swarm
particle swarm optimizer with local search for large scale global optimization. In:
2008 IEEE Congress on Evolutionary Computation, CEC 2008. https://doi.org/
10.1109/CEC.2008.4631320.
Zielinski, K., Laur, R., 2007. Stopping criteria for a constrained single-objective particle
swarm optimization algorithm. Inform. https://doi.org/10.31449/inf.v31i1.127.

11