J. Hydrol. Hydromech., 62, 2014, 1, 7–15
DOI: 10.2478/johh-2014-0008

Parameter estimation of soil hydraulic and thermal property functions forunsaturated porous media using the HYDRUS-2D code
Mohammad Nakhaei1, Jiří Šimůnek2*
1

Faculty of Geological Sciences, Kharazmi University, No.49, Mofateh Ave., Tehran, Iran. Tel.: +1-951-750-3157,
E-mail: nakhaei@khu.ac.ir
2
Department of Environmental Sciences, University of California, Riverside, USA.
*
Corresponding author. Tel.: +1-951-827-7854. E-mail: Jiri.Simunek@ucr.edu

Abstract: Knowledge of soil hydraulic and thermal properties is essential for studies involving the combined effects ofsoil temperature and water input on water flow and redistribution processes under field conditions. The objective of thisstudy was to estimate the parameters characterizing these properties from a transient water flow and heat transport field
experiment. Real-time sensors built by the authors were used to monitor soil temperatures at depths of 40, 80, 120, and160 cm during a 10-hour long ring infiltration experiment. Water temperatures and cumulative infiltration from a single
infiltration ring were monitored simultaneously. The soil hydraulic parameters (the saturated water content θ s, empiricalshape parameters α and n, and the saturated hydraulic conductivity Ks) and soil thermal conductivity parameters (coefficients b1, b2, and b3 in the thermal conductivity function) were estimated from cumulative infiltration and temperature
measurements by inversely solving a two-dimensional water flow and heat transport using HYDRUS-2D. Three scenarios with a different, sequentially decreasing number of optimized parameters were considered. In scenario 1, seven parameters (θ s, Ks, α, n, b1, b2, and b3) were included in the inverse problem. The results indicated that this scenario does
not provide a unique solution. In scenario 2, six parameters (Ks, α, n, b1, b2, and b3) were included in the inverse problem.
The results showed that this scenario also results in a non-unique solution. Only scenario 3, in which five parameters (α,n, b1, b2, and b3) were included in the inverse problem, provided a unique solution. The simulated soil temperatures andcumulative infiltration during the ring infiltration experiment compared reasonably well with their corresponding observed values.
Keywords: Infiltration; Heat transport; Parameter Estimation; Richards’ equation; HYDRUS-2D.
INTRODUCTION
The importance of the unsaturated zone as an integral part of
the hydrological cycle has long been recognized. The vadosezone plays an important role in many aspects of hydrology,
including infiltration, soil moisture storage, evaporation, plantwater uptake, groundwater recharge, runoff, and erosion. Simultaneous movement of water and heat in the vadose zone ofarid and/or semi-arid regions is of great interest in evaluatingwater and energy balance of subsurface environments in both
agricultural and engineering applications (e.g. Saito et al.,2006). Although it is widely recognized that the movement ofwater and heat is closely coupled and affect each other (e.g.
Nassar and Horton, 1992; Sakai et al., 2009), especially underarid conditions, their mutual interactions are rarely consideredin practical applications.
Unsaturated water flow through a porous medium is generally described using the Richards equation. The solution of thisequation requires knowledge of soil hydraulic functions, that is,
the retention curve, θ(h), and the hydraulic conductivity function, K(h). Currently there are many laboratory and field methods to determine these nonlinear functions (e.g. Dirksen, 1991;
Klute, 1986). The methods with minimum soil disturbance,
lowest time requirements, and lowest cost are preferred. Most
direct methods require restrictive initial and boundary conditions, which make measurements time consuming and expensive. Direct methods for determining the water retention curve
involve equilibrating soil at certain pressure heads or simultaneous measurements of water contents and pressure heads
during infiltration or drainage (e.g. Dane and Topp, 2002;
Dirksen, 2000).

When the inverse modeling approach is used, the unknown
hydraulic parameters are estimated by minimizing deviations
between observed variables and model-predicted output fortransient flow experiments (e.g. Kool and Parker, 1987; Kool etal., 1987). Collected transient flow data can be inversely analyzed using either analytical or numerical models (e.g. Šimůnek
and van Genuchten, 1996; Šimůnek et al., 1998). For the inverse problem solution, we encounter the question: under which
circumstances is the inverse problem well posed, utilizing theconcepts of identifiability, uniqueness, and stability? It has beenshown in the literature, that it may be possible to decrease nonuniqueness of optimized parameters by simultaneously considering multiple processes (e.g. Hopmans et al., 2002; Inoue etal., 2000; Mortensen et al., 2006; Šimůnek et al., 2002). Inoueet al. (2000) and Šimůnek et al. (2002) showed that by simultaneously analyzing data involving both water flow and solutetransport variables, they could significantly decrease the uncertainty in the optimized parameters. Similar conclusions weredrawn by Hopmans et al. (2002) and Mortensen et al. (2006)
who simultaneously analyzed data involving both water flow
and heat transport variables.
Analysis of water infiltration experimental data has becomea widely used practice for obtaining soil hydraulic properties.
The infiltration rate and its variation with time depends upon
the initial water content, as well as on the texture, structure, anduniformity (or layering sequence) of the soil profile (Reynolds
et al., 2002a, b). The single-ring infiltrometer is often used formeasuring cumulative infiltration, from which the infiltration
rate, and field-saturated hydraulic conductivity (Ks) are thenderived. The infiltration rate is initially large and decreases withtime until a near steady-state (gravity-induced) value is at-7

Mohammad Nakhaei, Jirka Šimůnek

tained. Factors that affect the infiltration rate have been divided
into categories of soil factors, surface properties, management
factors, and natural properties. Soil factors encompass soil
particle size, morphological, chemical, and soil water properties. Surface factors are those that affect the movement of water
through the air–soil interface such as soil crust. These factors
should be taken into account when parameterizing the infiltration model.
Considerable attention was given in the past to the estimation of soil hydraulic properties from ponded or tension infiltration experiments (e.g. Bohne et al., 1992; Russo et al., 1991;Šimůnek and van Genuchten, 1996). Russo et al. (1991) andŠimůnek and van Genuchten (1996) reported that the cumulative infiltration curve of ponded and tension infiltration experiments, respectively, does not provide enough information to
obtain a unique set of optimized parameters, and that some
additional information is needed in order to obtain unique solutions for a large number of parameters. Such additional information could include water contents and/or pressure heads (e.g.
Šimůnek and van Genuchten, 1996), solute concentrations (e.g.
Inoue et al., 2000), or temperatures (e.g. Mortensen et al., 2006)
measured at a certain depth under the soil surface.
In the past, relatively few studies have systematically evaluated the simultaneous estimation of soil hydraulic and heattransfer parameters, and therefore the effect of heat transport onwater flow has often been neglected. The coupled movement ofwater and heat in the subsurface is implemented in the HYDRUS codes (Šimůnek et al., 2006, 2008). The HYDRUS-1Dand HYDRUS-2D programs numerically solve the Richards
equation for saturated-unsaturated water flow and the convection-conduction equation for heat transport in one and two
dimensions, respectively. The heat transport equation involves
terms accounting for the transport due to conduction and convection with flowing water. The program may be used to analyze water and heat movement in unsaturated, partiallysaturated, or fully-saturated porous media. Examples in whichsoil thermal parameters have been estimated using HYDRUS
include studies by Hopmans et al. (2002), Mortensen et al.
(2006), Saito et al. (2007), and Sakai et al. (2009).
While it was shown by Russo et al. (1991) and several others(e.g. Šimůnek and van Genuchten, 1996) that cumulative infiltration data alone cannot be used to estimate more than two soil
hydraulic parameters due to problems with parameter uniqueness, the combined use of water flow and heat transport data inthe objective function of the inverse problem may significantly
reduce the uncertainty of the optimized parameters and enhancetheir identifiability. The first objective of this study was thus tocollect data during the ring infiltration experiment and to simultaneously measure cumulative infiltration and soil temperaturesin several depths under the center of the infiltration ring. The
second objective was to analyze the water flow and heattransport data using the parameter estimation approach in HYDRUS-2D to simultaneously estimate soil hydraulic and heattransport parameters.
MATERIAL AND METHODS
Study area
The study area is located on the campus of Kharazmi University in the city of Karaj, Alborz province, 48 kilometers
outside Tehran, Iran. The altitude is 1,297 m above sea level,
and the annual average precipitation is approximately 300 mm.
The annual average temperature is about 16°C. The geology isdominated by the Karaj tuff formations. Soils developed from8

this bedrock are sandy loams with a low content of clay. The
vegetation cover is sparse.

Fig. 1. Photo of the experiment and locations of thermal sensors.

Soil sampling
To determine the saturated hydraulic conductivity, Ks, aswell as the soil bulk density, ρ, the initial water content, and thegrain-size distribution curve (soil texture), two undisturbed soil
samples were collected using an auger from a depth of 1.0 m,
about 50 cm away from the center of the infiltration ring. The
soil samples were analyzed in the sedimentology laboratory of
Kharazmi University. Using the sieve analysis, the soil texture
was determined to be a silt loam (based on the USDA textural
triangle; e.g. Radcliffe and Šimůnek, 2010) with a low content
of clay (40% sand, 55% silt, and 5% clay). The soil could also
be classified as "Sandy Silt" according to the textural classification scheme for siliciclastics soil, modified from Shepard(1954), which subdivides the textural triangle into sandy clay,
silty clay, sandy silt, clayey silt, clayey sand, and silty sandclasses. The measured soil bulk density and the saturated watercontent were 1.24 g/cm3 and 0.43, respectively. Bulk density,
soil water content, and soil texture were measured using standard methods described by Klute (1986).
Single ring infiltration experiment
We used the following procedure to conduct the single-ring
infiltration experiment in the field. First, we cleared the sam-Parameter estimation of soil hydraulic and thermal property functions for unsaturated porous media using the HYDRUS-2D codepling area of surface residues that could impede insertion of thering. Second, using a sledge hammer and block of wood, we
inserted the single 40-cm diameter infiltration ring (Fig. 1)
down to a depth of 3 cm, keeping the ring sides vertical and thesoil disturbance to a minimum. Third, we consolidated the soilon the inside of the ring and sealed the space outside of the ringusing bentonite to prevent unwanted seepage. Finally, we instrumented the soil profile with 4 thermal sensors at depths of
40, 80, 120, and 160 cm (Fig. 2). Thermal sensors were installed from a trench wall, about 50 cm aside from the center ofthe ring, by horizontally drilling small diameter holes below thecenter of the ring. Drilled holes were filled with a wet fine silty
soil to guarantee a complete contact between soil and sensors.
The performance of the sensors was successfully comparedagainst the HACH SensION-156 Portable Multiparameter, a 5pin (pH/mV/temperature/conductivity/dissolved oxygen) probe.
The measurement period was from 8.00 to 18.00 hours on April
4th, 2012. Air temperature was 20°C and initial soil temperatures at observation depths were 17.5°C. The temperature of theinfiltration water was 61°C during the entire experiment. In our
calculations we assumed that the effect of air temperature variations during the experiment was negligible. The final steadystate infiltration rate was also measured and subsequently used
for calculating the saturated hydraulic conductivity of the soil
using the Wooding (1968) equation (cited from Radcliffe andŠimůnek, 2010, p. 106).

Soil hydraulic propertiesHYDRUS implements the soil hydraulic functions of van
Genuchten (1980), who used the statistical pore-size distribution model of Mualem (1976) to obtain a predictive equation
for the unsaturated hydraulic conductivity function in terms ofsoil water retention parameters. The expressions of van
Genuchten (1980) are given by Eq. (1) and Eq. (2):1  r n m
S e ( h)   [1 | h| ]s - r 1

2

K (h )  K s Sel [1- (1- Se1/ m ) m]

h0

,

(1)

,

(2)

h0

in which θ is the volumetric water content [-],θr and θs denote
the residual and saturated water content [-], respectively; Se is
the effective fluid saturation [-], Ks is the saturated hydraulicconductivity [LT–1], α [L–1], n [-], and m (=1-1/n) [-] are empirical shape parameters, and l is a pore-connectivity parameterassumed to be 0.5 [-] (Mualem, 1976).
Temperature dependence of the soil hydraulic functions
A scaling technique similar to one used to describe spatial
variability of soil hydraulic properties (Vogel et al., 1991) isused in HYDRUS to express the temperature dependence of thesoil hydraulic functions. Based on the capillary theory thatassumes that the influence of temperature on the soil waterpressure head can be quantitatively predicted from the influenceof temperature on the surface tension, Philip and de Vries(1957) derived the following equation:

dh h d ,=dT  dT(3)

where T is temperature [K] and σ is the surface tension at the
air-water interface [MT–2]. From (3) it follows that:

hT =

T
href =  h* href ,
 ref

(4)

where hT and href (σT and σref) are pressure heads (surface tensions) at temperature T and reference temperature Tref, respectively, and αh* is the temperature scaling factor for the pressure
head. Following Constantz (1982), the temperature dependenceof the hydraulic conductivity can be expressed as follows:

KT ( ) 
Fig. 2. Schematic of a single infiltration ring experiment. A - planeview, B - Frontal view, C-Lateral view.

The numerical model
The HYDRUS-2D software package (version 1.12 of HYDRUS (2D/3D); Šimůnek et al., 2006) was used to model thesoil water dynamics and heat transport assuming an axisymmetric transport domain.

 ref  T
Kref ( )   K* Kref ( ) , T  ref(5)where Kref and K denote hydraulic conductivities at the reference temperature ref and soil temperature , respectively, refand  (ref and ) represent the dynamic viscosity [ML–1T–1](density of soil water [ML–3]) at temperatures ref and , respectively, and K* is the temperature scaling factor for thehydraulic conductivity.

9

Mohammad Nakhaei, Jirka Šimůnek

Heat transport

Initial and boundary conditions

Neglecting the effects of water vapor diffusion, two- andthree-dimensional heat transport can be described as (Sophocleous, 1979):

The initial water content considered in the simulation (0.041
cm3/cm3) was measured in the field. A soil profile temperatureof 17.5°C, as measured in the field, was used as the initial condition for heat transport. All sides of the flow region were considered to be impervious, except for a small portion around theorigin at the ponded surface inside the ring infiltrometer wherea constant water content was imposed, and for the bottom of thesoil profile where a free drainage condition was imposed. The
volumetric water content of 0.43 was used as the water flowboundary condition under the single infiltration ring. A firsttype heat transport boundary condition was used for nodes
under the single infiltration ring with a constant temperature of
61°C, representing hot water infiltrating into the soil profile. A
third-type heat transport boundary condition was invoked at thebottom of the soil profile. No flux boundary conditions wereused for both water flow and heat transport along all otherboundary nodes.

C ( )

T
 
T 
T
,
=
 ij ( ) - Cw qi
t xi x j 
xi(6)where ij is the apparent thermal conductivity of the soil
[MLT–3K–1] (e.g. Wm–1K–1), and C and Cw are the volumetric heat capacities [ML–1T–2K–1] (e.g. Jm–3K–1) of the porous
medium and the liquid phase, respectively. Volumetric heat
capacity is defined as the product of the bulk density and thegravimetric heat capacity. The first term on the right-hand sideof (6) represents heat flow due to conduction, and the second
term accounts for heat being transported by flowing water. We
do not consider the transfer of latent heat by vapor movement.
According to de Vries (1963), the volumetric heat capacity canbe expressed as:

C ( )  Cn n  Co o  Cw w  Cg  g ,(7)

where  refers to a volumetric fraction [L3L–3], and subscripts n,o, g, and w represent solid phase, organic matter, gas phase andliquid phase, respectively.
Apparent thermal conductivity coefficient
The apparent thermal conductivity, ij combines thethermal conductivity  of the porous medium (solid plus
water) in the absence of flow, and the macro-dispersivity,
which is assumed to be a linear function of the velocity (deMarsily, 1986). In analogy with the dispersion coefficient forsolute transport, the apparent thermal conductivity ij isgiven by (Šimůnek and Suarez, 1993):

ij ( )  T Cw | q |  ij + (L - T )Cw

q j qi|q|+ 0 ( ) ij ,(8)where q is the absolute value of the Darcian fluid flux density[LT–1], ij is the Kronecker delta function, and L and  are thelongitudinal and transverse thermal dispersivities [L], respectively. The volumetric heat capacity of the liquid phase Cw is
included in the definition of the thermal conductivity in order tohave the dimensions of the thermal dispersivities in length units(de Marsily, 1986). The thermal conductivity,  accounts
for the tortuosity of the porous medium, and is described withthe simple equation (Chung and Horton, 1987):

0 ( ) = b1 +b2 +b3 0.5 ,(9)where b1, b2 and b3 are empirical parameters [MLT–3K–1] (e.g.
Wm–1K–1).
Numerical modelTransport domain
The axisymmetrical transport domain was 100 cm wide and200 cm deep. The domain was discretized into 441 nodes and
400 triangular finite elements. Simulations were carried out
over a period of 10 hours.

10

Inverse solution
The Levenberg-Marquardt optimization algorithm (Marquardt, 1963), in combination with the HYDRUS-2D numericalcode (Šimůnek et al., 2006), was used to inversely estimate
unsaturated soil hydraulic and heat transport parameters fromthe transient cumulative water flux from a single-ring infiltration experiment and soil temperatures measured at multiple
depths. An inverse problem requires the same type of information as a forward problem. It also requires initial estimates of
optimized parameters, positions of observation points, andmeasured data with corresponding measurement times andweighing factors.
A summary of various approaches that are routinely used ininverse modeling is provided by Hopmans and Šimůnek (1999)
and Vrugt et al. (2008). Optimized hydraulic and thermal parameters are determined by systematically minimizing thedifferences between observed and simulated state variables (e.g.
cumulative infiltration and temperatures). The sum of squaresof these differences is expressed using an objective function, ,which may be defined as:j myi n jj 1i 1 (  , y ) =  v j  wi , j  y *j  z , ti   y j  z , ti ,    ,2

(10)

where the right-hand side represents the residuals between
space–time variables that are measured (yj*) and modelpredicted (yj) using the soil hydraulic and thermal parameters of
the optimized parameter vector, . The first summation sums
the residual for all measurement types my (e.g. cumulative
infiltration and temperatures), whereas the variable nj in thesecond summation denotes the number of measurements for acertain measurement type j. In this paper, yj* represents cumulative infiltration from a single ring and temperature measurements at 4 different depths of the soil profile. Assuming that themeasurement errors within a measurement type are independentand uncorrelated, the weighting factor values for vj can be selected such that data types are weighted equally using the reciprocal of the measurement variance of measurement type j(Clausnitzer and Hopmans, 1995).
The Levenberg–Marquardt method, which is a gradientbased, nonlinear optimization technique, was used to minimize
the objective function (10). However, inverse problems areoften ill-posed, which results in optimized parameters being

Parameter estimation of soil hydraulic and thermal property functions for unsaturated porous media using the HYDRUS-2D code

unidentifiable and/or non-unique, producing divergent, unstable
results (Hopmans and Šimůnek, 1999). Because it is not apriori clear whether the inverse problem is sufficiently wellposed to allow parameter optimization when using cumulative
infiltration from a single-ring infiltration experiment and soil
temperatures measured at multiple depths, the parameter optimization procedure was repeated many times with different
initial estimates of the optimized parameters. The initial estimates of flow parameters (soil hydraulic parameters) weretaken from the ROSETTA program (Schaap et al., 2001).
The behavior of the inverse problem can be evaluated byplotting the value of the objective function, , against a pair of
optimized parameters to obtain the so-called response surfaces.
Each response surface is obtained by solving the flow andtransport equations, with appropriate initial and boundary conditions, for many possible combinations of a selected pair ofparameter values within a predetermined range, while keeping
the other parameters constant. Since only two parameters appear within a single response surface analysis, the behavior of in the various parameter planes is only an indication of theuniqueness of the solution in these planes and not in the fullparameter space. For example, other local minima may exist
elsewhere in the full parameter space but not appear in theselected cross-sectional planes (Šimůnek and van Genuchten,
1996).
The validity of estimated parameters can be assessed by anuncertainty analysis, which yields confidence intervals for theestimated parameters. This analysis assumes that uncertainty iscaused by measurement errors only, that is, that the model error
is zero, and that the inverse solution has converged to the global
minimum. The uncertainty analysis provides a means to compare confidence intervals among parameters; thereby indicating
which parameters should be measured or estimated independently.
We carried out three different scenarios with a different
number of optimized parameters to test the uniqueness of theinverse solution. In Scenario 1, we considered four parametersof the van Genuchten model s, n, and Ks) and three parameters of the thermal conductivity function (b1, b2, and b3) as
unknown parameters in the inverse solution. In Scenario 2, weconsidered three parameters of the van Genuchten model (, n,
and Ks) and three thermal conductivity parameters (b1, b2, andb3) as unknown parameters. Finally, in Scenario 3, we considered only the van Genuchten–Mualem parameters  and n, andthe empirical parameters b1, b2, and b3 of the thermal conductivity function as unknowns. The saturated water content, s, andthe saturated hydraulic conductivity, Ks, determined using
independent measurements, were considered as known parameters in both Scenario 2 and 3. The tortuosity parameter l was set
in all scenarios to a fixed value of 0.5 (Mualem, 1976). The
residual water content s was fixed to a value of 0.001.

Table 1 shows the values of soil hydraulic and thermal parameters that were used in all simulations when they were not
optimized. Table 2 then shows the initial estimate and final
values of the optimized soil hydraulic and thermal parameters,
R2, and SSQ for different optimization runs (different scenario).
In Scenario 1, since the optimized parameters were often very
different in different optimization runs with different initialvalues of optimized parameters, the inverse problem was clearly not well posed. In Scenario 2, although the inverse solutions
showed less variation in optimized parameters than in Scenario
1, the inverse solution was still non-unique. Only in Scenario 3,
the inverse solution for many different runs converged to thesame solution, indicating a unique solution. Table 2 lists R2
values for the best optimization runs that quantify the accuracyof regression. As all correlation values are very high, theirmagnitude is of little value for testing the model’s adequacy.

Fig. 3. Final simulated water contents (left) and temperatures(right) in the soil profile for Scenario 3.

RESULTS AND DISCUSSION
The steady-state infiltration rate from the ring was 0.00166
cm/s. Using the macroscopic capillary length of 8.3 cm for the
loamy soil (Radcliffe and Šimůnek, 2010), we calculated thesaturated hydraulic conductivity of 0.0011cm/s using Wooding's analytical solution (Wooding, 1968) and used this value as
an initial estimate in the inverse solution. Measured cumulative
infiltration and soil temperatures measured at four observation
depths were then used to optimize selected soil hydraulic andheat transport parameters of a silt loam soil profile.

Fig. 4. Simulated and observed cumulative fluxes (top) and temperatures at four depths (bottom) (blue - 40 cm, green - 80 cm, red- 120 cm, beige - 160 cm).

11

Mohammad Nakhaei, Jirka Šimůnek

Fig. 3 shows an example of the output simulated by the
HYDRUS-2D model for Scenario 3. It shows the final spatial
distribution of water contents and temperatures in the soil profile. Fig. 3 shows that at the end of the simulation the moisture
front reaches the horizontal distance of about 50 cm from thecentre of the ring, which is the distance of the trench fromwhich temperature sensors were installed. Notice that the moisture front is also visible in the trench in Fig. 1. This indicatesthat the trench starts affecting the flow field only at the very
end of the experiment and that the use of axisymmetrical domain in the analysis is acceptable. Fig. 4 compares simulated
and measured cumulative infiltrations and simulated and measured temperatures in the soil profile at depths of 40, 80, 120,
and 160 cm. Both temperatures and cumulative infiltrations
optimized with the model fitted the measured data well. Also,
the coefficient of determination (R2) (Table 2) indicated a goodfit. R2 for regression of predicted versus observed values isequal to 0.995.

Table 3 shows the upper and lower bounds (95% confidence
interval) of the optimized parameters for different optimization
scenarios and the final optimized soil hydraulic and heattransport parameters obtained in Scenario 3. Confidence intervals for optimized model parameters for a given level of significance were calculated from the variances of the estimatedmodel parameters. As shown in Table 2, the optimized soil
hydraulic and heat transport parameters have narrow confidence intervals only when the saturated water content and thesaturated hydraulic conductivity are assumed to be known(from independent measurements), and only the  and n parameters of the van Genuchten model are considered to be theunknown soil hydraulic parameters in the inverse problem.

Table 1. Values of soil hydraulic and thermal parameters that were
not optimized.
Parameters (unit)r(cm /cm )
33

s (cm /cm )
3

3

Parameter namevalue

Residual moisture content

0.001

Saturated moisture content

0.430

Saturated hydraulic conductivity

0.0011

Pore connectivity factor

0.500

n [-]

Solid fraction

0.600

o [-]

Organic matter fraction

0.001

L [cm]

Thermal dispersivity longitudinal

2.000

Ks (cm/s)l (-)

T [cm]

Thermal dispersivity transverse

0.200

Cn [J/M3/K]

Heat capacity of the solid phase

1920000

Co [J/M3/K]

Heat capacity of the organic matter

2510000

Heat capacity of water

4180000

3

Cw [J/M /K]

Fig. 5 shows a comparison of the optimized thermal conductivity function  for our silt loam soil with those for thethree main textural classes of clay, loam, and sand (Chung andHorton, 1987). The results indicate that the optimized thermal
conductivity function is relative similar to (between) those forloam and clay.

Fig. 5. Comparison of the optimized thermal conductivity function
() (denoted as "soil") with those for main textural classes clay,
loam, and sand (Chung and Horton, 1987).

12

Fig. 6. Uncertainty in optimized (h) (top) and log K(h) (cm/s)(bottom) functions in Scenario 3. Dotted lines represent the optimized functions and solid lines the 95% confidence intervals.

Fig. 7. A response surface in the -n parameter plane for Scenario
3. The white square indicates the minimum of the objective function  for a final set of parameters (Table 3).

Parameter estimation of soil hydraulic and thermal property functions for unsaturated porous media using the HYDRUS-2D code
Table 2. Optimized soil hydraulic and thermal parameters, R2, and SSQ for three scenarios (for optimization runs with the lowest value of
SSQ).
Heat transport parameters
Optimized value(kg cm s–3 K–1)

Water flow parameters:
Optimized values

R2

Objectivefunction

Scenario

Water flow parameters:
Initial estimates

s

1

0.43

0.035

2

#

0.075

2.0

0.00132

#

0.075

1.51

0.00122

88.3

10.1

47.6

0.995

0.016

3

#

0.048

1.60

#

#

0.075

1.51

#

43.3

29.6

106.0

0.998

0.007


(1/cm)n

Ks (cm/s)

s

1.5

0.00106

0.377



n

Ks
(cm/s)

b1

b2

b3

-

SSQ*

0.065

1.52

0.00118

24.82

67.64

96.76

0.995

0.014

(1/cm)*SSQ: Value of the objective function Φ being minimized during the parameter optimization process.
R2: Regression between observed and fitted values.
#: Parameter is not included in the parameter optimization process.
Table 3. The upper and lower bounds (95% confidence interval) of the optimized parameters for three scenarios. The final column gives
final parameter values for Scenario 3.
Parameters

Scenario 1
Lower
Upper
bound
bound

Scenario 2

Scenario 3

Lower bound

Upper
bound

Lower
bound

Upper boundFinal parameter values for
Scenario 3

#

#

0.43

s(cm3/cm3)

0.216

0.501

#

#

Ks (cm/s)

0.000230.010390.001190.00137##0.00132(1/cm)

0.009

0.102

0.066

0.084

0.072

0.078

0.075

n (-)

1.160

2.570

1.4

1.62

1.5

1.6

1.51

b1 (kg cm s–3 K–1)

0.071

78.90

28.7

122.1

30.3

70.1

43.3

–1

b2 (kg cm s K )0.0364

578.2

0.24

62

20.7

57.3

29.6

b3 (kg cm s–3 K–1)

1.295

425.7

10.6

147.1

78.5

140.3

106.0

–3

#: not included in the parameter optimization.

Uncertainties introduced into the optimized soil hydraulicfunctions by measurement errors, and other factors discussed
above, are shown in Fig. 6. When selecting the best estimates of
optimized parameters, we considered their 95% confidence
intervals, the magnitude of objective function and the R2. By
adjusting the weights of data points related to flow and temperature variables in the objective function, we attempted to givethe same weights to both measurement sets. Although not defining weights strictly based on measurement errors makes it
difficult to interpret the resulting optimized parameters, theirconfidence intervals, correlations, and, in general, their relationship with the true parameter values (Bard, 1974), this is acommon approach in many applications since the measurement
errors are often not known and thus the corresponding weights
are then not specified based on probabilistic assumptions.
However, weights still need to be selected with caution since an
improper selection of weights can influence not only the confidence regions of optimized parameters, but also the location of
the minimum of the objective function (Hollenbeck et al, 2000).
Even the robustness of the least-squares criterion for the estimation of model parameters has recently been questioned by
Finsterle and Najita (1998). They pointed out that the leastsquare criterion causes outliers to strongly influence the final
values of optimized parameters. Hence, outliers (e.g. individual
data points with large measurement errors, as is often the casewith field measurements) can introduce a significant bias in theestimated model parameters. Finsterle and Najita (1998) studied
several other more robust estimators with different error distributions that reduce the effect of outliers on the optimized
parameters. However, using these estimators in our analysis isbeyond the scope of this manuscript.

Finally, Fig. 7 shows a response surface in the -n parameter
plane for Scenario 3. The response surface was calculated on arectangular grid size of 50×50. While the  and n parameters
have been varied in the intervals of 0.05–0.1 and 1.46–1.56,
respectively, the values of fixed parameters were r,s, Ks=0.00132 cm/s, b1=43.3, b2=29.6, and b3=106 (kgcm s–3 K–1). The shape of the response surface indicates thatwhile the  parameter is relatively well defined in this interval,
the n parameter is more uncertain, confirming our results reported in Table 2.
CONCLUSIONSIn this study we used HYDRUS-2D to simulate infiltrationof relatively warm water from a single-ring infiltrometer into asilt loam soil profile and to estimate the soil thermal and hydraulics properties by inverse solution from collected experimental data. We showed that the cumulative infiltration andtemperatures collected at four sensors at four different depths
could provide the required information for simultaneous estimation of the soil thermal and soil hydraulic properties.
The results indicate that measurements of cumulative infiltration and temperature data would not necessarily result in
unique and reliable estimates of hydraulic parameters when thenumber of unknown soil thermal parameters is three and soil
hydraulic parameters is more than two. Prior information regarding the optimized parameters or additional measurements
are needed for the estimation problem to be well posed. Wheninfiltration data are used as input in the inversion procedure, theinclusion of prior information on measurable parameters, either
the saturated water content, s or the saturated hydraulic con-13

Mohammad Nakhaei, Jirka Šimůnek

ductivity, Ks, will enhance the likelihood of uniqueness andstability of the inverse solution in the estimation process. Using
available collected data, we cannot prove the uniqueness of ourestimated parameters without including s and Ks as known
parameters. Therefore, in a subsequent study we will be collecting additional data (especially pressure heads and water contents at different locations, as well as independently measured
retention and hydraulic conductivity functions) to be able tobetter evaluate results of our inverse analysis. We also plan toinclude nonreactive chemicals in the infiltrated water to study
the possibility of inverse modeling for additional estimation ofsolute transport properties.
REFERENCES

Bard, Y., 1974. Nonlinear parameter estimation. Academic
Press, New York, NY.
Bohne, K., Roth, C., Leij, F.J., van Genuchten, M.T., 1992.
Rapid method for estimating the unsaturated hydraulic conductivity from infiltration measurement. Soil Sci., 155, 4,
237–244.
Chung, S.O., Horton, R., 1987. Soil heat and water flow with apartial surface mulch. Water Resour. Res., 23, 12, 2175–
2186.
Clausnitzer, V., Hopmans, J.W., 1995. Non-linear parameter
estimation: LM_OPT. General-purpose optimization code
based on the Levenberg-Marquardt algorithm. Land, Air andWater Resources Paper No. 100032, University of California, Davis, CA.
Constantz, J., 1982. Temperature dependence of unsaturated
hydraulic conductivity of two soils. Soil Sci. Soc. Am. J.,
46, 3, 466–470.
Dane, J.J., Topp, G.C., (Eds.), 2002. Methods of Soil Analysis.
Part 4 - Physical Methods. Soil Sci. Soc. Am., Inc., Madison, WI.
de Marsily, G., 1986. Quantitative Hydrogeology. Academic
Press, London.
de Vries, D.A., 1963. The thermal properties of soils. In: vanWijk, R.W. (Ed.): Physics of Plant Environment, North Holland, Amsterdam, pp. 210–235.
Dirksen, C., 1991. Unsaturated hydraulic conductivity. In:
Smith, K.A., Mullins, C.E. (Eds.): Soil Analysis: Physical
Methods. Marcel Dekker, Inc., New York, pp. 209–269.
Dirksen, C., 2000. Unsaturated hydraulic conductivity. In:
Smith, K.A., Mullins, C.E. (Eds.): Soil Analysis: Physical
Methods, second ed., Marcel Dekker Inc., New York, pp.
183–237.
Finsterle, S., Najita, J., 1998. Robust estimation of hydrologic
model parameters. Water Resour. Res., 34, 2939–2947.
Hollenbeck, K., Šimůnek, J., van Genuchten, M.Th.., 2000.
RETCML: Incorporating maximum-likelihood estimation
principles in the hydraulic parameter estimation code RETC.
Comput. Geosci., 26, 319–327.
Hopmans, J.W., Šimůnek, J., 1999. Review of inverse estimation of soil hydraulic properties. In: van Genuchten, M.Th.,
et al. (Ed.): Characterization and Measurement of the Hydraulic Properties of Unsaturated Porous Media. Universityof California, Riverside, California, pp. 643–659.
Hopmans, J.W., Šimůnek, J., Bristow, K.L., 2002. Indirect
estimation of soil thermal properties and water flux fromheat pulse measurements: Geometry and dispersion effects.
Water Resour. Res., 38, 1, 7.1–7.14.
Inoue, M., Šimůnek, J., Shiozawa, S., Hopmans, J.W., 2000.
Estimation of soil hydraulic and solute transport parameters

14

from transient infiltration experiments. Adv. Water Resour.,23, 677–688.
Klute, A., (Ed.) 1986. Methods of Soil Analysis. Part 1, Second
edition. Agronomy Monograph, 9, 635–662, ASA and
SSSA, Madison, WI.,
Kool, J.B., Parker, J.C., 1987. Development and evaluation ofclosed-form expressions for hysteretic soil hydraulic properties. Water Resour. Res., 23, 1, 105–114.
Kool, J.B., Parker, J.C., van Genuchten, M.Th., 1987. Parameter estimation for unsaturated flow and transport models - A
review. J. Hydrol., 91, 255–293.
Nassar, I.N., Horton, R., 1992. Simultaneous transfer of heat,water, and solute in porous media: I. Theoretical development. Soil Sci. Soc. Am. J., 56, 1350–1356.
Marquardt, D.W., 1963. An algorithm for least-squares estimation of nonlinear parameters. Journal of the Society of Industrial and Applied Mathematics, 11, 2, 431–441.
Mortensen, A.P., Hopmans, J.W., Mori, Y., Šimůnek, J., 2006.
Multi-functional heat pulse probe measurements of coupledvadose zone flow and transport. Adv. Water Resour., 29, 2,
250–267.
Mualem, Y., 1976. A new model for predicting the hydraulicconductivity of unsaturated porous media. Water Resour.
Res., 12, 3, 513–522.
Philip, J.R., de Vries, D.A., 1957. Moisture movement in porous materials under temperature gradients. Trans. Am. Geophys. Union, 38, 222–232.
Radcliffe, D.E, Šimůnek, J., 2010. Soil Physics with HYDRUS.
CRC Press, Taylor & Francis Group, Boca Raton, FL,
ISBN: 978-1-4200-7380-5, pp. 373.
Reynolds, W.D., Elrick, D.E., Youngs, E.G., 2002a. Ring or
cylinder infiltrometers (vadose zone). In: Dane, J.H., Topp,
G.C. (Eds.): Methods of Soil Analysis. Part 4 - Physical
Methods. SSSA No. 5, Soil Science Society of America, Inc.
Madison, Wisconsin, pp. 818–820.
Reynolds, W.D., Elrick, D.E., Youngs, E.G., 2002b. Singlering and double or concentric-ring infiltrometers. In: Dane,J., Topp, G.C. (Eds): Methods of Soil Analysis. Part 4 Physical Methods. SSSA No. 5, Soil Science Society ofAmerica, Inc. Madison, WI, pp. 821–826.
Russo, D., Shani, U., Parker, J.C., 1991. Analyses of infiltration, events in relation to determining soil hydraulic properties by inverse problem methodology. Water Resour. Res.,27, 1361–1373.
Saito, H., Šimůnek, J., Mohanty. B., 2006. Numerical analysesof coupled water, vapor and heat transport in the vadosezone. Vadose Zone J., 5, 2, 784–800.
Saito, H., Šimůnek, J., Hopmans, J.W., Tuli, A., 2007. Numerical evaluation of the heat pulse probe for simultaneous estimation of water fluxes and soil hydraulic and thermal properties. Water Resour. Res., 43, W07408.
Sakai, M., Toride, N., Šimůnek, J., 2009. Water and vapor
movement with condensation and evaporation in a sandycolumn subject to temperature gradient. Soil Sci. Soc. Am.
J., 73, 3, 707–717.
Schaap, M. G., Leij, F. J., van Genuchten, M. Th., 2001. Rosetta: A computer program for estimating soil hydraulic parameters with hierarchical pedotransfer functions. J. Hydrol.,
251, 163–176.
Shepard, F. P., 1954. Nomenclature based on sand-silt-clay
ratios. Journal of Sedimentary Petrology, 24, 151–158.
Šimůnek, J., Suarez, D.L., 1993. The UNSATCHEM-2D codefor simulating two-dimensional variably saturated water
flow, heat transport, carbon dioxide production andtransport, and multicomponent solute transport with major

Parameter estimation of soil hydraulic and thermal property functions for unsaturated porous media using the HYDRUS-2D codeion equilibrium and kinetic chemistry. Version 1.1, Research
Report No. 128, U. S. Salinity Laboratory, USDA, ARS,
Riverside, California.
Šimůnek, J., van Genuchten, M.Th., 1996. Estimating unsaturated soil hydraulic properties from tension disc infiltrometerdata by numerical inversion. Water Resour. Res., 32, 9,2683–2696.
Šimůnek, J., Angulo-Jaramillo, R., Schaap, M., Vandervaere,
J.P., van Genuchten, M.Th., 1998. Using an inverse method
to estimate the hydraulic properties of crusted soils from tension disc infiltrometer data. Geoderma, 86, 1–2, 61–81.
Šimůnek, J., Jacques, D., Hopmans, J.W., Inoue, M., Flury, M.,
van Genuchten, M.Th., 2002. Solute Transport During Variably-Saturated Flow - Inverse Methods. In: Dane, J.H.,
Topp, G.C. (Eds.): Methods of Soil Analysis. Part 4 - Physical Methods. Third edition, Chapter 6.6, SSSA, Madison,
Wisconsin, pp. 1435–1449.
Šimůnek, J., Šejna, M., van Genuchten, M.Th., 2006. The hydrus software package for simulating two- and threedimensional movement of water, heat, and multiple solutes
in variably-saturated media. User Manual, Version 1.0, PC
Progress, Prague, Czech Republic.

Šimůnek, J., van Genuchten, M.Th., Šejna, M., 2008. Development and applications of the HYDRUS and STANMOD
software packages, and related codes. Vadose Zone J., 7, 2,
587–600.
Sophocleous, M., 1979. Analysis of water and heat flow inunsaturated-saturated porous media. Water Resour. Res., 15,
5, 1195–1206.
van Genuchten, M.Th., 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil
Sci. Soc. Am. J., 44, 892–898.
Vogel, T., Císlerová, M., Hopmans, J.W., 1991. Porous media
with linearly variable hydraulic properties. Water Resour.
Res., 27, 10, 2735–2741.
Vrugt, J.A., Stauffer, P.H., Wöhling, T., Robinson, B.A., Vesselinov, V., 2008. Inverse modeling of subsurface flow andtransport properties: A review with new developments. Vadose Zone J., 7, 2, 843–864.
Wooding, R.A., 1968. Steady infiltration from large shallow
circular pond. Water Resour. Res., 4, 1259–1273.
Received 11 June 2013
Accepted 22 November 201315