１

愛媛大学農学部紀要（Mem. Fac. Agr., Ehime Univ.）
５
７：１−３
６（２
０
１
２）

Inverse Modeling of Variably Saturated Subsurface
Water Flow in Isothermal/Non-isothermal Soil＊

Tomoki IZUMICONTENTS
Notations
List of Tables
List of Figures

4
2
4
4

1 Introduction
5
1.1 Background ………………………………………5
1.1.1 Worldwide water crisis………………………5
1.1.2 Status of water use and multifunctionalityof agriculture…………………………………5
1.2 Research Aim and Objectives …………………5
1.3 Structure of Thesis Framework …………………6
2 Literature Review
7
2.1 Governing Equations for Subsurface Water Flow………………………………………………7
2.1.1 Description of soil water movement
and numerical method ………………………7
2.1.2 Water movement in non-isothermal soil ……7
2.2 Soil Hydraulic Properties and Identification
Method ……………………………………………82.2.1 Representation of soil hydraulic properties………………………………………………82.2.2 Identification method ………………………8
3 Inverse Method to Identify Soil Hydraulic
Properties in Variably Saturated Subsurface
Water Flow Model ― "-based Form of Richards
Equation―9
3.1 Introduction ………………………………………93.2 Forward Problem …………………………………9
3.3 Inverse Problem …………………………………10
3.3.1 Parameterization ……………………………10
3.3.2 Solution procedure of IP …………………10
3.3.3 Optimization algorithm ……………………11
3.4 Validation ………………………………………11
3.4.1 Methods ……………………………………11
3.4.2 Results ………………………………………12
3.5 Conclusions ……………………………………14

Received March 15, 2012.
Accepted April 25, 2012.
＊
Lab. Water Resources Engineering

Inverse Method to Identify Unsaturated
Hydraulic Conductivity in Variably Saturated
Subsurface Water Flow Model in Non-isothermal
Soil ― !-based Form of Richards Equation―
15
4.1 Introduction ……………………………………15
4.2 Method …………………………………………15
4.2.1 Observation system…………………………15
4.2.2 Forward problem……………………………16
4.2.3 Parameterization of soil hydraulic properties………………………………………………16
4.2.4 Inverse problem and solution procedure …17
4.3 Validation ………………………………………17
4.3.1 Field observation and soil water retention
curve fitting…………………………………18
4.3.2 Results and discussion ……………………19
4.4 Conclusions ……………………………………21
5 Inverse Method to Identify Unsaturated
Hydraulic Conductivity in Variably Saturated
SubsurfaceWater Flow Model ―Mixed Form ofRichards Equation―
22
5.1 Introduction ……………………………………22
5.2 Forward Problem ………………………………22
5.2.1 Governing equation…………………………22
5.2.2 Parameterization of soil hydraulic properties………………………………………………22
5.2.3 Time-marching scheme ……………………23
5.3 Inverse Problem …………………………………23
5.3.1 Formulation of optimization problem ……23
5.3.2 Optimization algorithm ……………………24
5.4 Validation ………………………………………24
5.4.1 Observation system…………………………24
5.4.2 Identification of water retention property …25
5.4.3 Identification of RHC and validation of themethod………………………………………25
5.5 Conclusions ……………………………………26
6 Inverse Method to Identify Unsaturated
Hydraulic Conductivity in Variably Saturated
Subsurface Water Flow Model in Non-isothermal
Soil ―Mixed Form of Richards Equation―27
6.1 Introduction ……………………………………27
6.2 Governing Equations System …………………27
6.2.1 Water movement model ……………………27
6.2.2 Thermal transport model …………………27
6.2.3 Parameterization of soil hydraulic properties………………………………………………27
6.2.4 Numerical solution procedure ……………28
6.3 Parameter Identification Procedure ……………29

２

Tomoki IZUMI

6.3.1 Inverse problem ……………………………29
6.3.2 Optimization algorithm ……………………29
6.4 Validation ………………………………………29
6.4.1 Field observation and computational domain………………………………………………29
6.4.2 Identification of SWRC ……………………30
6.4.3 Identification of RHC ………………………31
6.4.4 Reproducibility of calibrated forward
simulation model for variably saturated
subsurface flow ……………………………31
6.5 Conclusions ……………………………………32
7 Summary and Conclusions
33
Acknowledgments
34
References
34

NOTATIONS
The following are the major symbols used in this thesis.
',
',',"
(
,
(
,
"
(
,
!,
!5
)
,
)
,
"
)
,
0
)
)
5
*,
*,*,"
"

= coefficient of piecewise cubic spline= coefficient of piecewise cubic spline= coefficient of piecewise cubic spline= coefficient of piecewise cubic spline= coefficient of piecewise cubic spline= coefficient of piecewise cubic spline= volumetric heat capacity= specific moisture capacity= coefficient of piecewise cubic spline= coefficient of piecewise cubic spline= coefficient of piecewise cubic spline= volumetric heat capacity of soil particle= volumetric heat capacity of water= coefficient of piecewise cubic spline= coefficient of piecewise cubic spline= coefficient of piecewise cubic spline= absolute error between observed and computed value")/.
= absolute error between observed and computed value with non-isothermal/free-form
parameterization (FFP) model"-2/= absolute error between observed and computed value with isothermal/FFP model"4+= absolute error between observed and computed value with non-isothermal/van Genuchten (VG) model#= gravitational acceleration!
= Hessian matrix+= hydraulic head"
= I × I unit matrix#= total least squared error$= unsaturated hydraulic conductivity$(*= admissible set of unknown parameter to be
identified$1= relative hydraulic conductivity" = relative hydraulic conductivity function over!
$1!,!
the ,th subdomain" = relative hydraulic conductivity function repre!
$1##%!
sented by FFP model" = relative hydraulic conductivity function repre!
$1'$!
sented by van Genuchten-Mualem model$2= saturated hydraulic conductivity$&= temperature factor of hydraulic conductivity#= vector of unknown parameter to be identified#/03= vector of optimal parameter-,= value of $1##% at node ,%= total number of available observed data.
= number of iteration step.4+
= parameter of VG model/= number of time level/4+= parameter of VG model&(*= admissible set of $$= vector of unknown parameter to be identified$/03= vector of optimal parameterInverse Modeling of Variably Saturated Subsurface Water Flow in Isothermal/Non-isothermal Soil

&
&/
&/
&6
&6
!
""4
"6
)
##
#4
#4,21
#42+4
#4
#4
'
(%
)

= flux on Neumann boundary= heat flux
= heat flux on Neumann boundary= water flux
= water flux on Neumann boundary= rainfall= effective saturation= specific storage= saturation= time domain= initial value of soil temperature= soil temperature= computed value of soil temperature= observed value of soil temperature= soil temperature= value of #4 on Dirichlet boundary= time= weighting factor= independent variable for a vertical axis
Greek letters!5.
= parameter of VG model"4= compressibility of soil
"6
= compressibility of water= Dirichlet boundaryΓ$= Dirichlet boundary for heat transportΓ$/= Dirichlet boundary for water movementΓ$6= Neumann boundaryΓ'= Neumann boundary for heat transportΓ'/= Neumann boundary for water movementΓ'6%= number of iteration#'= time step#= convergence criterion&= control variable in terms of search direction*= volumetric water content*= value of *on Dirichlet boundary!
= vector of unknown parameter to be identified,"= soil water retention function represented with*%%(!
FFP model," = soil water retention function represented with**&!
VG model*,21= computed value of *with non-isothermal/FFP
model
*042
= computed value of *with isothermal/FFP
model
*2+4
= observed value of **5.
= computed value of *with non-isothermal/VG
model*#= initial condition of **$,"at node $= value of *%%(!
," = soil water retention function over the $th*$!
subdomain*3= residual water content*4= saturated water content'= intrinsic permeability(= thermal conductivity(#= reference thermal conductivity)= dynamic viscosity

+
+3
+6
$
,
,#
,
,,21
,2+4
,5.
Ω

３

= water density= water density at reference temperature= water density= porosity= capillary head (on , ! 0) or pressure head(on ," 0)= initial value of ,= boundary value of ,on Dirichlet boundary= computed value of ,with FFP model= observed value of ,= computed value of ,with VG model= space domain

４

Tomoki IZUMILIST OF TABLES

LIST OF FIGURES

Common parameters ………………………………12
Parameters of Soil A ………………………………12
Parameters of Soil B ………………………………12
Identified parameter values ………………………19!"in simulationMinimum values of !!
optimizarion runs for different periods ……………21!"in simulation5.1 Minimum values of !!
optimization ………………………………………25!" …………………………32
6.1 Optimal values of !!

1.1
3.1
3.2
3.3
3.4
3.5
3.6
3.7

3.1
3.2
3.3
4.1
4.2

Structure of the main chapters ……………………6
Cubic splines ………………………………………10
Solution search algorithm …………………………11
Soil column ………………………………………12"'##% compared with "'&$ (Soil A) ………………13!##% compared with !&$ (Soil A) …………………13
Time-varying pressure heads (Soil A) ……………13
Absolute errors between computed and truepressure heads (Soil A) ……………………………13
3.8 "'##% compared with "'&$ (Soil B) ………………13
3.9 !##% compared with !&$ (Soil B) …………………13
3.10 Time-varying pressure heads (Soil B) ……………14
3.11 Absolute errors between computed and truepressure heads (Soil B) ……………………………14
4.1 Basic and optimal instruments for field observations……………………………………………………15
4.2 Free-form parameterization ………………………17
4.3 Full instrumentation ………………………………18
4.4 Observed data ……………………………………18
4.5 Correlation between suctions in daytime andnighttime …………………………………………18
4.6 Soil water retention property in the first desorption
period ………………………………………………19
4.7 Soil water retention property in the second
desorption period …………………………………19
4.8 Soil water retention property in the third desorption
period ………………………………………………19
4.9 Calibration result for the first desorption period …20
4.10 Calibration result for the second desorption period……………………………………………………20
4.11 Calibration result for the third desorption period …21
5.1 Free-form parameterization ………………………23
5.2 Observation system and discretization of targeteddomain ……………………………………………24
5.3 Observed data for validation ………………………24
5.4 Correlation between pressure heads in daytime andnighttime for the same volumetric water content …25
5.5 Soil water retention property………………………25
5.6 (a) Identified RHC and (b) Reproducibility of FP
solution for〈Data A〉 ……………………………26
5.7 (a) Identified RHC and (b) Reproducibility of FP
solution for〈Data B〉 ……………………………26
6.1 Free-form parameterization ………………………28
6.2 Flow chart for solving IP …………………………29
6.3 Schematic illustration of observation andcomputation ………………………………………30
6.4 Observed data ……………………………………30
6.5 Identified SWRC …………………………………30
6.6 Identified RHC ……………………………………31
6.7 Reproducibility of forward solution for water
movement …………………………………………31
6.8 Reproducibility of forward solution for heat transport……………………………………………………31
6.9 Identified RHC by isothermal model ……………32
6.10 Reproducibility of forward solution for water
movement in isothermal model ……………………32

Inverse Modeling of Variably Saturated Subsurface Water Flow in Isothermal/Non-isothermal Soil

CHAPTER 1
INTRODUCTION
1.1 Background1.1.1 Worldwide water crisis
The world is facing water crisis.
Water demand has been growing steadily in line withthe population increase and economic growth withurbanization and industrialization. According to the United
Nations Population Projections (medium variant of the2010 version)[56], the world population stood at some 6.9billion in 2010 and is expected to exceed 9.0 billion in
2050 and 10.0 billion by the end of this century. Though
Cohen (1995)[7] estimated the maximum allowable world
population as 7.0 billion based on the amount of
renewable freshwater resources, in fact, water resources
are running short because of the uneven distribution ofrenewable freshwater resources!in time and space. Using“Falkenmark indicator” or “water stress index” (Falkenmark, 1989[12]), people in the water stress state!when the
per-capita maximum available amount of renewablefreshwater resources falls below 1,700 cubic kilometers!
are numbered approximately 0.7 billion in 2005 and 2.0
billion in 2008 (UNDP, 2006[55]; Ministry of the Environment, 2011[40]).
Furthermore, the climate change adds momentum to thewater crisis. According to the Fourth Assessment Report(AR4) of the Intergovernmental Panel on Climate Change(IPCC)[24], the progress in global warming will expose
several hundred million people to increased water stress inthe coming years, and the increased frequency of droughts
and flooding is projected to adversely impact foodproduction. The rise in the global average temperature due
to climate change is feared to bring about various impacts
on water resources. The rapid urbanization is also
noticeable and not only demand for domestic andindustrial water is increased but also massive amounts ofwastes and effluents are generated in association withmass consumption which aggravates the environmentalload and deterioration of water quality (Ministry of theEnvironment, 2011[40]).
This crisis on both quality and quantity of water
resources indicates that human activities have reachedbeyond the environmental carrying capacity. Thus, we,human beings are required to wisely and efficiently
manage water resources based on the ideals of“sustainable development” propounded at the United
Nations World Commission on Environment andDevelopment in 1987.
1.1.2 Status of water use and multifunctionality ofagriculture
Driven by solar energy and gravitational force, water
globally circulates within the hydrosphere, the atmosphere,
surface water, subsurface water and plants. From theglobal hydrologic cycle, human beings withdraw freshwater like surface and subsurface water for agriculture,５

industry and household with the ratio of 70%, 20% and10%, respectively (FAO[13]). This tendency is closely thesame in Japan (Ministry of Land, Infrastructure, Transport
and Tourism, 2011[39]). Accordingly, it is impossible tosustainably plan or manage water use without considering
agriculture, which is the largest withdrawing sector.
Agriculture plays the role of not only its primaryfunction of food production but also the conservation ofnational lands, ecosystems and the environment throughcreating and maintaining a secondary nature. However,agriculture adversely impacts ecosystems, the global
environment and even agricultural environment itself. The
productivity-oriented agricide causes the reduction of soil
and land fertility and leads to desertification. Applicationof large amount of pesticide and fertilizer diminishes
water quality. Moreover, large scale irrigation in arid
areas increases the problem of salt accumulation. It is,therefore, of great importance to quantitatively assess theeffect of agricultural activities on the environment and tokeep it below the threshold limit of sustainable range.
1.2 Research Aim and Objectives
This thesis aims at developing a methodology toidentify the soil hydraulic properties with inverse
technique to contribute to sophisticated numericalmodeling of subsurface water flow phenomena that arebasic and representative in irrigation agriculture. It is alsoquite useful to accurately estimate or predict subsurface
water flow in substance transport problems dealing withthe assessment of the influence of pesticide or fertilizer
application and salt accumulation on water quality.
Generally, intensive and exhaustive investigations andanalyses are required for the estimation of flow fields
though it depends on the space scale of target area. To
reduce labor and cost, numerical analyses haveincreasingly attracted lots of attention because of the great
advances of computer technologies and capacities inrecent years. The numerical analyses are required toreproduce a physical phenomenon of interest moreaccurately, and their success depends on the model
structure given by the governing equations system andparameter identification which is a critical step inmodeling process.
A physically-based mathematical model governing
water flow in soil is represented by Richards equation(Richards, 1931[47]). Richards equation can be classified
into three different forms based on the decision variables:a pressure-head-based ("-based) form, a moisture-contentbased (!-based) form, and a mixed form. According toRichards equation, the water movement through soil isderived from the gradient of hydraulic head underisothermal assumption, while water movement is also
affected by soil temperature. It is, thus, necessary to deal
with a coupling problem of water movement and thermaltransport to represent of water movement near surface soil
where heat energy exchange is quite significant.
On the other hand, it is also very important to identify
model parameters included in the governing equations.

６

Tomoki IZUMIModel parameters in Richards equation are referred to as
the soil hydraulic properties, i.e. the soil water retention
curve and the unsaturated hydraulic conductivity.
However, an identification method has not yet been
established because these parameters are not directly
measurable and have to be determined from historicalobservations.
Based on the above premises, the following proposals
on inverse modeling of subsurface flow in isothermal andnon-isothermal soil are the main objectives of this thesis:(1) Inverse modeling for "-based form of Richards
equation(2) Inverse modeling for !-based form of Richards
equation(3) Inverse modeling for mixed form of Richards
equation
1.3 Structure of Thesis Framework
This thesis is organized into seven chapters including
three original research papers of the author[26]−[28].
The introductory chapter gives an overview of thecurrent status of water environment and discusses theimportance of the water resources management inagriculture, which is the largest water withdrawing sector.
The research aim and objectives are also presented,
emphasizing the necessity of sophisticated subsurface
water flow models to quantitatively assess the impacts ofagriculture on the environment.
In Chapter 2, recent and pertinent studies onmathematical and numerical models for subsurface waterflow in both isothermal and non-isothermal soil andrepresentation and identification of model parameters!
soil hydraulic properties!are reviewed.
In Chapter 3, an inverse modeling for "-based form ofRichards equation in variably saturated and isothermal soil
is proposed and discussed. To accurately describe
subsurface water flow, the soil hydraulic properties, whichis the unknown parameters to be identified, arerepresented by a free form parameterized function. As theparameterized function, a sequential piecewise cubic
spline function is used in lieu of the conventional fixed
form parameterized function. After an inverse problem issolved by simulation-optimization algorithm which
combines a numerical model to solve a forward problem
with a optimization method with the aid of theLevenberg-Marquardt method, the free form parameterizedsoil hydraulic properties are determined. The method
developed is validated through twin experiments.
In Chapter 4, an inverse modeling for !-based form ofRichards equation in unsaturated and non-isothermal soil
is proposed and discussed. Since water movement insurface soil is significantly affected by soil temperature,
the forward problem on subsurface water flow is based oncoupled water movement and heat transport equations.
Validation of the method proposed is examined throughapplying to upland soil based on field observations using
a proposed utilitarian observation system with simple
instrumentation.
In Chapter 5, an inverse modeling for mixed form ofRichards equation in variably saturated and isothermal soil
is proposed and discussed. The success of numericalmodeling depends on the model structure given by
governing equations. Since the mixed form Richards
equation is superior than the two other forms of Richards
equation in terms of the range of application to saturationand conservation of mass balance, the forward problem in
the inverse method is described as the mixed form
Richards equation. The method is validated based on thesame field observation in Chapter 4.
In Chapter 6, an inverse modeling for mixed form ofRichards equation in variably saturated and non-isothermal
soil is proposed and discussed. It is an extension of themethod proposed in Chapter 5 to non-isothermal
condition. The governing equation system is described asa coupled problem of mixed form Richards equation andheat conduction equation. From field observation,
validation of the inverse method is carried out.
And, Chapter 7 presents the summary and conclusionsof this thesis.
The matrix to illustrate the relation of the main chapters(Chapters 3 through 6) is shown in Figure 1.1.

Isothermal Model Non-isothermal Model"-based Form RE

Chapter 3!-based Form RE
Mixed Form RE

Chapter 4
Chapter 5

Chapter 6

RE: Richards Equation
Figure 1.1: Structure of the main chapters

Inverse Modeling of Variably Saturated Subsurface Water Flow in Isothermal/Non-isothermal Soil

CHAPTER 2
LITERATURE REVIEW
2.1 Governing Equations for Subsurface Water Flow
2.1.1 Description of soil water movement andnumerical method
Water movement in saturated porous media is describedas the Darcy’s law. The Darcy’s law was improved by
Buckingham in 1907 to present water flux in unsaturated
porous media (Jury and Horton, 2004[30]). Additionally, theDarcy-Buckingham law has been extended by Richards(1931)[47] to transient flow, combined with the water
conservation equation. Nowadays, the Richards equation
is mostly used as the governing equation for subsurface
water flow analysis.
Because of the highly nonlinear model parameters in
Richards equation, analytical solutions are not generally
possible except under very restricting assumptions
regarding the parameters (Hillel, 1998[16]). Instead, various
numerical methods such as the finite element method(FEM) have been developed and become a powerful
technique with the advancement of computer technologies.
Among those numerical methods, FEM is most frequently
used because it enables us to analyze flow regions having
complex geometric boundaries and arbitrary degrees ofheterogeneity and anisotropy. FEM, which was firstdevised as a procedure for structural analysis, has been
applied to steady and saturated seepage problems by
Zienkiewicz (1966)[64]. Neuman and Witherspoon (1971)[45]
have developed a finite element method to analyze
unsteady flow with a free surface, extending the techniqueto solve steady state problems proposed by Neuman andWitherspoon (1970)[44]. Moreover, FEM has been applied
to saturated-unsaturated flow analyses in porous media(Neuman, 1973[43]; Akai et al ., 1977[1]).
From the view point of decision variables, Richards
equation can be classified into following three forms:(1) the "-based form,(2) the !-based form and(3) the mixed form.
The "-based form of Richards equation has been mostcommonly used. This is because it can express all
possible range of saturation where " is continuous.
However, the disadvantage of "-based form is that it
requires very fine spatial and temporal discretization to
avoid computational instability and disproportionate massbalance since " is quite sensitive to change of ! in very
dry conditions (Hills et al ., 1989[17]). Various approaches
have been proposed to overcome these disadvantages.
Neuman (1973)[43] has proposed a mass lumping approachto improve numerical convergence and to eliminate
numerical oscillations. Milly (1985)[38] has presented massconserving schemes that used a modified definition of the７

storage term to satisfy a global mass balance.
The !-based form of Richards equation is derived by
Klute (1952)[32]. The advantages of the form are that (1)water mass is automatically conserved within thecomputational domain regardless of time and spatial step
size and (2) the hydraulic diffusivity!model parameter
included in the !-based form equation!does not varywith ! nearly as much as the unsaturated hydraulicconductivity varies with ". The disadvantage is that it
cannot be used for simulating flow in soil near or atsaturation since in that range the hydraulic diffusivity
becomes infinite (Hillel, 1998[16]; Hills et al ., 1989[17]).
The mixed form of Richards equation is an equation
whose decision variables are both " and !. Celia et al .
(1990)[6] have presented a numerical model which pertains
to the mixed form, perfectly conserving the mass using amodified Picard iteration method. After the work of Celia
et al . (1990)[6], the mixed form has been used often as thegoverning equation for subsurface water flow, and thecomputational efficiency of its numerical model has
improved (Huang et al ., 1996[19]). Since almost all thenumerical models proposed by Celia et al . (1990)[6] andsubsequently others deal with only unsaturated zone, someextended models for saturated model as well as
unsaturated model have been proposed (Takeuchi et al .,2008[54]; Zadeh, 2011[62]).
2.1.2 Water movement in non-isothermal soil
Richards equation expresses water movement throughsoil derived from the gradient of hydraulic head underisothermal condition while it is easily recognized thatwater movement is also forced by temperature gradient.
When water movement near surface soil where heat
energy exchange is quite significant is considered, coupledequation of water movement and heat transport is needed.
A theory on the simultaneous transfer of water and heathas formulated by Philip and de Vries (1957)[46] and deVries (1958)[9] (PDV model). There have been various
models based on PDV model. Kondo and Saigusa(1994)[33] had constructed a multi-layer soil model toestimate the evaporation from bare soil, incorporating aformula for vaporization in the soil pores in PDV model.
On the other hand, Milly (1982)[36] has converted theequations of PDV model to the "-based formulation inorder to extend the applicable range since PDV model isassociated with the !-based formulation. Milly model has
been widely applied to the studies on evaporation fromsoil surface (Milly, 1984[37]), simultaneous transfer ofwater and heat in desert soil (Scanlon and Milly, 1994[48]),and so on. Moreover, Moukalled and Saleh (2006)[41] had
proposed a numerical procedure to solve Milly model by
the finite volume method being inherently conservative
for the purpose of preventing the mass imbalance problem
reported in the literature when employing the "-based
formulation with other numerical methods.
Fujinawa (1995)[14] has presented a different formulation
coupling water movement and heat transport. An equation
derived from the generalized Darcy’s law and the８Tomoki IZUMIBoussinesq assumption and an advection-dispersion
equation are employed as governing equations for water
movement and heat transport, respectively.
2.2 Soil Hydraulic Properties and Identification Method
To solve the basic equations governing subsurface
water flow above mentioned, the soil hydraulicproperties!the soil water retention curve (SWRC) andunsaturated hydraulic conductivity (UHC)!are needed.
SWRC is defined as the relation between the volumetricwater content and the matric suction, and UHC as afunction of the saturation or suction.
2.2.1 Representation of soil hydraulic properties
The determination of UHC through laboratory or fieldmeasurement is generally laborious and time-consuming
and, thus, it has been attempted to estimate UHC from
SWRC which can be measured with relative ease. The
representative models are Brooks and Corey (BC) model(Brooks and Corey, 1966[5]), van Genuchten-Mualem (VG)
model (van Genuchten, 1980[59]; Mualem, 1976[42]) and
Kosugi’s lognormal pore-size distribution (LN) model(Kosugi, 1994[34]; 1996[35]). Among them, VG model is themost widely used for numerical analysis.
However, some drawbacks have been reported.
According to Iden and Durner (2007)[21], althoughparameterizations of UHC coupled to SWRC by either
empirical equations or capillary models are physicallymotivated, such approaches may give rise to practicalproblems, since (1) the capillary model itself may be only
approximately valid for the porous medium under study
and (2) any error of the employed model of SWRC
propagates into errors in the predicted UHC. Additionally,
Durner (1994)[10] has reported that the conventional
unimodal retention function cannot describe SWRC in soil
with heterogeneous pore systems, and proposed a flexible
SWRC function formed by superimposing unimodal
retention curve of VG model type. Moreover, Vogel et al .
(2001)[61] have pointed out that the UHC derived from VG
model can be quite sensitive to the shape of SWRC nearsaturation with unfavorable consequences for thecomputational performance in the numerical simulations,and presented a modified VG model. VG model could
also lead to wrong predictions of UHC under thecondition where air-entry effects are significant (Ippisch etal ., 2006[25]).
Another alternative approach is a free-form parameterization approach (e.g. Kastanek and Nielsen, 2001[31];Bitterlich et al ., 2004[3]). The approach describes the soil
hydraulic property with a sequential piecewise polynomial
function. The function is referred to as “the free form
function” in contrast with the function of conventionalmodel like VG model being “the fixed form function”.
The advantages of the approach is that it offers highflexibility because of no a priori assumption about the
specific shape of function and the number of parameters(Bitterlich et al ., 2004[3]).

2.2.2 Identification method
The soil hydraulic properties are determined by direct
measurement or indirect method. The direct measurement
methods are summarized by Japanese Association of
Groundwater Hydrology (2010)[29]. Though the direct
measurement of SWRC is possible with wide range ofpressure head, the direct measurement of UHC is difficult
because it is laborious and time-consuming (Hillel,
1998[16]).
Accordingly, the indirect method!inverse method!has
been focused on because of the advantage of optimization
techniques and computer power. The inverse method isbased on two methodologies: (a) the multistep outflow
method (Eching and Hopmans, 1993[11]; van Dam et al .,1994[58]) and (b) the evaporation method (Simunek et al .,1998[57]).
Bitterlich et al . (2004)[3] have proposed an inverse
method to estimate the soil hydraulic properties using freeform functions through column outflow experiments. Iden
and Durner (2007, 2008)[21],[22] have estimated the soil
hydraulic properties in the similar framework of inverse
modeling to the work of Bitterlich (2004)[3], based on themultistep outflow method and evaporation method,respectively.
Zhang et al . (2010)[63] have presented an inverse method
based on field experiments. In their method, FAO
approach (Allen et al ., 1998)[2] is embedded to calculate
the potential evapotranspiration in the forward solution
procedure, since measurement of actual evapotranspirationis generally difficult in field experiments.
"!

Inverse Modeling of Variably Saturated Subsurface Water Flow in Isothermal/Non-isothermal Soil

CHAPTER 3
INVERSE METHOD TO IDENTIFY SOIL
HYDRAULIC PROPERTIES IN VARIABLY
SATURATED SUBSURFACE WATER FLOW MODEL― !-BASED FORM OF RICHARDS EQUATION―
3.1 Introduction
The methodical way to numerically model andreproduce a real groundwater system is at the outset tosolve an inverse problem (IP) to estimate a priori
unknown parameters included in the groundwaterequations system, and subsequently to solve a forward
problem (FP) for obtaining solutions to the equations withthe parameters so estimated. A variety of well-developed
methods to solve FPs are currently available (e.g.,
Huyakorn and Pinder, 1983[20]). Since reliability of thesolutions to FP is contingent on the model parameters tobe estimated prior to these solutions, development of theperformance solver for IP is also an essential need.
Here, inverse modeling for Richards equation which isa standard mathematical model for saturated-unsaturated
groundwater flow in a porous medium is considered.
Conventionally, soil hydraulic properties (i.e., the relative
hydraulic conductivity "$ and the volumetric watercontent #), which are major parameters in the equation,
are described by the fixed-form functions expressed interms of the unknown variable (i.e., the pressure head $),which cover a possible global range of the variable. Mostof the earlier works on the solution of IP are based on theflow model with this type of parameterization in a fixed
functional form. For instance, Hsu and Liu (1990)[18], and"Constales and Kacur(2001)[8] proposed the methods toidentify a set of the functions which consists of theclosed-form equations by van Genuchten (1980)[59], called
VG model. Takeshita and Kohno (1993)[52] attempted to
identify the soil water retention curve relating # with $,based on the VG model. Si and Kachanoski (2000)[50]
proposed a method of the solution of IP based on anotherfixed functional model by Broadbridge and White(1988)[4]. Solution of IP based on this parameterization isindeed less laborious since the number of the parameters
to be estimated is limited. However, drawbacks arise fromless flexibility or small degree of freedom due to thespecific shape of the functions. There is every possibility
that assumptions made for formulating the functions (e.g.,assumed pore-size distribution) are invalid for the actual
soil hydraulic properties of interest. In use of the VG
model that is most commonly used, the conductivities
derived from this model can be quite sensitive to theshape of the soil water retention curve near saturation
with unfavorable consequences for the computational
performance in the solution of FP (Vogel et al ., 2001[61]).
Under particular conditions where effects of the air-entry
are significant, the VG model could lead to wrong
predictions of the conductivities (Ippisch et al ., 2006[25]).
As an alternative inverse modeling approach toovercome these drawbacks, Bitterich et al . (2004)[3] used９

the free-form parameterization in which particularpolynomial functions were employed to describe the soil
hydraulic properties within a small partitioned range ofthe entire pressure head domain, and in principle thenumber of parameters!the number of degrees of
freedom!could freely be selected. Two different
polynomial functions (quadratic B-splines and piecewise
cubic Hermite interpolation), expressed in terms of thepressure head , were employed to comparatively
investigate their performances in the solution of IP. The
validity of the proposed approach was examined based on
the data from the outflow experiments for soil columns.
Using a sequential programming method, the functions
best at fitting the soil hydraulic properties were identified
through minimizing the objective function which was
defined as a weighted sum of the squared errors between
measured (experimental) and numerically calculated
outflow rates from the bottom of the soil column.
However , no detailed algorithm for parameter
optimization was presented.
This chapter is associated with development of a fieldoriented approach for the inverse estimation of the soil
hydraulic properties which is based on the free-form
parameterization. Assuming the properties are dominantly
one-dimensional and independent of thermal transport in
the soil, i.e., considering an isothermal soil water
movement in a column, a methodology for parameter
estimation is developed and validated which would be
advantageous to the solution of IP. By assembly of thepiecewise cubic spline functions, each of "$ and # isexpressed as a continuous function of the pressure head.
To find the shapes of such functions best expressing headdependency of the soil hydraulic properties, a simulationoptimization algorithm with the aid of the LevenbergMarquardt method is developed which serves to iteratively
solve an optimization problem of minimizing errors
between the observed and computed values of the pressure
head in combination with the embedded simulation
module for forward solutions. The term “field-oriented ” isused to indicate that the errors are measured in terms of
the head rather than the outflow rate which is difficult to
be in-situ observed, and therefore solution of the IP canbe achieved having only recourse to the observed timeseries data of pressure head.
3.2 Forward Problem
The governing equation of one-dimensional saturatedunsaturated groundwater flow is expressed by theRichards equation as follows.
%
$
%
$
$"%
!
!&"#%" #! !!"!
"
%
%
%
%
&
&

(3.1)

where !& is the specific moisture capacity, #% the!
&%"the capillary head (on $"#) or
specific storage, $!
$!
&%" the hydraulic headpressure head (on $#), $!
!
$"$#&"$""!
,
unsaturated hydraulic conductivity, &the height defined as positive upward from the bottom of
the soil and %the time.

１
０

Tomoki IZUMIThe specific moisture capacity !5, the specific storage$3 and the unsaturated hydraulic conductivity #!
&"aredefined as follows.
! .$
!
&"#"
$
&"
%" .&
!5!
$# !
&'#"
#

(3.2)

$ #
&"
%%(!#3 "$
$3!
5"
$3

(3.3)

!#3#2!
&"!
&"#"
$
&"
%"
#!
$#3
!
&'#"
#

(3.4)
Figure 3.1: Cubic splines&"
where #3 is the saturated hydraulic conductivity, #2!
the relative hydraulic conductivity, $the volumetric watercontent, % the water density, ( the gravitational acceleration, #5 the coefficient of water compressibility and #3 thecoefficient of soil compressibility.
Eq.(3.1) is considered subject to the following initialand boundary conditions with the given or specified
values of &#, & (on Dirichlet boundary Γ ') and . (fluxon Neumann boundary Γ )).
&!
!
#"%&# in
1

Ω

(3.5)

&!
!
%& on
1
0"

Γ ' #+(3.6)'!#3#2 ) %. on
'
1

)
Γ #+

(3.7)

Given that the coefficients of the derivatives in Eq.(3.1)
are well parameterized or related to the unknown variable&, it can be numerically solved with respect to & by a
combined use of the standard Galerkin finite element
method and the implicit time-marching scheme (CrankNicolson scheme) which provides unconditionally stable
solutions. This forward solution model is embedded in theparameter optimization process described below, andserves to simulate the flow for the soil hydraulicproperties stepwise updated in the process.
3.3 Inverse Problem
3.3.1 Parameterization
For free-form parameterization, the individual relations
must be expressed by functions which are continuously
differentiable over the whole domain of interest. For thispurpose, as shown in Figure 3.1, the whole domain ispartitioned into ("-1) subdomains with "nodes ("denotes
the degrees of freedom of the parameterization), and eachof the functions for #2 and $ over the whole domain isexpressed by an assembly of piecewise cubic splineinterpolation functions each of which is locally definedover a confined subdomain bounded with two nodes.
Thus, the functions defined over *-th subdomain of&*"&*"$"&"and $*!
&"are described&*"$"!
&*!
(, #2!*!
as follows.

+
+
%
&!&*"
&!&*"
"'
! %*+"&*!
*!
&$ "(*+!
&!&*"
&*!
&*"$"
in (
&"
%"
#2!*!
(3.8)
$#
otherwise
#
$
$
%
&!&*"
&!&*"
"'
!%*$"&*!
*!
&$ "(*$!
&!&*"
&*!
&*"$"
in (
&"
$*!
%"
(3.9)
$#
otherwise
#

$
$
$
$
+
+
+
where %*+, &
*, '
*, (
*, %
*, &
*, '
* and (
* are the coeffi$&*&"" nodal
cients in the cubic splines and *!
number. The shapes of the functions to be determined are
therefore expressed as an assembly of the local functions
so defined.

&"&"%# #/!*!
#2((*!

(3.10)&"&"$((*!
%# $*!

(3.11)**&*"at&*"and $((*!
Hereinafter, the values of #2((*!
the node *are simply denoted by +* and $*, and their sets
by ! and !, respectively. Since, as well known, the soil
hydraulic properties monotonously increase with theincreasing head, Eqs.(3.10) and (3.11) must be identified
so that the following constraints are satisfied.
+*&+*"$,$*&$*"$(3.12)

3.3.2 Solution procedure of IP
At first, a set of decision variables " is defined torepresent both the unknown parameter sets ! and !.
!
!$
!
$
$
,%$
""
! -, %+,"%"(3.13)
!$
!
%
"$
$
$
,%"
""
# -, %$,!" !
Searching for the sets of the optimal solutions, "014, isthen achievable through minimizing the objectivefunction, defined as the total squares error integrated over"") andspace and time between the solution of FP (&-0/ !
0,3&), and therefore written as
the observation data (Inverse Modeling of Variably Saturated Subsurface Water Flow in Isothermal/Non-isothermal Soil

1"!
$0/
#"
"!
#236"
, #236, #%%),(3.14)

with$$ ! %
$ '(
'
"!
#"
)
* #"
%*$$

(3.15)$'*+20 !
!'*2*5#"#"'*!

(3.16)

where %), is the admissible sets of # and $ the total
number of observation data available in the space Ω#'.
In the process of minimizing Eq.(3.14) subject to Eq.
(3.12), the decision variables # are iteratively updated
while step by step solving FP with their assumed orpreviously estimated values. In this respect, identificationof the functions (Eqs.(3.10) and (3.11)) requires a sort ofsimulation-optimization technique. The entire procedure
for solving IP is illustrated in Figure 3.2.

where $is the iteration number, %a coefficient and "a
!
%
%
#!
!"
!"unit matrix. When % is equal to 0, '#!$"
reduces to the Gauss-Newton direction. On the other
hand, when %tends to infinity, '#!$" turns to the steepest
descent direction and size of '#!$" tends to zero.
Therefore, % is taken as 0 for an initial value, and if!!
"!
#!$"$""#!$""is not satisfied, then value of % isincreased and '#!$" is recomputed with Eq.(3.17) until!"!
#!$""is satisfied#!$"$""
the reduction condition "!
[51](Sun, 1994 ).
3.4 Validation
3.4.1 Methods
The validity of the methodology presently developed isprimarily examined with its application to a hypotheticalsoil column as shown in Figure 3.3. The column of 30 cm
high is assumed to be homogeneous, and bounded on its
bottom end by fully saturated groundwater. For theforward solution which is needed in solution searchsequence for parameter optimization, the spatial domain
being considered is discritized into 6 line-elements ofequal length with 7 nodes.
Observation of the pressure head is assumed to be madeat three different locations depicted by the black squaredots in Figure 3.3. Since actual observation data frommeasurements are unavailable, substitutes for them are
numerically generated. For this, the forward problem issolved with the same spatial discretization as in Figure
3.3, assuming that the soil hydraulic properties could bewell represented by the VG model written as
$

$

+7. %

+7." $'"$!!
$!
&-"%#$!!
&-"
#4(&!

&!&4
$
$
,7. +7.
&5!&4 !
"
$"!
#7.**
'"
&
!&
5
4
&(&!
'"
$&4"
,7. +7.
!
$"!
#7.**
"
'"
'"
$
&-!

Figure 3.2: Solution search algorithm

3.3.3 Optimization algorithm
Levenberg-Marquardt method, which is a modified
Gauss-Newton method, is adopted for the optimization
algorithm. With this method, the unknown parameters are
step by step updated with the following search sequence
through the iteration.
#!$"$"$#!$""'#!$"(3.17)

with
!$
&"!$"
'#$$!!
!!$""%
""
$"
!
*

$"
!
*

(
' (
'
&
(
(
-)
(
*$$
$

'
!!$"$%(3.18)(3.19)

１
１(3.20)(3.21)(3.22)

where &- is the effective saturation estimated by the VG
model, &5 the saturated water content, &4 the residual
water content, #7., +7. and ,7. parameters that shape the",7..
soil properties and have relation as +7. $$!$
Considering a one-way process of desorption, the wholedomain of interest is assumed to be initially saturated, andunsaturated due to quick drainage of the water absorbed insoils, i.e., due to quick fall of the phreatic surface fromthe top to the bottom of the column. This means that thesoil hydraulic properties under a condition repeating wet
and dry, i.e., affected by hysteresis, are not considered. In
addition, '$# cm and .$# cm/s are given to thebottom and top boundaries of the column, respectively.
Specification of no-flow through the top boundary iscorrect in reality because the desorption process could becompleted during a so short period of time that effects ofrainfall and/or evaporation can be neglected.
For comparative examinations, two different soil types(Soil A and Soil B) with different soil hydraulicproperties are considered which are characterized by lessand more sensitive changes of the relative hydraulic１
２

Tomoki IZUMITable 3.1: Common parameters
Parameter"/ [cm2/N]"2 [cm2/N]

2.0×10−4

Value

"0-

1.0

"*,2
#0-

0.0
#/

#*,2

#.

4.4×10−6

Table 3.2: Parameters of Soil A
Parameter
#/

Value

#.

0.581

!/ [cm/s]

7.0×10−4

0.801

VG model parameter!1) [cm−1]

0.0268#1)

3.249

Figure 3.3: Soil column
Table 3.3: Parameters of Soil B

conductivity near saturation, respectively. Soil A and Soil
B are those like volcanic ash soil (andosol ) and subsoil ofthe volcanic ash soil (loamy soil), respectively. Values ofthe parameters to be specified for simulation runs togenerate observation data as well as for simulationoptimization runs are summarized in Table 3.1 for thosecommon for both soils, Tables 3.2 and 3.3 for Soil A andSoil B, respectively.
In order to generate observation data for this desorption
process, it is enough in any case to simulate the flow overa period of 5 minutes during which the flow certainly
falls steady!timeindependent. For solving this forward
problem, the time is marched with the relatively small
increment of 0.1 s to reproduce the time-varying flow inthe column at high resolution. All of the time-dependent
discrete solutions of pressure head, obtained at selected
three observation points, are employed as observation dataas they are. This indicates that the generated pseudo
observation data are spatially selected part of the true
solutions to the flow problem under consideration, andtherefore the observation is made at intervals of 0.1 s
having no measurement errors and noises.
Taking into account the abovementioned characteristicsof the individual soils, the entire domain of the pressure
head is partitioned at even intervals into 6 subdomains forSoil A, and more finely, 10 subdomains for Soil B to findthe function with better fitness to the relative hydraulicconductivity which is quite sensitive near saturation.
3.4.2 Results
Figures 3.4 and 3.5 illustrate the identified functions of
the relative hydraulic conductivity and the volumetric$"$" and ###%!
, for Soil A inwater content, !.##%!
$"and #&$!
$"!
comparison with the true functions, .&$!
,
based on the VG model, respectively. Figure 3.6 shows
the timevarying pressure heads at three selected

Parameter
#/

Value

#.

0.218

!/ [cm/s]

4.5×10−3

0.760

VG model parameter!1) [cm ]#1)−10.01151.487

locations!of major concern in practical respects!for thesame soil type, which are part of the solutions obtained
from simulation practice with the identified functions forthe soil hydraulic properties, including those from pseudo
observation!true heads. Figure 3.7 illustrates the$(,+ !$,'/"of thosevariations of the absolute error "computed pressure heads from the observed (true) ones,
which helps closely investigating reproducibility(accuracy) of the solutions resulting from using theidentified functions. Correspondingly, the results for Soil
B are shown in the insert of Figures 3.8 through 3.11.
Function identifiability for Soil A is so high that$" andindistinguishable discrepancies between !.##%!
$"and #&$!
$"as well as between ###%!
$"can be!.&$!
found over the entire domain of pressure head (Figures
3.4 and 3.5). The function of relative hydraulicconductivity for Soil B is, however, perceptibly discrepant
within the limited range of the head near saturation,
though the function of volumetric water content is
identified favorably with its true one (Figures 3.8 and
3.9). This partially discrepant proneness in identificationof the conductivity function could be eliminated byaugmenting the number of subdomains for thecorresponding range. The forward solutions resulting fromusing the duly identified functions for the soil hydraulicproperties are satisfactory with good agreement to the true

Inverse Modeling of Variably Saturated Subsurface Water Flow in Isothermal/Non-isothermal Soil

１
３

Figure 3.4: !%!!# compared with !%$" (Soil A)

Figure 3.7: Absolute errors between computed and truepressure heads (Soil A)

Figure 3.5: !!!# compared with !$" (Soil A)

Figure 3.8: !%!!# compared with !%$" (Soil B)

Figure 3.6: Time-varying pressure heads (Soil A)

Figure 3.9: !!!# compared with !$" (Soil B)

１
４

Tomoki IZUMI

functions in a framework of the simulation-optimization.
Requiring field observations of the time-varying pressure
head alone, the approach is field-oriented, and therefore
less costs us to solve the problems. The approach has
been proved to be valid through its practical application toa hypothetical soil column, identifying the parameterswith well-performance. This suggests that the approachdeveloped could be a viable alternative to theconventional fixed functional form approaches.

Figure 3.10: Time-varying pressure heads (Soil B)

Figure 3.11: Absolute errors between computed and truepressure heads (Soil B)

solutions (Figures 3.6 and 3.10). Perceptible deviationfrom the true solutions!at most 1 mm and 5 mm for Soil
A and Soil B, respectively!occurs during a short periodof time just after the quick commencement of dewateringfrom the column, and dwindles away as time passes(Figures 3.7 and 3.11).
The results shown above prove that the present inverse
method enables us to well estimate the soil hydraulicproperties in a functional form, yielding unique and stable
inverse solutions even in a case of unrealistic
discontinuous flow dynamics like quick desorption.
3.5 Conclusions
In this chapter, an approach for identifying in afunctional form the unknown quantities!soil hydraulicproperties!parameterized in saturated-unsaturated groundwater flow model has been developed. The individual soil
hydraulic properties are interpolated by piecewise cubic
spline functions in a relation between the property and thepressure head. The inverse problems defined are solved as
the problems of determining the coefficients of theInverse Modeling of Variably Saturated Subsurface Water Flow in Isothermal/Non-isothermal Soil

CHAPTER 4
INVERSE METHOD TO IDENTIFYUNSATURATED HYDRAULIC CONDUCTIVITY IN
VARIABLY SATURATED SUBSURFACE WATERFLOW MODEL IN NON-ISOTHERMAL SOIL
― !-BASED FORM OF RICHARDS EQUATION―
4.1 Introduction
Because of increasing demand for numerical analysis
and importance of parameter identification being key tosuccess the analysis, Takeuchi et al . (2007)[53] and Izumi
et al . (2008)[28] proposed an inverse method with the freeform parameterization to estimate soil hydraulic properties
only from the observed time-series field data, and verified
its effectiveness through numerical tests considering
different combinations of the degree of freedom in theparameterization, logging interval and observation error.
Izumi et al . (2008)[28] was based on the governing
equation under isothermal assumption. The method
developed is, therefore, unsuited for the non-isothermal
situation where due to active heat energy exchange
through a soil surface the soil temperature distribution
becomes so non-uniform that it significantly affects water
movement. Generally, it is known that the soil
temperature influences the surface tension and viscosity ofwater. From the field observation conducted by theauthor, in fact, it is claimed that the temperaturedependency of the surface tension leads to a change in thewater flux, estimated by the formula in Kondo and
Saigusa (1994)[33], at the rate of a few percent, while thetemperature dependency of the viscosity to a change inthe saturated hydraulic conductivity at the rate of about
40% in the changeable range of soil temperature. Also,
Inoue et al . (2004)[23] showed that the air temperature
influenced readings of tensiometer, from the distinct
differences between the suction values observed in
daytime and nighttime at the same water content.
In this chapter, an extension of the previous work
(Izumi et al ., 2008[28]) to develop a general inverse
method applicable even for a non-isothermal seepage or
subsurface flow is presented. First, an in-situ observation
system with simple instrumentation is presented which
implements collection of the hydro-geological data!
volumetric water content, suction and soil temperature!
necessary for solving IP. After definition of FP, solution
procedure for the inverse problem in that the relative
hydraulic conductivity (RHC) is identified by use of thefree-form parameterization technique is described, andvalidity of the inverse method developed is examinedthrough its application to the soil of a real upland crop
field, or through calibration tests of the coupled water
movement and thermal conduction models, for thedesorption process in that considering the hysteretic
phenomenon is not needed. In addition, efficacy ofallowing for soil temperature variability!non-isothermal
condition!in the forward problem is verified.

１
５

4.2 Method
4.2.1 Observation system
A typical observation system for acquiring thenecessary data at a few different depths in a soil column
is shown in Figure 4.1. A set of three soil moistureprobes, one tensiometer and two thermometers, which are
drawn by black lines or painted black, is a basic
instrumentation or a minimum requirement for the
purpose of the observation. The remaining instruments,
depicted in gray, are optional. Two or three different
instruments should be paired and aligned in the same
depth.
The data of basic need, obtained from the basic
instrumentation, are put to use as follows. The volumetricwater contents and the soil temperatures at both top andbottom ends of a soil column, observed by the soil
moisture probes and the thermometers, respectively, areused as Dirichlet boundary values for the governing
equations to solve a FP. The volumetric water content
observed by the intermediate soil moisture probe is used
to evaluate the fitness of an unsaturated hydraulicconductivity function!actually a relative hydraulicconductivity function!in an optimization process. A pairof volumetric water content and suction data in the same
depth is used to identify a soil water retention property.
Supplement of one or two optional tensiometers serves
to improve reliability of the soil water retention property.
An additional thermometer in the intermediate depth
provides the benchmark data for verifying reproducibilityof forward solution.
The system is extremely simple. Thus, it is inexpensive
and can be buried in the ground without spending timeand effort, less disturbing the soil column of interest.

Figure 4.1: Basic and optimal instruments for field observations

１
６

Tomoki IZUMI

4.2.2 Forward problem
Assumed to be predominantly one-dimensional, theseepage flow in variably saturated and nonisothermal soil
can be described by coupled water mass (water
movement) and energy balance (thermal conduction)
equations (Kondo and Saigusa, 1994[33]). The mass balance
allowing for effects of soil temperature in soil, neglecting
vapor fluxes, is represented in terms of volumetric watercontent as follows.
)5

+
*5
+
(
$!
+
+
+
,

(4.1)

with
+
(
+
!)5". &3!)5$
+
+
,
,
+
* $
$
"( $$
+
( !5
+
*
". $$
+
&3

(4.3)

$$$2$3(4.5)*5 $!)5"(

$3$

%
)5,
'!
&3"

(4.2)

(4.4)

(4.6)

where )5 is the water density and assumed to beconstant here, ( the volumetric water content, +the time,, the height defined as positive upward, &3 the soil
temperature, * the suction, !5 the specific moisturecapacity, $ the unsaturated hydraulic conductivity, $2 therelative hydraulic conductivity, $3 the saturated hydraulicconductivity, %the intrinsic permeability, ' the viscosity
coefficient and ,the gravitational acceleration.
The energy balance equation assuming a local thermal
equilibrium between soil and water, to be coupled with
Eq.(4.1), is described as follows.
!1*1+&3 $ +
+
+ +
,

(4.7)

with
("
$!(3"
$!
'
'
!1!
3"(
5

(4.8)

+
&
*1 $!& 3
+
,

(4.9)
$

!
(
(&
&$&#"#(4.10)

where !1 is the volumetric heat capacity of soil, (3 the
3 and '
5 the volumetric heatsaturated water content, '
)
$
!
%
)
#
#$
[J/(m3・K)]) and thatcapacity of soil particles (
3
)
!
%
##$
# [J/(m・K)]), respectively and & theof water ('thermal conductivity of soil which is expressed in termsof the reference thermal conductivity &#.
Eqns.(4.1) and (4.7) are considered subject to thefollowing initial and boundary conditions.
(!
"
#"
$(#!
,"
,

in

Ω

(4.11)

"
#"
$&#!
,
,"
&3!

in

Ω

(4.12)

"
(!
"
$(!
,
+"
,
+"

on

Γ*
5

(4.13)

"
"
$&3!
,
+"
,
+"
&3!

on

Γ*
1

(4.14)

"
"
$*5!
,
+"
,
+"
*5!

on

Γ,
5

(4.15)

"
"
$*1!
,
+"
,
+"
*1!

on

Γ,
1

(4.16),"and &#!
,"are the initial value of thewhere (#!
volumetric water content and soil temperature, respec",+"the value of",+"and &3!
tively, Ω space domain, (!
the volumetric water content and soil temperature on the*
Dirichlet boundary, respectively, Γ *5 and Γ 1 the Dirichletboundary for the water and heat movement, respectively,"",+"and *1!
,+"the water and heat flux on the*5!
,
Neumann boundary, respectively and Γ ,5 and Γ 1 the
Neumann boundary for the water and heat movement,
respectively. Given that the coefficients of the derivatives
in Eqs.(4.1) and (4.7) are well parameterized or related tothe unknown variables (and &3, such a equations system
can be numerically solved with respect to (and &3 by the
combined use of the standard Galerkin finite element
method and the implicit time-marching scheme, i.e.,
Crank-Nicolson scheme. This forward solution model isembedded in the parameter optimization process described
below.
4.2.3 Parameterization of soil hydraulic properties
With the observed volumetric water content and suction
data collected from data loggers of soil moisture probes
and tensiometers which are in the same depth from theground surface, a scatter diagram relating volumetricwater content to suction is obtained. To represent in afunctional form the relation between them, or to obtain asoil water retention curve, VG model is employed which
is described as follows.
%/$

(!(2
$
$
)40 (40(3!(2 !
"
$"!
$40%%
*"

(4.17)

where %/ is the effective saturation, (2 the residual
water content and $40, (40 and )40 the unknown#)40.
parameters, the last two being related as (40 $$!$
To determine values of the parameters in Eq.(4.17) so as
to best fit the observed (- * relations, the least-squareapproach is used (Takeshita and Kohno, 1993[52]).
For functional representation of RHC, a free-form
parameterization method is used. For the free-form
parameterization, the relation must be expressed by afunction which is continuously differentiable over thewhole effective saturation domain of interest. For thispurpose, as shown in Figure 4.2, the effective saturation
domain is partitioned into (#-1) subdomains with #nodes(# denotes the degrees of freedom of the parameteriza%/" over the wholetion), and the function for $2++-!
!
(
"
"
!
(
"
%
%
/
2
/
3
'is expressed by a sequence ofdomain &piecewise cubic spline interpolation functions each ofwhich is locally defined over a confined subdomainbounded with two nodes. Thus the function defined over

Inverse Modeling of Variably Saturated Subsurface Water Flow in Isothermal/Non-isothermal Soil１
７

where $,/ is an admissible set of #, % the total numberof observed data available in space and time and 20(normally, taken as unity) a weighting factor.
In the process of minimizing Eq.(4.23) subject to Eq.
(4.21), the decision variables # are iteratively updated
while step by step solving FP with their assumed orpreviously estimated values. In this respect, identificationof the function (Eq.(4.18)) requires a sort ofsimulationoptimization technique.
Levenberg-Marquardt method, which is a modified
Gauss-Newton method, is adopted for the optimization
algorithm. With this method, the unknown parameters are
step by step updated with the following search sequence
through the iteration.
Figure 4.2: Free-form parameterization

the --th subdomain
described as follows.
$

of

'0"
'0#-"%,'0#-!
+, $7!-!
, is-$%(4.18)&"*'0!'0#-"'0!'0#-"
-!
-!
!(-")
'
$ "+-!
'0!'0#-"!
'0#-!
'0&+
'0#-"%,
$"
'0"
$7!-!
(4.19)$$!
#''0&'0#-"%,+0#-!
#- and +-, *- are coefficients in the cubic
where (-, )%%-%""a nodal number. Hereinafter, thesplines and -!
'0"at a node -are simply denoted by /values of $7!!&!
and its sets by #. Since, as well known, RHC
monotonously increases with the increasing saturation, Eq.
(4.18) must be identified so that the following constraints
are satisfied.

(4.20)

4.2.4 Inverse problem and solution procedure
The problem is to optimally decide a set of variables,#, defined as follows.
%%-%"*#$)/-!

(4.21)

Searching for the set of the optimal solutions, #569, isthen achievable through minimizing the objectivefunction, defined as the total least squares error integrated
over space and time between the solution of FP#") and the observed data (&5-8), and therefore,(&.53 !
written as
4#!
!
$32
#"
#569!
#!
#569"
#&$,/

with
(4.26)

$"
$"
!
!
'
'
,
,
0
0
$
/- '
/.
0$% '
%

%
!!$"$#

with/-%/-"%(4.25)!%'#!$"(#!$"$!!
!!$""%"""!%!
$% $7!-!
'0"
'0"

((*
7

#!$"%"$#!$""(#!$"(4.22)

with%(4.27)

where $is an iteration number, %a coefficient and "#" unit matrix. When % is equal to zero, (#!$"
the "reduces to the Gauss-Newton direction. On the other
hand, when %tends to infinity, (#!$" turns to the steepest
descent direction and size of (#!$" tends to zero. The
coefficient %is therefore taken as zero for an initial value,"#!
#!$""is not satisfied, then value of#!$"%""and if #!
% is increased and (#!$" is recomputed with Eq.(4.26)"#!
#!$""is satisfied (Sun, 1994[51]).
#!$"%""
until #!
4.3 Validation
Validity of the inverse method described above is
examined through its in-situ application to the test soil
(sandy soil) in an upland crop field in Miyoshi, Aichi
Prefecture. Primarily the examination is made on
performance of calibrating the numerical forward solution
model!obtained from coupled Eqns.(4.1) and (4.7)!
through identifying free-form parameterized RHC function
with the in-situ observation data, as well as on
reproducibility of the forward solutions obtained from themodel so calibrated. The validity is then verified incomparison with the method commonly used for
identifying RHC. In addition, efficacy of allowing for soil
temperature variability!non-isothermal condition!in theforward problem is examined in comparison with making
the isothermal assumption valid. The following are
preparatory descriptions for these examinations.
As a basis for comparison in parameterization for RHC,
the conventional method using VG model of a fixed-form,“VG method”, is considered in that RHC expressed interms of the effective saturation is given as follows.

%
&
$ %)
#"
,
#!
#"
*
0!
&0$%

(4.23)

1:1%!!
%!'0%#1:1"$'0%#&!
'0"
$7+)!
"

&0.53 !
$20!
!&05-8"#"#",0!

(4.24)

If the isothermal assumption is accepted or Eq.(4.1) is&(4.28)

１
８

Tomoki IZUMI

uncoupled from Eq.(4.7), the governing equation for FP is
simplified with independence of soil temperature as
follows.
%
%&
%
"
"!
%
%
&
'
%
"
%& "!#&"" !#&#
%
'
%
$ #
"
"" "#
%
" !&

#&

#"###$(4.29)(4.30)(4.31)(4.32)

where #$ takes a constant value at soil temperature of
15 ℃. RHC, ##, is estimated in the same manner as in
the non-isothermal case.
4.3.1 Field observation and soil water retention curve
fitting
The observation system with full instruments, as shownin Figure 4.3, was in-situ set up. The target domain is thesurface soil of 20 cm thick!from -10 cm to -30 cm!,where soil temperature varies perceptibly. Assumed
predominantly one-dimensional, the domain of a soilcolumn is divided into four equal elements with fivenodes. The instruments sense and automatically record atintervals of 10 min the variations of volumetric watercontent, suction and soil temperature at each depth of 10
cm, 20 cm and 30 cm below the ground level. To know
the response of soil moisture variation to occurrence ofrainfall events, an automatic rain gauge is placed close by
the system to record time-varying precipitation at thesame time interval.
The complete serial data of volumetric water content,
suction and soil temperature, observed during the periodof June 8, 2008 to July 9, 2008 and selected for thepresent examinations, are illustrated in Figure 4.4
including that of the rainfall. Here, a one-way process ofdesorption is considered for the estimation of the soil
hydraulic properties. The distinct desorption periods,
detectable at a glance from the suction variation in Figure
4.4, are those of June 12 to 19, June 24 to 28 and July 1
to 8. These are considered as mutually independent time
sequences to provide the different sets of the observed
serial data and thus to make more persuasive examinationsof the present propositions.
The examinations must be preceded by finding the soil
water retention curve Eq.(4.17), for each desorption
period, best fitting the scattered relations between
observed " and $. Here, "$ is taken as the value for$"!, and "# as the minimum value found in all theobserved serial data of ", including the data shown in
Figure 4.4. Thus, values of "$ and "# can directly be
identified to be 0.48 and 0.22, as common to threedifferent desorption periods, through analyzing theobserved serial data of $ and ", and of ", respectively.
Caution must be paid in identifying values of !%" and $%"by use of the least squares method. Theory predicts that $and " are uniquely correlated, or that, if under the same

Figure 4.3: Full instrumentation

Figure 4.4: Observed data

Figure 4.5: Correlation between suctions in daytime and nighttime

situation of the water content, the readings of $ must be
the same. In real observation, however, the daytime andnighttime readings of $ for the same value of " are atvariance, as exemplified in Figure 4.5, the former being
relatively reduced. This is a result of that the suction
sensors (UIZ-SMT, UIZIN Co., Ltd.) are positioned near
the mouth of a tensiometer tube and therefore prone to１
９Inverse Modeling of Variably Saturated Subsurface Water Flow in Isothermal/Non-isothermal Soilobserve the suctions affected by air temperature. For this,the observed #- $ relations over the nighttime of 20:00 to
8:00, less affected by air temperature, are exclusively
qualified as those to be fitted with a curve, discarding
those over the daytime. All the relations qualified and theresulting curves fitted with Eq.(4.17) are illustrated in
Figures 4.6 to 4.8 for the respective desorption periods.
The resultant identified values of "4. and '4. are listed inTable 4.1 including the value of #3$& that was measured
through laboratory tests for the soils sampled from thesite.

Figure 4.6: Soil water retention property in the first desorption
periodFigure 4.7: Soil water retention property in the second
desorption period

Figure 4.8: Soil water retention property in the third desorption period

Table 4.1: Identified parameter values
Desorption period#3$& [m/s]

6/12-6/19
6/24-6/28

4.5×107/1-7/8−4"4. [/m]

'4.

17.1

1.142

49.8

1.084

10.3

1.172

#3$& indicates #3 at 15 ℃.

4.3.2 Results and discussion
The results of the examinations concerning calibrationof the numerical model and solution reproducibility(accuracy) of the model calibrated are shown in the insert
of Figure X (X =4.9, 4.10, 4.11) for the respective
desorption periods. Note that for each all the results
except for #/31 in Figure X (c) are those for non-isothermal
condition. It should be also noted that for all (& in Eq.
(4.24) is equalized to unity. Figure X (a) comparatively
illustrates the RHC functions #2'') and #2*( which wereidentified by the present free-form parameterization andthe fixedform parameterization based on the VG model,
respectively. Including the time-varying #1+3 observed at
the intermediate depth of the soil column, Figure X (b)shows the corresponding #,10 and #*( computed as thesolutions of FP with the identified #2'') and #2*(,
respectively. In the lower half of the same, the absolute#,10 !#1+3# and !4. "##4.!#1+3# areerrors !,10 "#shown to demonstrate the time-varying difference insolution reproducibility between free- and fixed-form
parameterizations. Figure X (c) illustrated in the samemanner as in Figure X (b) is inserted to differentiate thenon-isothermal solution #,10 from the isothermal solution#/31 produced with the identified #2''). The daily periodic
variations of computed %3,10 and observed %31+3 at theintermediate depth are comparatively shown in Figure X(d), including the time-varying soil temperatures observed!
$) and bottom ()"!#!
%) of the soilat the top ()"!#
column, which are given as Dirichlet boundary values insolving FP. Minimum values of the objective function"!
!", which are reached in simulation-optimization
processes for identifying RHC function #2 and thus could
be a global indicator of solution reproducibility, aresummarized in Table 4.2, including those for the threedifferent combinations of free- or fixed-form parameterization and isothermal or non-isothermal setting for each of
the different three desorption periods.
From Figures 4.9(a), 4.10(a) and 4.11(a), it can readily$-"is relaxed
be seen that the shape of the function #2!
in the free-form parameterization, turning upward at alower level of effective saturation than in the conventional
fixed-form parameterization. From Figures 4.9(b), 4.10(b)and 4.11(b), it can be recognized that the model calibrated
with such a parameter identification technique reproduces
the forward solutions in better agreement with theobserved data. Table 4.2 demonstrates in a quantitative
sense that this is indeed true for any of desorption
periods. This could be strong evidence that the method
presently developed for solving IP surpasses the２０

Tomoki IZUMI

(a)

(a)

(b)

(b)

(a)

(c)

(a)

(d)

Figure 4.9: Calibration result for the first desorption period

Figure 4.10: Calibration result for the second desorption period

Inverse Modeling of Variably Saturated Subsurface Water Flow in Isothermal/Non-isothermal Soil

２
１

commonly used one which employs the well-defined VG
model for identification of RHC. From Figures 4.9(c),4.10(c), 4.11(c) and Table 4.2, it is also found that taking
into account the non-isothermal condition leads to theincrease in solution closeness to reality.
!"in simulationTable 4.2: Minimum values of !!
optimizarion runs for different periods
Desorption period
Non-isothermal/free-form

6/12-6/19 6/24-6/28

7/1-7/8

0.02390.003450.0144

Non-isothermal/fixed-form (VG) 0.0247

0.01188

0.03150.00597

0.0258Isothermal/free-form

0.0242

(a)

(b)

(c)

(d)
4.11: Calibration result for the third desorption period

4.4 Conclusions
A method for solving the inverse problem of variably
saturated seepage flow in non-isothermal soil has beendeveloped. A utilitarian in-situ observation system withsimple instrumentation has specially been contrived which
implements collection of the data necessary for solving
the inverse problem. The results of the examinations on
validity of the method reveals that the relative hydraulicconductivity of major parameter, described as a functionof the effective saturation, could successfully be identified
with a free-form continuous function formed by asequence of piecewise cubic spline functions whose
coefficients are found with the aid of the simulationoptimization technique. The results also show that thewater movement model provides the forward solutions ofhigh reproducibility, when coupled with thermal
conduction model and calibrated through identifying therelative hydraulic conductivity with such a freeformfunction in lieu of a fixed-form function. It is thus
concluded that the present approach of taking into accountthe effect of thermal conduction on water movement invariably saturated soil and employing the free-form
parameterization technique could be the viable alternative
to the conventional tedious and laborious approaches.
Additionally it should be mentioned that the inexpensive
observation system specially contrived enhances
usefulness of the approach.

２
２

Tomoki IZUMICHAPTER 5
INVERSE METHOD TO IDENTIFY UNSATURATEDHYDRAULIC CONDUCTIVITY IN VARIABLY
SATURATED SUBSURFACE WATER FLOW MODEL― MIXED FORM OF RICHARDS EQUATION―

mass conservation equation and the Darcy-Buckingham’slaw, and its mixed form is represented as follows(Huyakorn and Pinder, 1983[20]).
())", "" " )%)#! !!! #"
, *
)
)
)
)
(
(
'
'

(5.1)

with
5.1 Introduction
Richards equation can be expressed in three different
forms according to the state variables solved: the (-based
form, the &-based form, and the mixed form. Since the &based form is not applicable for the saturated zones where
the hydraulic diffusivity takes infinite values, the (-based
form which can express all possible ranges of saturation isoften used. However, Celia et al . (1990)[6] reported thatnumerical models which solve (-based form produce
significant errors in mass balance, and proposed analternative model which solves the mixed form perfectlyconserving the mass. After the work of Celia et al .
(1990)[6], the mixed form of Richards equation is often
used as the governing equation for seepage flow, and thecomputational efficiency of its numerical model isimproved (Huang et al ., 1996)[19]. Almost all thenumerical models proposed by Celia et al . (1990)[6], andsubsequently others, deal with only unsaturated zone,
while Takeuchi et al . (2008)[54] developed a numericalmodel which can solve the mixed form of Richards
equation for saturated zone as well as for unsaturated
zone.
Izumi et al . (2008[28], 2009[27]) have developed inverse
modeling for both the (-based form and &-based form ofRichards equation. In this chapter, an inverse modeling forthe mixed form using the free-form parameterization
technique is proposed. Firstly, the FP employing themixed form of Richards equation as the governing
equation is defined. Because the saturated hydraulicconductivity could be definitely known, the IP related tothe unsaturated hydraulic conductivity is reduced to a
problem of identifying the RHC. Secondly, RHC isrepresented by a free-form parameterized function which
is a sequence of piecewise cubic spline functions over thewhole effective saturation domain. Thirdly, the inverse
problem of minimizing errors between the observed andcomputed values of the pressure head is formulated basedon a simulation-optimization algorithm with the aid of theLevenberg-Marqurdt method to find optimal values of theparameters included in the RHC function. Finally, theinverse method developed is validated through applying it
for in-situ soil column in an upland crop field andcomparing the observed water movement with thecomputed one obtained from the forward simulation forthe desorption process in that considering the hystericphenomenon is not needed.
5.2 Forward Problem
5.2.1 Governing equation
The subsurface water flow in variably saturated soil canbe described by Richards equation which combines the$*"%$,""*#',*!

(5.2)!#!)!*(5.3)

where % is the porosity, ", the saturation, "* thespecific storage, ( the pressure head, ! the unsaturated#("("the hydraulic head, 'hydraulic conductivity, #!
the time, (the height defined as positive upward, ', thewater density assumed to be constant here, * thegravitational acceleration, $* and $, the compressibility
coefficients of soil and water, respectively, !) the relative
hydraulic conductivity and !* the saturated hydraulicconductivity.
Eq.(5.1) is subject to the following initial and boundary
conditions.
(!
!
#"
#(#!
("
(

in

!
(!
!
#(!
(
'"
(
'"

on

)
!
!!*!) # #&!
(
'"
)
(

Ω

(5.4)

Γ%
on

(5.5)
Γ&

(5.6)

("is the initial value of the pressure head inwhere (#!
!
('"the value of the pressure
the space domain Ω, (!
!
('"the waterhead on the Dirichlet boundary Γ % and &!
flux on the Neumannboundary Γ &.
Given that the coefficients of the derivatives in Eq.(5.1)
is well parameterized or related to the unknown variables(, such a equations system can be numerically solved
with respect to ( by the combined use of the standard
Galerkin finite element method and the fully implicit
time-marching scheme mentioned below. This forward
solution procedure is embedded in the parameter
optimization process.
5.2.2 Parameterization of soil hydraulic properties
To represent a constitutive relation between thevolumetric water content and pressure head, or to obtain asoil water retention curve, VG model is employed which
is described as follows.
"'#

&!&)
$
#
%+( $+(&*!&) !
"
$"!
#+($$
("

(5.7)

where "' is the effective saturation, & the volumetricwater content, &) the residual water content, &* thesaturated water content and #+(, $+( and %+( the unknown"%+(.
parameters, the last two being related as $+( #$!$
To determine values of the parameters in Eq.(5.7) so as tobest fit the observed &- ( relations, the least-square２
３

Inverse Modeling of Variably Saturated Subsurface Water Flow in Isothermal/Non-isothermal Soil

Picard iteration scheme can be written as&71"%!0"%!&71 &71"%!0&5 %1"%!0 !%1
"
#
&4
&4
1"%!0
&
%
% &
&
$1"%!0 #$! " "$1"%!0(5.11)
#!
#
# &
&
&
6
6
6
1"%!0"%

isAccording to Celia et al . (1990)[6], &7
expanded with respect to % by Taylor series expansion,
ignoring the higher order, as follows.
!
1"%!0
!
!
%1"%!0"%!%1"%!0 "&71"%!0"%#&71"%!0 "*&7!
!
*%!1"%!0 1"%!0"% 1"%!0!
%"(5.12)#&71"%!0 " 7!%#
where !7 is the specific water content. Substituting thetruncated series Eq.(5.12) into Eq.(5.11)

Figure 5.1: Free-form parameterization

approach is used.
For functional representation of RHC, the free-form
parameterization method is employed. For such aparameterization, the relation must be expressed by afunction which is continuously differentiable over thewhole effective saturation domain of interest. For this, as
shown in Figure 5.1, the effective saturation domain is!%) sub-domains with " nodes ("partitioned into ("denotes the degrees of freedom of the parameterization),&."over the whole domain
and the function for $4(()!
$4"
!&.!
$5"
#$
#%
&.!
)
*is expressed by a sequence ofpiecewise cubic spline functions each of which is locallydefined over a confined sub-domain bounded with two
nodes. Thus, the function defined over the -th sub-domain&."is described as follows.
&.#-"%*&.#-!
, $4!-!
of )"!%#$ $4!-!
&."&."$4(()!
-#%

(5.8)

#
&1"%!0 !&71
%1"%!0"%#%1"%!0 "&4 1"%!0 "! 7
&4
!7
1"%!0
1"%!0
1
&
%1"%!0!% % &&5 %&" " "$1"%!0
! 7
#
# &
#
&
&4
6
6
1"%!0
&
" $
(5.13)##&
6
When a soil domain of interest is saturated, &7 # const.
#%. Thus, the first term on the left-hand side of Eq.(5.1)
disappears, and the equation becomes as follows.
&%&&&5 #! "!$ ,#
&
&
&
6
6
4

Employing the modified Picard scheme, Eq.(5.14) canbe written in the time-marching form as follows.
&%1"%!0% &%1"%!0"%#%1"&4 " "$1"%!0
#
&
&5 &
6
6
& 1"%!0
" $
#
&
6

with
.
.
&
")
&.!&.#-"&.!&.#-"!'-."(-!
-!
'$ "*-.!
&.!&.#-"!
&.#-!
&.%)
&.#-"%*
#"
&."
$4!-!
(5.9)
$$
!
"&
&.%
&.#-"%*
).#-!
#
.
.
.
where '-., (-, )- and *
- are coefficients in the cubic%$-$""a nodal number. Hereinafter, thesplines and -!
&."at a node -are simply denoted by .-.
values of $4(()!
Since, as well known, the relative hydraulic conductivity
monotonously increases with the increasing saturation, Eq.
(5.8) must be identified so that the following constraints
are satisfied.

.-$.-"%(5.14)(5.15)

5.3 Inverse Problem
5.3.1 Formulation of optimization problem
Solving IP is defined as optimally deciding a set of
%$-$"
(
.-!
unknown variables, !#'. The objectivefunction is then defined as the total least squares error
integrated over space and time between the solution of FP!") and the observed data (%2+5), and therefore,(%,20 !
written as, !%$*1#!
#0/
!"
#!
!234"

(5.16)

with(5.10)

5.2.3 Time-marching scheme
The governing equation is discretized by a standard
Galerkin finite element method for space and a finite
difference method for time. As the time-marching scheme,
the modified Picard method proposed by Celia et al .
(1990)[6] is used. If the superscripts 0 and 1 stand for an
iteration step and a time level, respectively, then the

%

%
&
# $'
!"
+
#!
!"
(
/!
&/#%

(5.17)%/,20 !
#5/!
!%/2+5"!"!"+/!

(5.18)

where !236 is a set of ! optimal solutions, $*- anadmissible set of !, % the total number of observed data
available in space and time and 5/ (normally, taken as
unity) a weighting factor.

２
４

Tomoki IZUMIIn the process of minimizing Eq.(5.17) subject to Eq.
(5.10), the decision variables, #, are iteratively updated
while step by step solving FP with their assumed orpreviously estimated values. In this respect, identificationof the function (Eq.(5.8)) requires a sort of simulationoptimization technique.
5.3.2 Optimization algorithm
For the optimization algorithm to search for the set ofthe optimal solutions, #$%&, Levenberg-Marquardt method,which is a modified Gauss-Newton method, is adopted.
With this method, the unknown parameters are step bystep updated with the following search sequence throughthe iteration.
#!""#"$#!"""&#!""(5.19)

with
!#
%"!""
&#"$!!
!""#
""
""
""
!
!
$
$
$
$
(
(
"
'% $
'&
($# $

(5.20)

Figure 5.2: Observation system and discretization of targeteddomain##!!""$!

5.4 Validation
An in-situ application to the test soil (sandy soil) in an
upland crop field in Miyoshi, Aichi Prefecture, Japan, is
implemented to examine validity of the inverse method
described above. The validity is judged referring to
reproducibility of the forward solutions obtained from thecalibrated model.

between saturated and unsaturated situations withfrequency, and is discretized into four equal elements withfive nodes for numerical computations. The time-varying
pressure heads and volumetric water contents at the threedepths of -10 cm, -20 cm and -30 cm are automaticallylogged at intervals of 10 minutes.
A one-way process of desorption is considered for
estimation of the soil hydraulic properties to exclude thehysteric phenomenon, as in most of the preceding studies
associated with parameter identification. The data
observed during the no-rainfall period (or desorption
period) of July 1 to 9, 2008, used as basic data forcalibration of the numerical model as well as forvalidation of the inverse method presently developed, areillustrated in Figure 5.3. From Figure 5.3, it can readily
be seen that the pressure head involves periodic daily
variations while declining in time, having a distinct
difference from the volumetric water content. Even if thevolumetric water contents in daytime and nighttime are

5.4.1 Observation system
The data necessary for identification of RHC is
acquired by the observation system made up of different
three sets of instruments (tensiometer and soil moistureprobe), which are buried at different three depths in a soilcolumn (Figure 5.2). At individual depths, serial data ofpressure head and volumetric water content are observed
which are necessary to identify the soil water retention
property. The observed pressure heads at both top andbottom of the soil column are also used as Dirichletboundary values in solving FP based on the governing
equation. The pressure head observed at the intermediate
depth (-20 cm deep) is used to judge the fitness of
estimated unknown parameters (or RHC function) in an
optimization process.
The targeted domain is the surface soil of 20 cm thickbetween -10 cm and -30 cm deep, where soil alternates

Figure 5.3: Observed data for validation(5.21)

where "is an iteration number, #a coefficient and "#! unit matrix. When # is equal to zero, &#!""
the !
reduces to the Gauss-Newton direction. On the other
hand, when tends to infinity, &#!"" turns to the steepest
descent direction and size of &#!"" tends to zero.
Therefore, #is taken as zero for an initial value, and if!"!
#!"""is not satisfied, then value of # is"!
#!""#""
increased and &#!"" is recomputed with Eq.(5.20) until!"!
#!"""is satisfied (Sun, 1994[51]).
"!
#!""#""

２
５

Inverse Modeling of Variably Saturated Subsurface Water Flow in Isothermal/Non-isothermal Soil

Eq.(5.7) can be best fit to the observed relations between$ and #. Note that in this curvefitting the temperaturedependent #- $ relations observed during the daytime of
6:00 to 18:00 are not considered. Eventually the values of"0+ and %0+ are identified to be 2.87 m−1 and 1.29,
respectively, to have a best fitting curve as shown in
Figure 5.5.
5.4.3

Figure 5.4: Correlation between pressure heads in daytime andnighttime for the same volumetric water content

the same, therefore, the pressure heads are different asshown in Figure 5.4. This result is against the theory that
there should be a one-to-one correspondence between thepressure head and the volumetric water content. As Inoueet al . (2004)[23] reported, this is because the sensing partof this type of a tensiometer (UIZ-SMT, UIZIN Co., Ltd.)is put near the top of the tube, and therefore the readingsof the pressure head are significantly affected by air
temperature. Thus, as described later, identification of
RHC function is comparatively done considering anotherdata series of the pressure head in which temperaturedependency of the pressure heads in the daytime of 6:00
to 18:00 is removed by means of the linear interpolation
method.
5.4.2 Identification of water retention property
Identifying RHC function must be preceded by fixing
the form of the soil water retention curve expressed by
Eq.(5.7). The saturated water content #/ is defined as thevalue of #in case that $is equal to zero, and the residual
water content #. is defined as the minimum value in all
the observed serial data of #. Thus, straightforwardly
from the data observed, #/ and #. can be identified to be
0.45 and 0.22, respectively. By use of the least squares
method proposed by Takeshita and Kohno (1993)[52], theremaining parameters α and %0+ are determined so that

Identification of RHC and validation of themethodAs mentioned above, the time-varying pressure headobtained from a field observation is significantly affectedby air temperature. Different two data sets of time-varying
pressure head are thus considered for identifying RHC (orcalibrating the model) by free-form parameterization andexamining solution reproducibility (accuracy) of the modelso calibrated. One is the observed data as they are, andthe other is the corrected data which can be obtained fromfiltering out daytime temperature-dependent noises fromthe observed data, or from assuming that the pressure
head varies linearly from the observed head at 6:00 to thatat 18:00. These are referred to as〈Data A〉and〈Data B〉,
respectively.
!
&"$#!% m/s from the
The value of #/ is given as %result of the laboratory tests for the soils sampled fromthe site. In application of Eq.(5.18), all &$ are equalized
to unity.
The results for〈Data A〉and〈Data B〉are shown in
Figures 5.6 and 5.7, respectively, illustrating the identified
RHC function #.''( and the solution reproducibility. The
solution reproducibility is examined through comparing
the computed $*-, (or the forward solution with theidentified #.''() with the non-corrected (observed) or
corrected $-)/, and measured in terms of the absolute$*-, !$-)/$. Values of the objective function
error !*-, #$"!
!", which are reached in the simulation-optimization
processes for identifying RHC function and thus could bea global indicator of solution reproducibility, are listed inTable 5.1.
From Figure 5.6 (b), it can be seen that the model
calibrated for〈Data A〉produces the forward solutions inpartially good agreement with the original observed data.
From Figure 5.7 (b), on the other hand, it is recognizedthat the model calibrated for〈Data B〉reproduces thecorrected observed data in closer agreement. Table 5.1
also demonstrates the same in a quantitative sense. It isthus concluded that the method presently developed could
be a viable tool for finding RHC in a functional form,
fulfilling its function with better performance if airtemperature dependency of the pressure head variation isappropriately eliminated.

!"in simulation-optimization
Table 5.1: Minimum values of "!
Data typeFigure 5.5: Soil water retention property!""!

〈Data A〉

1.88

〈Data B〉

0.737

２
６

Tomoki IZUMI

(a)

(a)

(b)

(b)

Figure 5.6: (a) Identified RHC and (b) Reproducibility of FP
solution for〈Data A〉

Figure 5.7: (a) Identified RHC and (b) Reproducibility of FP
solution for〈Data B〉

5.5 Conclusions
A method for solving the inverse problem associatedwith variably saturated seepage flow has been developed.
The mixed form Richards equation which provides massconservative solutions to saturated-unsaturated seepage
flow problem is considered as the governing equation
which describes the forward problem. Examinations forvalidation of the method developed show that the relative
hydraulic conductivity which is a major hydraulicproperty of the governing equation system, expressed by afree-form sequential piecewise cubic spline function,
could successfully be identified through determining thevalues of the coefficients of the function with the aid of
the simulation-optimization technique. The results also
show that the forward solution model with identified
relative hydraulic conductivity could reproduce actual oneway desorption process with higher accuracy when airtemperature dependency of the observed pressure headvariation is appropriately eliminated.

Inverse Modeling of Variably Saturated Subsurface Water Flow in Isothermal/Non-isothermal Soil

CHAPTER 6
INVERSE METHOD TO IDENTIFY UNSATURATEDHYDRAULIC CONDUCTIVITY IN VARIABLY
SATURATED SUBSURFACE WATER FLOW
MODEL IN NON-ISOTHERMAL SOIL
―MIXED FORM OF RICHARDS EQUATION―
6.1 Introduction
Richards equation expresses the water movement
through soil derived from the gradient of hydraulic head
under isothermal assumption. However, water movement
is also affected by soil temperature. Thus, this equation is
not suitable to represent water movement, especially, nearsurface soil where heat energy exchange is quite
significant. As extended models, coupling models for
simultaneous transfer of water and heat based on the workof Philip and de Vries (1957)[46] have been proposed (e.g.,
Milly, 1982[36]; Kondo and Saigusa, 1994[33]). Fujinawa(1995, 2010)[14],[15] proposed another coupling model. Whilethese models are related to the '-based form or %-based
form of Richards equation coupled with heat transport
equation, few coupling models of the mixed form ofRichards equation and the heat transport equation has
been studied. In this chapter, thus, an inverse modelingfor the mixed form of Richards equation in nonisothermal soil is purposed. Firstly, the governing
equations system (i.e., FP) is described, and the soil
hydraulic properties which are model parameters includedin the system are parameterized. The RHC which is amajor unknown parameter to be identified in this study is
described by a free-form parameterized function which isa sequence of piecewise cubic spline functions over thewhole effective saturation domain. For the representationof SWRC, VG model is employed due to being timeproven. Secondly, the IP is defined as minimizing errors
between the observed and computed values of the pressure
head based on a simulation-optimization algorithm withthe aid of the Levenberg-Marqurdt method to determine
the functional shape of RHC. Finally, the applicability of
the inverse method developed is assessed through in-situ
soil column experiments in terms of reproducibility forobserved water movement during the desorption process
in that considering the hysteric phenomenon is notneeded.
6.2 Governing Equations System
6.2.1 Water movement model
For the water movement, the mixed form of Richards
equation derived from the generalized Darcy’s law and theBoussinesq assumption (Huyakorn and Pinder, 1983[20]) is
employed to obtain the mass-conservative numericalsolutions. Considering the dependency of density andviscosity of water on soil temperature to couple withthermal transport, the equation in onedimensional verticaland saturated-unsaturated flow which in the liquid phasehas considerable magnitude, i.e., neglecting the vapor
fluxes, is described as follows;

２
７

'
(
( &*!&#/ "# # (#($! !!"! &"/ .
"" (6.1)
(
(
(
&(
)
)
(
(
with
"."#"/"#.$&-)!

(6.2)

#+"
"*!
$."
".
"$"-!

(6.3)

'
")$'")
&-)

(6.4)

&$

where # is the porosity, #/ the saturation, #. thespecific storage, ' the pressure head, " the unsaturated
hydraulic conductivity, & the hydraulic head, (the time, )the height defined as positive upward, &* the water
density at the soil temperature $., &- the reference water
density at the reference soil temperature $-, ) thegravitational acceleration, ". and "/ the compressibility
coefficients of soil and water, respectively, "- the relative
hydraulic conductivity, "* the correction-factor functionof soil temperature, ". the saturated hydraulic conductivity, #+ the effective saturation and 'the water pressure.
6.2.2 Thermal transport model
The heat flux due to the water movement in soil issmaller than the heat conduction by the solid soil and thus
can be neglected. Accordingly, the heat conduction
equation assuming a local thermal equilibrium between
soil and water is employed for the thermal transport, anddescribed based on Kondo and Saigusa (1994)[33] asfollows;
(!
!&$."
(
(
$! !!$ $."
(
(
(
(
)
)

(6.5)

with
$!#"
!, $!
%
%
."%
/
$

!
(
%&
$$$#"#

(6.6)
(6.7)

where !, is the volumetric heat capacity of soil, % the. and %/ the volumetric heatvolumetric water content, %!
%)#$#) [J/(m3・K)]) and thatcapacity of soil particles ($!
%##$#) [J/(m3・K)]), respectively and $ theof water ('thermal conductivity of soil expressed in terms of thereference thermal conductivity $# and volumetric watercontent.
6.2.3 Parameterization of soil hydraulic properties
Soil hydraulic properties to be identified are theunsaturated hydraulic conductivity and soil water retention
curve.
The unsaturated hydraulic conductivity is described as
the product of three variables shown in Eq.(6.3). The
correction-factor function of soil temperature and saturated
hydraulic conductivity are represented as follows.

２
８

Tomoki IZUMI

&.
&*
(.,$
"/#
&.
"* #

(6.8)
(6.9)

where &. and &* are the dynamic viscosity coefficientat temperature $. and $/, respectively and $the intrinsic
permeability. Because the dynamic viscosity is thefunction of the temperature and the saturated hydraulicconductivity can be determined through laboratory
experiments by definition, the parameterization of RHC isneeded.
To represent RHC, a sequential piecewise cubic spline
function is employed. This is referred to as a free-form
approach indicating that no a priori shape of RHC isassumed except for their monotonicity (Bitterlich et al .,2004[3]). In the free-form approach, as shown in Figure
6.1, the effective saturation domain of interest is!%) subdomains with ! nodes (!
partitioned into (!
denotes the degrees of freedom of the parameterization),#+"over whole domain '
'."
#+!
,
and the function for ".!
'/"(is expressed by a summation of piecewise cubic#+!
spline interpolation functions which is locally defined overa confined subdomain bounded with two nodes. Thus, the#+"with ".#*!
#+"defined over the *thfunction for ".!
#+#*, #+#*"%(is described as follows;
subdomain '!
!%#! ".!*!
#+"
#+"
".!
*#%

(6.10)

with&"'#+!#+#*"#+!#+#*"
*!
*!
!%*"&
'
$ "(*!
#+!#+#*"!
#+#*!
#+%'
#+#*"%(
#"
#+"
".!*!
(6.11)$$!
"##+#*"%(#+%'+#*!
#* and (*, '* are coefficients in the cubic
where %*, &%$*$!"a nodal number. Hereinafter, thesplines and *!
#+"at a node *are simply denoted by +*.
values of ".!
Since, as well known, RHC monotonously increases with
the increasing saturation, Eq.(6.10) must be identified sothat the following constraints are satisfied.

+*$+*"%(6.12)

SWRC is a constitutive relation between the volumetricwater content and pressure head. Various models havebeen proposed to represent SWRC. Among them, VG
model has been more frequently used for the simulationof water movement in soil, and thus, VG model is
employed for the representation of SWRC in this study.
Additionally, to account for the effect of soil temperaturewhen calculating ' using VG model, the value of *obtained from Eq.(6.1) is multiplied by two temperaturecorrection coefficients. One is the ratio of the surface
tensions at the soil temperature of interest and a reference
temperature, and another is that of water density. This
results in the following;
*."
#
#+!

%
'!'.
-0, ,0, #
'/!'.
"
!
%"!
#0,&&
*."

with
%
,0, #%!
-0,) (**.# . *)* (.

(6.14)(6.15)

where '. is the residual water content, '/ the saturated
water content, #0,, ,0, and -0, the unknown parameters,*. the pressure head at a reference soil temperature $. and)* and ). the surface tension at soil temperature $/ and areference temperature $., respectively.
6.2.4 Numerical solution procedure
After discretization with the combination of thestandard Galerkin finite element method for space and thefinite difference method for time, Eqs.(6.1) and (6.5) are
subjected to the following initial and boundary conditions
and numerically solved with the iterative partitioned
method;
*!
!
$"
#*$!
0"
0

in

Ω

(6.16)

!
$"
#$$!
0
0"
$/!

in

Ω

(6.17)

!
*!
!
#*!
0
/"
0
/"

on

Γ(
1

(6.18)

!
!
#$/!
0
/"
0
/"
$/!

on

Γ(
-

(6.19)

+
!
)!
0
/"
!
!"/".
#.1!
0
/" on Γ )
1
+
0
+
!
$/!
0
/"
!
!%
#.-!
0
/" on Γ )+0

Figure 6.1: Free-form parameterization(6.13)(6.20)(6.21)

0"and $$!
0"are the initial value of thewhere *$!
pressure head and soil temperature, respectively, Ω the!
!
0/"and $/!
0/"the values of thespace domain, *!
pressure head and soil temperature on the Dirichlet(boundary, respectively, Γ (1 and Γ - the Dirichlet boundary
for the water movement and heat transport, respectively,!
!
0/" and .-!
0/" the water and heat flux on.1!
)
Neumann boundary, respectively and Γ )1 and Γ - theInverse Modeling of Variably Saturated Subsurface Water Flow in Isothermal/Non-isothermal Soil

Neumann boundary for the water movement and heattransport, respectively.
6.3 Parameter Identification Procedure
The unknown parameters to be identified are the soil
hydraulic properties, i.e., RHC and SWRC, as describedabove. If the parameters can be determined with certainty
in advance, they should be treated as known parameters
because the number of identified parameters should be
kept as small as possible. RHC, in general, cannot be
determined by the direct measurement while SWRC canbe obtained from the time-series data of the pressure headand volumetric water content with relative ease. Hence,
RHC is treated as the unknown parameter and is identified
with use of inverse technique.
6.3.1 Inverse problem
To solve IP is to optimally decide the unknown
parameters, #$(&!#%&%!, which minimize theobjective function defined as the total least squares error
integrated over space and time between the solution of FP#") and the observed data ($,&.), and thereby IP is($',* !
defined as follows;
+"!
!
$*)
#"
#,-/!
"!
#,-/"
#&#%(

(6.22)$)',* !
$*)!
!$),&."
#"
#"
%
)!

２
９

(6.24)

where #,-/ is the set of optimal solutions, #%( anadmissible set of #,-/ and #, $ the total number ofobserved data available in space and time and *)(normally, taken as unity) a weighting factor.
The process of solving IP is shown in Figure 6.2. The
decision variables, #, are iteratively modified or updated
while step by step solving FP with their assumed orpreviously estimated values. In this respect, identificationof RHC function (Eq.(6.10)) requires a sort of simulationoptimization technique.
6.3.2 Optimization algorithm
Levenberg-Marquardt method, which is a modified
Gauss-Newton method combining the Gauss-Newton
method and the gradient method, is employed for the
optimization algorithm to search for the set of optimal
solution with following search sequence through theiteration (Sun, 1994[51]);#!""#"$#!"""(#!""(6.25)
with
!#
'"!""
(#!""$!!
!!"""#
""
""
""
!
!
%
%
%
%
)
)
"
(& %
('
)$# %

(6.26)$#!!""$!

with#$ ! $
$ #)
%
"!
#"
*
) #"
$)$#

(6.23)(6.27)

where "is an iteration number, # a coefficient which
controls search strategy between the Gauss-Newton and#!unit matrix.
the steepest decent direction and "the !
6.4 Validation
Compared with the inverse modeling under isothermal
assumption proposed by Izumi et al . (2010)[26], theapplicability of this inverse modeling is assessed in termsof reproducibility for observed water movement in test
soil. The time-series observed data is obtained through insitu experiments (sandy soil) in Matsuyama, Ehime
Prefecture.

Figure 6.2: Flow chart for solving IP

6.4.1 Field observation and computational domain
The observation system and computational domain forparameter identification are illustrated in Figure 6.3.
The observation system consists of three sets of
instruments!tensiometer (UIZ-SMT), soil moisture probe(UIZ-SM-2X) and thermometer (TMC20-HD)!and
automatically records time-series data of the pressure
head, volumetric water content and soil temperature atintervals of 10 minutes. To reduce the influence of direct
solar radiation on the observed pressure heads, the sensorsof tensiometers attached above ground are shielded bywhite-colored box with slits for ventilating air. Each set of
instruments is buried at three different depths in a soilcolumn: -10 cm, -20 cm and -30 cm. A pair of pressure
head and volumetric water content data at the same depth
is used to determine SWRC. The observed values of thepressure head at -20 cm deep is used to determine the３０

Tomoki IZUMINovember 25 to December 2, 2010!and used forparameter identification and validation of the inverse
modeling, is shown in Figure 6.4.

function shape of RHC (Eq.(6.10)) and to assess thefitness of estimated RHC function in the simulationoptimization runs. The pressure head and soil temperatureobserved at both top and bottom (-10 cm and -30 cmdeep) of test soil are utilized as Dirichlet boundary values
for solving FP to intentionally avoid measuring flux likeevapotranspiration whose measurement is generally
difficult in the field experiments. The time-series soil
temperatures at -20 cm deep are used as the benchmark
data for confirming reproducibility of forward solution.
The computational domain therefore is the surface soilof 20 cm thick between -10 cm and -30 cm deep, where
the saturation and soil temperature vary significantly withmeteorological conditions, and is divided into four equal
elements with five nodes for numerical computations.
Like in most earlier works on parameter identification,a one-way process of desorption is considered for
estimation of the soil hydraulic properties to exclude thehysteric phenomenon in this study. Hence, the data series
observed during no-rainfall or desorption period!

6.4.2 Identification of SWRC
SWRC is defined as the approximate curve for the
scatter plots relating volumetric water content to pressure
head, and should be determined before estimating RHC.
To obtain "- # relations at a reference soil temperature
expressed by Eq.(6.13), van Genuchten et al . (1991)[60]
proposed the RETC program which is an estimation
method for unknown parameters included in SWRC
function with use of the nonlinear least-squares approach.
Generally, it is very important to input an adequate
estimate of the initial value of unknown parameters forsuch nonlinear fitting. Since the RETC program does nothave the function of automatically estimating the initial
parameters, Seki (2007)[49] developed the alterativealgorithm called SWRC fit. This nonlinear fitting software
not only can automatically determine the initial estimateof unknown parameters, but also can accomplish theparameter identification for the five SWRC models
including VG model. Accordingly, SWRC fit is employed
to identify the reference "- #relations.
For the implementation, observed data of " and # arecorrected for a reference temperature based on the soil
temperature observed at the same depth, using thecoefficient of thermal expansion and Eq.(6.15 ) ,respectively. Additionally, the saturated water content "#is defined as the value of "in case that #is equal to zero
while the residual water content "" as the minimum value
in all observed data of ". Thus, "# and "" are determined
to be 0.30 and 0.05 from all observed data series
including them shown in Figure 6.4, respectively. As thereference soil temperature, the average value of observed
data is adopted which equals to 12.4 degrees Celsius.
Conclusively, a best fitting curve for "- #relation at thereference soil temperature is obtained with the value of!$! and !$! being 28.9 m−1 and 1.43, respectively, as
shown in Figure 6.5.

Figure 6.4: Observed data

Figure 6.5: Identified SWRC

Figure 6.3: Schematic illustration of observation and computation

Inverse Modeling of Variably Saturated Subsurface Water Flow in Isothermal/Non-isothermal Soil

6.4.3 Identification of RHC
Firstly, the number of identified unknown parameters isdetermined. In the free form approach, the large numberof parameters can be set for the functional form of RHC
to have high flexibility or degree of freedom. However,the large number of parameters cause the increasing
uncertainty of identified parameters, making IP difficult tosolve. For this, the identified parameter number should be
kept to the minimum in general. Since the optimal numberof parameters is approximately from seven to nine
according to Bitterlich et al . (2004)[3] and Iden andDurner (2007)[21], the number of unknown parameters is
set to be 10.
Secondly, the manner of division on the definition
domain!effective saturation domain!is determined. The
functional shape of RHC generally has steeper gradientnear the saturated zone. In the context of accuracy, therange with steeper gradient should be divided into finer
subdomain. Therefore the effective saturation range from
0.8 to 1.0 where RHC function rapidly changes ispartitioned into sub-range in the manner of geometricseries.
Eventually, RHC function is identified as shown in
Figure 6.6 after the saturated hydraulic conductivity at the!
$&"$#!% m/s
reference temperature is determined to be %through laboratory experiment.

３
１

Figure 6.7: Reproducibility of forward solution for water
movement

Figure 6.8: Reproducibility of forward solution for heattransport

Figure 6.6: Identified RHC

6.4.4 Reproducibility of calibrated forward simulation model for variably saturated subsurface
flow
The applicability of the inverse modeling presently
developed is assessed in terms of reproducibility for thetime-series observed data by calibrated forward simulation
model through comparing the observed and computed
values. The result of reproducibility for water movement
and heat transport are shown in Figures 6.7 and 6.8,
respectively. In the lower half of Figure 6.7, the absolute"(*) !"*'+#are also shown to demonstrate theerrors #time-varying difference in solution reproducibility. From
the result, it is found that both water movement and heattransport are reproduced with high accuracy.
For the purpose of comparison, the preceding inverse
modeling proposed by Izumi et al . (2010)[26], which did
not consider the influence of soil temperature on water
movement, is applied to the same observed data. In
Figures 6.9 and 6.10, the results for identification of RHC
and the reproducibility of forward solution for water
movement are shown in a similar trend of Figures 6.6 and
6.7, respectively. Hereinafter, the proposed model in thischapter is referred to as the non-isothermal model and themodel proposed by Izumi et al . (2010)[26] is referred to as
the isothermal model. In addition, the values of the!"for both models are summarized
objective function !!
in Table 6.1 to supply superiority of the non-isothermal
model with a quantitative underpinning. The comparison
results indicate that consideration of temperaturedependency is beneficial to simulate the water flow near
the soil surface, and thus it concludes that the inverse
model proposed provide a high performance for themodeling of subsurface water flow.

３
２

Tomoki IZUMI

approach is employed. The inverse problem to identify the
relative hydraulic conductivity function is defined andsolved for the practical field experiments with in-situ soilcolumn near the surface soil to evaluate the performanceof proposed method. From the examination andcomparison of the preceding method which did notconsider the influence of soil temperature on the water
movement, the advantage and appropriateness of thismethod is shown and the importance of consideration for
the temperature-dependency in the water movement near
the soil surface is emphasized.

Figure 6.9: Identified RHC by isothermal model

Figure 6.10: Reproducibility of forward solution for water
movement in isothermal model!"
Table 6.1: Optimal values of !!
Model

!"
!!

Non-isothermal model

1.091

Isothermal model

2.362

6.5 Conclusions
An inverse modeling for variably saturated and nonisothermal subsurface water flow is developed. To
consider the water movement depending on the soil
temperature, a couple of mixed form of Richards equation
with heat conduction equation is employed as thegoverning equation which describes the forward problem.
The soil hydraulic properties which are the model
parameters included in the governing equation are the soil
water retention curve and unsaturated (or relative)
hydraulic conductivity. Since the former is determinedthrough experiments with relative ease, the latter is
unknown parameter in this study. For the functional formof the relative hydraulic conductivity, the free-form

Inverse Modeling of Variably Saturated Subsurface Water Flow in Isothermal/Non-isothermal Soil

CHAPTER 7
SUMMARY AND CONCLUSIONSIn this thesis, inverse modeling of variably saturated
subsurface water flow in isothermal/non-isothermal soil isstudied, considering the following three forms of Richards
equation which describe subsurface water flow.
(1) "-based form Richards equation(2) !-based form Richards equation(3) Mixed form Richards equation
In Chapter 3, an inverse method to identify a function
form of the soil hydraulic properties, i.e., the water
retention curve and unsaturated hydraulic conductivity, invariably saturated subsurface water flow model describedas the "-based form Richards equation is developed. The
soil hydraulic properties are represented by a free form
function composed of sequential piecewise cubic spline
functions. The functions are used as an alternative to theconventional fixed form functions in order to increase
accuracy of forward problems (numerical analyses orsimulations). The inverse problem is formulated into the
optimization problem minimizing the objective function.
The objective function is defined as the total squared error
integrated over space and time between the computed andobserved pressure head, which is the decision variable inthe "-based form Richards equation. It is then solved by asimulation-optimization algorithm with the aid of theLevenberg-Marquardt method. Validity of the method is
examined through twin experiments applying it to two
different soil types which are representative to thesensitivity of the relative hydraulic conductivity nearsaturation. The results show that the soil hydraulicproperties could be successfully identified with a freeform function, and therefore, the approach could be a
viable alternative to the conventional fixed form function
approaches.
In Chapter 4, an inverse method to identify theunsaturated hydraulic conductivity in variably saturated
subsurface water flow model in non-isothermal soil
described as coupled equation of the !-based form
Richards equation and heat conduction equation is
proposed . An in-situ observation with simple
instrumentation set-up collects the hydro-geological datanecessary in solving the inverse problem. For theparameterization of the soil hydraulic properties, a freeform function formed by sequential piecewise cubic spline
functions is used for the unsaturated hydraulicconductivity and van Genuchten model is used for the soil
water retention curve. Using the same manner as that
described in Chapter 3, an inverse problem is defined forthe objective function in terms of the decision variable !
in the !-based form Richards equation which is used todescribe the forward problem for water movement. The

３
３

validity confirmation is carried out by applying themethod to three serial data during desorption in field soil.
From the results, the following issues are concluded:(1) the relative hydraulic conductivity could
successfully identified by a free form functionbe(2) the present approach of reflecting the effect ofthermal conduction on water movement andemploying the free form parameterization techniquecould be a practical alternative to the conventional
approaches, and(3) the specially-contrived and inexpensive observation
system enriches utility of the approach.
In Chapter 5, an inverse method to identify theunsaturated hydraulic conductivity in variably saturated
subsurface water flow model described as the mixed form
Richards equation is proposed. As the basic equation
governing the forward problem, the mixed form Richards
equation is considered which is conservative in terms ofmass balance. The unsaturated hydraulic conductivity, themodel parameter to be identified in the mixed form
Richards equation, is represented by a free form function
composed of sequential piecewise cubic spline functions.
An inverse problem to determine the function form of theunsaturated hydraulic conductivity is defined as the
optimization problem to minimize errors between theobserved and computed pressure heads, and solved withthe use of a simulation-optimization technique. The
developed method is applied to an in-situ soil column andits validity is examined in terms of reproducibility ofdesorption process in the soil. The results indicate that acoupled model of water movement and heat transportshould be considered.
In Chapter 6, an inverse method to identify theunsaturated hydraulic conductivity in variably saturated
subsurface water flow model in non-isothermal soil
described as coupled equation of the mixed form Richards
equation and heat conduction equation is discussed. The
present method is an extension of the method presented in
Chapter 5 to consider water movement in surface soil
depending on soil temperature, and thus, is developed
using the same method, except for the governing equation.
Validity of the method is tested and confirmed throughpractical application to an in-situ soil column and throughcomparison with the previous method.

３
４

Tomoki IZUMIACKNOWLEDGMENTS
The author would like to express his sincerest
appreciation and deepest gratitude to Dr. Toshihiko
KAWACHI, Professor of Water Resources Engineering,
Division of Environmental Science and Technology,
Graduate School of Agriculture, Kyoto University, andChairperson of the Dissertation Committee, for his keen
guidance, valuable suggestions, constructive criticisms,
continuous support and encouragement in the completionof this dissertation. Author’s gratitude to him cannot be
fully acknowledged here.
The author would also like to extend his sincere
appreciation to the members of his Dissertation
Committee: Professor Dr. Akira MURAKAMI andProfessor Dr. Shigeto KAWASHIMA, for their valuable
discussions and advice in reviewing and improving thisthesis.
The author wishes to express special thanks to Dr.
Masayuki FUJIHARA, Professor of Water Resources
Engineering, Department of Rural Engineering, Faculty ofAgriculture, Ehime University, and Dr. Junichiro
TAKEUCHI, Instructor of Water Resources Engineering,
Division of Environmental Science and Technology,
Graduate School of Agriculture, Kyoto University, fortheir warm support, valuable comments and technicalassistance over the research period.
Heartfelt thanks are also due to Associate Professor Dr.
Koichi UNAMI and Lecturer Dr. Shigeya MAEDA of
Water Resources Engineering, Division of Environmental
Science and Technology, Graduate School of Agriculture,
Kyoto University, for their fruitful comments andtechnical assistance.
The author would also like to express his appreciation
to the researchers and professors of the Department of
Rural Engineering, Faculty of Agriculture, Ehime
University, and Associate Professor Dr. Kunihiko
HAMAGAMI of Iwate University, for their hearty
support, cooperation and useful recommendation.
The author also likes to thank all members of Water
Resources Engineering Laboratory of both Ehime andKyoto University, the former and current students and Ms.
Masumi HASHIMOTO, secretary of the Laboratory of
Kyoto University, for their assistance, cooperation andinspiration in all respects of the study period.
The author is grateful to his family for theirwholehearted support and is dedicating this thesis to his
parents, Akihiro and Keiko. Finally, the author also likes
to thank “Aizen-myoo”, the principal Buddhist statue inKongoshin-in Temple, and Kobo Daishi, for theirprotection.

REFERENCES[1] Akai, K., Y. Ohnishi and M. Nishigaki ; Finite element
analysis of saturated-unsaturated seepage in soil ,Proceedings of the Japan Society of Civil Engineers,(264), pp. 87−96, 1973 (in Japanese)[2] Allen, R. G., L. S. Pereira, D. Raes and M. Smith ; Cropevapotranspiration: Guidelines for computing crop waterrequirements, Irrigation and Drainage Paper, (56), FAO,
Rome, 300p 1998
[3] Bitterlich, S., W. Durner, S. C. Iden and P. Knabner ;Inverse estimation of the unsaturated soil hydraulicproperties from column outflow experiments using freeform parameterizations, Vadose Zone Journal , 3(3), pp.
971−981, 2004
[4] Broadbridge, P. and I. White ; Constant rate rainfall
infiltration: a versatile nonlinear model. 1. analytic
solution, Water Resources Research, 24(1), pp. 145−154,
1988[5] Brooks, R. H. and A. T. Corey ; Properties of porous
media affecting fluid flow, Journal of the Irrigation andDrainage Division, Proceedings of the American Society of
Civil Engineers, 92(IR2), pp. 61−88, 1966
[6] Celia, M. A., E. T. Bouloutas and R. L. Zarba ; A generalmass-conservative numerical solution for the unsaturated
flow equation, Water Resources Research, 26(7), pp. 1483
−1496, 1990
[7] Cohen, J. E. ; How Many People Can the Earth Support? ,Norton, 532p, 1995
[8] Constales, D. and J. Ka!
cur ; Determination of soil
parameters via the solution of inverse problems ininfiltration, Computational Geosciences, 5(1), pp. 25−46,
2001[9] de Vries, D. A. ; Simultaneous transfer of heat andmoisture in porous media, Transactions of American
Geophysical Union, 39(5), pp. 909−916, 1958
[10] Durner, W. ; Hydraulic conductivity estimation for soils
with heterogeneous pore structure, Water Resources
Research, 30(2), pp. 211−223, 1994
[11] Eching, S. O. and J. W. Hopmans ; Optimization of
hydraulic functions from transient outflow and soil waterpressure data, Soil Science Society of America Journal ,57(5), pp. 1167−1175, 1993[12] Falkenmark, M., J. Lundqvist and C. Widstrand ; Macroscale water scarcity requires micro-scale approaches,
Natural Resources Forum, 13(4), pp. 258−267, 1989
[13] FOA ; AQUASTAT Online Database, 2011
http://www.fao.org/nr/water/aquastat/main/index.stm[14] Fujinawa, K. ; Coupled analysis of heat transport andwater flow with hysteresis in saturated-unsaturated porous
media, Journal of Groundwater Hydrology, 37(3), pp.
175−192, 1995 (in Japanese)[15] Fujinawa, K. ; Environmental Groundwater Science,
Kyoritsu Publishing, 354p, 2010 (in Japanese)[16] Hillel, D. ; Environmental Soil Physics, Academic Press,
771p, 1998[17] Hills, R. G., I. Porro, D. B. Hudson and P. J. Wierenga ;Modeling one-dimensional infiltration into very dry soils:
1. Model development and evaluation, Water Resources
Research, 25(6), pp. 1259−1269, 1989
[18] Hsu, R. and C.L. Liu ; Use of parameter estimation in
determining soil hydraulic properties from unsaturated

Inverse Modeling of Variably Saturated Subsurface Water Flow in Isothermal/Non-isothermal Soil

transient flow, Water Resources Management, 4(1), pp. 1−19, 1990[19] Huang, K., B. P. Mohanty and M. Th. van Genuchten ; A
new convergence criterion for the modified picard
iteration method to solve the variably saturated flow
equation, Journal of Hydrology, 178(1−4), pp. 69−91,
1996[20] Huyakorn, P. S. and G. F. Pinder ; Computational
Methods in Subsurface Flow, Academic Press, 473p,
1983[21] Iden, S. C. and W. Durner ; Free-form estimation of theunsaturated soil hydraulic properties by inverse modelingusing global optimization, Water Resources Research,
43(7), W07451, 2007[22] Iden, S.C. andW. Durner ; Free-form estimation of soil
hydraulic properties using Wind’s method, European
Journal of Soil Science, 59(6), pp. 1228−1240, 2008[23] Inoue, M., T. Morii, T. Nishimura and H. Fujimaki ; Insitu determination of unsaturated hydraulic conductivity
using profile probe, Transactions of the Japanese Society
Irrigation, Drainage and Rural Enginnering, (231), pp.
39−45, 2004 (in Japanese)[24] IPCC ; Climate Change 2007: Synthesis Report, Contribution ofWorking Groups I, II and III to the Fourth
Assessment Report of the Intergovernmental Panel on
Climate Change, 104p, 2007[25] Ippisch, O., H. J. Vogel and P. Bastian ; Validity limits
for the van Genuchten-Mualem model and implications
for parameter estimation and numerical simulation,
Advances in Water Resources, 29(12), pp. 1780−1789,2006[26] Izumi, T., M. Fujihara, J. Takeuchi and T. Kawachi ;Inverse modeling of massconservative numerical model
for variably saturated seepage flow, Journal of Rainwater
Catchment Systems, 17(1), pp. 11−16, 2011[27] Izumi, T., J. Takeuchi, T. Kawachi and M. Fujihara ; Aninverse method to estimate unsaturated hydraulicconductivity in seepage flow in non-isothermal soil,
Transactions of the Japanese Society Irrigation, Drainage
and Rural Enginnering, (264), pp. 35−42, 2009[28] Izumi, T., J. Takeuchi, T. Kawachi, K. Unami and S.
Maeda ; An inverse method to estimate soil hydraulicproperties in saturated-unsaturated groundwater flow
model, Journal of Rainwater Catchment Systems, 13(2),pp. 23−28, 2008[29] Japanese Association of Groundwater Hydrology ;Simulation of Groundwater Flow and Solute Transport,
Gihodo Publishing, 232p, 2010 (in Japanese)[30] Jury, W.A. and R. Horton ; Soil Physics, John Wiley &Sons, 370p, 2004[31] Kastanek, F. J. and D. R. Nielsen ; Description of soil
water characteristics using cubic spline interpolation, Soil
Science Society of America Journal , 65(2), pp. 279−283,
2001[32] Klute, A. ; A numerical method for solving the flow
equation for water in unsaturated materials, Soil Science,73(2), pp. 105−116, 1952[33] Kondo, J. and N. Saigusa ; Modeling the evaporationfrom bare soil with a formula for vaporization in the soil
pores, Journal of Meteorological Society of Japan, 72,
pp. 413−421, 1994 (in Japanese)

３５[34] Kosugi, K. ; Three-parameter lognormal distribution
model for soil water retention, Water Resources
Research, 30(4), pp. 879−890, 1994[35] Kosugi, K. ; Lognormal distribution model forunsaturated soil hydraulic properties, Water Resources
Research, 32(9), pp. 2697−2703, 1996[36] Milly, P. C. D. ; Moisture and heat transport in hysteretic,
inhomogeneous porous media: A matric head-based
formulation and a numerical model, Water Resources
Research, 18(3), pp. 489−498, 1982[37] Milly, P. C. D. ; A similation analysis of thermal effects
on evaporation from soil, Water Resources Research,
20(8), pp. 1087−1098, 1984[38] Milly, P. C. D. ; A mass-conservative procedure for timestepping in models of unsaturated flow, Advances in
Water Resources, 8(3), pp. 32−36, 1985
[39] Ministry of Land, Infrastructure, Transport and Tourism ;Water Resources in Japan, 209p, 2011 (in Japanese)[40] Ministry of the Environment ; 2010 Annual Report on theEnvironment, the Sound Material-cycle Society and the
Biodiversity in Japan : Our Responsibility andCommitment to Preserve the Earth −Challenge 25−,
169p, 2011[41] Moukalled, F. and Y. Saleh ; Heat and mass transport inmoist soil, part I. formulation and testing, Numerical Heat
Transfer, Part B , 49, pp. 467−486, 2006[42] Mualem, Y. ; New model for predicting the hydraulicconductivity of unsaturated porous media, Water
Resources Research, 12(3), pp. 513−522, 1976[43] Neuman, S. P. ; Saturated-unsaturated seepage by finite
elements , Journal of the Hydraulics Division ,Proceedings of the American Society of Civil Engineers,
99(HY12), pp. 2233−2250, 1973[44] Neuman, S. P. and P. A. Witherspoon ; Finite element
method of analyzing steady seepage with a free surface,
Water Resources Research, 6(3), pp. 889−897, 1970[45] Neuman, S. P. and P. A. Witherspoon ; Analysis of
nonsteady flow with a free surface using the finite
element method, Water Resources Research, 7(3), pp.
611−623, 1971[46] Philip, J. R. and D. A. de Vries ; Moisture movement inporous materials under temperature gradients, Transactions of American Geophysical Union, 38(2), pp. 222−
232, 1957[47] Richards, L. A. ; Capillary conduction of liquids throughporous mediums, Journal of Applied Physics, 1(5), pp.
318−333, 1931
[48] Scanlon, B. R. and P. C. D. Milly ; Water and heat fluxes
in desert soils: 2. Numerical simulations, Water Resources
Research, 30(3), pp. 721−733, 1994
[49] Seki, K. ; SWRC fit −A nonlinear fitting program with awater retention curve for soils having unimodal andbimodal pore structure, Hydrology and Earth System
Sciences Discussions, 4(1), pp. 407−437, 2007[50] Si, B.C. and R.G. Kachanoski ; Estimating soil hydraulicproperties during constant flux infiltration: Inverse
procedures, Soil Science Society of America Journal ,64(2), pp. 439−449, 2000
[51] Sun, N. Z. ; Inverse Problems in Groundwater Modeling,
Kluwer Academic Publishers, 337p, 1994
[52] Takeshita, Y. and I. Kohno ; A method to pedict

３
６

Tomoki IZUMI

hydraulic properties for unsaturated soils and its
application to observed data, Ground Engineering;
Journal of Chugoku Branch, JGS , 11(1), pp. 95−113,
1993 (in Japanese)[53] Takeuchi, J., T. Izumi and T. Kawachi ; An inverse
method for unsaturated hydraulic conductivity based on
in-situ observation, Proceedings of the 15th Annual
Congress of Japan Rainwater Catchment Systems
Association, 15, pp. 5−8, 2007 (in Japanese)[54] Takeuchi, J., T. Kawachi and T. Izumi ; Explicit scheme
for mixed form of Richards equation, Proceedings of theJapanese Society of Irrigation, Drainage and RuralEngineering, pp. 322−323, 2008 (in Japanese)[55] United Nations Development Programme (UNDP) ;Human Development Report 2006: Beyond scarcity;
power, poverty and the global water crisis, 422p, 2006
[56] United Nations Population Division ; World Population
Prospects: the 2010 Revision, 2010
http://esa.un.org/unpd/wpp/index.htm[57] !
Simunek, J., O. Wendroth and M. Th. van Genuchten ;Parameter estimation analysis of the evaporation methodfor determining soil hydraulic properties, Soil Science
Society of America Journal , 62(4), pp. 894−905, 1998
[58] van Dam, J. C., J. N. M. Stricker and P. Droogers ;Inverse method to determine soil hydraulic functions frommultistep outflow experiments, Soil Science Society ofAmerica Journal , 58(3), pp. 647−652, 1994
[59] van Genuchten, M. Th. ; A closed-form equation forpredicting the hydraulic conductivity of unsaturated soils,
Soil Science Society of America Journal , 44(5), pp. 892−898, 1980
[60] van Genuchten, M. Th., F. J. Leij and S. R. Yates ; The
RETC code for quantifying the hydraulic functions of
unsaturated soils, Technical Report EPA/600/2-91/065,
US Environmental Protection Agency, 85p, 1991
[61] Vogel, T., M. Th. van Genuchten and M. Cislerova ;Effect of the shape of the soil hydraulic functions nearsaturation on variably-saturated flow predictions ,Advances in Water Resources, 24(2), pp. 133−144, 2001[62] Zadeh, K. S. ; A mass-conservative switching algorithm
for modeling fluid flow in variably saturated porous
media, Journal of Computational Physics, 230(3), pp.
664−679, 2011[63] Zhang, K., I. G. Burns, D. J. Greenwood, J. P. Hammondand P. J. White ; Developing a reliable strategy to infer
the effective soil hydraulic properties from fieldevaporation experiments for agro-hydrological models,
AgriculturalWater Management, 97(3), pp. 399−409,
2010[64] Zienkiewicz, O. C., P. Mayer and Y. K. Cheung ;Solution of anisotropic seepage by finite elements,
Journal of the Engineering Mechanics Division,
Proceedings of the American Society of Civil Engineers,
92(EM1), pp. 111−120, 1966!