Advances in Water Resources 29 (2006) 250–267
www.elsevier.com/locate/advwatres

Multi-functional heat pulse probe measurements of coupledvadose zone ﬂow and transport
Annette P. Mortensen a,1, Jan W. Hopmans b,*, Yasushi Mori c, Jirka Šimůnek dabGeological Institute, Copenhagen University, Øster Voldgade 10, 1350 Copenhagen, Denmark
Hydrology Program, Department of Land, Air and Water Resources, University of California,
123 Veihmeyer Hall, Davis, CA 95616, USAcFaculty of Life and Environmental Science, Shimane University, Matsue 690-8504, Japand
Department of Environmental Sciences, University of California, Riverside, CA 92521, USA
Received 1 March 2005; accepted 14 March 2005
Available online 1 July 2005

Abstract
Simultaneous measurement of coupled water, heat, and solute transport in unsaturated porous media is made possible with themulti-functional heat pulse probe (MFHPP). The probe combines a heat pulse technique for estimating soil heat properties, waterﬂux, and water content with a Wenner array measurement of bulk soil electrical conductivity (ECbulk). To evaluate the MFHPP, we
conducted controlled steady-state ﬂow experiments in a sand column for a wide range of water saturations, ﬂow velocities, and solute concentrations. Flow and transport processes were monitored continuously using the MFHPP. Experimental data were analyzed
by inverse modeling of simultaneous water, heat, and solute transport using an adapted HYDRUS-2D model. Various optimization
scenarios yielded simultaneous estimation of thermal, solute, and hydraulic parameters and variables, including thermal conductivity, volumetric water content, water ﬂux, and thermal and solute dispersivities. We conclude that the MFHPP holds great promise as
an excellent instrument for the continuous monitoring and characterization of the vadose zone.
 2005 Elsevier Ltd. All rights reserved.
Keywords: Heat pulse probe; Unsaturated ﬂow; Heat transport; Solute transport; Vadose zone properties; Water content; Water ﬂux; Dispersion

1. Introduction
Multi-functional measurements of vadose zone processes have received increased attention over the past
few years. New sensors are being developed, by which
several already well-known measurement techniques
are combined into a single device. Examples of multifunctional instruments include a combined TDR andtensiometer probe [44], a combined TDR and heat pulse
probe [28], and a combined heat pulse probe with aWenner array [25], as used in the current study. Also*

Corresponding author. Tel.: +1 530 752 3060; fax: +1 530 752 5262.
E-mail address: jwhopmans@ucdavis.edu (J.W. Hopmans).
1
Now at Niras, Sortemosevej 2, 3450 Allerød, Denmark.
0309-1708/$ - see front matter  2005 Elsevier Ltd. All rights reserved.
doi:10.1016/j.advwatres.2005.03.017

multi-functional tracing techniques have been developedthat apply several tracers with diﬀerent properties simultaneously (e.g., [27]). The overall motivation for thedevelopment of these multi-functional techniques is toachieve an improved characterization of ﬂow and transport processes. Several beneﬁts are achieved by combining measurements. First, by measuring several
parameters at the same time and place, the coupling ofrelated transport properties are determined in concert,
thereby allowing examination of the nature of theirinterdependency, such as for the coupled transport ofwater and solute, and water and heat. Second, by using
the same instrument for various measurements withinapproximately the same measurement volume at about
the same time, the need to interpolate diﬀerent

A.P. Mortensen et al. / Advances in Water Resources 29 (2006) 250–267

measurement types in space and time is largely eliminated. Thirdly, simultaneous analysis of ﬂow and transport using combined soil measurements of watercontent, temperature, and solute concentration, decreases parameter uncertainty [24,37,39]. Thus, using
multi-functional measurement techniques allows determination of interdependent soil properties and processes, providing an improved understanding of
coupled ﬂow and transport.
The presented multi-functional heat pulse probe(MFHPP) originates from the dual-probe heat pulse(DPHP) method developed by Campbell et al. [11]. By
inducing a short heat pulse from one sensor needle
and measuring the temperature response at a second
sensor, the soil thermal properties (i.e., heat capacity,
C; thermal conductivity, k0; and thermal diﬀusivity, j)
and the water content, h, were estimated. This method
has been tested in both laboratory settings (e.g. [2,4–7]) and in ﬁeld soils [16,40]. Over the past few years,
the DPHP probe has been reﬁned and was developed
into various multi-sensor probes that are capable ofsimultaneous measuring a suite of soil properties. The
so-called thermo-TDR was developed by combining
the DPHP probe with TDR technology, to achieve aprobe that in addition to soil thermal properties also
estimates soil solute concentration from the simultaneous measurement of the bulk soil electrical conductivity (ECbulk) [28,30,32,34]. The potential of measuring
ECbulk using a modiﬁed heat pulse probe was shownby Bristow et al. [8], by including two extra sensor needles in the DPHP probe to create a so-called four-electrode Wenner array. Since ECbulk measurements are
dependent on both solute concentrations and water contents, ECbulk is an integral variable characterizing bothwater ﬂow and solute transport that can be beneﬁciallyused to simultaneously estimate soil hydraulic and solute transport parameters [21,37].
By inclusion of an extra temperature sensor to measure both the temperature responses at the upstream
and downstream end of the heater sensor, Ren et al.
[33] demonstrated the added beneﬁt of water ﬂux estimation. The resulting three rod heat pulse probe was
successfully used to measure water ﬂuxes in the range
of 1–3 m d1, but signiﬁcantly underestimated larger
ﬂuxes [33]. In their modeling and sensitivity study, Hopmans et al. [19] suggested that this underestimation
could be caused by the neglection of thermal dispersion
at these larger ﬂow velocities, and they showed that
more accurate water ﬂux estimates could be obtained,
if an extra transverse thermistor needle was includedin the HPP to take into account thermal dispersion.
The multi-functional heat pulse probe (MFHPP) of
the current study combines this type of heat sensor witha four-needle Wenner array. As demonstrated by Moriet al. [25,26], the MFHPP allow estimation of the soil
thermal properties (i.e., heat conduction, thermal con-251

ductivity, and thermal diﬀusivity), simultaneously withthe soil water properties (water content and water ﬂux).
Furthermore, the MFHPP estimates ECbulk from Wenner array measurements, from which the soilÕs soluteconcentration and water content can be determined.
In this study, the MFHPP technique was applied in
ﬂow and transport column experiments for the simultaneous and coupled measurement of water, heat, and solute transport variables and properties. Controlled
steady-state ﬂow experiments through variably-saturated sand were conducted for a range of water ﬂux,
water saturation, and solute concentration values. Heat
pulse and ECbulk measurements were analyzed with inverse modeling using a modiﬁed version of the HYDRUS-2D code [19,36]. The main objective of thepresented study was to evaluate the MFHPP as a means
to fully characterize coupled soil water ﬂow, heat, andsolute transport properties and processes.
2. Theory
Soil thermal and hydraulic properties were evaluatedfrom analysis of temperature responses of the thermistors of the MFHPP, solving for heat transport by bothconduction and convection. Solution of solute transport
is needed for interpretation of bulk soil electrical conductivity (ECbulk) measurements. The presented analysis
will solve the coupled heat, water, and solute transport
equations using HYDRUS-2D [36]. Although we willassume that both water and solute transport is onedimensional in our experiments, the two-dimensional
form of the heat ﬂow equation is required to accountfor heat transport in the vertical and lateral directions
between the sensors of the MFHPP [19].
2.1. Heat transport
Under the assumption of instantaneous heat transfer
between the solid, liquid, and gas phases in homogenous
porous media with constant uniform vertical water ﬂow,
heat transport is described by [3,19,36,38]:



 

oT
o
oT
o
oT
oT
¼
jxx
jzz
þ
 Vh
ð1Þ
ot
ox
ox
oz
oz
oz
where T is temperature (K), t is time (s), z is vertical position (m), and jzz and jxx denote the eﬀective thermal
diﬀusivities (m2 s1) in the z and x directions, respectively, deﬁned by:jzz ðhÞ ¼

k0 ðhÞ þ kd;L ðqw;z Þ;C bulk ðhÞ

jxx ðhÞ ¼

k0 ðhÞ þ kd;T ðqw;z ÞC bulk ðhÞ
ð2Þwhere k0 is the bulk soil thermal conductivity(W m1 K1), or the so-called stagnant thermal conductivity. kd,L = bLCwqw,z and kd,T = bTCwqw,z denote the252

A.P. Mortensen et al. / Advances in Water Resources 29 (2006) 250–267

longitudinal and transversal thermal dispersion coeﬃcient, respectively, where bL and bT are the longitudinal
and transverse heat dispersivity (m), and the water ﬂux
density, qw,z, is uniform and parallel to the z-axis.
Neglecting the heat capacity for air, the soil volumetricheat capacity, Cbulk (J m3 K1), can be described by[10]:
C bulk ¼ C s ð1  /Þ þ C w hð3Þwhere C = qc (J m3 K1) is the heat capacity and thesubscripts bulk, s, and w denote the bulk soil, solid,
and water phases, respectively, q is the density (kg m3),
c is the speciﬁc heat (J kg1 K1), h is the volumetricwater content (m3 m3), and / is the porosity(m3 m3). In analogy with solute transport, velocity
variations within the pore spaces will cause mixing ofpore water and thereby create dispersion-like spreadingof the temperature ﬁeld. By including hydrodynamic
dispersion on heat transport, the eﬀective thermal conductivity in the longitudinal direction can be written
as [19]:keff;zz ¼ k0 þ kd;Lð4ÞBecause of its order-of-magnitude smaller value, we
ignore the eﬀect of transverse dispersion. Although heatmixing occurs through all phases, whereas solute mixing
is limited to the water phase, values for heat and solute
dispersivities should be related [13]. The bulk soil thermal conductivity is relatively large, thus dispersive effects on keﬀ are only expected for large watervelocities. Hopmans et al. [19] deﬁned the KJJ-number
to evaluate the contribution of hydrodynamic thermal
dispersion to the eﬀective conductivity term:
KJJ ¼

kd;L
k0

C w qw;z hC w vw¼
C bulkC bulkð7Þ

describing heat ﬂow by the moving liquid phase, relative
to the stationary bulk porous medium [33]. The heatvelocity, Vh, lags behind the water front velocity vw,
since heat is assumed to be instantaneously transferred
between the solid, liquid, and gas phases giving thermal
homogeneity.
2.2. Water ﬂow
Steady-state variably-saturated water ﬂow in the z
direction is described by the Darcy ﬂux, qw,z (m s1):


ohm
þ1
ð8Þ
qw;z ¼ KðhÞoz
where K(hm) is the hydraulic conductivity function(m s1) and hm is the soil water matric head (m). The
two hydraulic relations, the soil water retention curve,h(hm), and the unsaturated hydraulic conductivity function, K(h), that are needed to solve Eq. (8), are here described by the van Genuchten and the van Genuchten–
Mualem relationships, respectively [43]:hðhm Þ  hr
n m
¼ ½1 þ ðajhm jÞ 
hs  hrKðhÞ ¼ K S 0.5 ½1  ð1  S 1=m Þm 2

S e ðhm Þ ¼

s e

e

ð9aÞ
ð9bÞ

where Se is the eﬀective saturation, hs and hr (m3 m3)
denote saturated and residual water contents, respectively, a (m1) and n are empirical constants,
m = 1  1/n, and Ks denotes the saturated hydraulic conductivity (m s1).
2.3. Solute transport

ð5Þ

At increasing pore-water velocities, hydrodynamic effects may become more important to heat transport,so that the need to include dispersion in the bulk thermal diﬀusivity will depend on the magnitude of the
KJJ-number. In their sensitivity analysis, Hopmans etal. [19] showed that no thermal dispersion correction isneeded if KJJ < 1.
The bulk soil thermal conductivity, k0, is a functionof mineral type, the geometrical arrangement of various
phases, and the water content [14]. Chung and Horton[12] proposed a three-parameter polynomial expression
to describe the water content dependence of soil thermal
conductivity. However, we found that a linear expression was adequate for the soil moisture range of ourexperiments, or:
k0 ¼ b0 þ b1 h

Vh ¼

ð6Þ

where b0 and b1 are empirical constants (W m1 K1).
In Eq. (1), Vh denotes the convective heat pulse velocity:

Assuming single domain transport for a conservative
tracer with steady-state water ﬂow, one-dimensional solute transport of our experiments is described by the conventional advection–dispersion equation:


oC
o
oC
ohC
h
¼
hD
ð10Þ
 vw
ot
oz
oz
oz
where C is the solute concentration (kg m3), D is theeﬀective solute dispersion coeﬃcient (m2 s1), andvw = qw,z/h is the average pore-water velocity (m s1).
For unsaturated porous media, hydrodynamic dispersion depends on both the water velocity and the watercontent [41]. Assuming that molecular diﬀusion is negligible, the hydrodynamic dispersion coeﬃcient reduces
to:
D ¼ aL ðhÞvw

ð11Þ

where aL is the longitudinal solute dispersivity (m).
Nützmann et al. [29] examined the dependency of dispersivity on the water content in unsaturated porous

A.P. Mortensen et al. / Advances in Water Resources 29 (2006) 250–267

media, and derived a linear relation between dispersivity
and relative water velocity ﬂuctuations. Their experimental work demonstrated that the relationship between
water velocity variations and water content was well described by a power function, to yield the following relation for solute dispersivity:aL ðhÞ ¼ a1 ha2

ð12Þ

where a1 (m) and a2 are constants, to be determined inthis study.

3. Materials and methods
Flow column experiments were performed to examine
steady-state ﬂow and transport for a range of water saturation, water velocity, and solute concentration values.
A newly-developed multi-functional heat pulse probe(MFHPP) by Mori et al. [25,26] was used to monitor
coupled heat, water, and solute transport. Inverse
modeling techniques were applied for analysis of theexperimental data, allowing simultaneous estimation ofsoil thermal properties (heat conductivity, heat dispersivity, and heat diﬀusivity) and solute properties (soluteconcentration and solute dispersion) that are coupledthrough the soil water properties of water content andwater ﬂux. A summary of this new MFHPP device will
be presented, followed by a detailed description of theexperimental column design and inverse modeling
analysis.
The MFHPP is constructed from six sensors (26 mm
long and 1.27-mm O.D. stainless steel needles), which
are incorporated into a 27-mm radius probe withapproximately 6 mm distance between the sensors(Fig. 1). The six sensors include a single central heater

253

sensor, four thermistor sensors located around the heater, and a four-needle Wenner array. Heat pulse experiments are conducted by generating an 8-s heat pulse at
the central heater sensor, corresponding to a heat ﬂux
of about 60 W m1, after which the temperature responses at the surrounding thermistor sensors aremeasured. A detailed description of the design, manufacturing, and implementation of the MFHPP was presented in [25]. Since the heat dissipation depends on thesoilÕs thermal properties, analysis of the temperature
signals provides estimates for the soilÕs thermal diﬀusivity and conductivity, water content, and water ﬂux
density. Various analytical solutions have been developed for heat pulse analysis (see for example [25,33,45]) to successfully estimate these thermal properties
and ﬂow variables. In this study, we take it a step further, and will employ inverse numerical modeling techniques for the simultaneous analysis of coupled ﬂow
and transport in an experimental column, as describedin detail hereafter.
3.1. Calibration of the MFHPP
Temperature signals were found to be highly sensitive
to sensor spacings, so that calibration of each thermistor
probe was necessary to achieve satisfactory measurements [22]. Analytical solutions of the heat transport
equation were used to calibrate the sensor distance, r,
between each thermistor and the heater sensor. The
solution describes the temperature increase, DT (K),for an inﬁnite line heat source in a homogeneous andisotropic medium after a heat pulse of duration t0 (s)at a distance r (m) at time t > t0 [6,14,22]:  2 
q0
r2
r
DT ðr; tÞ ¼
Ei
 Ei
; t > t0
4pCj
4jðt  t0 Þ
4jt
ð13Þ

Fig. 1. Schematic of the MFHPP design, with heat pulse measurement
using the central heater sensor and the four surrounding thermistors(A) and EC measurements using the horizontal four-needle Wennerarray (B).

where q 0 is the energy input per unit length of heater
per unit time (W m1), and Ei(x) denotes the exponential integral operator with argument x [1]. By applying a heat pulse in a media with known heat capacity,
C, the distance between the heater and the thermistor,r, can be estimated by ﬁtting the measured temperature
response to this analytical solution. Traditionally, anagar–agar solution was used for this calibration since
the heat capacity for water is known. Alternatively, insitu calibration was suggested by Mori et al. [25], usinga priori information on the heat capacity for the fully
saturated media, thereby accounting for changing
needle-spacings after insertion into the soil. We followed this approach in our study as well, and repeatedit ten times, during the total experimental period, for
no-ﬂow conditions at full saturation. Measured temperature signals were ﬁtted to Eq. (13) using Solver, Excel
[46], by ﬁtting the sensor spacing, r, and the thermal

254

A.P. Mortensen et al. / Advances in Water Resources 29 (2006) 250–267

neglecting the soil surface conductance of the sandy soil,
ECbulk is given by:

Table 1
Physical properties of washed Tottori Dune sandBulk density
Saturated hydraulic conductivity
Saturated water content
Residual water content
van Genuchten parametervan Genuchten parameter
Speciﬁc heatqb (kg m3)
Ks (m s1)hs (m3 m3)
hr (m3 m3)
a (cm1)ncs (J kg1 K1)

1.63
2.4 · 1040.371
0.05350.028812.18
795.0

ECbulk ¼ c1 ECw h2 þ c2 ECw h

where c1 and c2 are empirical parameters. Traditionally,
the ﬁrst step in Wenner array calibration is to determine
the relationship between measured probe-speciﬁc electrical resistance and ECbulk (mS cm1), using a range of
CaCl2 solution concentrations. However, we measured
probe-speciﬁc electrical resistance directly to determine
the relationship between ECbulk and ECw, both expressed in V V1, as resistance was measured as the voltage drop relative to the voltage across a reference
resistor. For each single probe, ECbulk was measured
at various known chloride concentrations for a rangeof soil water content values (Fig. 2). Subsequently, solute concentration C (mol L1) was computed from
ECw, using a linear calibration relationship that relates
probe-speciﬁc electrical resistance to solute concentration. The ﬁnal result was four diﬀerent calibration equations for each of the Wenner arrays, with RMSE-values
ranging between 21 and 58 V V1, corresponding withchloride concentration uncertainties of 0.0005–
0.0017 mol L1, respectively. We note that the experimental range of water content values varied between
probes. Speciﬁcally, MFHPP5 (top of experimental col-diﬀusivity, j, using a priori information on the heat
capacity, C, for the sand (Table 1). Using 10 datasets,
average r-values of each sensor distance and thermal
diﬀusivity were estimated with an average coeﬃcientof variation of 0.34% and 1.30%, respectively. For the
total of four MFHPPÕs, average sensor spacing was
5.91 mm with r-values ranging between 5.52 mm and6.14 mm.
The four horizontal sensors of the MFHPP (Fig. 1B)
were used as a four-electrode Wenner array sensor forbulk soil electrical conductivity (ECbulk) measurements[21,25]. The ECbulk value depends on the solution electrical conductivity, ECw (mS cm1), the soil surface conductivity, ECs (mS cm1), soil water content, h, soil bulk
density, q, temperature, T, and sensor geometry. Assuming the Rhoades et al. [35] relationship to be valid and
500

500

MFHPP2

MFHPP3

C=0.01 M

400

C=0.01 M

400

C=0.03 M

C=0.03 M

ECbulk [V V-1]

C=0.06 M

C=0.06 M

C=0.08 M

300

C=0.08 M

300

C=0.1 M

C=0.1 M

Model c1=1.68 c2=-0.026

Model c1=1.86 c2=-0.11

200

200

100

100

0

0
0

0.1

0.2

0.3

0.4

0

500

500
C=0.01 M

400

0.3

0.4

0.3

0.4

C=0.03 M

C=0.06 M

C=0.06 M

C=0.08 M

300

0.2

C=0.01 M

400

C=0.03 M

0.1
MFHPP5

MFHPP4ECbulk [V V-1]

ð14Þ

C=0.08 M

300

C=0.1 M

C=0.1 M

Model c1=1.23 c2=0.16

200

200

100

100

0

Model c1=1.46 c2=0.022

0
0

0.1

0.2

Water content θ [m3 m-3]

0.3

0.4

0

0.1

0.2

Water content θ [m3 m-3]

Fig. 2. Calibration of the four-needle Wenner array for the four MFHPPs measuring ECbulk as a function of water content and solute concentration.
The model given in Eq. (14) was ﬁtted to the data for each probe with the parameter values c1 and c2.

A.P. Mortensen et al. / Advances in Water Resources 29 (2006) 250–267

255

Fig. 3. Schematic of the experimental ﬂow column. The relative positions of the MFHPPÕs and tensiometers are presented in the cross-sections(right): A–A at 4 cm depth, B–B at 10 cm depth, C–C at 16 cm depth, and D–D at 22 cm depth.

Table 2Flux and suction boundary values for the saturated and unsaturated experiments
Experiment

Flux Q (mL min1)

Velocity qw,z (m d1)

Boundary pressure h (m)–I. Saturated with no ﬂow

0.0

0.0

II. Saturated with ﬂow

63
46
32
15
3.3

17.84
13.39
9.37
4.37
1.05

III. Unsaturated with no ﬂow

0.0

0.0

IV. Unsaturated with ﬂow

20
10
5
2.5
0.6

6.25
3.12
1.56
0.78
0.19

C1 (mol L1)

C2 (mol L1)

–

–

0.1
0.01
0.01
0.01
0.01

0.01
0.1
0.1
0.1
0.1

0.2

–

–

0.2
0.2
0.2
0.2
0.2

0.1
0.1
0.1
0.06
0.1

0.01
0.01
0.06
0.01
0.01

0.0
0.07
0.14
0.265
0.21

Several tracer experiments were performed for each ﬂow experiment; however the ones used for the inverse modeling are presented with initialconcentration C1 and pulse concentration C2.

umn) and MFHPP2 (bottom of column, Fig. 3) corresponded with the largest and smallest water content
ranges, respectively. Rather than having independent
water content values available for the calibration, values
were estimated from MFHPP measurements for theunsaturated experiments with known chloride concentrations, to be discussed later (Experiments III and IV
in Table 2).

3.2. Experimental setup
Steady-state water ﬂow experiments were conductedin a 7.94-cm inner diameter and 30 cm long Plexiglas
column packed with 28 cm of Tottori Dune sand (Fig.
3). This Tottori Dune sand was extensively studied inprevious investigations [21,25,26], and was selected forthis study since its unsaturated hydraulic properties

256

A.P. Mortensen et al. / Advances in Water Resources 29 (2006) 250–267

allowed for a relatively large range of measurable waterﬂux density values across a wide water content range.
The sand was carefully washed before use with a detergent solution to remove ﬁne clay and silt materials that
could potentially clog the porous membrane at the bottom of the column. The retention and unsaturated conductivity characteristics of the washed sand were
measured using multi-step outﬂow experiments, independently [25], whereas the saturated hydraulic conductivity was determined from constant head permeametermeasurements. Optimized van Genuchten parameters,
as deﬁned in Eqs. (9), were estimated using the SFOPT
code [42]. Optimized hydraulic parameters along withthe measured saturated hydraulic conductivity are presented in Table 1, whereas the combined soil hydraulicfunctions are shown in Fig. 4A. The speciﬁc heat of

1.E-04

0.6

Matric head hm [m]1.E-050.5
1.E-060.4

1.E-070.3

1.E-08
1.E-090.2
1.E-10

0.1

Hydraulic conductivity K (h) [m s-1]

1.E-03A 0.7

1.E-11
1.E-12
0.4

0
0

0.1

0.2

0.3

Water content θ [m3 m-3]
B

0.0
0

Column depth [cm]

4

0.1

0.2

0.3

0.4
qw,z=0qw,z=0.19

MFHPP5qw,z=0.78
qw,z=1.56

10

MFHPP4qw,z=3.12qw,z=6.25

16

MFHPP322
MFHPP228

Fig. 4. (A) Retention curve and hydraulic conductivity function
measured for the Tottori Dune sand by multi-step outﬂow experiment.
(B) Calculated steady-state distributions of the volumetric watercontent in the 28 cm sand column for the diﬀerent inﬁltration rates qw,z(Table 2). The calculations are based on measured soil hydrauliccharacteristics (Table 1) and a lower boundary suction of 20 cm. The
four MFHPPs were located at 4, 10, 16, and 22 cm depth.

the sand, cs (Table 1), was measured by diﬀerential scanning calorimetry as described in [25]. Prior to the ﬂow
experiments, the sand was dry-packed in the column ata bulk density of 1.63 g cm3 to a total height of
28 cm. The same sand pack was used for all presented
experiments to eliminate diﬀerences due to changing
heterogeneities. In order to attain full initial saturation
prior to each experiment, the sand column was ﬂushedwith several pore volumes of CO2 before saturationfrom below with de-aired water.
Flow through the sand column was controlled
through application of a range of steady-state water ﬂux
values, qw,z, between 0 and 17.84 m d1, with corresponding variable suctions at the bottom boundary.
Assuming spatial uniformity, these steady-state water
ﬂuxes were achieved by an artiﬁcial 12-needle rainmaking device. A microtube pump (Masterﬂex L/S variablespeed digital drives, Cole-Parmer Instrument Company,
Vernon Hills, IL) provided for the necessary ﬂow,
whereas a six-piece manifold (Fig. 3) connecting to 12
syringe needles provided for a uniform water-sprinkling
application to the sand surface. The bottom of the sandcolumn included a metal screen, on which a 20-lm porous nylon membrane (Osmonics R22SP14225, GE
Osmonics Labstore, Minnetonka, MN) was glued, witha protective wet strength ﬁlter (Whatman Qualitative
Wetstrength Filter Paper 114, Whatman Inc., Clifton,
NJ) in between the screen and the nylon. The membrane
assured hydraulic continuity between the column andthe drainage outlet. The magnitude of the appliedboundary suction was controlled by adjusting the outﬂow height of a tube connected to the bottom outlet.
Tracer solutions were mixed prior to the experimentsby adding CaCl2 to distilled water, achieving a range ofconcentrations of 0.01, 0.03, 0.06, 0.08, and 0.1 mol L1.
The tracer solutions were allowed to equilibrate with thetemperature of the laboratory, which was maintained
constant within the interval of 18–20 C. The column
was covered with a black fabric to protect it from light,
thereby minimizing bacterial growth in the sand.
The various column experiments were monitoredwith four MFHPPÕs and four miniature tensiometers(Fig. 3). All instruments were inserted into the Plexiglas
column before packing the sand to ensure good andevenly contact. The MFHPPÕs were inserted such thatthe four-needle Wenner array was positioned horizontally. The vertical distance between the MFHPPÕs was
6 cm, with an oﬀset of a 60 angle in between, to minimize ﬂow disturbance along the vertical direction. A
CR10 data logger and three AM416 multiplexers(Campbell Scientiﬁc) were used to control operation of
the MFHPPÕs, with separate datalogger programs forthe temperature and EC measurements. A series of heatpulse measurements was repeated every 15 min. Each 8-sheat pulse was followed by a 120-s period of temperature measurements at the four thermistor sensors.

A.P. Mortensen et al. / Advances in Water Resources 29 (2006) 250–267

Immediately before the heating, the background temperature of each thermistor sensor was measured and used
for calculating the exact temperature increase at eachsensor. Separately, during the tracer experiments, ECbulk
was measured continuously at each MFHPP with a time
resolution of 10 s.
The soil water matrix potential was measured withminiature tensiometers with a 1 cm long and 0.635-cm
O.D. ceramic cup [15]. Four tensiometers were inserted
approximately 2 cm into the column, at the same depth
as the MFHPPÕs. For the top and bottom MFHPPÕs, thetensiometers were positioned at an 180 angle with the
MFHPP, whereas for the two middle MFHPPÕs, theirangle with the tensiometers was 150 (Fig. 3). The tensiometers were read every 10 s by 1-psi transducers using aseparate CR21X datalogger (Campbell Scientiﬁc). However, as also pointed out in [25], h values are extremely
sensitive to small measurement errors in hm for thiscoarse-textured sandy material (see Fig. 4A). Therefore,
we could not use the tensiometric data in the
optimizations.
3.3. Experimental matrix
Four diﬀerent categories of experiments were conducted on the sand column: I. Saturated experiments
without ﬂow; II. Saturated experiments with steadystate ﬂow; III. Unsaturated experiments without ﬂow;
and IV. Unsaturated experiments with steady-state ﬂow(Table 2). Five diﬀerent saturated ﬂow experiments wereconducted, at which the outﬂow tube was adjusted tomaintain full saturation of the sand column, while preventing water ponding. Also ﬁve diﬀerent unsaturated
steady-state ﬂow experiments were performed, all witha constant suction of 20 cm at the lower boundary.
For these experiments, the sand was at near-saturation
at the bottom of the column and water content decreased upwards with values being a function of the applied ﬂux. The various water content proﬁles arepresented in Fig. 4B, using the soil hydraulic function
parameters of Table 1.
Heat pulse measurements were conducted for eachsingle experiment for estimation of thermal and hydraulic properties with a minimum of two measurements foreach ﬂow condition and solute concentration. Tracer
experiments were conducted for each of the ﬂow experiments by applying CaCl2 as solute tracer using diﬀerent
combinations of the ﬁve tracer solutions. The tracer
experiments were performed after initial establishmentof a constant initial solute concentration, C1 across thecolumn, switching instantaneously to the invading solution (C2) at time t = 0, while maintaining the same ﬂow
rate. This was achieved by connecting two containers
with their respective concentrations to the same pump,and simply changing from one to the other concentration by turning a valve (Fig. 3). The transport of the sol-

257

utes was measured with the four MFHPPÕs providing
tracer breakthrough curves at their four locations in
the column. For each ﬂow condition, a minimum of
three tracer experiments were conducted by using diﬀerent combinations of the ﬁve tracer solutions. The experiments used for the inverse modeling are listed inTable 2.
3.4. Inverse modeling
Inverse modeling techniques has successfully been
applied for a wide range of vadose zone problems [20].
Heat pulse signals have traditionally been analyzed using
analytical solutions, both without and with convective
heat transport [25,26,33,45]. The opportunities of using
inverse modeling techniques for heat pulse analysis were
described by Hopmans et al. [19]. Most importantly,
inverse modeling allows for the simultaneous estimationof coupled ﬂow and transport parameters. Furthermore,both the number and the type of measurements to be
included in the optimization can be varied.
We used the HYDRUS-2D ﬁnite element code by Šimůnek et al. [36] for analysis of the MFHPP experiments. The code solves simultaneously water ﬂow withheat and solute transport in two spatial dimensions (xand z). The code solves numerically RichardsÕ equation
for saturated–unsaturated water ﬂow and convection–
dispersion equations for heat and solute transport. It
considers heat transport by conduction, convection,and dispersion, and solute transport by Fickian-based
advection–dispersion. HYDRUS-2D uses the Levenberg–Marquardt (LM) algorithm for parameter optimization by inverse modeling. The experiments were
simulated for a 3.97 cm (horizontal) by 28 cm (vertical)
transport domain, corresponding to the radius and theheight of the sand column. A ﬁnite element grid was
generated so that nodes were positioned at exact locations of the thermistor and EC probe for all four
MFHPPÕs. For the purpose of the numerical analysis,
all sensors were assumed to extend inﬁnitely in the y
direction, perpendicular to the x–z plane of the simulated domain, thereby ignoring possible three-dimensional eﬀects. The simulated x–z domain consisted of
9855 elements, with a ﬁne nodal spacing of 0.34 mm
near the sensor needles, increasing to a coarser resolution of 4.59 mm near the domain boundaries. Depending on the type of experiment, top and bottom ﬂow
boundary conditions were deﬁned by water potential
or constant ﬂux values and solute transport by a thirdtype boundary condition. A zero-order source term
was added to Eq. (1), deﬁning a non-zero heating rate
for the nodes representing the heater needles equal to
the speciﬁc heat ﬂux measured in the experiment. On
all lateral boundaries of the sand domain, a zero water,
solute and heat ﬂux was deﬁned. Thermal properties
were deﬁned as presented in Eqs. (12), (3), and (6),258

A.P. Mortensen et al. / Advances in Water Resources 29 (2006) 250–267

and hydraulic properties as in Eq. (9). Although we will
demonstrate that hr can be estimated concurrent with allother soil hydraulic parameters, we chose to ﬁx it to the
independently-estimated value of 0.0535 (Table 1) for
most optimizations, since it will generally lead to moreunique solutions. Temperature eﬀects on soil hydraulicproperties [17] and thermal properties [18] were ignored.
Moreover, distillation eﬀects causing enhanced heattransport by latent heat through vaporization and subsequent condensation [9] was not considered at this time.
The dependency of ECbulk on solute concentration andwater content was included in the model by incorporating the four calibration equations found for each single
Wenner array.
For the heat pulse experiments, the objective function
U to be minimized during the parameter estimation wasdeﬁned as followsUðpHP Þ ¼ W 1

N1X2½T  ðx; ti Þ  T ðx; ti ; pHP Þ

ð15Þ

ature peak. For the saturated experiments, pHP included
the bulk soil volumetric heat capacity, Cbulk, and thermal conductivity, k0 and also water ﬂux density, qw,z,and heat dispersivity, bL, for the ﬂow experiments. For
the unsaturated experiments, the parameter optimization vector pHP was extended to also include volumetricwater content, h (for the no-ﬂow experiments optimized
through Cbulk). For the multiple probe optimizations,the objective function included data from all four
MFHPPÕs (in total 1776 data points) with the same
weights provided as for single MFHPP optimization.
Since the multiple probe optimizations provided column-average parameter values, the pHP for the unsaturated experiments included functional relationships
for thermal conductivity (b0 and b1) and unsaturated hydraulic properties functions (a, n, Ks, and hr)(Table 3).
For the tracer experiments, the objective function to
be minimized was deﬁned byi¼1

where the right hand size represents the residuals between measured temperatures, T*, and corresponding
predicted temperatures, T. The vector x denotes the spatial coordinate of each measurement i, N1 the total number of temperature measurements, W1 is the weightassociated with a particular measurement data point,
and the vector pHP contains the optimized parameters.
Two diﬀerent optimization techniques were applied
depending on whether the temperature signals from asingle MFHPP or from all four MFHPPÕs were used
to optimize parameters. The number and type of optimized parameters depend on experiment type (Table3). For the single MFHPP optimizations, the objectivefunction contained temperature data (total of 444) for
each of the four thermistors sensors, with measurements
each second for a duration of 120 s after heating started.
Improved ﬁts to the temperature response data wereachieved by introducing a value of 5 to the weighting
factor (W1) of the 16 measurements around the temper-UðpEC Þ ¼ W 2

N2X½ECbulk ðx; ti Þ  ECbulk ðx; ti ; pEC Þ

2

ð16Þ

i¼1

where ECbulk and ECbulk denote the measured and predicted bulk soil EC, respectively, N2 deﬁnes the numberof EC measurements, W2 is the weight associated with aparticular measurement data point, and pEC contains the
optimized parameters. For both saturated and unsaturated experiments, the objective function contained EC
measurements from all four MFHPPs. The number
and type of optimized parameters were however diﬀerent for the two experimental groups (Table 3). For saturated ﬂow, optimized parameters included water ﬂux
density, qw,z, and solute dispersivity, aL. For the unsaturated tracer experiments, pEC contained the functionaldescription of the solute dispersivity (i.e. parameters a1
and a2 as deﬁned in Eq. (12)), the hydraulic properties(a, n, and Ks) and the water ﬂux density, qw,z.
Finally a coupling of both heat pulse and ECbulk data
was tested by combining Eqs. (15) and (16) into a single

Table 3
Optimized parameters and variables for options A, B, C, and D: soil volumetric heat capacity, Cbulk; thermal conductivity, k0; thermal conductivity
parameters b0 and b1, water ﬂux, qw,z; water content, h; residual water content, hr; retention parameters a and n; saturated hydraulic conductivity Ks;
thermal dispersivity, bL; solute dispersivity, aL; solute dispersivity parameters a1 and a2
Experiment

Saturatedwith no ﬂow

Saturatedwith ﬂow

Unsaturated
with no ﬂow

Unsaturated
with ﬂow

A. Heat pulse from single MFHPP

k0 Cbulk

k0 bL
qw,z

k0
h

k0 bL
qw,z hB. Heat pulse from multiple MFHPPÕs

k0 Cbulk

k0
qw,z

b0 b1
a n hrb0 b1
a n hr Ks qw,z

C. Solute from multiple MFHPPÕs

–

aL
qw,z

–

a1 a2
a n Ks qw,z

D. Heat pulse and solute from multiple MFHPPÕs

–

k0 aL
qw,z

–

b0 b1
a1 a2
a n Ks qw,z

A.P. Mortensen et al. / Advances in Water Resources 29 (2006) 250–267

objective function to simultaneously estimate all ﬂow
and transport parameters and variables of Table 3.
The objective function included both temperature andECbulk measurements for all four MFHPPÕs. Internal
weighting was accomplished by including the numberof measurements and variance in W1 and W2, thereby
achieving equal weight between the diﬀerent data types,
as recommended in [20].

4. Results and discussion
Measured temperature diﬀerences of all four thermistors for each MFHPP device were used to estimate bothsoil thermal and hydraulic parameters. As an example,
we present a comparison of measured with optimized
temperature signals for unsaturated experiments atzero-ﬂow and at steady-state ﬂow of qw,z = 6.25 m d1
in Fig. 5. In the absence of water ﬂow, heat is transported by conduction only, resulting in symmetricalspreading of the heat pulse around the heater needle
to all four surrounding thermistors equally, with slightvariations caused by variable sensor spacing and soil
heterogeneities. In contrast, asymmetrical temperature
responses occur for the steady-state ﬂow water experiments with heat convection, resulting in temperature differences between the upstream, downstream, andtransverse thermistor locations [19]. However, in either
case, the temperature response will increase withdecreasing soil water content. EC measurements wereanalyzed to estimate solute dispersivity and water ﬂux,
with and without simultaneous temperature response
measurements. Hereafter, we will discuss the optimization results of four diﬀerent options, A through D, of

259

which the experimental settings are listed in the ﬁrst column of Table 3.
4.1. Thermal and hydraulic properties from single
MFHPP measurements
The single probe optimizations of both unsaturated
zero-ﬂow and steady-state experiments included estimations of the volumetric water content, h. These estimatedwater content values are compared with the corresponding simulated water content values (from hydraulicparameters, Table 1) in Fig. 6, with a RMSE of
0.019 m3 m3. The relative water content errors, as deﬁned by (hest  htrue)/htrue, were on average smaller than
±10% between all MFHPPÕs, and were unbiased. In our
simulations with HYDRUS-2D, temperature data for
all four thermistors were included to estimate averagewater content. Yet, signiﬁcant diﬀerences in water contents with vertical position might occur, because of thelarge sensitivity of the sand soil moisture to the matrichead (Fig. 4A). In part, this explains slightly larger errors in the intermediate water content range of Fig. 6,
as diﬀerences in h between the upstream and downstream thermistors could be as large as 0.03 m3 m3for some ﬂow experiments.
Mori et al. [25,26] also used the MFHPP techniquefor water content estimation and applied the analyticalsolution of heat transport equation using the horizontal
thermistor only, assuming that this thermistor is least affected by convective transport. Although they foundexcellent agreement for water ﬂuxes smaller than
0.5 m d1 (RMSE of 0.0056 m3 m3), estimated water

1.4
Upstream
Downstream
Transverse
Zero-flux
Optimized

Temperature difference [C]

1.2
1.0
0.8
0.6
0.4

Simulated water content [m3 m-3]

0.4

0.3

0.2

0.1

0.2

Steady-state flow
No flow

0.0
0

20

40

60

80

100

120

Time [s]

0.0
0.0

θr

0.1

0.2

0.3

θs

0.4

Estimated water content from MFHPP [m3 m-3]

Fig. 5. Example of measured and estimated temperature responses at
upstream, downstream, and transverse thermistor sensors for unsaturated experiments with no ﬂow (h = 0.147; qwz = 0 m d1) and steadystate ﬂow (h = 0.213; qwz = 6.25 m d1).

Fig. 6. Volumetric water contents estimated from single MFHPP
optimizations (option A) compared with the simulated water contents(calculated from retention characteristics in Table 2).

260

A.P. Mortensen et al. / Advances in Water Resources 29 (2006) 250–267

content values were signiﬁcantly overestimated for higher water ﬂuxes, as their analytical solution neglected
convective heat transport. Also Kluitenberg and Heitman [23] examined the inﬂuence of water ﬂux density
and thermistor orientation on water content estimation
errors using analytical solutions, and concluded thatwater content errors increased with increasing water
ﬂux.
Simultaneously with soil water content, the soil thermal conductivity, k0, was estimated from both saturated
and unsaturated experiments for both zero-ﬂux and aseries of steady-state ﬂow conditions (Table 2). Estimated thermal conductivities for all experiments combined are presented in Fig. 7 as a function of the soil
water content, as estimated from the MFHPP measurements. The linear relation of Eq. (6) was ﬁtted to thedata, thereby characterizing the water content dependency of thermal conductivity of the Tottori sand as afunction of water content, resulting in values ofb0 = 1.13 and b1 = 1.68 (W m1 K1), with a RMSE
value of 0.003 W m1 K1. Our data compared remarkably well with the independent sandy soil data of Hopmans and Dane [17], as was also found from previous
work by Mori et al. [25,26], who used the MFHPP technique in combination with analytical solutions. Yet,Mori et al. [26] concluded that the MFHPP measurements overestimated the thermal conductivity for water
ﬂuxes larger than 2 m d1. Their analytical solutions
necessitated the ﬁtting to the horizontal thermistor measurements only, thereby largely neglecting the inﬂuenceof convection. In our case, we used all four thermistors
simultaneously in the numerical model, allowing for the

coupling between thermal transport and water ﬂow.
Therefore, the ﬁnal outcome in Fig. 7 was not aﬀectedby the magnitude of the water ﬂux density. Note that
we do not attempt to estimate the thermal conductivities
for water content values smaller than about 0.1, since
this is below our measurement range.
In addition to water content and thermal properties,
the single probe optimizations also yielded estimatedwater ﬂux densities. Estimated water ﬂuxes are compared with the true ﬂuxes, as computed from volumetric drainage ﬂow rates in Fig. 8. More accurate ﬂux
estimates were generally obtained for the unsaturated
ﬂow experiments (RMSE of 0.49 m d1) as comparedto the saturated experiments (RMSE of 1.94 m d1),
with the mostly higher ﬂow rates. Relative ﬂux errors,
as deﬁned by (qwz,est  qwz,true)/qwz,true, were calculated
for each MFHPP for both the saturated and unsaturated ﬂow. For the saturated ﬂow experiments, water
ﬂuxes were generally underestimated, with relative
errors between 10% and 20%; however, the errors of
MFHPP2 at the column bottom were signiﬁcantly larger(13–26%). For the unsaturated steady-state ﬂow experiments, the MFHPP optimizations generally underestimated the ﬂow at high water ﬂuxes as well, but
relative errors were positive for the low ﬂow rate experiments, with MFHPP2 showing the largest relative
errors. In contrast to the ﬂux rate-independent errors
for the saturated experiments, the relative ﬂux error increased as water ﬂux decreased for the unsaturated
experiments. We note that the lower water ﬂuxes correspond with smaller water content values in the column(Fig. 4B).

2.0

1.6

1:1 line15

1.2

True flux [m d-1]

Thermal conductivity λ [Wm-1 K -1]

200.8
MFHPP saturated with no flow

10

MFHPP saturated with flow
MFHPP unsaturated with no flow

0.4

5

MFHPP unsaturated with flow
. 0+b1θλ=bSaturated with flow

Hopmans and Dane [1986]

0.0
0.0

Unsaturated with flow

0.1

0.2

0.3

0.4

Estimated water content from MFHPP [m3 m-3]

0
0

5

10

15

20

-1

Estimated flux from MFHPP [m d ]Fig. 7. Thermal conductivities estimated from single MFHPP optimizations (option A) as a function of the estimated water content fromsingle MFHPP.

Fig. 8. Water ﬂuxes estimated from single MFHPP optimizations(option A) compared with the true ﬂuxes applied as top boundary.

–
NO
3.27
1.70
0.94
0.19
–
1.134
0.822
1.11
1.16
1.15
1.04
1.05
qw,z = 6.25 m d1
qw,z = 3.12 m d1
qw,z = 1.56 m d1
qw,z = 0.78 m d1
qw,z = 0.19 m d1
qw,z = 0 m d1
–
15.05
10.52
7.79
3.72
0.90
–
1.80
1.82
1.75
1.78
1.80
1.79
1.80
Single probe/multi-step outﬂow
qw,z = 17.84 m d1
qw,z = 13.39 m d1
qw,z = 9.37 m d1
qw,z = 4.37 m d1
qw,z = 1.05 m d1
qw,z = 0 m d1

For comparison, thermal properties from single probe optimization and hydraulic properties from multi-step outﬂow experiments are included in the ﬁrst row.
NO: not optimized. Solutions did not converge.

1.68
2.77
1.97
1.63
1.63
3.28
1.91

b1 (W mK )b0 (W mqw,z (m d )K )1

k0 (W m261

2.4E42.5E4
1.8E4
1.6E4
1.9E4
1.4E4–0.05350.000240.00150.000180.03950.024
0.0650
12.18
17.93
19.93
13.32
13.04
13.55
19.23
0.02880.02920.03110.03050.03010.02900.0278Ks (m s1)hr (m3 m3)na (cm1)
K )Hydraulic properties

1
1
1
1

Thermal properties

1
1

Unsaturated
experiments
Hydraulic properties
Thermal properties

For the multiple probe optimizations the temperature
diﬀerences for all four MFHPPÕs, for a total of 16 temperature signals, were included into the objective function. For the saturated experiments, this resulted incolumn-average values for the saturated thermal conductivity and water ﬂux (Table 3). For the unsaturated
experiments, the simultaneous use of all MFHPP data
provided optimized parameter values for the water content dependent thermal conductivity function in Eq. (6),a, n and Ks values of the soil hydraulic functions (Eq.
(9)), and the water ﬂux density, qw,z (Table 3).
The results of the multiple probe optimizations is presented in Table 4 and Fig. 9, whereas the optimized
water ﬂux density values are included in Fig. 11. For
the saturated conditions, the optimized thermal conductivity and water ﬂux values agreed with corresponding
average values for the single probe optimizations. For
the unsaturated experiments, each steady-state ﬂow creates water content proﬁles that are determined by waterﬂux density (Fig. 4B), so that optimization results apply
to each speciﬁc water content range only. Speciﬁcally, atSaturated experiments

4.2. Thermal and hydraulic properties from multiple
MFHPP measurements

Table 4
Option B. Optimized thermal and hydraulic properties for multiple probe optimizations of unsaturated experiments

Various other studies also concluded that HPP measurements underestimated water ﬂux at high ﬂow rates[19,25,26,33]. However, both studies of Mori et al.
[25,26] point out the diﬃculty in estimating water ﬂuxes
if values are smaller than 0.1 m d1. It was suggested
that the limitation of the MFHPP in the low water ﬂux
range is controlled by the temperature resolution of thethermistors (0.01 C).
The underestimation in the high water ﬂux range wasattributed by Hopmans et al. [19] to thermal dispersion.
Therefore, in the single probe optimizations, we included the longitudinal heat dispersivity (Eq. (4)), bL,
as an additional estimation parameter in both saturated
and unsaturated ﬂow experiments. However, the estimated values were generally small, with average valuesof 0.004 cm (CV = 1.7 · 105) for saturated experiments
and 0.006 cm (CV = 9.75 · 105) for unsaturated experiments. These low values suggest that thermal dispersion
is not important in our sand, as is to be expected withaverage KJJ-numbers (Eq. (5)) of 0.01 (CV =2.3 · 104) and 0.008 (CV = 3.1 · 104) for the saturated and unsaturated experiments, respectively. These
values are considerably lower than the threshold valueof KJJ = 1, for which heat dispersion becomes important [19]. Optimizations that excluded thermal dispersion conﬁrmed that this process was not importantfor our experiments. Moreover, it is expected that dispersivity values will be small if determined from measurements over such small travel distances (6 mm)
between the heater and the upstream and downstream
thermistors.

qw,z (m d1)

A.P. Mortensen et al. / Advances in Water Resources 29 (2006) 250–267

262

A.P. Mortensen et al. / Advances in Water Resources 29 (2006) 250–267
2.0

Thermal conductivity λ [W m-1 K -1]

A
1.6

1.2

All single data

0.8

q=6.25 m d-1
q=3.12 m d-1
q=1.56 m d-1
0.4

q=0.78 m d-1
q=0.19 m d-1
q=0

m d-1

0.0
0

0.1

0.2

0.3

0.4

0.6

B

Pressure head hm [m]

0.5

0.4

0.3

Multi-step data0.2

q=6.25 m d-1

the case with qw,z = 0.19 m d1 (black dashed line in
Fig. 9A) is an exception. Table 4 and Fig 9B presents
the optimized soil water retention parameters for the different steady ﬂux rates, and compare them with independent values, as obtained from the multi-step
outﬂow experiments. As expected, the most accurate
optimized retention parameters across the experimental
water content range were determined from the low waterﬂux experiments. Although not shown, when optimizations were carried out with a ﬁxed residual water content
value of hr = 0.0535 m3 m3, a much improved description of the retention curve for the high ﬂux steady-stateexperiments was obtained.
Using multiple datasets provided several advantages
for characterization of thermal and hydraulic properties.
First, by selecting the most optimal experimental conditions, a single experiment, with MFHPPÕs installed at
diﬀerent depths that combined cover a wide range ofwater contents, may provide suﬃcient information for
a complete description of the soil thermal conductivity
and soil hydraulic functions. In contrast, to obtain thesame information from a single probe, many measurements at diﬀerent water contents would be necessary.
Second, the single probe optimizations would requiresoil water potential measurements at the MFHPP locations, to obtain soil water retention curve parameters.
The multiple probe approach only needs measured values for the bottom boundary head (hm = 20 cm forour experiments) for successful optimization.

q=3.12 m d-1
q=1.56 m d-1

4.3. Tracer experiments

q=0.78 m d-1

0.1

-1

q=0.19 m d

m d-1

q=0

0.0
0

0 .1

0.2

0.3

0.4

Water content θ [m3 m-3]

Fig. 9. Thermal and hydraulic properties from multiple MFHPP
optimization of unsaturated experiments (option B). Based on
estimated parameters (Table 4) the thermal conductivity function (A)and retention curve (B) is presented for the water content range
applicable for each ﬂow rate.

high water ﬂuxes, the water content range is limited. For
example, for a steady-state ﬂux of qw,z = 6.25 m d1, theexperimental water content range is limited to only thatbetween 0.371 and 0.25. Alternatively, the experimental
water content range is between 0.11 and 0.37 m3 m3 forthe applied steady-state ﬂux of 0.19 m d1, thereby providing a much complete and accurate description of thefunctional relationship across the entire water content
range.
The optimized thermal conductivity parameters forthe unsaturated experiments are compared with the values found from the single MFHPP optimization (optionA) in Table 4 and Fig. 9A. Although, the multiple probe
optimizations agree well with the single probe results,

For the tracer experiments, the objective function
consists of ECbulk values, as measured by the four-electrode Wenner array of the four MFHPPÕs for each of the
10 steady-state water ﬂuxes (duplicate for each ﬂux density). Since the ECbulk contains information on both solute concentration and volumetric water content,
breakthrough measurements were used to simultaneously estimate both soil hydraulic and solute transport parameters [21,37]. The individual calibration
equations of Eq. (14) for each MFHPP were includedin the HYDRUS-2D code, so that calculated solute concentrations and water contents could be used to evaluate
ECbulk, and then used in the optimization process. The
dependency of water content on ECbulk resulted in distinctly diﬀerent breakthrough curves for both the saturated (Fig. 10A) and unsaturated (Fig. 10B) ﬂow
experiments, with ECbulk diﬀerences at equal concentrations between probes (such as the oﬀset) attributed toindividual calibration equations. When comparing theresults, we also notice the slightly larger noise in concentration values for the unsaturated experiments, eventhough the average velocity was much larger for the saturated experiments. More disturbing are the unstableconcentration values found along some of the saturated

A.P. Mortensen et al. / Advances in Water Resources 29 (2006) 250–267
20.0

A 500
MFHPP5450

MFHPP4400350

MFHPP3
MFHPP2

300

Model

17.5

15.0
-1

]

250

True velocity [m dECbulk [V V -1]

263

200
150
100
50
0
0

2

4

6

8

10

12.5

10.0

7.5

Saturated, EC (OptionC)

5.0Time [min]Saturated, MFHPP (Option B)

B 300

2.5

MFHPP5

Unsaturated, MFHPP (OptionB)

MFHPP4250

0.0

MFHPP3ECbulk [V V -1]

Unsaturated, EC (Option C)

0.0

MFHPP2

200

2.5

5.0

7.5

10.0

12.5

15.0

17.5

20.0

-1

Estimated velocity [m d ]ModelFig. 11. Water ﬂuxes estimated from multiple probe ECbulk optimizations (option C) compared with the true water ﬂuxes applied as top
boundary (values presented in Table 5). Estimated water ﬂuxes frommultiple probe heat pulse optimizations (option B, Table 4) are
included for comparison.

150
100
50
0
0

20

40

60

80
100
Time [min]

120

140

160

Fig. 10. Example of measured and estimated ECbulk breakthroughcurves for (A) saturated ﬂow at qwz = 17.84 m d1 and (B) unsaturated
ﬂow at qwz = 0.78 m d1. Measured ECbulk is both a function of theprobe and the water content, which is described by individual
calibration equations.

breakthrough curves at high ECbulk values resulting in‘‘bumps’’ on the breakthrough curves for reasons thatare unclear. Alike disturbances in the break throughcurves were measured consistently for all MFHPPÕs at
ECbulk values larger than about 350 mS cm1, indicatingslightly unstable Wenner array measurements at thesehigh concentrations. Finally, the unsaturated breakthrough curves showed signiﬁcant tailing, which has
been observed in general for unsaturated breakthroughcurves (e.g. [31]).
Optimized parameters for the saturated and unsaturated experiments are presented in Table 5. The opti-

mized water ﬂux values are compared in Fig. 11 withthe true water ﬂuxes and with those estimated with option B (i.e., multiple MFHPPÕs heat pulse signals).
RMSE values were 0.458 and 0.357 m d1 for the saturated and unsaturated experiments, respectively. Generally, the tracer optimizations gave better ﬂux estimates
for the saturated experiments, especially for the highﬂux experiments for which options A and B underestimated the water ﬂux. This is not surprising as solute
breakthrough is largely controlled by convection, especially for the relatively high water ﬂuxes of this experiment. For the unsaturated experiments, very similar
optimization results were found for the tracer experiments and the heat pulse experiments.
In addition to water ﬂuxes, solute dispersivity values
were also optimized (Table 5). For the saturated experiments, solute dispersivity was optimized directly withestimated values between 0.306 and 0.580 cm. For theunsaturated experiments, the water content dependenceof solute dispersivity was optimized through the two

Table 5
Option C. Optimized solute and hydraulic properties from multiple probe optimizations of saturated and unsaturated tracer experiments
Saturated
experiments

SolutepropertiesaL (cm)1

qw,z = 17.84 m d
qw,z = 13.39 m d1
qw,z = 9.37 m d1
qw,z = 4.37 m d1
qw,z = 1.05 m d1

0.580
0.360
0.312
0.306
0.480

Hydraulicproperties
qw,z (m d1)
18.49
12.70
8.99
4.16
1.00

Unsaturated
experiments
1

qw,z = 6.25 m d
qw,z = 3.12 m d1
qw,z = 1.56 m d1
qw,z = 0.78 m d1
qw,z = 0.19 m d1

Solute properties

Hydraulic properties

a1 (cm)a2

a (cm1)n

Ks (m s1)qw,z (m d1)

0.06430.000886
0.000223
0.00205
0.000635

1.73
3.17
3.47
2.17
2.59

0.02870.02930.03210.03270.0308

6.64
6.98
14.54
17.13
7.48

3.5E42.2E4
1.8E4
0.9E4
2.7E4

6.94
3.33
1.84
0.97
0.21

264

A.P. Mortensen et al. / Advances in Water Resources 29 (2006) 250–267

4.4. Coupled water, heat, and solute transport
optimizations5
q=6.25 m d-1

A
Solute dispersivity αL [cm]q=3.12 m d-14

q=1.56 m d-1
q=0.78 m d-1
q=0.18 m d-1

3

Nützmann et al. (2002)
Toride et al. (2003)

2

1

0
0

0.1

0.2

0.3

0.4

5

B

q=6.25 m d-1

Solute dispersivity α L [cm]q=3.12 m d-14

q=1.56 m d-1
q=0.78 m d-1
q=0.18 m d-1

3

Nützmann et al. (2002)
Toride et al. (2003)

2

1

0
0

0.1

0.2

0.3

0.4

Water content θ [m3 m-3]

Fig. 12. Solute dispersivity aL as a function of water content (Eq. (12))
based on estimated parameters a1 and a2 from (A) multiple optimization of for unsaturated ﬂow experiments (option C) and (B) multiple
optimization of coupled heat and solute ﬂow experiments (option D).
Corresponding estimated parameter values are presented in Tables 5
and 6.

parameters a1 and a2, as deﬁned in Eq. (12), with theparameter values and the functional forms presented
in Table 5 and Fig. 12A. As for option B, the mostwidely applicable results are expected for the lower ﬂux
experiments (qw = 0.18 m d1) for which the water content range was the largest. Interestingly, for these optimizations, our results most closely approach those of
Toride et al. [41]. The optimized parameter values compared reasonably well with those found by Nützmann
et al. [29] for a coarse-textured sand with a1 =0.0517 cm and a2 = 2.25 (Fig. 12A).
Also the soil hydraulic property parameters a, n, andKs were optimized from the unsaturated tracer experiment results, simultaneously with the dispersivity
parameters (Table 5). As expected, the best results wereobtained for the steady-state ﬂow experiments with thelowest ﬂuxes that correspond with the largest water content range. However, in general we conclude that the
ECbulk data alone was less suitable for estimation of
the hydraulic properties.

As will be demonstrated hereafter, the presented
MFHPP technique allows simultaneous measurementof parameters and variables characterizing coupledwater, heat, and solute transport. The option B results
demonstrated that coupled water ﬂow and heat transport was characterized from optimization of the heatpulse measurements, allowing simultaneous estimationof both thermal and hydraulic properties. Likewise, coupled water ﬂow and solute transport parameters were
estimated from the ECbulk measurements, providing
information for simultaneous estimation of solute transport and hydraulic properties. In this last section, we
will demonstrate that coupled water, heat, and soluteprocesses can be analyzed simultaneously using a combination of temperature and ECbulk measurements in thesame objective function.
Optimized parameters for both saturated and unsaturated ﬂow experiments are presented in Table 6. We
note that as for option B, the highest ﬂux optimizations
did not converge. Several diﬀerences were noticed whenusing coupled heat and solute data for the inverse analysis instead of heat (option B, Table 4) or solute (option
C, Table 5) data alone. First, whereas excellent waterﬂux estimates were obtained for saturated ﬂow in option
C, the coupled approach estimates were not as good,though closer than the underestimation of option B.
This was expected. However, water ﬂux errors resulted
in slightly larger estimated values for solute dispersivity(Table 6). For the unsaturated experiments, the diﬀerences in ﬂux estimates between options B and C wereminor, hence the coupled estimates (option D) wereequally accurate. Coupling the temperature with the
breakthrough data resulted in an excellent descriptionof the soil water retention characteristics, especially forthe low ﬂux density experiments. The improved soil
water retention optimization also aﬀected the optimized
solute dispersivity functions (compare Fig. 12A and B),
because of the inﬂuence of water content on ECbulk.
However, we note that the same general soil moisture effects on solute dispersivity were determined. Concurrent
optimized values for the thermal conductivity parameters were similar or better than using option B, andcompared very well with the single MFHPP optimization results of option A.
Although the coupled optimizations were successful
in many ways, the corresponding ﬁtting to the solute
breakthrough curves was generally disappointing. We
believe it to be mostly caused by the relative high uncertainty of the ECbulk measurements. First, as pointed out
in Section 3.1, the number of data pairs for the Wennerarray calibrations were limited, and was diﬀerent between probes. Second, we hypothesize that the Wennerarray sensor is unstable in the range with the highest

NO
3.41
1.73
0.91
0.25
NO
1.75
1.97
1.60
1.99

NO
0.0051
0.018
0.016
0.338

NO
3.85
2.22
2.31
0.86

NO
0.03080.03100.03020.0298

NO
16.27
14.29
12.97
11.21

NO
2.7E42.4E4
1.9E43.0E4

5. Summary and conclusions

qw,z = 6.25 m d1
qw,z = 3.12 m d1
qw,z = 1.56 m d1
qw,z = 0.78 m d1
qw,z = 0.19 m d1
NO: not optimized. Solutions did not converge.

16.12
11.85
8.33
3.99
0.98
0.688
0.633
0.516
0.434
0.509
qw,z = 17.84 m d1
qw,z = 13.39 m d1
qw,z = 9.37 m d1
qw,z = 4.37 m d1
qw,z = 1.05 m d1

1.87
1.86
1.84
1.83
1.79

qw,z (m d1)aL (cm)k0 (W m1 K1)

265

water content and solute concentration values, corresponding with the largest experimental range of ECbulk
values.

NO
1.18
1.08
1.44
1.07

Ks (m s1)na (cm1)a1 (cm)b1 (W m1 K1)b0 (W m1 K1)

Thermal properties
Unsaturated
experiments
Hydraulicproperties
Soluteproperties
Thermalproperties
Saturated
experiments

Table 6
Option D. Optimized parameters from coupled heat pulse and tracer experiments for saturated and unsaturated conditions

Soluteproperties

a2

Hydraulic properties

qw,z (m d1)

A.P. Mortensen et al. / Advances in Water Resources 29 (2006) 250–267

The new MFHPP technique combined with inverse
modeling allowed for simultaneous estimation of coupled thermal, hydraulic, and solute properties. Single
probe inverse optimization allowed simultaneous estimation of thermal characteristics, such as thermal conductivity and heat dispersion, soil hydraulic properties,
water ﬂux density and volumetric water content. Estimated parameters were generally in good agreementwith independently-measured values. Thermal dispersivity was found to be insigniﬁcant, conﬁrming that heatdispersion was not important in our experiments. The
main advantage of multiple probe optimizations is thatthe water content dependency of thermal and hydraulicrelationships could be determined simultaneously fromthe coupled column experiments, provided that a wide
range in experimental water content values can beachieved.
The Wenner array of the MFHPP provided for bulk
soil EC measurements during tracer breakthrough.
Since the EC measurements include information on both
solute concentration and soil water content, the coupled
estimation of solute dispersivity and hydraulic parameters was possible. For the unsaturated tracer experiments, we conclude that the solute dispersivityÕs
dependence on water content was described by a powerfunction. However, we also found that the Wenner arrayis not performing accurately at full saturation, in therange of high solute concentrations.
The combination of MFHPP measurements with inverse numerical analyses is shown to be a promising
method for the simultaneous analysis of coupled water,
heat, and solute ﬂow processes in variably-saturated
porous media. The probeÕs ability to simultaneously
measure thermal, hydraulic, and solute properties within
the same sample volume provides for a new powerfultool for vadose zone monitoring and characterization.
Moreover, we show that new insights can be learnt fromthe analysis of coupled ﬂow and transport measurements. Although the current study examines steady-state
ﬂow only, the MFHPP performs equally well for transient ﬂow conditions.

Acknowledgments
The authors wish to thank Dr. Keith Bristow of
CSIRO Land and Water, Australia for allowing us to
use their prototype of the four-needle heat pulse
probe. We also thank Dr. Gerard Kluitenberg for his

266

A.P. Mortensen et al. / Advances in Water Resources 29 (2006) 250–267

comments and encouragement. APM acknowledges
funding from VILLUM KANN RASMUSSEN andDONGÕs Jubilæumslegat. YM acknowledges funding
from Japan Society for the Promotion of Science(2001–2003). The work was also supported in part by
SAHRA (Sustainability of semi-Arid Hydrology andRiparian Areas) under the STC Program of the National
Science Foundation, Agreement No. EAR-9876800.
References[1] Abramovitz M, Stegun I. Handbook of mathematical functions. New York: Dover Publications; 1972 [p. 231].
[2] Basinger JM, Kluitenberg GJ, Ham JM, Frank JM, Barnes PL,
Kirkham MB. Laboratory evaluation of the dual-probe heat pulse
method for measuring soil water content. Vadose Zone J 2003;
2:389–99.
[3] Bear J. Dynamics of ﬂuids in porous media. New York: Elsevier;1972.
[4] Bilskie JR, Horton R, Bristow KL. Test of a dual-probe heatpulse method for determining thermal properties of porous
materials. Soil Sci 1998;163:346–55.
[5] Bristow KL, Campbell GS, Calissendorﬀ K. Test of a heat-pulse
probe for measuring changes in soil water content. Soil Sci Am J
1993;57:930–4.
[6] Bristow KL, Kluitenberg GJ, Horton R. Measurement of soil
thermal properties with a dual-probe heat-pulse technique. Soil
Sci Soc Am J 1994;58:1288–94.
[7] Bristow KL. Measurement of thermal properties and watercontent of unsaturated sandy soil using dual-probe heat-pulse
probes. Agric Forest Meteorol 1998;89:75–84.
[8] Bristow KL, Kluitenberg GJ, Goding CJ, Fitzgerald TS. A small
multi-needle probe for measuring soil thermal properties, watercontent and electrical conductivity. Comput Electron Agric 2001;
31:265–80.
[9] Cass A, Campbell GS, Jones TL. Enhancement of thermal vapor
diﬀusion in soil. Soil Sci Am J 1984;48:25–32.
[10] Campbell GS. Soil physics with BASIC—transport models forsoil-plant systems. New York: Elsevier; 1985.
[11] Campbell GS, Calissendorﬀ C, Williams JH. Probe for measuring
soil speciﬁc heat using a heat-pulse method. Soil Sci Am J 1991;
55:291–3.
[12] Chung S-O, Horton R. Soil heat and water ﬂow with a partial
surface mulch. Water Resour Res 1987;23:2175–86.
[13] de Marsily G. Quantitative hydrogeology: groundwater hydrology
for engineers. San Diego: Academic; 1986.
[14] de Vries DA. Thermal properties of soils. In: van Wijk WR,
editor. Physics of plant environment. New York: North-Holland; 1963.
[15] Eching SO, Hopmans JW. Optimization of hydraulic functions
from transient outﬂow and soil water pressure data. Soil Sci Soc
Am J 1993;57:1167–75.
[16] Heitman JL, Basinger JM, Kluitenberg GJ, Ham JM, Frank JM,
Barnes PL. Field evaluation of the dual-probe heat pulse method
for measuring soil water content. Vadose Zone J 2003;2:552–60.
[17] Hopmans JW, Dane JH. Thermal conductivity of two porous
media as a function of water content, temperature, and density.
Soil Sci 1986;142:187–95.
[18] Hopmans JW, Dane JH. Temperature dependency of soil
hydraulic properties. Soil Sci Am J 1986;50:4–9.
[19] Hopmans JW, Šimůnek J, Bristow KL. Indirect estimation of soil
thermal properties and water ﬂux using heat pulse probe
measurement: geometry and dispersion eﬀects. Water Resour
Res 2002;38(1). doi:10.1029/2000WR000071.

[20] Hopmans JW, Šimunek J, Romano N, Durner W. Simultaneous
determination of water transmission and retention properties—
inverse modeling of transient water ﬂow. In: Topp GC, Dane JH,
editors. Methods of soil analysis, Part 4. Physical methods. Soil
Science Society of America Book Series, no. 5.
[21] Inoue M, Šimůnek J, Shiozawa S, Hopmans JW. Simultaneous
estimation of soil hydraulic and solute transport parameters fromtransient inﬁltration experiments. Adv Water Resour 2000;23:
677–88.
[22] Kluitenberg GJ, Ham JM, Bristow KL. Error analysis of the heatpulse method for measuring soil volumetric heat capacity. Soil Sci
Am J 1993;57:1444–51.
[23] Kluitenberg GJ, Heitman JL. Eﬀect of forced convection on soil
water content measurement with then dual-probe heat-pulse
method. In: Smiles D, Raats PAC, Warrick A, editors. Heatand mass transfer in natural environment, the Philip volume.
AGU, geophysical monograph series, no. 129.
[24] Mishra S, Parker JC. Parameter estimation for coupled unsaturated ﬂow and transport. Water Resour Res 1989;25:385–96.
[25] Mori Y, Hopmans JW, Mortensen AP, Kluitenberg GJ. Multifunctional heat pulse probe for the simultaneous measurement ofsoil water content, solute concentration, and heat transportparameters. Vadose Zone J 2003;2:561–71.
[26] Mori Y, Hopmans JW, Mortensen AP, Kluitenberg GJ. Estimation of vadose zone water ﬂux from multi-functional heat pulse
probe measurement. Soil Sci Soc Am J 2005;69(3):599–606.
[27] Mortensen AP, Jensen KH, Nilsson B, Juhler RK. Multiple
tracing experiments in unsaturated fractured clayey till. Vadose
Zone J 2004;3:634–44.
[28] Noborio K, McInnes KJ, Heilman JL. Measurements of soil
water content, heat capacity, and thermal conductivity with asingle TDR probe. Soil Sci 1996;161:22–8.
[29] Nützmann G, Maciejewski S, Joswig K. Estimation of water
saturation dependence of dispersion in unsaturated porous media:
experiments and modelling analysis. Adv Water Res 2002;25:
565–76.
[30] Ochsner TE, Horton R, Ren T. Simultaneous water content, airﬁlled porosity, and bulk density measurements with thermo-time
domain reﬂectometry. Soil Sci Am J 2001;65:1618–22.
[31] Pot V, Šimůnek J, Benoit P, Coquet Y, Yra A, Martı́nez-Cordón
M-J. Impact of rainfall intensity on the transport of two
herbicides in undisturbed grassed ﬁlter strip soil cores. J Contam
Hydrol, in press.
[32] Ren T, Noborio K, Horton R. Measuring soil water content,
electrical conductivity, and thermal properties with a thermo-time
domain reﬂectometry probe. Soil Sci Soc Am J 1999;63:450–7.
[33] Ren T, Kluitenberg GJ, Horton R. Determining soil water ﬂux
and pore water velocity by a heat pulse technique. Soil Sci Am J
2000;64:552–60.
[34] Ren T, Ochsner TE, Horton R. Development of thermo-time
domain reﬂectometry for vadose zone measurement. Vadose Zone
J 2003;2:544–51.
[35] Rhoades JD, Raats PAC, Prather RJ. Eﬀects of liquid-phaseelectrical conductivity, water content, and surface conductivity on
bulk soil electrical conductivity. Soil Sci Soc Am J 1976;40:651–5.
[36] Šimůnek J, Šejna M, van Genuchten MT. The HYDRUS-2D
software package for simulating two-dimensional movement ofwater, heat, and multiple solutes in variably-saturated media,
version 2.0, Rep IGCWMC-TPS-53, 251 pp. Int Ground Water
Model. Cent, Colo Sch of Mines, Golden; 1999.
[37] Šimůnek J, Jacques D, Hopmans JW, Inoue M, Flury M, vanGenuchten MTh. Solute transport during variably-saturated
ﬂow—inverse methods. In: Dane JH, Topp GC, editors. Methodsof soil analysis, Part 1. Physical methods. 3rd ed. Madison,
WI: SSSA; 2002. p. 1435–49 [Chapter 6.6].
[38] Sophocleous M. Analysis of water and heat ﬂow in unsaturated–saturated porous media. Water Resour Res 1979;15:1195–206.

A.P. Mortensen et al. / Advances in Water Resources 29 (2006) 250–267[39] Sun N-Z, Yeh WW-G. Coupled inverse problems in groundwatermodeling. 1. Sensitivity analysis and parameter identiﬁcation.
Water Resour Res 1990;26(6):2507–25.
[40] Tarara JM, Ham JM. Measuring soil water content in thelaboratory and ﬁeld with dual-probe heat-capacity sensors. Agron
J 1997;89:535–42.
[41] Toride N, Inoue M, Leij FJ. Hydrodynamic dispersion in an
unsaturated dune sand. Soil Sci Soc Am J 2003;67:703–12.
[42] Tuli A, Denton MA, Hopmans JW, Harter T, MacIntyre JL.
Multi-step outﬂow experiment: from soil preparation to parameter estimation. Land, Air and Water Resources Report No.
100037, University of California, Davis, CA 95616, USA, 2001.

267[43] van Genuchten MTh. A closed-form equation for predicting thehydraulic conductivity of unsaturated soils. Soil Sci Soc Am J
1980;44:892–8.
[44] Vaz CMP, Hopmans JW, Macedo A, Bassoi LH, Wildenschild D.
Soil water retention measurements using a combined tensiometercoiled time domain reﬂectometry probe. Soil Sci Am J
2002;66:1752–9.
[45] Wang Q, Ochsner TE, Horton R. Mathematical analysis of heatpulse signals for soil water ﬂux determination. Water Resour Res
2002;38(6). doi:10.1029/2001WR001089.
[46] Wraith JM, Or D. Nonlinear parameter estimation using spreadsheet software. J Nat Resour Life Sci Edu 1998;27:13–9.