Soil Science Society of America Journal 65:320-323 (2001)
© 2001 Soil Science Society of America
DIVISION S-1-SOIL PHYSICS
Improved evaluation procedure for heat-pulse soil water flux density method
G.J. Kluitenberga and
A.W. Warrickb
a Dep. of Agronomy, Kansas State Univ., Manhattan, KS 66506
b Dep. of Soil, Water and Environ. Sci., Univ. of Arizona, Tucson, AZ 85721
Corresponding author (gjk{at}ksu.edu)
 |
ABSTRACT
|
|---|
Analytical solutions of the heat equation must be evaluated in order to implement the heat-pulse method for measuring soil water flux density. We developed an improved procedure for evaluating the integrals in these solutions by recognizing that they can be reduced to a single function, W, known as the well function for leaky aquifers. The evaluation procedure was improved further by developing an efficient method to approximate W. This method, which involves summing the first few terms of an infinite series, also provides a simple means of determining the approximation error. For a wide range of input parameters, at most two terms of the series are needed to approximate W with error less than 10-4 in absolute value. Thus, numerical integration is not required in order to implement the heat-pulse method for measuring soil water flux density. Our results regarding the evaluation of W are relevant to other problems in soil science and groundwater hydrology in which the function W appears.
|
INTRODUCTION
|
|---|
REN ET AL. (2000) proposed a new method for measuring soil water flux density in which a heat tracer was used to quantify the magnitude of convective heat transfer resulting from soil water movement. Their method utilized analytical solutions of the heat equation to describe temperature changes that occur upstream and downstream of a line heat source following the emission of a heat pulse. Unfortunately, these solutions contain integrals that must be evaluated using numerical integration techniques. Insofar as one of the integrals is improper, specialized numerical integration procedures are required.
In this note, we show that the integrals in the aforementioned solutions can be reduced to a single integral, known in the groundwater hydrology literature as the well function for leaky aquifers (Hantush, 1964; Hantush and Jacob, 1955; Neumann and Witherspoon, 1969). This simplifies the task of evaluating the analytical solutions. Also, we introduce an efficient method, based on the results of Bruggeman (1999), for evaluating the well function for leaky aquifers. Evaluation of the analytical solutions is simplified further by using this method. Example calculations are included to illustrate how the improved evaluation procedure simplifies implementation of the method of Ren et al. (2000).
|
Theory
|
|---|
Problem of Interest
Ren et al. (2000) based their method on a solution to the heat equation
 | (1) |
in which T is temperature (°C), t is time (s), 
is the thermal diffusivity (m2 s-1), and x and y are space coordinates (m). This equation describes coupled conduction and convection in a porous medium through which water moves at constant velocity in the x direction. The heat pulse velocity, V, is related to the water flux density, J, by the expression V = J(
c)
/(c) where (c) and (c) are the volumetric heat capacities of water and the soil–water–air system, respectively (Ren et al., 2000). The heat pulse velocity (m s-1) is taken to be a positive quantity for flow in the positive x direction.
Ren et al. (2000) considered an infinite line source, coinciding with (x, y) = (0,0) and normal to the x–y plane, which is heated at the rate q'(W m-1) during the time interval 0 < t 
t0. For this heat source, the solution of Eq. [1] is (their Eq. [7])
 | (2) |
where 
= 4
T/q', and is the thermal conductivity (W m-1 °C-1). This solution gives dimensionless temperatures downstream (x > 0) or upstream (x < 0) from the line source.
We note here that Eq. [2] is identical in form to the solution for a two-dimensional solute plume in a porous medium through which water moves at constant velocity in one direction. In this case, the solution describes the solute concentration in the porous medium after a pulse of solute is released from an infinite line source (cf. Fetter, 1999, Eq. [2.31]). Thus, results presented hereafter apply to the solute plume problem as well as the heat transfer problem.
Transformation of Solution
The substitution = (x2 + y2)/4s transforms Eq. [2] into
 | (3) |
where v = Vx/2, 
= V
/2, 
= 
/4t, and ' = /4
. But the integrals in Eq. [3] are identical to the well function for leaky aquifers, defined as (Hantush, 1964, p. 321)
 | (4) |
Thus, Eq. [3] can be written in the reduced form
 | (5) |
Series Approximation of W(u,ß)
Consider the Maclaurin series of the function exp(-
) for > 0 (Kaplan, 1984, p. 431)
 | (6) |
with remainder term
 | (7) |
in which 0 < 1 < . Equation [7] indicates that the Maclaurin series can be used to approximate exp(-) with an error of exactly Rn if truncated after term n, provided 1 is known. Although there is no convenient means of determining 1, an upper limit for the absolute value of Rn can be obtained. Because is positive, we know that
 | (8) |
Thus, an upper limit for the magnitude of Rn is
 | (9) |
Comparing Eq. [6] and [9] shows that the upper limit for |Rn| is no greater than the absolute value of the first truncated term of the series.
Upon substituting ß2/4z for in Eq. [6] and using the result in Eq. [4], we get
 | (10) |
Then, by substituting ß2/4z for in Eq. [9] and using the result in Eq. [10], we obtain the inequality
 | (11) |
in which Em+1(u) is an exponential integral (see Eq. [5.1.4] of Gautschi and Cahill, 1972). The finite series in Eq. [11] was obtained by performing term-by-term integration of the series in Eq. [10]. Rearrangement of the integral in Eq. [11] eventually leads to the result
 | (12) |
which indicates that the error in truncating the alternating series is no greater in absolute value than the absolute value of the first truncated term. This leads to the approximation
 | (13) |
which is a truncated form of the infinite series given by Bruggeman (1999)(Eq. [16], p. 878). The recurrence relation (Gautschi and Cahill, 1972, p. 229)
 | (14) |
is useful for the evaluation of Eq. [13].
|
Results and Discussion
|
|---|
Numerical Evaluation of W(u,ß)
We have shown that the integrals in Eq. [2] can be expressed in terms of the well function for leaky aquifers. If W(u,ß) is evaluated by numerical integration of Eq. [4], the computational task is now simplified inasmuch as a single function must be evaluated. We note here that values of W(u,ß) tabulated by Hantush (1964) provide a useful means of checking numerical integration schemes for select values of u and ß. Nevertheless, numerical integration of Eq. [4] still requires some care because it is an improper integral and the error analysis can be vague.
Series approximations for W(u,ß) are available (Hantush and Jacob, 1955), but they are too cumbersome for routine use. Equation [13] forms the basis of a simpler approach. The series in Eq. [13] converges rapidly for u > ß/2. An alternate form
 | (15) |
follows from Eq. [13] with use of the identity (Hantush, 1964, p.321) W(u,ß) = 2K0(ß) - W(ß2/4u,ß), where K0 is the zero-order modified Bessel function of the second kind. The series in Eq. [15] converges rapidly for u < ß/2 and converges at the same rate as the series in Eq. [13] along the line u = ß/2. Thus, by summing a finite number of the terms in Eq. [13] or [15], W(u,ß) can be approximated with known accuracy.
Table 1 shows the number of terms that must be used in Eq. [13] or [15] to approximate W(u,ß) with error no greater than 10-4 in absolute value. In performing the calculations for Table 1, E1 was evaluated by using Eq. [5.1.53] and [5.1.56] of Gautschi and Cahill (1972); K0 was evaluated by using Eq. [9.8.1], [9.8.2], [9.8.5], and [9.8.6] of Olver (1972). The results indicate that W(u,ß) can be approximated with excellent accuracy by using only zero, one, or two terms in Eq. [13] or [15] throughout much of the u-ß domain (Table 1). More than two terms are required only near the line u = ß/2 for large values of ß.
View this table:
[in this window]
[in a new window]
|
Table 1. Number of terms required to approximate W(u,ß) with error no greater than 10-4 in absolute value. Entries above those in parentheses (upper right portion of table) refer to the number of terms required in Eq. [15]. Entries below those in parentheses (lower left portion of table) refer to the number of terms required in Eq. [13]. For entries in parentheses, Eq. [13] and [15] both require the same number of terms
|
|
Indiscriminate use of Eq. [13] or [15] with a large number of terms generally is not recommended, because Eq. [14] suffers from numerical instability under conditions described by Gautschi and Cahill (1972)(p. 233–234, ex. 4–6). They indicate that the recurrence relation can be used safely in the forward direction (increasing m) for arguments less than about five. Thus, numerical instability becomes an issue only when using Eq. [14] to evaluate Eq. [13] with u > 5 or Eq. [15] with ß2/4u > 5. Fortunately, this limits potential instability problems to a small region defined by 10 < ß < 100 and 5 < u < 50, a region not likely to be encountered in practice. If calculations must be performed in this region, numerical instability can be avoided by using Eq. [14] in the backward direction (decreasing m) after obtaining a starting value with Eq. [5.1.52] of Gautschi and Cahill (1972).
Example Calculations
The following example illustrates how Eq. [5] can be evaluated using the procedure described in the previous section. We also illustrate the procedure with the solution given by Ren et al. (2000) for the temperature difference (
T) between positions xd and xu, arbitrary positions downstream and upstream from the heat source, respectively. A generalized version of their Eq. [18] (y 
0) can be written as
 | (16) |
where = 4T/q' and the subscripts "d" and "u" indicate values of 
, , , and ' in which xd and xu are substituted for x, respectively. For the special case where x = xd = |xu|, Eq. [16] reduces to
 | (17) |
with a positive quantity. Substituting tm (time at which reaches a maximum) for t in and ' changes Eq. [16] and [17] into transformed versions of the maximum dimensionless temperature difference expression of Ren et al. (2000). See their Eq. [21].
Example calculations were performed with xd = 6.0 mm, xu = -6.0 mm, = 6 x 10-7 m2 s-1, and t0 = 15 s, values similar to those used by Ren et al. (2000). In addition, we set V = 2 x 10-5 m s-1 ( = = 0.1), which approximates the lowest heat pulse velocity used by Ren et al. (2000). The temperature was calculated for a position directly downstream from the heater by using Eq. [5] with y = 0. Equation [17] with y = 0 was used to calculate . With y = 0, is the temperature difference between positions directly downstream and directly upstream from the heater. Results are given for the time interval 0 < t 60 s (Table 2). Equation [13] was used for all approximations of W(,) and W(',) in this example, because all values of and ' were greater than /2. The series in Eq. [13] was truncated, so that approximation errors were no greater than 10-4 in absolute value. The values of resulting from these calculations are identical to those used to plot the curve for V = 2 x 10-5 m s-1 in Fig. 2 of Ren et al. (2000). The results (Table 2) show that at most two terms of the series in Eq. [13] were needed to approximate W(,) and W(',) to within the specified level of accuracy.
View this table:
[in this window]
[in a new window]
|
Table 2. Example calculation of (downstream position) and using y = 0, xd = 6.0 mm, xu = -6.0 mm, = 6 x 10-7 m2 s-1, t0 = 15 s, and V = 2 x 10-5 m s-1 ( = 0.1). The numbers in parentheses following W(,) and W(',) indicate the number of terms used in Eq. [13] to calculate these quantities with error no greater than 10-4 in absolute value. Also shown are the errors (%) in and for = 0.04, 0.1, and 0.4 when the calculations were performed using only E1() and E1(') to approximate W(,) and W(', ), respectively. Errors in and are identical when expressed on a percentage basis
|
|
Hantush (1964) suggested that the approximation W(u,ß) 
E1(u) is useful for u > 2ß. In this example, errors in and did not exceed 1% (Table 2) when W(,) and W(',) were approximated with E1() and E1('), respectively. Furthermore, the maximum value of , a quantity used in the method of Ren et al. (2000), was calculated with error less than 0.3%. The error resulting from the approximation W(u,ß) E1(u) increases rapidly for larger heat-pulse velocities, but also diminishes rapidly for smaller heat-pulse velocities. The previous calculations were repeated with V = 8 x 10-6 m s-1 ( = 0.04) and V = 8 x 10-5 m s-1 ( = 0.4) to illustrate this (Table 2), but note that the condition u > 2ß is violated for t > 18 s when = 0.4. Errors in the maximum value of are 0.05% for = 0.04 but 4.5% for = 0.4.
|
Summary
|
|---|
We have shown that the integrals in the analytical solutions of Ren et al. (2000) can be expressed in terms of W(u,ß), the well function for leaky aquifers. This simplifies evaluation of the analytical solutions inasmuch as only a single function must be evaluated. We also have presented an efficient method for evaluating W(u,ß), based on Bruggeman's (1999) series representation of W(u,ß). Use of this method further simplifies evaluation of the solutions of Ren et al. (2000). The method for approximating W(u,ß) involves use of Eq. [13] for u > ß/2 and Eq. [15] for u < ß/2. We have shown that W(u,ß) can be approximated with excellent accuracy by summing only a few terms of the series in Eq. [13] and [15]. In addition, we have shown that the magnitude of the error incurred by using the truncated series in Eq. [13] and [15] is no greater than the absolute value of the first truncated term. Thus, our method allows W(u,ß) to be approximated with known accuracy.
Finally, we reiterate that our results regarding the evaluation of W(u,ß) are relevant to other problems in soil science and groundwater hydrology. The function W(u,ß) is encountered when solving the convection–dispersion equation for a two-dimensional contaminant plume (Fetter, 1999) and when solving for the transient drawdown near a well pumping from a leaky, confined aquifer of infinite extent (Hantush, 1964; Hantush and Jacob, 1955; Neumann and Witherspoon, 1969).
|
NOTES
|
|---|
Contribution no. 00-331-J from the Kansas Agric. Exp. Stn., Manhattan, KS. Research supported by Western Regional Research Project W-188.
Received for publication March 20, 2000.
|
REFERENCES
|
|---|
- Bruggeman, G.A. 1999. Analytical solutions of geohydrological problems. Developments in Water Science 46. Elsevier, Amsterdam.
- Fetter, C.W. 1999. Contaminant hydrology. 2nd ed. Prentice Hall, Upper Saddle River, NJ.
- Gautschi, W., and W.F. Cahill. 1972. Exponential integral and related functions. p. 227–254. In M. Abramowitz and I.A. Stegun (ed.) Handbook of mathematical functions. Dover, New York.
- Hantush, M.S. 1964. Hydraulics of wells. Adv. Hydrosci. 1:281–432.
- Hantush, M.S., and C.E. Jacob. 1955. Nonsteady radial flow in an infinite leaky aquifer. Trans. Am. Geophys. Union. 36:95–100.
- Kaplan, W. 1984. Advanced calculus. 3rd ed. Addison-Wesley, Reading, MA.
- Neumann, S.P., and P.A. Witherspoon. 1969. Applicability of current theories of flow in leaky aquifers. Water Resour. Res. 4:817–829.
- Olver, F.W.J. 1972. Bessel functions of integral order. p. 355–433. In M. Abramowitz and I.A. Stegun (ed.) Handbook of mathematical functions. Dover, New York.
- Ren, T., G.J. Kluitenberg, and R. Horton. 2000. Determining soil water flux and pore water velocity by a heat pulse technique. Soil Sci. Soc. Am. J. 64:552–560.[Abstract/Free Full Text]
This article has been cited by other articles:
|
|
|
|
 
J. Gao, T. Ren, and Y. Gong
Correcting Wall Flow Effect Improves the Heat-Pulse Technique for Determining Water Flux in Saturated Soils
Soil Sci. Soc. Am. J.,
March 29, 2006;
70(3):
711 - 717.
[Abstract]
[Full Text]
[PDF]
|
|
|
|
|
|
|
T. E. Ochsner, R. Horton, G. J. Kluitenberg, and Q. Wang
Evaluation of the Heat Pulse Ratio Method for Measuring Soil Water Flux
Soil Sci. Soc. Am. J.,
May 6, 2005;
69(3):
757 - 765.
[Abstract]
[Full Text]
[PDF]
|
|
|
|
|
|
|
Y. Mori, J. W. Hopmans, A. P. Mortensen, and G. J. Kluitenberg
Estimation of Vadose Zone Water Flux from Multi-Functional Heat Pulse Probe Measurements
Soil Sci. Soc. Am. J.,
April 11, 2005;
69(3):
599 - 606.
[Abstract]
[Full Text]
[PDF]
|
|
|
|
|
|
|
G. W. Gee, G. W. Gee, Z. F. Zhang, and A. L. Ward
A Modified Vadose Zone Fluxmeter with Solution Collection Capability
Vadose Zone J.,
November 1, 2003;
2(4):
627 - 632.
[Abstract]
[Full Text]
[PDF]
|
|
|
TOP
ABSTRACT
INTRODUCTION
