VOL. 13, NO. 12, JUNE 2018ISSN 1819-6608

ARPN Journal of Engineering and Applied Sciences©2006-2018 Asian Research Publishing Network (ARPN). All rights reserved.

www.arpnjournals.com

GLOBAL SENSITIVITY ANALYSIS, INVERSE MODELLINON
SOIL WATER INFILTRATION
Goh Eng Giap1, Kosuke Noborio2 and Asmadi Ali1
1

School of Ocean Engineering, Universiti Malaysia Terengganu, Kuala Terengganu, Terengganu, Malaysia
2
School of Agriculture, Meiji University, Higashimita, Tamaku, Kawasaki, Japan
E-Mail: sunnygoh@gmail.com

ABSTRACT
Richards’ equation was used to govern water infiltration into the unsaturated soil. Haverkamp’s hydraulicproperties were used to govern the matric pressure head and hydraulic conductivity relations. Sobol’ variance-based
method as a global sensitivity analysis tool was coded to study sensitivity and uncertainty values of input parameters. The
same input parameters uncertainty values used in Sobol’ variance-based method could be subjected to an inverse method tosearch for a set of globally optimized input parameters. Therefore, a simple Genetic Algorithm, i.e. an inverse method, wascoded and used to improve the current framework for the uncertainty analysis. A simple significant digit approximation
method was developed, applied on the input parameters, and the generated uncertainty parameter values were found tocause 15.4 % of error, which could be reduced if the number of input parameter significant digit used in the simulation
were increased. The simple Genetic Algorithm was able to search for a set of globally optimized input parameters at a low
population size, i.e. 50, on three cases tested and also was found able to handle noise data. Data points at the lower andupper plains of the water infiltration front are necessary for the inverse method.
Keywords: Richards’ equation sobol’, variance-based method, simple genetic algorithm, significant digit approximation.

1. INTRODUCTION
Traditionally, a sensitivity analysis was carriedout by changing a single input parameter while keeping
other parameters at constant values. It is better known to
the statistical community as one-at-a-time (OAT)
procedure. A better method is to simultaneously changing
more than one parameter at a time known as global
sensitivity analysis (GSA). It is because the simultaneous
variation of few input parameters at a time could result insensitivity coefficient that is different from adding up theentire varied individual OAT from base case sensitivity
coefficients.
GSA tools can be ranged from qualitative
screening method to quantitative technique (Sobol’ 1990;
Morris 1991; Campolongo et al. 2007).Sobol’ variancebased method is one of the GSA tools. It is useful forsensitivity analysis and also for uncertainty analysis. In theformer, it is used to identify the most and the leastsensitive parameters. While in the latter, it is used todetermine the greatest influence parameter for parameterprioritization and the least influential parameter forparameter fixing (Saltelli et al. 2008). Some models were
subjected to this method are SWAT (Nossent et al. 2011),
inhalation dose model (Avagliano and Parrella 2009),
CERES-EGC (Drouet et al. 2011) and SLAMM 5 (ChuAgor et al. 2011). Moreover, the variance-based method
has been listed in the US Environmental Protection
Agency list of attributes as preferred sensitivity analysis
method (EPA 2009).
The range of input parameters uncertainty used in
Sobol’ variance-based analysis could be further processedby an inverse method to determine an optimized set of
input parameters. In order to obtain globally optimized setof input parameters, a simple Genetic Algorithm (SGA)
was used for this purpose. SGA was introduced byHolland (1975). It is a class of evolutionary algorithms

that is a nature inspired computational method. Varietiesof GA methods were used in estimating hydraulicproperties in unsaturated soils. For instance, SGA used by
Schneider et al. (2013) and micro-GA used by Ines andDroogers (2002) and Shin et al. (2012). Inverse method
other than GA, e.g. Levenberg-Marquardt (L-M)
algorithm, was used by Kelleners et al. (2005) to predicthydraulic properties. Similarly, the L-M algorithm wasalso used by Mollerup et al. (2008). SGA is a global
search method, whereas L-M is a local search method. Itmeans that L-M could potentially found a local minimumthat is not a global minimum found by SGA. Global
minimum means the prediction is closer to experimental or
benchmark data than local minimum prediction.
Therefore, in the current work, it was limited to SGA.
Numerical simulation of water infiltration is not
only allowed prediction of water movement in theunsaturated zone, it has significant implication forcontaminant movement (Šimůnek and Bradford 2008).
Generally, Richards’ equation (Richards 1931)is used togoverned water infiltration in variably saturated soils. It
has been applied in various disciplines such as hydrology,
meteorology, agronomy, and environmental protection(Pachepsky et al. 2003). Mollerup et al. (2008) coupled
Philip’s semi-analytical solution of Richards’ equation
with an inverse method, i.e. L-M algorithm. In the current
studies, Philip’s semi-analytical solution from Haverkamp
et al. (1977) was used as benchmark data to verify model
simulation for SGA method. The published data fromHaverkamp is one of the resources available for modelers
to verify their newly developed simulation source code,e.g. Noborio et al. (1996), Fayer and Jones (1990), Fayer(2000) and Goh and Noborio (2013).Therefore, the current
studies are built on the research for the chosen case thathas been a common case study for modelers, and hence,
the finding from the present work would be a useful3895

VOL. 13, NO. 12, JUNE 2018ISSN 1819-6608

ARPN Journal of Engineering and Applied Sciences©2006-2018 Asian Research Publishing Network (ARPN). All rights reserved.

www.arpnjournals.com
reference for other researchers with interest in developing
model and couple with GSA and SGA methods on
Richards’ equation.
The aim of these studies is to manage uncertainty
input parameter values. Firstly, the water infiltration
model was subjected to Sobol’ variance-based method tostudy global sensitivity analysis. Secondly, the model was
run through SGA. While the former manages to determine
important parameters and also to identify influential
parameter for prioritization, the latter extends uncertainty
analysis to search for a globally optimum set of input
parameters that is best described the water infiltration
front. The first objective was achieved by three cases.
Case A was based on significant digit approximation thatwas discussed in Goh and Noborio (2013), and Cases B
and C were based on 10 and 20 % variation of input
parameters,
respectively.
Thesignificant
digit
approximation was to determine the effect of uncertainty
in model output due to the uncertainty values of input
parameters established from the published data. The model
output in all case studies is volumetric water content
value.In the SGA analysis, the Gaussian noise was
imposed on the benchmark data and the inverse method
was used to predict the volumetric water content, under
the uncertainty values of input parameters in Cases A, B,
and C. The extreme variation of input parameter
uncertainty values was also examined.
2. MATERIALS AND METHODS
2.1 The governing equation for water flow inunsaturated soil
The θL −based form used in the simulation
model is as following:
∂θL
∂

=

∂

∂

[ K

∂ψm ∂θL
⃑]
− Kk
∂θL
∂

(1)

where: θL is the volumetric water content(m3/m3); t is the time of simulation (s); z indicates thevertical distance of simulation (m); K is hydraulicconductivity of the medium (m/s);ψm is the matric⃑ is the vector unit with a value ofpressure head (m);kpositive one when it is vertically downwards. The 𝜕θL ⁄𝜕tterm indicates partial derivative of θL with respect to time,
while 𝜕ψm ⁄𝜕θL term is partial derivative of ψm withrespect to θL . The 𝜕θL ⁄𝜕z is the partial derivative of θLwith respect to vertical distance. The initial condition
volumetric water content value, in Table 1, indicates theinitial spatial distribution of water content. The value wasused along with upper and lower boundary conditions for
simulation. The upper boundary was fixed to a constantvolumetric water content given in Table-1. The lower
volumetric water content was set equal to adjacent cell.
The Equation.1, in general, is known as Richards’
equation (Richards 1931). The equation was numericallyapproximated using finite-difference method, and its
algebra was coded in FORTRAN 2008.
2.2 The constitutive functions of matric pressure headand hydraulic conductivity
In this study, the constitutive functions ofHaverkamp et al. (1977) were used:ψm = −

K=K

−

exp [

+ −

ψm

α θ −θ
θL −θ

− ]

β

(2)
(3)

where: , , A and B are fitting parameters; θ(m3/m3) is residual volumetric water content; θ (m3/m3) issaturated volumetric water content; and K (m/s) is
saturated hydraulic conductivity.
2.3 Numerical experiment, initial and boundary
conditions
Infiltration of Yolo light clay was used as ournumerical experiment, and the coefficient values (basecase) were listed in Table-1. Initial condition for θL was
0.2376 m3/m3. The Lower boundary was set permeable towater. The upper boundary was set at 0.495 m3/m3. Afterconsidering the mass balance ratio (MBR) and the numberof iteration, the simulation parameters such as time-step,
simulation time, convergence value (CV), and spatial
discretization size were set as 500 s, 105 s, 10-12 m3/m3 and1 cm, respectively.

3896

VOL. 13, NO. 12, JUNE 2018ISSN 1819-6608

ARPN Journal of Engineering and Applied Sciences©2006-2018 Asian Research Publishing Network (ARPN). All rights reserved.

www.arpnjournals.com
Table-1. Goes approximately here.
Parameter

Base Case

Case A

Case B

Case C739

738.5 – 739.499

739 – 812.9

739 – 886.8

3

0.124

0.1235 – 0.124499

0.124 – 0.1364

0.124 – 0.1488

θ (m /m )0.495

0.495 – 0.495499

0.495 – 0.5445

0.495 – 0.594

4

3.5 – 4.499

4 – 4.4

4.0 – 4.8

A

124.6

124.55 – 124.6499

124.6 – 137.06

124.6 – 149.52

1.77

K (cm/hr)

4.428x10-2

1.765 – 1.77499
4.4275x10-2 –4.428499x10-2
0.23755 –0.2376499

1.77 – 1.947
4.428x10-2 –4.8708x10-2

1.77 – 2.124
4.428x10-2 –5.3136x10-2

0.2376 – 0.26136

0.2376 – 0.28512

0.495

0.495

0.495

0.495

500

500

500

500

1

1

1

1

θ (m /m )
3
3

3

B

θL initial cond.
(m3/m3)
θL upper bound.
(m3/m3)∆t (s)

0.2376∆z (cm)

2.4 Sobol’ variance-based method as global sensitivity
analysis
Model output (Y) is a function of input
parameters (X , X , … . , X ). In a variance-based method,
total unconditional variance, V Y , is given by Sobol’(1990):V Y = ∑ V + ∑ ∑ >𝑖 V + ⋯ + V…

(4)where ∑ V is the sum of partial variances that
include individual effects of each input parameter,∑ ∑ > V , all the partial variances of two input
parameters interaction, and so on.
Dividing equation. 4 by total unconditionalvariance, V Y , the equation becomes,∑ S + ∑ ∑ >𝑖 S + ⋯ + S

…

=

(5)

where the first term on the left side of theequation S = V ⁄V Y is first order indices, S =V ⁄V Y is second order indices for the second term, andso on. The indices are scaled between 0 and 1. If there isno interaction between input parameters, ∑ S would beequal to unity, and if ∑ S is less than unity, it indicatesthe presence of interaction effects, which could be any
combination of indices that is greater than first order.
The quasi-Monte Carlo estimators compiled by
Saltelli et al. (2010) allowed estimation of first ordersensitivity index (S ) and total effect index (STi ) for input
parameter i. As an example, three input parameters model
has a total sensitivity index on input parameter that isgiven by:
ST = S + S

ST = S + S

+S

+S

+S

+S

(6)
(7)

ST = S + S

+S

+S

(8)

Since the quasi-Monte Carlo estimators canestimates S and STi , the difference between the two
estimators, for instance (ST − S ) would remove single
effect value and only higher than first order sensitivity
indices are left, which is usually used to indicates(4)input parameters.
interaction effects between
The first order sensitivity index and total effectindex were estimated by the following quasi-Monte Carloestimators (Saltelli et al. 2010):
S =
STi =

−

/ N ∑N
m=

/N ∑N
m=

−

−

m
i

m

/N ∑N
m=

/ N ∑N
m=

m

m

m

m
i

(9)

−fo

(10)

−fo

where y m , y m , and y m
are model outputs.
iEach of these model outputs was given by a set of input
parameter values by matrix of N × k. The k columns werereferred to the number of input parameters. Sobol’ quasirandom sequences were used to generate quasi-random
numbers for k-dimensional input parameters, where the kvalue used in this study was eight. The quasi-random
number generated was ranged between zero and unity, andit was scaled with respect to the uncertainty values inTable-1 for Cases A, B, and C. For each row of thematrix, e.g. matrix input parameter values ofy m , was acomplete set of input parameter values for an execution ofa model simulation. The number of model execution
depended on the N rows, where the N value used in thisstudy was fifteen thousands. Nossent et al. (2011) had
shown that an N value of twelve thousands for twenty-six
input parameters was sufficient to obtain reliable3897

VOL. 13, NO. 12, JUNE 2018ISSN 1819-6608

ARPN Journal of Engineering and Applied Sciences©2006-2018 Asian Research Publishing Network (ARPN). All rights reserved.

www.arpnjournals.com
estimation of sensitivity indices. Hence, the matrix input
parameter values of y m and y m would require 2N model
executions. The matrix input parameter values of y m
wasi
referred to parameteri, and it has N rows for eachparameter. Since there were k input parameters that must
be tested, and the model executions would be kN. The
matrix input parameter values of y m
for parameter i was
iobtained by using all the columns of matrix input
parameter values ofy m , except the column of parameteriwas taken from the matrix input parameter values of y m .
This step was repeated based on the number of input
parameter as indicated by the value k. Therefore, in thecurrent studies, the total model executions that werecarried out for each case (Cases A, B, and C) was given byN(2+k)=15,000(8+2)=150,000. All the cases combined
would require 450,000 model executions. The fo wasgiven by

m
/N ∑N
m= y

.

2.5 Simple genetic algorithm (SGA) as inverse method
The uncertainty values of input parameters wereused to generate a random population of initial solutions.
The solution is referring to a set of input parameter values
for use in model simulation. The random population was
generated using Sobol’ sequence to spread widely between
individual solutions in order to encompass every possiblecombination of input parameter values that could best
describe the model output. Vrugt et al. (2008) used Latin
hypercube sampling to generate the initial population.
Each solution of the population is known as chromosome,which is consisting of input parameter values that wereencoded into binary string of 1 and 0. Two initialchromosomes (parent) were randomly selected to go
through crossover process in order to generate two new
chromosomes (offspring). In this study, a single crossover
point was used (Figure-1).

Figure-1. Goes approximately here.
The offspring was decoded into a real value ofeach input parameter values, and it was then used as input
for the simulation and evaluated with an objectivefunction. The objective function values from the parent
and offspring were compared, and if offspring objectivefunction value is better than the parent, the offspring genes
will overwrite the parent genes. Similarly, under theprocess of mutation, two parents were selected, andrandom genes in the parent were selected for the mutation(Figure-2). They were then decoded and used as input forthe simulation and then evaluated with the objectivefunction. Genes from offspring with better objectivefunction value would replace the parent genetic values.

Figure-2. Goes approximately here.
The crossover and mutation processes would
continue based upon the preset number of population size.
A complete iteration of crossover and mutation would
signify a complete creation of a new generation where thesize of the initial population and the new generation are
the same.
3 RESULTS AND DISCUSSIONS
3.1 Simulation outcomes of Sobol’ variance-based
methodThere are various forms of an uncertainty of input
parameter distribution, for instance, normal, uniform,
triangular and lognormal. In this study, the uniform
distribution was used; the distribution was based on thesignificant digit approximation method (Case A), 10 and
20 % increment for Cases B and C, respectively, as in
Table-1.The uniform distribution used in the current
studies was similar to those used in Saltelli et al. (2004),
Vrugt et al. (2008), Drouet et al. (2011), Yang (2011), andYounes et al. (2013). The significant digit approximation
for Case A was developed from the base case in Table-1.
For instance, the residual volumetric water content, i.e.
0.124 m3/m3 as the base case, could be approximated bylower and upper bounds of 0.1235 and 0.124499 m3/m3,
respectively. Vice versa, the lower and upper bounds
could be reduced to three significant digits; it would returnto the original value as before. Hence, the lower and upperbounds were the minimum and maximum uncertainty
limits resulted from the given significant digit of an input
parameter. If the base case was to increase its significant
digit to 0.1240 m3/m3 the uncertainty range for theparameter would be bounded between 0.12395 and0.1240499 m3/m3, respectively, which has a lower
uncertainty range than the former example. All theparameters in Case A were developed from the base casevalues with the same method. Hence, the uncertainty
values of input parameters distribution on model output
due to the effect of the significant digit approximation wasable to be investigated.

3898

VOL. 13, NO. 12, JUNE 2018ISSN 1819-6608

ARPN Journal of Engineering and Applied Sciences©2006-2018 Asian Research Publishing Network (ARPN). All rights reserved.

www.arpnjournals.com
Table-2. Goes approximately here.
Parameter

θ
θ

A

B

K

θL initial cond.

Case A𝐢

0

Case B
0.00001𝐢

0.0008900.000020.22980

Case C
0.00130𝐢

00.0014300.0000500.00006

0.26600

0.37800

0.46200

0.152900.279140.517700.58400

0.01630

0.023800.00904

0.0268400.000100.003220.007680.00027

0.00855

0.18580

0.21300

0.30900

0.38900

0.204410.3317000.000020.009780.01460

0.005780.016490.000170.000190.19100

0.19400

0.480300.48563𝐢

Using the quasi-Monte Carlo estimators from
Equations. 9 and 10, the first order sensitivity index (S )and total effect index (STi ) were estimated for every input
parameter (Table-2). The model variability was calculated
by simply multiplying the first order sensitivity index andtotal effect index with one hundred. The results for Case A
showed that had the highest S that 51.77 % of the model
variability was due to this factor alone.It was followed byθ andB at 22.98 and 18.58 %, respectively. A minor
percentage was accounted by θL initialcond. at 0.017 %.
Insignificant contributions were found from parameters ,θ , A and K . The combined model parameters ,θ andBwere found responsible for 93.33 % of all themodel variability. By taking − ∑8= S , only 6.67 % was
due to the interaction between parameters for a
combination of indices higher than first order. The specificcontribution of each input parameter on interaction effect
was investigated through the use of STi − S . Theinteraction effect calculation was found similarly andmainly contributed by and then followed by θ andB at
6.63, 3.62 and 2.72 %, respectively. It is important to notethat the combined value of the three parameters was
greater than 6.67 %, as given by − ∑8= S , that was
found in the previous calculation. It can be justified by thatsome second and higher order sensitivity indices wererepeatedly calculated and accounted for as shown in theexample of Equations. 6-8. In addition, similar to S , theSTi for parameters , θ , A, K and θL initialcond. were
found to be insignificant.
While some input parameters were found to be
more significant than the others, it should be noted that
this was basedon the test case (base case) of theuncertainty values of input parameters that were given in
Table 1. Under different uncertainty values of input
parameters, the results of S and STi would be different. To
demonstrate our claim, the base case input parameters in
Table-1 were modified to have an increment of 10 % toform a new uniform distribution for each parameter, i.e.

𝐢𝐢Case B. In Table-2, the parameter
was no longer
providingdominant
influenceon
modelvariability.Parameters demonstrated significant influenceswereθ , B and θL initialcond. at 37.8, 30.9 and 19.1 % of
the model variability, respectively. Some characteristics
were remain unchanged, for instance, STi was greater thanS , and also parameters with greater S resulted in greater
STi . Since the results in Case B, as in Table-2, were basedon 10 % increase in value for all input parameters, inshort, it could be termed as a sensitivity analysis.
However, the results in Case A were based upon different
percentage variations on each input parameter andtherefore, it should be termed as uncertainty analysis. The
word uncertainty signifies that each time the sensitivity
indices are to be calculatedin different cases with different
uncertainty values for each input parameter. Hence, in
uncertainty analysis, the importance of parameters would
be differed from one case to the others. The consistency ofsensitivity analysis results was further tested in Case C
that all input parameter values were increased by 20 %.
The results of Case C showed that the model variability
was dominated by θL initialcond. , B, and θ indecreasing order. The parameter importance ranking
appeared to be different from those from Case B. The
discrepancy could be explained, refer to Figures 3(b) and3(c), that the interception point of three S input
parameters (θ , θL initialcond. , B) was located below
and above the horizontal line of 0.185 m depth
corresponding to Cases B and C. A general comparison
between first order sensitivities of Cases A, B, and C(Figures 3(a)-(c)) showed that there were similarities ofsensitivity index curvature for Cases B and C andirrespective of percentage increment applied, while Case
A was distinctively different. In a general perspective,
sensitivity analysis is a particular case of uncertainty
analysis, and it happened only when all parameters werevaried with an equal percentage.

3899

VOL. 13, NO. 12, JUNE 2018ISSN 1819-6608

ARPN Journal of Engineering and Applied Sciences©2006-2018 Asian Research Publishing Network (ARPN). All rights reserved.

www.arpnjournals.com

(a)

(c)

(b)

0.185 m0.185 m0.185 mFigure-3(a)-(c). Go approximately here.
In Case A, refer to Figure-3(a), the parameter
was the most important parameter where almost 100 % ofmodel variability could be explained by the parameter
alone, except at the near surface layer and the depth of
0.185 m. In Case B, Figure-3(b), θ was the single mostimportant parameter at near surface depth. It was gradually
decreased and followed by a steady increase of B and alsoa slow minor increase of θL initialcond. at adeeper depth
between 0.05 and 0.15 m. After 0.15 m depth, the increaseof θL initialcond. was exponential and the decreaseofθ was steep with increasing depth, while B gradually
increase to a peak before it decreases to zero after water
infiltration front. The upper depth volumetric watercontent was largely influenced by θ , and at the lower
depth it was entirely governed by θL initialcond. . A
similar observation was found in Case C and only that theincreased and decreased of parameters were found to takeplace earlier than Case B.
The simulation output was also processed forother statistical measures, Figures 4-6. The figure thatindicate the difference between mean and median (Figure4-6(d)) was found in close trend to that skewness inFigures 4-6(c). The mean value in Figures 4-6(a) appeared
to have surrounded by boundaries from 95 percentile
confident intervals, as inner boundaries, and Minimum,
Maximum values as outer boundaries. Case C has greater
variability of volumetric water content than Case B, andalso the variability of Case B was found higher than Case
A; as shown in Figures 4-6(b). The figure for standard
deviation has a similar trend as those observed in MaxMin and standard error of the mean. The standard error ofmean data was not shown here. At 0.185 m depth of CaseA, there was sudden decreased of uncertainty as indicatedin standard deviation, Max-Min and the 95 percentile
confident intervals. It was rather an uncommon event as it
was not evident in Cases B and C. At a distance before theend of water infiltration front, a peak of kurtosis andskewness was observed in all cases, except in Case A thepeak of skewness was rather subtle.At deeper depth ofwater infiltration front, i.e. some distances somewhere in
the middle before the end of water infiltration front andafter the rest of 0.4 m depth, kurtosis remained subtle, i.e.
approximate -1.2, signifying low and flat peak as shown inthe distribution curve in Figure-4(e). For Cases B and C,
the low and flat peak happened after the depth of 0.3 m.
Skewness and kurtosis were found to have reached its
peak after the peak of standard deviation, which suggestthat the distribution curves were initially horizontally
spread (as indicated by standard deviation) before skewed
right (or positive), and raised high frequency (as indicatedby kurtosis). The positive (or right) and negative (or left)
skewness values appeared to be rather intermittent, but it
became zero (no effect of skewness) at deeper depth after
the water infiltration front. A normal distribution was
rarely found, because the skewness and kurtosis, inFigures4-6(c), were consistently deviated from zero (orperfect normal distribution curve) as shown in Figures46(e). In addition, an apparent relation between kurtosis andsensitivity index was observed. In regions, i.e. depths,
where only single input parameter has dominant
influenced on simulation output, kurtosis value tended toapproximate equal or less than -1.0. Since the negative
kurtosis value indicates the appearance of small and flat
peak, a normal distribution with a unique peak could not
be found in these regions. Such observation was
predominantly observed in all depths of Case A, except atnear surface and areas surrounding 0.185 m depth;
whereas in Cases B and C, predominant regions weresimilarly observed at near surface and depth after water
infiltration front. The generated low and flat peak
distributions should be defined as uniform distribution
where
was found primarily responsible for Case A,
while θ and θL initialcond. were found responsible fornear surface and depth after water infiltration front,
respectively, for corresponding to Cases B and C.

3900

VOL. 13, NO. 12, JUNE 2018ISSN 1819-6608

ARPN Journal of Engineering and Applied Sciences©2006-2018 Asian Research Publishing Network (ARPN). All rights reserved.

www.arpnjournals.com
Volumetric water content(m3/m3)

(a)

0.20
0.00

0.30

0.40

0.50

(b) Std. Dev., or (c) Skewness, or (d)Max-Min0.00 0.03 0.06

Kurtosis

-2 -1

0

1

2

MeanMedian-4.E-03 1.E-03(e) Frequency
00.005 m0.045 m0.05
Depth (m)5500

0.085 m0.10

0.125 m0.15
0.165 m0.20

0.205 m0.245 m0.25

MeanMin, Max95 percentile

0.30

Std. Dev.
Max-Min

MeanMedian

Kurtosis
Skewness

0.285 mFigure-4(a-e). Go approximately here.

(a)Volumetric water content(m3/m3)

0.20
0.00

0.30

0.40

0.50

(b) Std. Dev., or (c) Skewness, or(d)

0.0 0.1 0.2 0.3-5.E-04 1.E-02

Max-Min

Kurtosis-2 -1 0 1 2 3

MeanMedian01500 3000
0.005 m0.05
Depth (m)(e) Frequency

0.045 m0.085 m0.10

0.125 m0.15
0.165 m0.20
0.25

0.205 m0.245 mMeanMin, Max

0.30

95 percentile

Std. Dev.
Max-Min

Kurtosis
Skewness

MeanMedian

0.285 mFigure-5. Go approximately here.

3901

VOL. 13, NO. 12, JUNE 2018ISSN 1819-6608

ARPN Journal of Engineering and Applied Sciences©2006-2018 Asian Research Publishing Network (ARPN). All rights reserved.

www.arpnjournals.com(a)Volumetric water content(m3/m3)

0.20
0.00

0.30

0.40

0.50

(b) Std. Dev., or (c) Skewness, or (d)MeanMedian(e) Frequency

0.0 0.1 0.2 0.3

0.01

0

Max-Min

Kurtosis-2 -1 0 1 2 3-0.01

0.005 m0.045 m0.05
Depth (m)1500 3000

0.085 m0.10
0.125 m0.15
0.165 m0.20
0.205 m

0.25
0.30

Mean
Min, Max95 percentile

0.245 mStd. Dev.
Max-Min

Kurtosis
Skewness

MeanMedian

0.285 mFigure-6. Go approximately here.
A high standard deviation and Max-Min values ofCase C than Case B were not unusual given that 20 %increment of input parameter values of the former case,than the latter case of only 10 % increment. The
corresponding largest values of Max-Min were found to be
0.245 and 0.202 m3/m3 that were equivalent to 51.6 and45.7 % of error from the base case water content. The
significant digit approximation method in Case A resulted
in 0.06 m3/m3 that can be translated into 15.4 % of error
from base case water content.
3.2 Simulation outcomes of simple genetic algorithm(SGA)
In SGA, there are three basic input data required,
which were population size, mutation number, andgeneration number. In the first stage of the studies, thegeneration number was set at 100 while mutation number
was varied between 0 and 90, and population size wasstudied between 20 and 300. The input parameters of CaseC were used at this stage because it has the largest
uncertainty values of input parameters compared to theother two cases. The simulation results were compared toPhilip’s semi-analytical solution from Haverkamp et al.

(1977), and from now onis known as Philip (H). The
objective function of the model was to minimize the sumof the absolute difference between simulated and Philip(H), from now onis known as OF. Each chromosome (orindividual) in the population would be assigneda particular
OF value based on its simulation output. The primaryobjective of the SGA is to minimize the OF value. The
results showed that the OF value reduced to a minimum
value after a number of generations (Figures 7(a)-(b)), i.e.
approximately 10 generations. The OF value reached theminimum value at later generation as the population size
increases from 50 to 300, as shown in Figures 7(b)-(d). In
addition, it was found that neither extremely high nor low
mutation number was required to obtain the lowest OF
value at each population size. In the case of 20 population
size, mutation number was found to be an optimum at 18.
The mutation sites needed to obtain small OF values
reduced as the population size increased from 50 to 300.
At population size of 50, a mutation number of 18 (Figure7(b)) was found sufficient to result in the smallest OF
value, whereas at 200 and 300 population sizes the thelowest OF values were achieved respectively at 9 (Figures
7(c)) and 4 (Figure-7(d)) mutation numbers.

3902

VOL. 13, NO. 12, JUNE 2018ISSN 1819-6608

ARPN Journal of Engineering and Applied Sciences©2006-2018 Asian Research Publishing Network (ARPN). All rights reserved.

www.arpnjournals.com

Figure-7. Goes approximately here.
Volumetric water content (m3/m3)
0.10
0.00
0.05

Depth (m)Moreover, the global search for lower OF values
indicates that at higher population size the minimum OF
value would reduce from 0.1395 to 0.1123 at eachconsecutive population size from 20 to 300, respectively.
The water infiltration front data points from Philip(H) was
plotted in Figure-8. They were inverse searched with input
parameter uncertainty range of Case C, which was statedearlier to have the most wide uncertainty range that its
success in inverse search should presumably be applicable
to Cases A and B. The Philip(H) was inverse searchedwith different population sizes from 20 to 300. The resultsof the inverse search by population sizes of 50 and 300
were showed in Figure-8. A similar trend of water content
was observed from the two different population sizes fromthe surface of infiltration model until the depth of 0.25 m.
However, a significant discrepancy was observed between
the predictions of both population sizes at water
infiltration depth after 0.25 m. Although both population
sizes of 50 and 300 have showed approximately equal
prediction of water content at the near surface, their values
were also significantly deviated from the Philip (H), after
0.25 m infiltration depth.

0.10

0.30

0.50

Philip(H) semianalytical solution
p=50,m=18
p=300,m=4

0.15
0.20
0.25
0.30

Figure-8. Go approximately here.

3903

VOL. 13, NO. 12, JUNE 2018ISSN 1819-6608

ARPN Journal of Engineering and Applied Sciences©2006-2018 Asian Research Publishing Network (ARPN). All rights reserved.

www.arpnjournals.com
Volumetric water content (m3/m3)
0.20
0.00

0.30

0.40

0.50

200 data points

Depth (m)0.05

p=20,m=18
p=50,m=18

0.10
0.15
0.20
0.25
0.30

Figure-9. Go approximately here.
Apparently, the result up to the current stage, as
in Figure-8 indicates that SGA was not able to reproduce agood result at two end tails of water infiltration front. A
hypothesis was established that the insufficient data at
both ends were the result of unabling the SGA to searchproperly and subsequently result in failure to represent thecharacteristic of water infiltration adequately. Hence, theVolumetric water content(m3/m3)

(a)

0.20
0.00

0.40

(b)

0.60

hypothesis was tested by increasing the data points from18 to 200 and re-run the SGA. The result in Figure-9
shows that the population size of 50 was already sufficient
to obtain a good agreement between numerical data andsimulation. Similarly, at higher population sizes, e.g. 100500, were able to reproduce a good fitting result. It
suggests that sufficient water content data at the upper andlower plains of water infiltration front were necessary, andthey dictate the accuracy of SGA prediction. Without suchinformation, the SGA would not be able to predict aleveling off of the saturation water content and the residual
water content.
In addition, SGA was subjected to artificial noise
to test its ability to operate on non-smooth water content
data. Three different Gaussian noises were imposedin eachcase of A, B, and C, which is equivalent to a total of nine
independent investigations. An example of Gaussian noise
application on data has been demonstrated inBalasubramanian and Schwartz (2002). SGA was foundable to generate smooth curve for each type of noise foreach case. Figure-10 shows the selected few cases of thewater infiltration front where a reasonably good fitting
result was achieved for each case. It was noticeable thatnoise data do not impede the function of SGA as aninverse method, but the noise data were rather influence
the curvature shape of the water infiltration front.

Volumetric water content(m3/m3)

0.20

0.40

0.60

(c)
0.20

Volumetric water content(m3/m3)0.40

0.60

Depth (m)

0.05
0.10
0.15
0.20
0.25
0.30

Simulated
Noise-Type1

Simulated
Noise-Type2

Simulated
Noise-Type3

Figure-10. Goes approximately here.
Additional tests on the robustness of the SGA
were also carried out. Apart from the input parameter
range in Cases A, B, and C, random variation of input
parameter was tested. Again, the population size of 50 was
found sufficient to reproduce the water infiltration front.
At extreme case when all the lower boundaries of the input
parameters uncertainty values were to reduce to zero,
except the saturated water content parameter, the SGA wasfound able to reproduce a reasonably good result at theexpense of higher population size, i.e. 300. When theupper limit of each input parameter was increased to few
hundreds of percentage, on top of the previous casecondition, a population size of 500 was needed toreproduce a reasonable good result.
4. CONCLUSIONS
Saltelli and Annoni (2010) suggested that allmodels should be tested by global sensitivity analysis inorder to avoid perfunctory sensitivity analysis, i.e. local
sensitivity analysis, due to an incomplete investigation of3904

VOL. 13, NO. 12, JUNE 2018ISSN 1819-6608

ARPN Journal of Engineering and Applied Sciences©2006-2018 Asian Research Publishing Network (ARPN). All rights reserved.

www.arpnjournals.com
input spaces. Identification of relevant parameters fromthose of least importance is necessary. The uncertainty
values of input parameters could be subjected to the SGA
to find a set of globally optimized input parameters values
that was best described the water infiltration model. Thosemodels with large number of input parameters may preferto first screening out uninfluential parameters using
elementary effect method (EEM) before processed by
Sobol’ variance-based method and SGA. Hence, thecomputational time would tremendously reduce under
reduced number of input parameters.
In Sobol’ variance-based method, the equal
percentage variation on input parameters was showed tolead to sensitivity analysis (Cases B and C), while
different percentages variation of each input parameter
would result in uncertainty analysis (Case A). It is simply
because those in different case studies withdifferentpercentage variations on each input parameter
would result in different first order and total effect indicesvalues. In addition, the simulated output was significantly
influenced by saturated volumetric water content at nearsurface depth; and at deeper depth after the water
infiltration front, the output was solely influenced byinitial water content used in the simulation. It thus
suggested that the accuracy of water content prediction inthese regions would depend on a narrower range of
uncertainty values of the related parameters. However, in asituation when a particular parameter, e.g. in Case A,
has a greater values of uncertainty than other parameters,
the previous observation does not apply. A simple casestudy of Case A has demonstrated that significant digit
approximation that was developed based on existing
significant digit of published data was found responsiblefor 15.4 % of error. Such error could be reduced byproviding a greater number of significant digits on the
input parameters values.
The SGA was found capable of searching for
globally optimized input parameter values. A population
size of 50 was found successful of searching andreproduced water infiltration front in all cases tested,
provided that sufficient data at the upper and lower plainsof the water infiltration front were provided. SGA wasable to handle a random range of input parameter values
and also noise data. At extreme cases when extremely low
and high values were used as search space for input
parameters, the population size needed to perform theinverse method of SGA would be increased.

Campolongo F., Cariboni J., Saltelli A. 2007. An effective
screening design for sensitivity analysis of large models.
Environmental Modelling & Softw. 22(10): 1509-1518.
doi:http://dx.doi.org/10.1016/j.envsoft.2006.10.004.
Chu-Agor M.L., Muñoz-Carpena R., Kiker G.,
Emanuelsson A., Linkov I. 2011. Exploring vulnerabilityof coastal habitats to sea level rise through global
sensitivity and uncertainty analyses. Environmental
Modelling&
Softw.
26(5):
593-604.
doi:http://dx.doi.org/10.1016/j.envsoft.2010.12.003.
Drouet J.L., Capian N., Fiorelli J.L., Blanfort V.,
Capitaine M., Duretz S., Gabrielle B., Martin R., Lardy R.,
Cellier P., Soussana J.F. 2011. Sensitivity analysis formodels of greenhouse gas emissions at farm level. Casestudy of N2O emissions simulated by the CERES-EGC
model. Environmental Pollut. 159(11): 3156-3161.
doi:http://dx.doi.org/10.1016/j.envpol.2011.01.019.
EPA: 2009. Guidance on the Development, Evaluation,
and Application of Environmental Models, Technicalreport, Office of the Science Advisor, Council for
Regulatory Environmental Modeling, March 2009,
EPA/100/K-09/003,
http://www.epa.gov/crem/library/cred_guidance_0309.pdf,
accessed 2 June 2014 In.
Fayer M.J. 2000. UNSAT-H Version 3.0: Unsaturated Soil
Water and Heat Flow Model - Theory, User Manual, andExamples. In., vol. PNNL-13249. Pacific Northwest
Laboratory, Richland, Washington.
Fayer M.J., Jones T.L. 1990. UNSAT-H Version 2.0:
Unsaturated soil water and heat flow model. In., vol. PNL6779. Pacific Northwest Laboratory, Richland,
Washington.
Goh E.G., Noborio K. 2013. Sensitivity Analysis on the
Infiltration of Water into Unsaturated Soil. In: Proceedingsof Soil Moisture Workshop 2013, Hiroshima University
Tokyo Office in Campus Innovation Center. pp. 66-68.
Haverkamp R., Vauclin M., Touma J., Wierenga P.J.,
Vachaud G. 1977. A Comparison of Numerical Simulation
Models for One-Dimensional Infiltration1. Soil Sci. Soc.
Am.
J.
41(2):
285-294.
doi:10.2136/sssaj1977.03615995004100020024x.

REFERENCES
Avagliano S., Parrella L. 2009. Managing Uncertainty in
Risk-Based Corrective Action Design: Global Sensitivity
Analysis of Contaminant Fate and Exposure Models Usedin the Dose Assessment. Environ. Model Assess. 14(1):47-57. doi:10.1007/s10666-008-9163-5.
Balasubramanian M., Schwartz E.L. 2002. The Isomap
Algorithm and Topological Stability. Science 295(5552),
7. doi:10.1126/science.295.5552.7a.

Holland J.H. 1975. Adaptation in Natural and Artificial
Systems. University of Michigan Press.
Ines A.V.M., Droogers P. 2002. Inverse modelling inestimating soil hydraulic functions: a Genetic Algorithm
approach. Hydrol. Earth Syst. Sci. 6(1): 49-66.
doi:10.5194/hess-6-49-2002.
Kelleners T.J., Soppe R.W.O., AyarsJ.E., Šimůnek J.,
Skaggs T.H. 2005. Inverse Analysis of Upward Water

3905

VOL. 13, NO. 12, JUNE 2018ISSN 1819-6608

ARPN Journal of Engineering and Applied Sciences©2006-2018 Asian Research Publishing Network (ARPN). All rights reserved.

www.arpnjournals.com
Flow in a Groundwater Table Lysimeter. Vadose Zone J.
4(3): 558-572. doi:10.2136/vzj2004.0118.
Mollerup M., Hansen S., Petersen C., Kjaersgaard J.H.
2008. A MATLAB program for estimation of unsaturated
hydraulic soil parameters using an infiltrometer technique.
Computers&
Geosciences.
34(8):
861-875.
doi:http://dx.doi.org/10.1016/j.cageo.2007.12.002.
Morris M.D. 1991. Factorial sampling plans forpreliminary computational experiments. Technometrics.
33(2): 161-174. doi:10.2307/1269043.
Noborio K., McInnes K.J., Heilman J.L. 1996. TwoDimensional Model for Water, Heat, and Solute Transport
in Furrow-Irrigated Soil: II. Field Evaluation. Soil Sci.
Soc.
Am.
J.
60(4):
1010-1021.
doi:10.2136/sssaj1996.03615995006000040008x.
Nossent J., Elsen P., Bauwens W. 2011. Sobol’ sensitivity
analysis of a complex environmental model.
Environmental Modelling & Softw. 26(12): 1515-1525.
doi:http://dx.doi.org/10.1016/j.envsoft.2011.08.010.
Pachepsky Y., Timlin D., Rawls W. 2003. Generalized
Richards' equation to simulate water transport inunsaturated soils. J. of Hydrology 272(1-4): 3-13.
doi:http://dx.doi.org/10.1016/S0022-1694(02)00251-2.

layer-specific soil moisture assimilation scheme. Water
ResourcesRes.
48(6):
W06529.
doi:10.1029/2010WR009581.
Šimůnek J., Bradford S.A. 2008. Vadose zone modeling:
Introduction and importance. Vadose Zone J. 7(2): 581586.
Sobol I.M.: 1990. On sensitivity estimates for nonlinear
mathematical models. Matematicheskoe Modelirovanie.
112-118.
Vrugt J.A., Stauffer P.H., Wöhling T., Robinson B.A.,
Vesselinov V.V. 2008. Inverse Modeling of Subsurface
Flow and Transport Properties: A Review with New
Developments. Vadose Zone J. 7(2): 843-864.
doi:10.2136/vzj2007.0078.
Yang J. 2011 Convergence and uncertainty analyses inMonte-Carlo based sensitivity analysis. Environmental
Modelling&
Softw.
26(4):
444-457.
doi:http://dx.doi.org/10.1016/j.envsoft.2010.10.007.
Younes A., Mara T.A., Fajraoui N., Lehmann F., Belfort
B., Beydoun H. 2013 Use of Global Sensitivity Analysis
to Help Assess Unsaturated Soil Hydraulic Parameters.
Vadose Zone J. 12(1): doi:10.2136/vzj2011.0150.

Richards L.A. 1931. Capillary Conduction of Liquids
through Porous Mediums. J. of Appl. Phys. 1(5): 318-333.
doi:http://dx.doi.org/10.1063/1.1745010.
Saltelli A., Annoni P. 2010. How to avoid a perfunctory
sensitivity analysis. Environmental Modelling & Softw.
25(12):
1508-1517.
doi:http://dx.doi.org/10.1016/j.envsoft.2010.04.012.
Saltelli A., Annoni P., Azzini I., Campolongo F., Ratto M.,
Tarantola S. 2010. Variance based sensitivity analysis ofmodel output. Design and estimator for the total sensitivity
index. Computer Phys. Communications. 181(2): 259-270.
doi:http://dx.doi.org/10.1016/j.cpc.2009.09.018.
Saltelli A., Ratto M., Andres T., Campolongo F., Cariboni
J., Gatelli D., Saisana M., Tarantola S. 2008. VarianceBased Methods. In: Global Sensitivity Analysis. The
Primer. pp. 155-182. John Wiley & Sons, Ltd.
Saltelli A., Tarantola S., Campolongo F., Ratto M. 2004.
The Screening Exercise. In: Sensitivity Analysis in
Practice. pp. 91-108. John Wiley & Sons, Ltd.
Schneider S., Jacques D., Mallants D. 2013. Inverse
modelling with a genetic algorithm to derive hydraulicproperties of a multi-layered forest soil. Soil Res. 51(5):372-389. doi:http://dx.doi.org/10.1071/SR13144.
Shin Y., Mohanty B.P., Ines A.V.M. 2012. Soil hydraulicproperties in one-dimensional layered soil profile using3906