WATER RESOURCES RESEARCH, VOL. 38, NO. 1, 10.1029/2000WR000071, 2002

Indirect estimation of soil thermal properties and water flux using
heat pulse probe measurements: Geometry and dispersion effects
Jan W. Hopmans
Department of Land, Air and Water Resources, University of California, Davis, California, USA

Jirka Šimunek
George E. Brown Jr. Salinity Laboratory, Agricultural Research Service, U.S. Department of Agriculture, Riverside,
California, USA

Keith L. Bristow
CSIRO Land and Water, Townsville, Queensland, Australia
Received 6 November 2000; revised 24 August 2001; accepted 24 September 2001; published 29 January 2002.

[1] Traditionally, analytical solutions for heat transport in soils have been used in combination
with heat pulse probe (HPP) measurements to estimate soil thermal properties. Although theanalytical method has resulted in accurate estimation of soil thermal properties, we suggest thatparameter estimation using inverse modeling (IM) provides new and unique opportunities for soil
thermal characterization. Moreover, we show that the IM approach provides accurate estimationof soil water flux density in both unsaturated and saturated soil conditions for a wider range ofwater velocities than originally thought possible. Specifically, we show that accurate soil watervelocity is obtained, simultaneously with soil thermal properties, if heat dispersion is included in
the heat transport equation. The requirement for including heat dispersivity depends on the valueof the newly defined dimensionless Keith Jirka Jan (KJJ) number, which is equal to the ratio ofthermal dispersion to thermal conductivity. For example, when KJJ > 1, ignoring thermal
dispersivity leads to errors in the water flux density which can exceed 10%. By including thermal
dispersivity, water flow velocities were accurately determined for water flux densities ranging
from 1.0 to >10 m d1. We also demonstrate the general application of inverse modeling toestimate soil thermal properties and their functional dependence on volumetric water content in a
separate numerical experiment. We suggest that inverse modeling of HPP temperature data may
allow simultaneous estimation of soil water retention (when combined with matric potential
measurements) and unsaturated hydraulic conductivity (through water flux estimation) from
simple laboratory experiments.
INDEX TERMS: 1866 Hydrology: Soil moisture; 1875
Hydrology: Unsaturated zone; 1894 Hydrology: Instruments and techniques; KEYWORDS: Soilwater flow; soil heat flow; inverse modeling; dispersivity

1. Introduction[2] In the past decade, soil thermal research has advanced development of the heat pulse technique for estimation of soil thermal
properties, i.e., soil thermal conductivity (l), volumetric heat
capacity (C), and soil thermal diffusivity (k). Apart from needing
to characterize the soil’s physical properties, knowledge of the soil
thermal properties is required for accurate prediction of soil temperature and its influence on seed emergence and crop growth [e.g.,
Gilmore and Rogers, 1958], soil water retention and unsaturated
hydraulic conductivity [Hopmans and Dane, 1986b], and soil water
vapor flow [Nassar and Horton, 1992] in coupled water and heattransport. It is expected that interest in the soil thermal regime will
become increasingly important as both fundamental and appliedresearch questions are posed regarding the water movement inradioactive waste repositories [Tsang and Pruess, 1987] and the fate
and transport of high vapor pressure liquids (e.g., solvents) andmicrobes (e.g., bacteria and viruses). In addition, many chemical
processes in soils are temperature-dependent, and their functioning

Copyright 2002 by the American Geophysical Union.
0043-1397/02/2000WR000071$09.00

in contaminant transport requires better prediction of soil temperature regimes.
[3] The soil thermal properties can be quickly and conveniently
measured using heat pulse probes [Campbell et al., 1991; Bristow etal., 1993, 1994a, 1994b]. The dual probe consists of thin needle-likeheater and temperature probes (1-mm outer diameter), which aremounted in parallel with a 6-mm separation distance. A heat pulse isapplied to the linear heater, and the temperature response is recordedat the sensor probe. Campbell et al. [1991] presented the analyticalsolution for the temperature rise of an instantaneous heat pulse,
allowing estimation of the soil’s volumetric heat capacity andvolumetric water content if some additional data are available.
Subsequently, Bristow et al. [1993] extended the probe to threeneedles to correct for drift in background temperature. Furthermore,they emphasized the relatively large sensitivity of water content
estimations to error in the heater-sensor spacing (r). Continued
development of radial heat transport theory resulted in application of
the heat pulse probe (HPP) method to simultaneously estimate thesoil’s volumetric heat capacity and thermal conductivity (and hence
thermal diffusivity) as well as water content, using a short-duration
heat pulse from an infinite line source [Bristow, 1998]. Error
analyses by Kluitenberg et al. [1993] for the instantaneous heatpulse and by Kluitenberg et al. [1995] for the short-duration heat

7-1

7-2

HOPMANS ET AL.: ESTIMATION OF SOIL THERMAL PROPERTIES AND WATER FLUX

pulse technique highlighted the importance of accurate r and timeto-maximum (tmax) temperature rise measurements on the accuracyof the thermal and water content estimations. It was determined thatrigid needles are required to minimize changes in mutual probe
positions while inserting the heat pulse probe into soils. Typical heatpulse strengths used for making HPP measurements cause a temperature rise of 1C at the sensing probe, so that a 5% precision in heat
capacity estimation requires temperature measurements with
0.05C precision [Kluitenberg et al., 1993].
[4] Further developments have led to the simultaneous measurement of soil thermal properties, water content, and electrical conductivity using TDR time domain reflectometry (TDR) combinedwith the HPP [Noborio et al., 1996; Ren et al., 1999]. Bristow et al.
[2001] showed that a simple modification of the dual probe with anadditional two needles for bulk soil electrical conductivity measurements provides an alternative measure of the soil’s electrical conductivity. In a subsequent study, Ren et al. [2000] reported on thepossibility of using a three-needle heat pulse probe to indirectly
estimate water flux density from temperature responses, measured
upstream and downstream of the heat source. Their experimental
results, using the maximum difference between the upstream anddownstream temperature signals (MDTD), indicated that such analysis can be successful for flux density values larger than 105 m s1(0.864 m d1). However, their results also suggested that the method
is limited at higher fluxes (larger than 2.4 m d1), because of asystematic overestimation of MDTD by 10% or more at these higher
water velocities. The discrepancy may be the result of ignoring thefinite heater geometry in the heat analysis, which assumes infinitely
small point sources and sensors with negligible diameter. Although
not important for zero-water flow conditions [Kluitenberg et al.,
1995], it could be an important consideration when using the HPP for
estimation of water flux. Moreover, the finite physical size of theneedles may distort the water flow field between the needles, thereby
affecting effective flow of water and heat between the heat sourceand sensor. Although the proposed flux density measurement of Ren
et al. [2000] requires a priori knowledge of the thermal properties,
these can be obtained separately from zero-water flux heat probe
measurements, using the method of Bristow et al. [1994a]. Also, Ren
et al. [2000] suggested that the limit of sensitivity could be loweredto 106 m s1 (0.0864 m d1) if the temperature measurement
precision is 0.01C. It is expected that even lower water velocities
can be measured for unsaturated soils because of the reduction in thebulk soil heat capacity with decreasing water content.
[5] In this paper, we hypothesize that difficulties in using the HPP
for water flux density measurements is partly caused by the omissionof dispersive heat transport in the heat flow equation. As solute
dispersion is generally important for solute transport in soils, thermal
dispersion should also be considered for heat transport. For solutetransport, molecular diffusion and hydrodynamic dispersion combined make up the soil’s dispersion coefficient. By analogy, theeffective thermal conductivity (leff) is the sum of the bulk soil
thermal conductivity with stationary fluids (l0) and a kinematic
thermal dispersion term (ld), resulting from the heterogeneity ofwater velocities within and between water-filled soil pores, to be
defined later. It is expected that the influence of the additional
dispersion term on heat transport in soils is most significant onlyat high water velocities, since it is only then that the thermal
dispersion term is of the same order of magnitude or larger thanthe conductive term.
[6] If indeed successful, the multineedle heat pulse techniquemight be an essential new tool for general flow and transport studies
in both the laboratory and the field because of its simple, accurate,

Figure 1. Schematic showing the heat pulse probe (HPP). The
distance L = 30 mm in this particular case.
and versatile design. Moreover, it is becoming increasingly clear thatmultifunctioning sensors are critical for accurate in situ soil physicalmeasurements [Hopmans et al., 1999] because of inherent soil
spatial variability. In particular, we also note the promise of using
inverse methodologies [Hopmans et al., 2002] to infer in situ soil
hydraulic, thermal, and solute transport characteristics and thepossibilities of combining inverse modeling with multineedle heatpulse probe measurements.
[7] The objective of this paper is therefore threefold. First, we
provide a sensitivity analysis to illustrate the effect of the finite heatsource needle on the temperature response for a range of water flux
densities in both saturated and unsaturated soil conditions. We also
use these data to demonstrate the effect of the finite probe
geometry on inverse modeling estimation of water flux and thermal
diffusivity. Second, in an additional sensitivity analysis we show
the influence of thermal dispersion on heat transport and on theestimation of thermal properties and water flow density usingparameter estimation by inverse modeling. The results suggest analternative experimental design for the HPP that will provide moreaccurate water flux density measurements, even at high watervelocities. Finally, the third objective is to evaluate the inverse
modeling method for accurate estimation of soil water content,water flux density, and the water content dependence of soil
thermal properties from HPP measurements.

2. Materials and Methods
2.1. Heat Pulse Probe (HPP)[8] The theory assumes that the soil around the HPP is homogeneous and isotropic and that contact resistance between the soil andthe probe needles is negligible. Moreover, it is assumed that thesensor needle has an infinitely small heat capacity and large thermal
conductivity, so that temperature measurements are instantaneous.
Typical heat pulse strengths and heat pulse durations are 480 J m1
and 8s, respectively, so that the total quantity of heat released is 60
J m1 s1. The assumed distance between the centers of the heaterand sensor needles (r) is 6 mm, whereas the outside diameter of bothneedles is 1 mm. This basic structure depicted in Figure 1 with aneedle length (L) of 30 mm agrees with the recommended dimensions as outlined by Bristow et al. [1994a, 2001]. The diagram in
Figure 2 shows the physical dimensions of the simulated twodimensional soil domain with the three-needle heat pulse probe(heater with upstream (U) and downstream (D) thermocouple) andincludes an additional thermocouple location (T), which is transverse to the principal water flow direction. Later, we will suggest that

7-3

HOPMANS ET AL.: ESTIMATION OF SOIL THERMAL PROPERTIES AND WATER FLUX[11] In addition, the HPP method [Bristow et al., 1994a] allows
easy estimation of the thermal diffusivity k (m2 s1) from thethermal conductivity l0 (W m1 K1) and volumetric heat capacity
Cbulk (J m3 K1) using
ð4Þ

k ¼ l 0 =Cbulk :2.3. Heat Transport Equation[12] The general two-dimensional heat transport equation for ahomogeneous porous medium can be written as [Bear, 1972;
Sophocleous, 1979; Šimunek et al., 1999]
Cbulk ðqÞFigure 2. Simulated soil domain showing front view of heater
needle (H), downstream (D), upstream (U), and transverse (T)sensors. Water flow direction is downstream.
additional temperature information from this transverse sensor
needle will allow estimation of high water flow velocities. The
HPP extends perpendicularly to the illustrated (x-z) plane of the soil
domain, so that the HPP needles are represented by their crosssectional circular planes within the vertical plane.
2.2.



@T
@
@T
@T
¼
 Cw qw;iLeff ;ij ðqÞ
;
@t
@xi
@xj
@xi

where the subcript ij denotes the Einstein notation, T istemperature, and Leff,ij is a second-order tensor of the effective
thermal conductivity. The first term on the right-hand side
represents heat flow due to conduction, whereas the second term
accounts for convective heat transport by water flow as
controlled by the water flux density qw,i. This formulation
assumes instantaneous heat transfer between phases, so that at
any point in the bulk soil the solid, liquid, and gas temperatures
are identical. For an isotropic porous medium the effective
thermal conductivity can be written as [De Marsily, 1986; Bear,1972]:
Leff ;ij ðqÞ ¼ l 0 ðqÞdij þ dij bT Cw jqw j þ ðbL  bT ÞCw

Thermal Properties[9] The soil volumetric heat capacity, Cbulk = (rc)bulk (J m3
K1), can be determined from the sum of the heat capacities of theindividual constituents according to [de Vries, 1963]
Cbulk ¼ ðrcÞs ð1  fÞ þ ðrcÞw q þ ðrcÞa ðf  qÞ;

ð1Þ

where r is the density (kg m3), c is the specific heat (J kg1 K1),
f is the porosity (m3 m3), q is the volumetric water content (m3m3), and the subscripts s, w, and a indicate the soil’s solid phase,water, and air, respectively. If the heat capacity of air is ignored andthe solid phase includes both mineral and organic matter fractions,setting Cw = (rc)w and Cs = (rc)s allows (1) to be rewritten as[Campbell, 1985]
Cbulk ¼ Cs ð1  fÞ þ Cw q:

ð2Þ

It is clear from (2) that the volumetric water content can becomputed directly if Cbulk is measured using the HPP method andvolumetric heat capacity values of the solid phase and water are
known.
[10] The bulk soil thermal conductivity (l0) is a function ofmineral type and geometrical arrangement of the various phases, aswell as the water content [de Vries, 1963], and hence is soil
specific. Although various empirical relationships have been suggested (including those by Campbell [1985], a three-parameter
polynomial expression can be used for the water content range
beyond the initial sharp increase in thermal conductivity with watercontent. This relationship is [Chung and Horton, 1987]
l 0 ¼ b 0 þ b 1 q þ b 2 q0:5 ;where b0, b1, and b2 are empirical constants.

ð5Þ

qw;j qw;i
jqw j

ð6Þ

to include the effect of hydrodynamic dispersion on heattransport. In analogy with solute transport, velocity variations
within the water-filled pore spaces cause mixing of pore waters
with different velocities, thereby causing dispersion-like properties of the temperature field. In (6), jqwj denotes the magnitudeof the Darcy flux density, dij is the Kronecker delta (dij = 1 if
i = j and dij = 0 if i 6¼ j), and bL and bT are the longitudinal andtransverse heat dispersivities (m), respectively. In principle, theirvalues should be of similar magnitude to the dispersivities
defined for solute transport [De Marsily, 1986], but there iscurrently little information about their relationship. For waterflow conditions where flow is uniform and parallel to the z axis,(6) can be simplified toleff ;zz ¼ l 0 þ l d;Lleff ;xx ¼ l 0 þ l d;T ;ð7Þ

where ld,L = bLCwqw,z and ld,T = bTCwqw,z denote thelongitudinal and transversal thermal dispersion coefficients,
respectively. If the water flux density is zero, the convective
term disappears and leff equals the bulk soil thermal
conductivity l0. While applying (7), it must be realized thatthe bulk soil thermal conductivity is relatively large so that
dispersive effects on leff may be small except for large watervelocity values. To quantify the contribution of hydrodynamic
dispersion, relative to the thermal conductivity (l0), we define
the dimensionless Keith Jirka Jan (KJJ) number:
KJJ ¼ l d;L l 0 ;

ð8Þ

ð3Þ
which is equal to one when the thermal dispersion is equal tothe bulk soil thermal conductivity (l0).

7-4

HOPMANS ET AL.: ESTIMATION OF SOIL THERMAL PROPERTIES AND WATER FLUX[13] For uniform vertical water flow conditions the simplified
heat flow equation, after dividing the left-and right-hand side of (5)by Cbulk, can be written as



 

@T
@
@T
@
@T
@T
¼
kxx
þ
kzz
 Vh
;
@t
@x
@x
@z
@z
@z

ð9Þ

where kxx and kzz denote the effective thermal diffusivity in the x
and z direction, respectively, defined by 
l 0 ðqÞ þ l d;T qw;zkxx ðqÞ ¼
Cbulk ðqÞ 
l 0 ðqÞ þ l d;L qw;zkzz ðqÞ ¼
Cbulk ðqÞð10Þ

and where
Vh ¼

Cw qw;z qCw Vw
¼
Cbulk
Cbulk

ð11Þ

denotes the convective heat pulse velocity, describing heat flow by
the moving liquid phase, relative to the stationary bulk porous
medium [Ren et al., 2000]. Since (qCw/Cbulk) < 1.0 (from (2)), it
follows that Vh < Vw, where Vw denotes the average pore watervelocity. As pointed out by Bear [1972], the magnitude of
convective relative to conductive heat transport can be summarized
by the dimensionless Peclet number (Pe)Vw qCw d Vh d
Pe ¼
¼
;
leff
Keff

ð12Þ

where d denotes the characteristic length of the porous medium(assumed to be equal to the mean grain size). The need to include
dispersion in the bulk soil’s thermal diffusivity will depend on themagnitude of the KJJ and Pe numbers, as hydrodynamic effects onheat transport may become significant at increasing pore watervelocities.
2.4.

Inverse Modeling[14] Inverse modeling (IM) has been applied successfully
mostly to the parameter estimation of soil water retention andunsaturated hydraulic conductivity functions. An extensive review
of inverse modeling regarding parameter optimization is presented
by Hopmans et al. [2002]. Using numerical simulations, soil
hydraulic parameters can be determined for transient flow, so that
IM can accommodate much more flexible experimental conditions
than required for steady state methods for which analytical solutions are available. In the analysis to follow, we show how IM canbe applied to coupled heat and water transport. Advantages of IM
for analysis of the HPP measurements are the following: (1)
Various parameters and/or variables can be estimated simultaneously (e.g., heat transport parameters and pore water velocity); (2)IM can be applied to coupled flow of heat and water, allowing
optimization of both heat and water flow parameters; and (3) the
number and type of measurements to achieve a well-posed solution
is unlimited.
[15] The objective function  to be minimized during theparameter estimation process may be defined as
fðp; T Þ ¼ w1nX

½T *ðx; ti Þ  T ðx; ti ; pÞ 2i¼1

þw2 ½T *ud ðt*max Þ  Tud ðt*max ; pÞ 2 ;ð13Þ

Figure 3. (a) Soil hydraulic functions (adapted from Inoue et al.
[2000]) and (b) thermal conductivity function (adapted from
Chung and Horton [1987]) for the Tottori sand.
where the right-hand side represents the residuals between theexperimental (T*) and corresponding predicted (T) temperatures.
The vector x denotes the spatial coordinate of each measurement
i, whereas the vector p contains the optimized parameters. The
first term sums the residuals for all temperature measurements(n), whereas the second term provides the temperature difference
between the upstream (u) and downstream (d) temperature
measurement at time t* = tmax
* , where tmax
* is the time of the
maximum temperature at the downstream sensor. Weighting
factor values for w1 and w2 were set equal to one. The
Levenberg-Marquardt (LM) method was used to minimized (13),thereby including computed confidence intervals for the optimized parameters, assuming zero model error. For the sensitivity
analysis discussed in this paper, the ‘‘experimental’’ temperature
data are provided by simulated data, using forward simulations
with the same flow and transport code as used for the IM(HYDRUS-2D). The results of these types of studies will help inselecting optimum experimental conditions so that IM can beapplied most effectively by yielding accurate and uniqueparameter values.
2.5. HYDRUS-2D Simulation Model[16] The model environment selected for the sensitivity analysisof the HPP is the HYDRUS-2D code [Šimunek et al., 1999], which
incorporates transient simultaneous water and heat flow in two
spatial dimensions (x and z) combined with LM optimization
algorithm for parameter estimation using IM. The domain size
was fixed to 0.025 m (horizontal) by 0.1 m (vertical) with a nodal
spacing of 0.5 mm for the structured mesh used for the point heater
needle (total of 10,000 elements) and a nodal spacing varying
between 0.2 mm (near needles) and 3 mm (near domain boundaries) for the unstructured mesh (total of 1700 elements) used for theHOPMANS ET AL.: ESTIMATION OF SOIL THERMAL PROPERTIES AND WATER FLUX

7-5

Table 1. Soil Hydraulic and Thermal Properties Used in This Study
Parameter

Units
3

qr
qs
a
n
Ks
l
bL
bT
d
b0
b1
b2
Cw
Csa
a

3

m m
m3 m3m1
dimensionlessm s1
dimensionlessmm
mm
dimensionlessdimensionlessdimensionless
J m3 K1
J m3 K1

Value
0.027
0.310
13.8
2.01
0.00045
1.16
0.002210.1 bL
0.36
0.228
2.4064.909
4.18 106
1.92 106

Reference
Inoue et al. [2000]

Šimunek et al. [1999], Bear [1972]
Shiozawa and Inoue [2001]
Chung and Horton [1987]

Cs = (rc)s, with cs = 717 J kg1 K1 and rs = 2650 kg m3.

other simulations that considered a finite size heater needle. The
numerical analysis assumes that the temperature and needle probes
extend infinitely in the y direction, perpendicular to the (x-z) plane of
the simulated soil domain (Figure 2), thereby ignoring possible
three-dimensional effects that may occur for finite HPP needles. The
soil hydraulic properties for the Tottori sand used in the study arepresented in Figure 3a [Inoue et al., 2000], whereas the mean grain
size (d) was taken from Shiozawa and Inoue [2001]. The functional
dependence of thermal conductivity with water content (Figure 3b)
for this sandy soil was taken from Chung and Horton [1987].
[17] To simplify the analysis, either steady state vertical waterflow or no water flow was assumed. Hence water flux is governedby the Darcy flux qw and is written for the respective z and x
directions as



@hm
þ1
qw;z ¼ K ðhm Þ
@z



@hm
qw;x ¼ K ðhm Þ
;
@x

ð14Þ

where K (m s1) denotes the unsaturated hydraulic conductivity
and hm (m) denotes the soil water matric potential. The flux in the x
direction (transverse direction) is included for completeness, sincesome lateral flow will occur around the finite heater needle. Zeroflux is imposed on the lateral boundaries of the soil domain(parallel with the main water flow direction), whereas the imposed
steady state vertical water flux is controlled by water potential
values at the horizontal boundaries of the simulated soil domain.
The soil hydraulic properties are described by the van Genuchten(soil water retention) and van Genuchten-Mualem (unsaturated
hydraulic conductivity) relationships [Van Genuchten, 1980]:
q  qrm¼ ½1 þ ðjh m jÞn
qs  qr

ð15aÞ

h
m i2
K ¼ Ks Sel 1  1  Se1=m
;

ð15bÞ

Se ¼

where Se is the effective saturation, qs, qr, and q are saturated,
residual, and actual water contents, respectively, a, n, and l areconstants, m = 1  1/n, and Ks denotes the saturated hydraulicconductivity. Parameters for both the soil hydraulic and thermal
properties are summarized in Table 1.
[18] Although HYDRUS-2D solves for both transient heat andwater flow, their coupling is limited to the mutual occurrence ofwater flux density and water content in both transport equations.

Simultaneous solution of both transport equations is done by firstsolving the water flow equation and by subsequent substitution of
the computed flux density and water content values into the heatflow equation for each time step. However, temperature effects on
soil hydraulic properties [Hopmans and Dane, 1986a] and thermal
properties [Hopmans and Dane, 1986b] were ignored. Moreover,distillation effects causing enhanced heat transport by latent heatthrough vaporization and subsequent condensation [Cass et al.,1984] were also not considered. Although we expect that temperature effects are minor for the experimental temperature regimes
presented herein, subsequent studies will fully account for thetemperature influence on flow and transport.
[19] Several small modifications of HYDRUS-2D were implemented to make numerical heat transport experiments possible.
The public version of HYDRUS-2D can handle either Dirichlet(temperature) or Cauchy (heat flux) boundary conditions, but it has
no provisions for a heat source without a corresponding water flux.
The heat pulse probe boundary condition was implemented byadding a zero-order source term to (5), which was active only forthe nodes representing the heat pulse probe during heat pulse
application. In simulations for which a single node represented theheat pulse probe (case P), the corresponding heat flux at this single
node was 57.7 W m1. When simulating the finite heat source
experiments (case F), several nodes were used to describe thegeometry of the heat pulse probe, with the heat source strength
distributed uniformly across these nodes.
2.6. Numerical Experiments[20] Corresponding to objectives 1 through 3, we present three
numerical heat transport experiments (experiments I, II, and III).
The first experiment (I) analyzes the influence of the finite probe
geometry of the heater probe on the optimized soil thermal andwater flux density values. To achieve this, we compare optimization
results with a line heater needle of infinitely small diameter (case P)
with those using a finite size 1-mm-diameter heater probe (case F).
In the second experiment (II) we compare results of inverse
modeling simulations with and without the dispersion term in (7).
We show results for both unsaturated and saturated flow conditions
and for a range in water flux density (qw) and thermal dispersivity(b) values. We show the contribution of the dispersive component
to be a function of the KJJ and Pe numbers, with increasing
significance as the water flux density becomes larger. Finally, inthe third experiment (III) we demonstrate how the functional
dependence of the soil thermal properties with water content and7-6

HOPMANS ET AL.: ESTIMATION OF SOIL THERMAL PROPERTIES AND WATER FLUX

Figure 4. (a) Simulated spatial temperature distributions, 15 and
30 s after heating, for qw values of 0, 1, and 10 m d1 (b) Simulated
temporal distributions at the upstream (U) and downstream (D)temperature sensors for qw values of 0, 1, and 10 m d1 (no thermal
dispersion).

pore water velocity can be estimated from a single experiment,using inverse modeling. The outcome of this experiment willshow that the functional dependence of bulk soil thermal conductivity and bulk soil heat capacity with water content can be
estimated using IM.
[21] Parameter optimizations were carried out with a single setof initial estimates, in which we evaluated effects of needle
geometry (experiment I) and thermal dispersion (experiment II)
on parameter estimation of thermal and water flow properties.
Although, in general, we recommend using various combinationsof different initial estimates, a single optimization was justified
because of the relative closeness of the initial parameter values
to their true values. Multiple sets of initial parameter estimates
were used in experiment III, with optimization results reported
only for the smallest value of the objective function. The initial
estimate for qw was twice its true value; initial estiamtes for q
were 0.33 and 0.28 m3 m3 for the saturated and unsaturated
problems, respectively, whereas the initial value of l0 was set to
2.5 W m1 K1. For cases with variable bL its initial value was
set to twice its true value as was done for the flux density
estimate.

3. Results and Discussion
3.1.

Analysis of Convective Heat Transport[22] To improve our intuitive understanding of the contribution of diffusive and convective heat transport during HPP
measurements, we show in Figure 4a the spatial temperature

distribution around the heater (H) for two simulation times,
following a typical 8-s heating period (57.7 W m1) at location
H. For illustrative purposes the origin is defined at the heater
location, and the spatial coordinate is along the z direction only.
For this example the soil is saturated (qs = 0.31 m3 m3), andtemperature distributions are presented for water flux density(qw) values of 0.0 (no flow), 1.0, and 10.0 m d1. The finite
dimensions of the heater source (1 mm) were accounted for in
the flow calculations. In addition, Figure 4b shows the corresponding temperature values as a function of time at both thedownstream (D) and upstream (U) thermocouple locations. Notethat for both flow velocities used here, the KJJ number is equal
to zero since dispersion was not considered.
[23] Focusing on the case without heat convection first (noflow), the symmetrical temperature distribution around the heater
location in Figure 4a clearly shows the dissipation of the temperature signal around the heater. Heat transport is by conduction only,
causing identical pulse-type temperature signals at thermocouple
locations D and U (no flow (Figure 4b)). When including convection, the temperature signal generated at position H will movein the direction of D with a convective heat pulse velocity Vh of15.95 m d1, as calculated from (11) for a qw value of 10 m d1(or averagepore water velocity Vw of 32.3 m d1). This will movethe heat pulse 0.0055 m after 30 s. As is evident by visualinspection of Figure 4a, the temperature signals for the case with
qw = 1.0 m d1 is about equal to the no-flow case. Hence it isexpected that accurate water flow velocity measurements with the
HPP are unlikely for these or lower water flow velocities. Clearly,
the sensitivity of the water flux density measurements will dependon the temperature measurement accuracy, but as Ren et al. [2000]
concluded, the lower detection limit of water velocity is 0.06 m
d1 for a thermocouple measurement accuracy of 0.01 K. For
steady state water flow, with uniform water content distribution
along its pathway, the shape of the temperature signal should beapproximately independent of flow rate but translated over a
distance determined by the convective heat pulse velocity Vh.
Inspection of the temperature signal as a function of time (Figure4b) also shows that the downstream peak temperatures will
increase with increasing water flow velocities. This is caused by
the slightly smaller travel time of the faster heat pulse (notedifferences in tmax between water flow velocities, as determinedfrom its position relative to the dotted vertical line), thereby
reducing lateral heat loss by conduction. As is expected, thetemperature peak of the downstream thermocouple is always
higher than for the upstream location, with the difference between
thermocouple locations increasing as the flow velocity increases.
The results of Figure 4 also clearly illustrate that incorrect thermal
property estimations with HPP measurements are obtained forwater fluxes >1.0 m d1, unless the water flux is taken into
consideration.
[24] Figure 5 shows the effect of dispersivity on the temperature
signal at the upstream and downstream needle, using a qw value of
1 m d1. To document the heat movement in the transverse
direction, we also included an additional measurement point 6mm away from the heat source in the horizontal direction, i.e.,
perpendicular to the direction of water flow (transverse thermocouple T (Figure 2). As is obvious from (5) and (7), dispersion has
an effect similar to thermal conduction on the temperature distribution. However, while thermal conduction acts symmetrically inall directions, dispersion manifests itself preferentially in the
direction of water flow (since it is usually assumed that bT =0.1bL). For the case with no flow the temperature profile will be7-7

HOPMANS ET AL.: ESTIMATION OF SOIL THERMAL PROPERTIES AND WATER FLUX

Figure 5. Temperature as a function of time for upstream, downstream, and transverse thermocouple locations fordispersivity values of (a) 0.0, (b) 0.001, (c) 0.01, and (d) 0.1 m (qw = 1.0 m d1).
identical for all three measurement positions, as demonstrated in
Figure 4a.
[25] Water flow (without dispersion) results in an asymmetrical temperature profile, with a larger temperature peak at thedownstream location than at the upstream location and with thetransverse thermocouple recording intermediate temperatures(Figure 5a). Larger dispersivities will cause higher temperaturepeaks (Tmax) and earlier arrival of the temperature front (tmaxdecreases) at both downstream and upstream locations becauseof the increased effective thermal conductivity (Figures 5b
through 5d). As the effective conduction (including the dispersion contribution) becomes dominant over convection (smaller
Pe numbers), the temperature difference at the peak between the
upstream and downstream tends to disappear (Figure 5d),thereby reducing the flux information. This result is in agreement with Ren et al. [2000], who found increasing discrepancies
between measured and predicted temperature fields for increasing flux densities within their range of water flux densities
between 2 and 6
105 m s1 (1.73 – 5.18 m d1). Table 2
summarizes numerically results presented in Figure 5 asobtained with different values of dispersivity for the same water

flux of 1 m d1. Note that the results span conditions for heattransport dominated by conduction (KJJ values between 0 and
0.022) to conditions where the dispersion component is about
twice as large as conduction (KJJ = 2.2). We find that while
increasing the dispersivity (at this water flux of 1 m d1) fromzero to 0.1 m, the corresponding temperature difference (Tmax)
between the upstream and downstream locations decreases bymore than a half. Consequently, using the predicted temperature
difference (Tmax) assuming zero dispersion will underestimate
the water flux density [see Ren et al., 2000, Figure 9], increasingly more so for soils with larger dispersivities or withincreasing water velocities.
[26] As dispersion acts preferentially in the direction of water
flow, less heat is available for conduction in the lateral flow
directions (perpendicular to water flow), thereby significantly
reducing the peak temperature at the transverse thermocouple
location (T). For a dispersivity of only 0.01 m the peak temperature
at the transverse location is already smaller than at the upstream
location (Figure 5c). Thus dispersion causes differences in temperature between longitudinal (upstream and downstream thermocouple locations) and transverse locations (transverse thermocouple

Table 2. Influence of Longitudinal Dispersivity bL on Maximum Temperatures at the Downstream (TD), Upstream (TU), and Transverse(TT) Locations at Time tmax on the Temperature Difference (Tmax) Between the Upstream and Downstream Locations at Time tmax andon the KJJ NumberbL, mtmax,a sTD, C

TU, CTT, CTmax, C

KJJ,b dimensionless

0
0.001
0.01
0.05
0.10

15.13515.13513.666
10.6709.452

20.60120.606
20.64120.74720.80620.52920.53520.57820.70320.77520.56020.55520.50820.35520.271

0.072
0.071
0.063
0.044
0.031

0
0.022
0.22
1.1
2.2

a
b

Here tmax is time of the maximum temperature at downstream location.
Here qw = 1 m d1.

7-8

HOPMANS ET AL.: ESTIMATION OF SOIL THERMAL PROPERTIES AND WATER FLUX

Water flux density distribution along transect AA0 of

Figure 7. Comparison of estimated with ‘‘true’’ flux density
values for the case with no thermal dispersion (experiment I).

needle). We therefore suggest that inclusion of a third transverse
thermocouple needle in the heat pulse sensor (see Figure 2) isneeded to extend the range of accurate water velocity measurements to higher flow velocities.

near the top and bottom of the heater needle, leading to reduced
water flow velocities in these regions. An example of the vertical
distribution of water flow rate (qw) for saturated conditions (qw = 1
m d1) as calculated numerically across the AA0 section (Figure 2)
is given in Figure 6. Velocities decreased as the moving water
approached the top of the needle but increased as water moved
around the needle and then decreased again immediately below theneedle. The maximum water velocity was 3 times larger than theminimum velocity (Figure 6).
[28] A second effect may be caused by the decreased distancefor the heat pulse to move before arriving at the downstream sensor
needle. Specifically, for a needle diameter of 1.0 mm and a centerto-center distance (r) between the heater and downstream needlesof 6.0 mm the finite diameter of the heater needle reduces the total
travel distance of the heat front from 6.0 to 5.5 mm, or 8%, ascompared to a 6-mm travel distance assuming that the heater
needle is inifinitely thin.
[29] In order to evaluate these two effects of needle geometry on
estimation of thermal and water flow properties we carried out anumerical study, in which we considered the two following cases.
In the first case we simplified the problem by assuming a point
heater source (case P), thus neglecting effects of geometry, i.e.,

Figure 6.
Figure 2.

3.2. Experiment I: Effect of Needle Geometry on IM
of Thermal and Water Flow Properties[27] We assume no significant influence of the geometry of thetemperature needle of the HPP on temperature measurements, as
we expect that the thermocouple will sense an integrated value
around the sensor needle. We thus assume that the higher temperature arriving at the upstream end of the temperature needle due toheat pulse will be balanced by the lower temperature at thedownstream end of the needle. The finite geometry of the heater
probe can, however, affect the temperature response of the sensor
probe in two ways. First, for conditions of unidirectional waterflow in the z direction the physical presence of the heater probe
will distort the flow field of water around the heater needle, thereby
affecting the convective heat transport component in (5). Simulation results, indeed, showed the presence of such flow distortion.
Specifically, water flow converged in the x direction around theheater needle to cause increased local flow velocities but diverged

Table 3. Influence of Probe Geometry on Parameter Optimization Results With No Dispersiona
Case

P
F
P
F
P
F
P
F
P
F
P
F
P
F
P
F
a
b
c

Peb

True Values
1qw, m d1

qw, m dq,dimensionlessl0,
W m1 K1

0

0.310

2.215

0

0

0.155

1.788

0

0.1

0.310

2.215

0.0008

0.5

0.310

2.215

0.004

1.0

0.310

2.215

0.008

1.0

0.155

1.788

0.01

5.0

0.310

2.215

0.04

10.0

0.310

2.215

0.079

ErrorcOptimized Values

0.09890.09730.5005
0.497
1.016
0.994
1.013
1.001
5.031
4.994
9.995
9.997

P represents the point case; F represents the finite case. No dispersion is bL = 0.
Cw = 48.38 W d m3 K1, and d = 0.00036 m.
Error = [( yopt  ytrue)/ytrue]100%, y = qw , q, or l0.

q,dimensionlessl0,
W m1 K1

qw, %

q, %

l0, %

0.29960.310
0.14740.155
0.30350.310
0.29990.309
0.30360.311
0.14950.15480.30310.308
0.3029
0.310

2.286
2.215
1.845
1.778
2.308
2.217
2.287
2.210
2.307
2.195
1.857
1.787
2.294
2.216
2.279
2.216

...
...
...
...
1.1
2.7
0.1
0.6
1.6
0.6
0.3
0.0
0.6
0.1
0.1
0.0

3.4
0.0
4.9
0.0
2.1
0.0
3.3
0.3
2.1
0.3
3.5
0.0
2.2
0.6
2.3
0.0

3.2
0.0
3.2
0.6
4.2
0.1
3.3
0.2
4.2
0.9
3.9
0.0
3.6
0.0
2.9
0.0

1.0

1.0

1.0

10.0

10.0

1.0

1.0

1.0

0.00640.00380.0025

0.064

0.025

0.00950.00770.0026

0.155

0.155

0.155

0.310

0.310

0.310

0.310

0.310

0.310

0.1

0.01

0.001

0.01

0.001

0.10

0.05

0.01

0.001

bL, m

True Values

1.788

1.788

1.788

2.215

2.215

2.215

2.215

2.215

2.215

6.62

2.27

1.83

7.05

2.70

7.05

4.63

2.70

2.26

l0, W m1 leff, W m1
K1
K1

2.7

0.27

0.027

2.2

0.22

2.2

1.1

0.22

0.022

KJJ,
dimensionless

1.192
0.973

0.984
1.156

0.884
0.994

10.38
13.33

9.135
13.57

0.602
1.151

0.720
1.133

0.904
1.037

0.989
1.042

qw ,
m d1
0.304
0.334
0.310
0.259
0.338
0.315
0.144
0.388
0.309
0.076
0.400
0.313
0.259
0.384
0.310
0.0012
0.511
0.310
0.103
0.152
0.156
0.149
0.225
0.155
0.0002
0.144
0.153

q,
dimensionless

1.84
1.797

1.56
1.786

1.80
1.784

2.04
2.215

2.31
2.212

2.36
2.196

2.18
2.220

2.17
2.196

2.13
2.212

5.36

1.81

2.01

6.87

2.46

4.33

3.35

2.46

2.24

l0, W m1 leff, W m1
K1
K1

Optimized Values

0.097
0.099

0.012
0.01010.00960.00110.011
0.01001

0.00210.001010

0.10
0.10120.056
0.04970.013
0.01000.00440.00105

bL, m

b

a

Cw = 48.38 W d m3 K1, and d = 0.00036 m.
Error = [(yopt  ytrue)/ytrue]100%.
c
For case a, assume bL = 0; optimize qw, q, and leff(=l0). For case b, optimize bL, qw, q, and l0. For case c, fix qw to its ‘‘true’’ value; optimize bL, q, and l0.

1.00.0077

1ac
1b
1c
2a
2b
2c
3a
3b
3c
4a
4b
4c
5a
5b
5c
6a
6b
6c
7a
7b
7c
8a
8b
8c
9a
9b
9c

qw ,
q,
m d1 dimensionless

Pea

Case

Table 4a. Influence of Dispersion on Parameter Optimization Results Using Information From the Upstream and Downstream Locations

19.2
2.7

1.6
15.6

11.6
0.6

3.9
33.3

8.6
35.7

39.8
15.1

28.0
13.3

9.6
3.7

1.1
4.2

qw, %

1.9
7.7
0.0
16.5
9.0
1.6
53.5
25.2
0.3
75.5
29.0
1.0
16.5
23.9
0.0
99.6
41.1
0.0
33.5
1.9
0.6
3.9
45.2
0.0
99.9
7.1
1.3

q, %

2.9
0.5

12.8
0.1

0.6
0.2

7.9
0.0

0.1

6.5
1.1

1.6
0.2

2.0
1.1

3.8
0.1

l0, %

Errorsb

19.0

20.3

9.8

2.5

8.9

38.6

27.6

8.9

0.9

leff, %

3.0
1.1

20.0
0.1

4.0
10.0

10.0
0.1

110.0
1.0

0.0
1.2

12.0
0.6

3.0
0.8

340
5.0

bL, %

HOPMANS ET AL.: ESTIMATION OF SOIL THERMAL PROPERTIES AND WATER FLUX

7-9

1.0

1.0

1.0

10.0

10.0

1.0

1.0

1.0

0.00640.00770.0038

0.064

0.038

0.00950.00770.0026

3

1

0.155

0.155

0.155

0.310

0.310

0.310

0.310

0.310

0.310

q,
dimensionless

0.1

0.01

0.001

0.01

0.001

0.10

0.05

0.01

0.001

bL, m

True Values

1.788

1.788

1.788

2.215

2.215

2.215

2.215

2.215

2.215

l0, W m1
K1

6.62

2.27

1.83

7.05

2.70

7.05

4.63

2.70

2.26

leff, W m1
K1

2.7

0.27

0.027

2.2

0.22

2.2

1.1

0.22

0.022

KJJ,
dimensionlessb

Cw = 48.38 W d m K , and L = 0.00036 m.
Error = [(yopt  ytrue)/ytrue]100%.
c
For case d, fix qw to its ‘‘true’’ value; optimize bL, q, and l0. For case e, optimize bL, qw, q, and l0.

a1.00.0077

1dc
1e
2d
2e
3d
3e
4d
4e
5d
5e
6d
6e
7d
7e
8d
8e
9d
9e

qw ,
m d1

Pea

Case

0.999

0.998

0.999

9.996

9.994

0.998

1.000

0.999

0.999

qw ,
m d1
0.310
0.310
0.311
0.310
0.310
0.310
0.310
0.310
0.310
0.310
0.310
0.310
0.155
0.155
0.155
0.155
0.155
0.155

q,
dimensionless
2.213
2.215
2.224
2.215
2.213
2.215
2.215
2.216
2.215
2.214
2.216
2.216
1.788
1.788
1.788
1.787
1.788
1.788

l0, W m1
K1

Optimized Values

0.001020.000991
0.01040.01000.04980.05000.10080.10020.001004
0.000999
0.10000.01000
0.00100.00100.01000.0100
0.100
0.100

bL, m

0.1

0.2

0.1

0.0

0.0

0.2

0.0

0.1

0.1

qw, %

0.0
0.0
0.3
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0

q, %

0.1
0.0
0.4
0.0
0.1
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0

l0, %

Errorsb

2.0
0.9
4.0
0.0
0.4
0.0
0.8
0.2
0.4
0.1
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0

bL, %

Table 4b. Influence of Dispersion on Parameter Optimization Results Using Temperature Information From the Transverse Location in Addition to Upstream and Downstream Locations

7 - 10
HOPMANS ET AL.: ESTIMATION OF SOIL THERMAL PROPERTIES AND WATER FLUXHOPMANS ET AL.: ESTIMATION OF SOIL THERMAL PROPERTIES AND WATER FLUX

Figure 8. Comparison of estimated with ‘‘true’’ flux density
values without (solid circles) and with (open circles) inclusion of
the transverse (T) temperature data in the objective function(experiment II).

flow distortion and decreased distance between probes. In thesecond case, by considering a finite geometry source (case F)
we took geometry effects fully into account. The sensitivity
analysis involved the following procedural steps. First, a forward simulation was conducted with the HYDRUS-2D code,with specific values for q, steady state vertical water flux
density (qw), a finite geometry source, and the hydraulic andthermal properties of Table 1. The resulting simulated T(z,t)
data for both the upstream and downstream thermocouple
sensor (case F) were assumed to be the ‘‘true’’ temperature
values (as measured with the 1.0-mm-diameter heater needle).
In the subsequent inverse modeling (IM) step these true data
were used in the objective functions as the measured values,
whereas q, qw, and leff were optimized for both the point (P)and finite (F) case. The simulation results are given in Table 3.
The estimation of these three parameters allows direct computation of the thermal diffusivity (k) and convective heat pulse
velocity (Vh), provided the bulk volumetric heat capacity (Cbulk)
is known. Hence it is assumed that bulk soil density (rb) andsoil specific heat (cs) values are available. For mineral soils(contribution of organic matter content to bulk soil heat
capacity can be ignored) the cs is usually between 0.7 and
0.9 J kg1 K1 and varies with mineral type and quartz
content [de Vries, 1963]. Although ignored, there is also sometemperature dependence on cs [Kay and Goit, 1975; Hopmans
and Dane, 1986a].
[30] The geometry effect may be considered negligible if the
optimized heat and flow parameters are equal to their true values.
The error in the last three columns of Table 3 is expressed relative
to the true value. When considering all results combined, ignoring
the effects of the true geometry of the heat source on heat transport
resulted in a slight underestimation of the water content values by
3% and an overestimation of the thermal conductivities by 3 –4%. For no-flow conditions the error analysis by Kluitenberg et al.
[1995] concluded that errors in the thermal properties were <1%using a ratio of heater radius to sensor needle spacing of 0.08 and aheating time of 8 s. Therefore we tentatively conclude that theconvective heat transport component increases the uncertainty ofthe soil thermal parameters. However, if the finite probe geometry
is included in the heat transport simulations, we conclude thatpredictions of water flux density, water content, and thermal
conductivity are excellent, with estimation errors <1% for qw >7 - 11

0.5 m d1. Moreover, even though the water flux estimation error
is close to 3% for the smallest water velocity (qw = 0.1 m d1), theprediction errors of water content and thermal conductivity arepractically zero. We must note that the finite probe geometry
simulations assumed that all HPP needles have an infinitely smallheat capacity and large thermal conductivity, so that the ‘‘true’’temperature distribution is not affected by the thermal properties ofthe HPP. The real thermal properties of the high-conductive epoxy
within the needles though are about equal to those of a wet soil,
implicating that future modeling studies should more accurately
represent the thermal properties of the HPP. The estimated qw
values are compared with the corresponding true values in Figure
7. We note that the optimized values in Table 3 were obtained forboth unsaturated and saturated soil conditions with water flux
density values ranging between 0.1 and 10.0 m d1, whileassuming zero thermal dispersion.
[31] Taking into account the physical size of the heater needle isalso important for optimizing experimental conditions. For example, heat pulse strength must be controlled so that maximum soil
temperatures will not enhance water flow as may occur by vapor
transport. For equal source strengths, because of its larger surface
area, the finite heat source will cause a smaller increase in local soil
temperature than the point source. For example, on the basis ofnumerical simulations, maximum soil temperature at the sourcesoil interface for the 57.7 W m1 heat pulse reached 35C for the
point source, while it was 29C near the finite source (1-mm
diameter). Hence one can safely use larger energy sources, thereby
improving the resolution at the sensor needle, when using larger
heater needles. For example, if the temperature range during the

Figure 9. (a) Temperatures measured with heat pulse probes atfive locations within the soil column and (b) temperatures along
the z axis at four different times for experiment III.

7 - 12

HOPMANS ET AL.: ESTIMATION OF SOIL THERMAL PROPERTIES AND WATER FLUX

b2

4.195
4.900
(0.0443)
4.865(0.0406)

experiment must be below 35C, a stronger source than 57.7 W
m1 can be used when considering the actual dimension of theheater needle.

c

b

2.105
0.260
0.05

1.919
0.193
0.10

1.766
0.150
0.15

1.651
0.123
0.20

1.544
0.105
0.25
a

The 95% confidence intervals are indicated within parentheses. (True values of b0 = 0.228, of b1 = 2.406, and b2 = 4.909).

106
1.920

Cs, J mK1

l0, W mK1q,dimensionless
Position
z, m
Case

a

1.550
(0.001)1.617(0.0044)1.720(0.0051)
1.886(0.0018)
2.056(0.0047)0.110(0.0001)
0.128(0.0006)
0.160(0.0007)
0.196(0.0002)
0.265(0.0007)l0, W m1
K1q,dimensionless

Estimated

3
1

True Values

Table 5. Estimation of Thermal Properties Using Inverse Modeling of a Simulated Column Experimenta

Cs, J m3
K1

1.920 106(242.2)

0.332
0.230
(0.0089)
0.237(0.0082)b0, W m1
K1

b1

1.628
2.395(0.0533)
2.354(0.0489)

3.3. Experiment II: Effect of Thermal Dispersion on
Parameter Estimation of Thermal and Water Flow Properties[32] To investigate the influence of thermal dispersion on soil
heat transport and on indirectly estimated heat transport and waterflow parameters, forward simulations were conducted for a range ofwater flux density, water content, and dispersivity values, as
indicated by the true values in Tables 4a and 4b. Cases 1, 2, and4 from Tables 4a and 4b, as well as the case with zero dispersion, areshown in Figure 5. Initial estimates of optimized parameters wereselected as in experiment I. We note that the dispersion contribution
to the total effective thermal conductivity increases as the flux
density (qw) and dispersivity (bL) increase and the water content (q)decreases. The increasing dispersion effect is described by increasing values of the KJJ number. In the second IM step the resulting
temperatures at the upstream and downstream locations were used
as input to the inverse model, from which the water flux qw, thewater content q, and either the effective thermal conductivity leff(zero thermal dispersion (case a)) or the bulk soil thermal conductivity l0 and the thermal dispersivity bL (case b) were estimated.
[33] It is rather difficult to interpret the optimization results when
information from only the upstream and donwstream locations isused. In case a, optimization results tend to underestimate water
fluxes and significantly underestimate water contents and effective
thermal conductivities with increasing true dispersivity. If thermal
dispersion is neglected (case a), the bulk soil thermal conductivity is
increasingly overestimated with increasing dispersivity. Optimization results tend to slightly improve when l0 and bL are estimatedsimultaneously (case b). When dispersion is included in the inverse
analysis, dispersivity is estimated with a reasonable precision, but
the water fluxes, as well as the water contents, tend to be overestimated. Temperature signatures at the upstream and downstream
locations apparently do not contain enough information for successful simultaneous optimization of these two sets of parameters,i.e., either qw, q, and leff in case a or qw, q, l0, and bL in case b.
[34] The difficulty of the simultaneous estimation of flux qw anddispersivity bL stems from the fact that these two processes tend tohave opposite effects on the difference in temperatures at the
upstream and downstream locations. Increasing flux tends toincrease this temperature difference, while increasing dispersivity
tends to decrease it. At the same time, increasing dispersion has thesame effect as increasing heat conductivity on the arrival of thetemperature peak at these two locations. Optimization results
improve significantly when water flux qw is known. In this case
(case c) the water content q, the bulk soil thermal conductivity l0, andthe heat dispersivity bL were estimated with reasonable precision.
[35] As hypothesized above, valuable information about fluxand dispersivity can be collected when temperature measurements
at the upstream and downstream locations are complemented withtemperature measurements at the transverse location of the heater
needle. This is demonstrated in Table 4b where temperatures
measured at three locations (upstream, downstream, and transverse) were used to estimate simultaneously q, l0, and bL for case dor additionally qw for case e. In either case, leff values are
determined from optimized bL and qw values using (7) with aknown Cw value. Inspection of the results in Table 4b reveals thatboth parameter sets were successfully back calculated using
inverse analysis with zero estimation errors for all but one case.
A comparison of estimated with true water flux density values forcases c (without transverse temperature measurememts) and e (withHOPMANS ET AL.: ESTIMATION OF SOIL THERMAL PROPERTIES AND WATER FLUX

Figure 10. True and fitted thermal conductivity functions for thecolumn experiment III. Thermal conductivity function was fitted tothe estimated data values of case a.
transverse temperature measurements) is presented in Figure 8. We
note that the number of displayed data points is nine (see Table 4b)
for each case, however, most of these coalesce because of nearly
identical values. These results suggest that information from all
three sensor needles can be used advantageously in estimating
simultaneously the water flux and the thermal dispersivity, especially at the higher water fluxes.
3.4. Experiment III: Parameter Optimizationof Thermal Properties[36] To demonstrate the application of IM for estimation of thefunctional dependence of the soil thermal properties with watercontent, a numerical column experiment was conducted. We
suggest that a similar laboratory experiment may be used for
estimating a soil’s thermal properties, as determined by its thermal
conductivity and bulk soil heat capacity.
[37] The simulated soil column was 30 cm long, with water
being in static hydraulic equilibrium with the water table locatedat the base of the column. Hence the soil water matric potential
was zero at the base of the column and linearly decreased to 30
cm at the soil surface. A no-flow condition was selected since it
can be achieved relatively easily in a laboratory. However, wesuggest that the same experimental setup can be used for steadystate flow, allowing simultaneous estimation of soil thermal andhydraulic function parameters. Soil hydraulic and thermal properties used in this experiment are identical to those used in theprevious experiments and are summarized in Table 1. Although
the matric head range for this experiment may appear small, thecorresponding water content range for the Tottori sand (seeFigure 3a) was between 0.31 (base of column) and 0.11 (topof column). The temperature responses of five vertically placedheat pulse probes were simulated at 5, 10, 15, 20, and 25 cm
below the soil surface. Using the vertical HPP arrangement, thecentral heater needle was placed at the exact depth, whereas thetwo sensor needles where either above or below the heater. An 8s heat pulse of 57.7 W m1 was applied to each probe, andtemperature responses were recorded at the sensor needles duringa subsequent 120-s monitoring period. The temperature information was used in the inverse analysis by including T(z,t) data inthe objective function of (13). Because water content values
differed between probe locations (from 0.105 to 0.260 m3
m3), the temperature signatures were different for each depth(Figure 9), with increasing temperature amplitudes for the nearsurface sensors because of the corresponding decrease in soil
water and hence bulk soil heat capacity.

7 - 13[38] We conducted two inverse modeling analyses with thesimulated temperature data (Table 5). First, the temperature signatures of each HPP were analyzed separately (case a), from whichwater content and thermal conductivity data were estimated foreach location. The resulting soil water retention (matric potential is
known since the flow experiment is in hydraulic equilibrium) andthermal conductivity data can be directly used to fit soil water
retention and thermal conductivity function parameters. Results of
the optimized thermal conductivity data are presented in Table 5
and in Figure 10 (diamonds and dashed line). We note that watercontent values were slightly overestimated, leading to a small
underestimation (2%) of thermal conductivity across the experimental water content range. This small estimation error was likely
caused by the underlying assumption that each HPP measurement
represented uniform water content for the soil between the HPP
needles. However, because of the vertical arrangement of the HPPs
in a soil column that is in hydraulic equilibrium with the watertable, the soil water content between the heater and downstream
sensor needles must increase downward. Vertical placement of theheat probe was chosen, since it allows simultaneous estimation ofvertical flux density, either downward (infiltration) or upward(evaporation). Hence we would expect an even more accurate
estimation of the thermal conductivity function if the HPP were
installed horizontally. The IM method used does not allow for thesimultaneous estimation of the bulk soil heat capacity and thewater content, since these two variables are fully correlated by thethermal diffusivity term in (9).
[39] In the second estimation procedure (cases b and c) thetemperature signatures at all five locations were analyzed simultaneously so that parameters of the thermal conductivity function were
estimated directly, with (case b) or without (case c) additional
estimation of the bulk soil heat capacity of the solid phase (Cs).
Simultaneous estimation of both soil thermal properties is now
possible, since the water content is assumed to be known a priorifrom the soil water retention curve. In both cases the estimatedthermal conductivity function parameters were practically identicalto their true values, and the function coincides with the true functionin Figure 10 (solid line). Also, the heat capacity of the mineral phase
could be estimated accurately, once the water content was known.
[40] In summary, this third experiment (III) has shown that thefunctional dependence of the bulk soil thermal conductivity andbulk soil heat capacity with water content can be estimated from a
simple laboratory experiment when integrated with inverse modeling. Moreover, if the bulk soil heat capacity is known, thesuggested experimental setup allows simultaneous estimation ofsoil water retention (when combined with matric potential measurements) and unsaturated hydraulic conductivity (through water
flux estimation) parameters.

4. Summary and Conclusions[41] Traditional data analysis of the heat pulse probe is based on
the analytical solution of the heat transport equation in a uniform
static or steady state water flow field. Although such analysis
proved to be a powerful tool in estimating soil thermal properties,some of the underlying assumptions may not be valid in typicallaboratory and/or field situations. In this study, we used a numerical solution of the heat transport equation to evaluate the effects ofthe physical size of the heater needle on the water flow field andthe temperature response for a range of water flux densities. We
conclude that error estimates of water flux and soil thermal
properties can be reduced from 3% to near zero if the physical7 - 14

HOPMANS ET AL.: ESTIMATION OF SOIL THERMAL PROPERTIES AND WATER FLUX

size of the heater needle is included in the numerical solution of theheat transport equation.
[42] The influence of thermal dispersion on heat transport andestimated thermal properties and water flow densities was analyzed
using parameter estimation by inverse modeling. On the basis of
our analyses we suggested an explanation for increasing discrepancies between measured and predicted temperature fields forincreasing flux densities as reported by Ren et al. [2000]. Since
increasing flux tends to increase the difference in temperatures
between the upstream and downstream measurement locations,while increasing dispersivity causes decreasing temperature differences, the flux density cannot be accurately estimated without
considering dispersive heat flux. Thus we conclude that flux
density can be estimated only when thermal dispersivity is known
and vice versa. Since the analytical solution of Ren et al. [2000] did
not take thermal dispersion into account, their predicted water
fluxes were significantly underestimated. Although it is intuitively
expected that solute and thermal dispersivity values should be
similar, few if any studies have been conducted to determine
thermal dispersivity in unsaturated soils.
[43] To avoid the need for prior knowledge of dispersivity, wesuggest using an additional sensor located in the transverse
direction of the line heat source. From an additional sensitivity
analysis we evaluated the feasibility of such an approach for theestimation of thermal properties and water flow densities using
inverse modeling. We showed that temperatures collected at three
sensors (two parallel and one perpendicular to water flow) could
provide the required information for simultaneous and accurate
estimation of water flux, heat dispersivity, and soil thermal properties. If heat dispersivity values are of similar magnitude as thedispersivity defined for solute transport (as suggested by De
Marsily [1986]), the inverse modeling method presented herein
may be a valuable tool for estimation of solute transport dispersivity. Further experimental evidence, however, is needed tosupport this conclusion.
[44] Finally, we demonstrated the general application of inverse
modeling to estimate the soil thermal properties and their functional dependence on volumetric water content. In a subsequent
study we will present additional results, regarding the potential
application of inverse modeling to estimate solute transport properties in addition to the thermal and flow properties when combinedwith multineedle heat pulse probe measurements.

References
Bear, J., Dynamics of Fluids in Porous Media, Dover, Mineola, New York,1972.
Bristow, K. L., Measurement of thermal properties and water content of
unsaturated sandy soil using dual-probe heat-pulse probes, Agric. For.
Meteorol., 89, 75 – 84, 1998.
Bristow, K. L., G. S. Campbell, and C. Calissendorff, Test of a heat-pulse
probe for measuring changes in soil water content, Soil Sci. Soc. Am. J.,
57, 930 – 934, 1993.
Bristow, K. L., G. J. Kluitenberg, and R. Horton, Measurement of soil
thermal properties with a dual-probe heat-pulse technique, Soil Sci.
Soc. Am. J., 58, 1288 – 1294, 1994a.
Bristow, K. L., R. D. White, and G. J. Kluitenberg, Comparison of single
and dual probes for measuring soil thermal properties with transientheating, Aust. J. Soil Res., 32, 447 – 464, 1994b.
Bristow, K. L., G. J. Kluitenberg, C. J. Goding, and T. S. Fitzgerald, A
small multi-needle probe for measuring soil thermal properties, watercontent and electrical conductivity, Comput. Electron. Agric., 31, 265 –
280, 2001.
Campbell, G. S., Soil Physics with BASIC — Transport models for soilplant systems, 150 pp., Elsevier Sci., New York, 1985.

Campbell, G. S., C. Calissendorff, and J. H. Williams, Probe for measuring
soil specific heat using a heat pulse method, Soil Sci. Soc. Am. J., 55,
291 – 293, 1991.
Cass, A., G. S. Campbell, and T. L. Jones, Enhancement of thermal vapor
diffusion in soil, Soil Sci. Soc. Am. J., 48, 25 – 32, 1984.
Chung, S.-O., and R. Horton, Soil heat and water flow with a partial surface
mulch, Water Resour. Res., 23, 2175 – 2186, 1987.
De Marsily, G., Quantitative Hydrogeology: Groundwater Hydrology forEngineers, Academic, San Diego, Calif., 1986.
de Vries, D. A., Thermal properties of soils, in Physics of Plant Environment, edited by W. R. van Wijk, pp. 210 – 235, New York, North-Holland, 1963.
Gilmore, E., and J. S. Rogers, Heat units as a method of measuring maturity
in corn, Agron. J., 50, 611 – 615, 1958.
Hopmans, J. W., and J. H. Dane, Thermal conductivity of two porous media
as a function of water content, temperature, and density, Soil Sci., 142,
187 – 195, 1986a.
Hopmans, J. W., and J. H. Dane, Temperature dependence of soil hydraulicproperties, Soil Sci. Soc. Am. J., 50, 4 – 9, 1986b.
Hopmans, J. W., J. M. N. Hendrickx, and J. S. Selker, Emerging measurement techniques for vadose zone characterization, in Vadose Zone Hydrology: Cutting Across Disciplines, edited by M. B. Parlange and J. W.
Hopmans, pp. 279 – 316, Oxford Univ. Press, New York, 1999.
Hopmans, J. W., J. Šimunek, N. Romano, and W. Durner, Simultaneous
determination of water transmission and retention properties — Inverse
modeling of transient water flow, in Methods of Soil Analysis, Part I, 3rded., Monogr. 9, edited by G. C. Topp and J. H. Dane, in press, Am. Soc.
of Agron., Madison, Wis., 2002.
Inoue, M., J. Šimunek, S. Shiozawa, and J. W. Hopmans, Simultaneous
estimation of soil hydraulic and solute transport parameters from transient
infiltration experiments, Adv. Water Resour., 23, 677 – 688, 2000.
Kay, B. D., and J. B. Goit, Temperature dependent specific heats of dry soilmaterials, Can. Geotech. J., 12, 209 – 212, 1975.
Kluitenberg, G. J., J. M. Ham, and K. L. Bristow, Error analysis of the heatpulse method for measuring soil volumetric heat capacity, Soil Sci. Soc.
Am. J., 57, 1444 – 1451, 1993.
Kluitenberg, G. J., K. L. Bristow, and B. S. Das, Error analysis of the heatpulse method for measuring soil heat capacity, diffusivity and conductivity, Soil Sci. Soc. Am. J., 59, 719 – 726, 1995.
Nassar, I. N., and R. Horton, Simultaneous transfer of heat, water, andsolute in porous media, I, Theoretical development, Soil Sci. Soc. Am.
J., 56, 1350 – 1356, 1992.
Noborio, K., K. J. McInnes, and J. L. Heilman, Measurements of soil watercontent, heat capacity, and thermal conductivity with a single TDR probe,
Soil Sci., 161, 22 – 28, 1996.
Ren, T., K. Noborio, and R. Horton, Measuring soil water content, electricalconductivity and thermal properties with a thermo-TDR probe, Soil Sci.
Soc. Am. J., 63, 450 – 457, 1999.
Ren, T., G. J. Kluitenberg, and R. Horton, Determining soil water flux andpore water velocity by a heat pulse technique, Soil Sci. Soc. Am. J., 64,
552 – 560, 2000.
Shiozawa, S., and M. Inoue, Solute dispersion in unsaturated unit-gradientwater flow monitored in soil columns, Soil Sci. Soc. Am. J., in press,
2001.
Šimunek, J., M. Šejna, and M. T. van Genuchten, The HYDRUS-2D software package for simulating two-dimensional movement of water, heat,and multiple solutes in variably-saturated media, version 2.0, Rep.
IGWMC-TPS-53, 251 pp., Int. Ground Water Model. Cent., Colo. Sch.
of Mines, Golden, 1999.
Sophocleous, M., Analysis of water and heat flow in unsaturated-saturated
porous media, Water Resour. Res., 15, 1195 – 1206, 1979.
Tsang, Y. W., and K. Pruess, A study of thermally induced convection nears
a high-level nuclear waste repository in partially saturated fractured tuff,
Water Resour. Res., 23, 1958 – 1966, 1987.
Van Genuchten, M. T., A closed-form equation for predicting the hydraulicconductivity of unsaturated soils, Soil Sci. Soc. Am. J., 44, 892 – 898,
1980.
















































































































































K. L. Bristow, CSIRO Land and Water, PMB Aitkenvale, Townsville,
Queensland 4814, Australia. (Keith.Bristow@csiro.au)
J. W. Hopmans, Department of Land, Air and Water Resources, 123
Veihmeyer Hall, Davis, CA 95616-8628, USA. (jwhopmans@ucdavis.edu)
J. Šimunek George, E. Brown Jr. Salinity Laboratory, Agricultural
Research Service, U.S. Department of Agriculture, 450 West Springs Road,
Riverside, CA 92507-4617, USA. (JSIMUNEK@ussl.ars.usda.gov)