Soil Science Society of America Journal 64:1219-1225 (2000)
© 2000 Soil Science Society of America
DIVISION S-1-SOIL PHYSICS
A Thermal Wave Approach for Heat Transfer in a Nonuniform Soil
Mostafa A. Karam
GenCorp Aerojet, 1100 W. Hollyvale St., Azusa, CA 91702 USA
mostafa.karam{at}aerojet.com
 |
ABSTRACT
|
|---|
A novel thermal wave model is developed for studying heat transfer in nonuniform soils. The model, which is based on the wave-like characteristics of periodic heat flow, discretizes the soil into sublayers. The total thermal reflection and transmission coefficients at the interfaces of those sublayers are formulated and then employed to construct the soil temperature and heat flux profiles from their values at the air–soil interface. Numerical simulations showed that the thermal wave model could predict the predictions of two exact analytic models for the temperature profiles and their derivatives.
|
INTRODUCTION
|
|---|
IT HAS BEEN SHOWN that periodic heat flow in a medium can be described by the propagation of thermal waves in such a medium. Those thermal waves have characteristics similar to the characteristics of other types of waves (e.g., electromagnetic waves and acoustic waves). Among these characteristics are the reflection and transmission at interfaces separating two media, scattering from thermal anomalies, and interference. There has been considerable interest in exploiting the characteristics of thermal waves in detecting, analyzing, and imaging thermal and nonthermal features in a variety of materials (Bennett and Patty, 1982; Rosencwaig et al., 1985; Tang and Araki, 1996). In this study, two characteristics of thermal waves in layered media are explored and then used to study heat transfer in nonuniform soils. This study is the first step in our attempt to answer the following questions: can we apply techniques developed for studying wave propagation to study heat transfer, and how accurate are the results stemming from applying those techniques? Answering these questions may lead to an alternative efficient technique for studying heat transfer in soils, or it may enhance some of the existing techniques for studying heat transfer.
To achieve the objective of this study, the thermal reflection and transmission coefficients characterizing thermal waves at the interfaces of layered soils are formulated and then used to create a thermal wave model for constructing temperature profiles in nonuniform soils. The model is tested and validated using the predictions of two exact analytic models. One of those models has been developed for soils having thermal properties described by exponential profiles (Massman, 1993), and the other model has been developed for soils having thermal properties described by linear profiles (Nerpin and Chudnovskii, 1984). The mathematical formulations describing those two analytic models are usually given in terms of Kelvin's functions of different orders. To ensure that those formulations are correct and their predictions are exact, they are rederived, and programmed on a digital computer using FORTRAN language.
|
Theory and methods
|
|---|
Thermal Wave Characteristics in Layered Soils
Let us explore the thermal wave characteristics of periodic heat transfer in layered soils and derive two of those characteristics that will be used to study heat transfer in nonuniform soils. To do so, we consider a layered soil occupying the lower-half space and having (N + 1) layers (Fig. 1)
. Any arbitrary layer
within the layered soil has thermal properties (e.g., thermal conductivity 
n, and heat capacity cn) that are taken to be uniform, time independent, and different from the corresponding properties of other layers. The one-dimensional equation governing the temperature Tn (z,t) within the arbitrary layer n can be written as
,
 | (1) |
where T is temperature (°C), t is time (s), z is an arbitrary depth below the surface within a layer, and d is a depth below the surface at a layer interface.
In the quasi steady-state portion of the solution of Eq. [1], the initial conditions are not required and time domain problem can be extended to become -
< t < . This permits writing the daily temperature in the following Fourier series.
 | (2) |
In Eq. [2], 
m = m, m = 1,2, ..., is the daily radial frequency, m is the order of temperature harmonics, and j is the complex number (
- 1). In addition, Re( ) is the real part operator, and Tn(z,m) is the complex frequency-dependent temperature. For the sake of brevity, this study will focus on the first harmonics only
, and the harmonic subscript m will be dropped. Substituting Eq. [2] into Eq. [1] yields the differential equation governing Tn(z,).
 | (3) |
where
 | (4) |
Similar to electromagnetic wave equations, Eq. [3] can be considered a thermal wave equation governing time harmonic temperature Tn(z,) propagating within the nth layer with propagation constant 
n. Accordingly, the solution of Eq. [3] for Tn(z,) can be written as
 | (5) |
In addition, the heat flux 
n(z,) associated with Tn(z,) is
 | (6) |
In Eq. [5], Bn is the amplitude of the upward thermal wave at the lower boundary of the nth layer
. An is the amplitude of downward thermal wave at the upper boundary of the layer
. By analogy to electromagnetic wave theory (Karam, 1996), the thermal reflection coefficient inside the nth layer at its lower boundary
can be obtained from the ratio of amplitudes of the upward and downward thermal waves at such a boundary. Moreover, the thermal transmission coefficient can be obtained from the ratio of the amplitudes of forward waves across such a boundary. The physical meaning of the reflection and transmission coefficients can be found in any textbook dealing with electromagnetic waves or optics, such as the text written by Wait (1975).
Thermal Reflection Coefficient
To obtain the thermal reflection coefficient inside the nth layer at the lower boundary of such a layer
, we apply the continuity of both temperature and heat flux across such a boundary. The continuity of temperature at
can be obtained by using Eq. [5] to get
 | (7) |
The continuity of heat flux at z = -dn+1, can be obtained by using Eq. [6] to obtain
 | (8) |
In Eq. [7] and [8],
is the height of the nth layer, and hn+1 is the height of the (n + 1)th layer. Dividing Eq. [7] by [8] yields
 | (9) |
In the above Rn is the thermal reflection coefficient inside the nth layer at its lower boundary
.
 | (10) |
As for Rn+1, it is the reflection coefficient inside the (n + 1)th layer at the lower boundary of such a layer
. The relation governing Rn+1 can be obtained from Eq. [10] through replacing n with n + 1. A further reduction for Eq. [9] leads to
 | (11) |
Then by grouping the coefficients of Rn+1 with each other and using Eq. [4], we simplify Eq. [11] to
 | (12) |
where
 | (13) |
In the above 
n,n+1 is the unbounded thermal reflection coefficient inside the nth layer at its lower boundary
. The unbounded thermal reflection coefficient is similar to the Fresnel reflection coefficient in electromagnetic wave theory (Wait, 1975; Karam, 1996). Both coefficients describe the reflection of waves at an interface, as if the layer beneath such an interface has an infinite height. This can be easily envisioned by setting
in Eq. [12]. In addition, both coefficients arise from discontinuities in properties of media surrounding the interface. For instance, if the layers surrounding an interface have similar thermal properties, the value of the unbounded thermal reflection coefficient reduces to zero. On the other hand, setting
in Eq. [13] produces what is introduced by Van Wijk and Derksen (1963) under the name "an auxiliary variable" (Van Wijk and Derksen, 1963, Eq. [6.9]).
Now let us examine the value of the reflection coefficient (Eq. [12]) at the lower boundary of the Nth layer
. Since the layer beneath such a boundary has an infinite height
, setting
into Eq. [12] yields
 | (14) |
This value for the reflection coefficient confirms the definition introduced earlier for the unbounded reflection coefficient.
Thermal Transmission Coefficient
To obtain the thermal transmission coefficient at the lower boundary of the nth layer
we should determine the relations between the amplitudes of the forward waves across such a boundary. To do so, we use Eq. [10] into Eq. [7] to get
 | (15) |
Then from Eq. [15] we can relate An+1 to An as
 | (16) |
which can be reduced by using Eq. [12]
 | (17) |
In Eq [17], In,n+1 is the thermal transmission coefficient from the nth layer to the (n + 1)th layer across the boundary
with
 | (18) |
In Eq. [18], 
n,n+1 is the unbounded thermal transmission coefficient from the nth layer to the (n + 1)th layer.
 | (19) |
The above unbounded transmission coefficient represents the transmission coefficient across the lower boundary of the nth layer as if the layer beneath such a boundary has an infinite height. This can be justified by letting hn+1 in Eq. [18] go to .
Thermal Wave Model for a Nonuniform Soil
Here the thermal reflection Eq. [12] and transmission Eq. [18] coefficients are used to create a thermal wave model for constructing temperature profiles in nonuniform soils. Soil nonuniformity arises from the dependence of soil thermal properties on soil water content and soil bulk density. Both water content and bulk density vary strongly with soil depth near the soil–air interface because of evaporation, drainage, transpiration by plants, tillage, and soil formation process. At soil depths far below the soil–air interface, the dependence of water content and bulk density, and hence the thermal properties, on the above parameters decrease until it vanishes (Novak and Black, 1985; Novak, 1986; Karam, 1999). At those soil depths, soil thermal properties are almost uniform. Thus, in creating the thermal wave model for a nonuniform soil, the soil is divided into two zones: an upper zone extending across soil depths where variations in thermal properties occur, and a lower zone extending across soil depths where thermal properties are uniform. The upper zone is discretized into N layers. The heights of those layers are taken to be small enough to replace thermal properties within each layer by their average values. The lower zone is taken as a semi-infinite layer with uniform thermal properties and is labeled as the (N + 1)th layer.
In the next step in constructing the thermal wave model, we use Eq. [14] to calculate the thermal reflection coefficients RN at the interface separating the two zones
. Then we set
in Eq. [12] and substituting for RN from Eq. [14] to obtain RN-1. Furthermore, we use RN-1 in Eq. [12] after setting
to obtain Rn-2, and so on, until we obtain R1 at the lower boundary of the upper layer
. At this point, the problem of heat transfer in the upper zone is reduced to the problem of heat transfer in a uniform soil layer with two boundaries, a lower boundary located at
and an upper boundary located at
. At the lower boundary
, the thermal waves are described by the thermal reflection coefficient R1, and at the upper boundary
, the thermal waves are controlled by the boundary conditions at the soil–air interface. There are two types of such boundary conditions reported in the literature. The use of each type of those boundary conditions is discussed below.
The first type of boundary condition is based on the continuity of heat fluxes. To apply this type of boundary conditions we set
in Eq. [8], equate its left-hand side to the heat flux at the soil–air interface air(), and use Eq. [10] to get
 | (20) |
The second type of boundary conditions is based on a linearized form of energy balance equation (Novak, 1993).
 | (21) |
where ß is the known constant temperature coefficient and f is the known function.
Taking the Fourier transform of the above equation yields
 | (22) |
In Eq. [22], f() is the Fourier transform of f(t).
 | (23) |
Then by using Eq. [5], [10], [22] we get
 | (24) |
At this point the value of A1 is obtained. It can be used in Eq. [17] after setting
to obtain A2. The latter will be used in Eq. [17] after setting
to obtain A3, and so on, until the values of all Ans are obtained. Incorporating those values in Eq. [10] yields the values of Bns. The values of both Ans and Bns will be used in Eq. [5] and [6] to construct frequency-dependent temperature and heat flux profiles within the upper zone of the soil. Furthermore, An+1 can be used in the following relation to construct the temperature profile in the lower zone.
 | (25) |
The above equation is recovered from Eq. [5] through setting
. Equating BN+1 to zero in Eq. [25] is required to insure finite values for the thermal waves as z 
- .
Testing and Validating the Thermal Wave Model
In this section, the thermal wave model is tested and validated using two exact analytic models describing periodic heat transfer in two types of nonuniform soils. One soil type has thermal properties characterized by exponential profiles, and the other type has thermal properties characterized by linear profiles.
Exact Analytic Model for Soils with Exponential Profiles
In this type of model, the thermal properties in the upper zone of the soil are taken to increase or decrease exponentially with depth.
 | (26) |
The thermal properties of the lower zone are taken to be uniform and equal to their corresponding values at the interface separating the two zones.
The exact analytic model describing the temperature profile in the upper zone is (Hildebrand, 1976; Massman, 1993).
 | (27) |
where
 | (28) |
P and Q are two unknown constants to be determined. In addition, berµ(
), beiµ(), kerµ(), keiµ() are Kelvin functions of order µ (Hildebrand, 1976; Abramowitz and Stegun, 1972). In the lower zone, the temperature is described by a decaying exponential function similar to that given in Eq. [25].
Now we consider two exponential profiles: an increasing exponential profile described by:
; and a decreasing exponential profile described by:
, and
, and
. The values of 
, and c considered here are selected to give
. The Kelvin functions associated with such a value of µ can be expressed in terms of exponential functions (Massman, 1993). Exploiting such a property of the Kelvin function reduces Eq. [27] to [see Appendix]
 | (29) |
In Eq. [29], which is new and has not been published before, a and b are unknowns to be determined from the boundary conditions.
To use Eq. [29] for evaluating the performance of the thermal wave model, we write the frequency-dependent temperature amplitude and its first derivative as
 | (30) |
In the above |Tn
| and (z,) are the magnitude and phase of the frequency-dependent temperature amplitude. In addition |T'n| and d(z,) are the magnitude and phase of the first derivative of the frequency-dependent temperature amplitude. The predictions of the thermal model for the temperature amplitude represented by
, and for the phase represented by -(z,) are compared with the corresponding exact values predicted by Eq. [29] in Fig. 2 and 3 , respectively. Moreover, the predictions of the thermal wave model for temperature derivative represented by
, and for the phase difference |(z,) - d(z,)| are compared with the corresponding predications based on Eq. [29] in Fig. 4 and 5
. In calculating the thermal wave model predictions, the interface separating the upper and lower zones is located at
, and the upper zone is divided into N sublayers. We started with
and increased the number of sublayers until all the predictions of the thermal wave model converge to the corresponding predictions of the analytic model (Eq. [29]). We found that the convergence for the magnitude of first derivative occurs at
. As for the magnitude of the temperature profiles and its phase, and for the difference between the phase of the temperature amplitude and phase of the first derivative the convergence occurs at
. Such a higher value for the number of sublayers is attributed to the higher values of the gradient (rate of change with respect soil depth) of the thermal properties. Accounting for the scattering characteristics of the thermal waves could reduce number of sublayers at which the predictions of the thermal wave model converge to the corresponding predictions of the analytic model. This point will be investigated in the future.
We would like to mention that comparing Fig. 2 through 5 with the corresponding figures in Massman (1993) indicates that the line drawings in Fig. 2 through 5 in Massman's paper are misplaced. In addition, the right curve of Fig. 2 and 4 in Massman's paper should be revised.
Exact Analytic Model for Soils with Linear Profiles
For this type of analytic model, the thermal properties of the soil upper zone are taken to decrease or increase linearly with depth.
 | (31) |
In the lower zone the thermal properties are taken to be uniform and equal to their corresponding values at the interface separating the two zones of the soil volume.
The analytic model describing the temperature profile in the upper zone is described in terms of Kelvin function of order zero (Hildebrand, 1976; Nerpin and Chudnovskii, 1984)
 | (32) |
The algorithm for evaluating the Kelvin functions of order zero and their derivatives is taken from Section 9.11 of Abramowitz and Stegun (1972). The analytic model describing the temperature profile in the lower zone is given in terms of a decaying exponential function similar to that given in Eq. [25].
To use the predictions of Eq. [32] in evaluating the performance of the thermal wave model, two linear profiles are considered: an increasing linear profile characterized by
,
; and a decreasing linear profile characterized by
. In addition, the soil interface separating the upper and lower regions is assumed to be located at
. Furthermore, the maximum value of the first harmonic of time-dependent temperature at the soil–air interface is equated to unity, and it is assumed to occur at noon.
The predictions of the thermal wave model for the first harmonic of time-dependent temperature are compared with the corresponding predictions of Eq. [32] in Fig. 6 and 7
. Both predictions are calculated at two different times of the day: 0600 and 1200 h, for the increasing linear profile (Fig. 6) and for the decreasing linear profile (Fig. 7). In calculating the predictions of thermal wave model, the upper zone is discretized into 10 sublayers. From Fig. 6 and 7, it is obvious that predictions of the thermal wave model for the first harmonic of time-dependent temperature based on 10 sublayers are in excellent agreement with the corresponding predictions of analytic model described by Eq. [32].
View larger version (19K):
[in this window]
[in a new window]
|
Fig. 6 The time-dependent temperature for a soil having thermal properties described by the increasing linear profile of Eq. [31]
|
|
View larger version (19K):
[in this window]
[in a new window]
|
Fig. 7 The time-dependent temperature for a soil having thermal properties described by the decreasing linear profile of Eq. [31]
|
|
|
Discussion
|
|---|
The thermal wave characteristics of periodic heat transfer in layered soils were investigated and used to derive the thermal reflection (Eq. [12]) and transmission (Eq. [18]) coefficients across the interfaces of the sublayers. Those coefficients, which are new and have not been published before, were used to create a thermal wave model for constructing the temperature profiles in nonuniform soils. The model discretizes the nonuniform soils into sublayers, and instead of applying the boundary conditions at all interfaces of sublayers, the model uses the thermal reflection and transmission coefficients to construct the temperature and heat flux profiles. Using the reflection and transmission coefficients simplifies the mathematical analysis. On the other hand, if the boundary conditions are applied at all interfaces of sublayers, the mathematical formulations will be too cumbersome to be handled, as indicated by Van Wijk and Derksen (1963)(Sec. 6-19).
Numerical simulations were performed to test the model predictions against the predictions of two exact analytic models describing heat transfer in two different types of nonuniform soils, a soil type having thermal properties described with exponential profiles and a soil type having thermal properties described with linear profiles. The numerical simulations showed that the model predictions could reach the predictions of the exact analytic models. On the other hand, since the thermal wave model is not restricted to certain profiles for heat capacity or thermal conductivity, as in the case of analytic models, the thermal wave model has a wider domain of applicability than those models.
Two features of the thermal wave model could be used to enhance the performance of numerical models (e.g., those based on finite difference schemes). The first feature is the implication of an infinite thickness layer with constant thermal properties, namely, the lower zone of the soil. This represents an advantage over finite difference schemes. In most finite difference schemes one would need to place a lower boundary condition at the interface separating the two zones of the soil or at some finite depth because there is no way to represent a semi-infinite layer in those schemes. However, using Eq. [25] as a lower boundary condition in finite difference schemes could enable those schemes to handle an infinite thickness bottom layer with constant thermal properties. The other feature is incorporating the thermal wave approach into the finite difference schemes, which could reduce numerical phase error without increasing the order and hence without increasing the processing time or memory requirement. This feature has been discussed by Nehrbass et al. (1998) and there is no need to discuss it here.
Received for publication May 3, 1999.
|
Appendix
|
|---|
Analytic Model for a Special Exponential Profile
Here the analytic model associated with
will be recovered from the general analytic model described by Eq. [27]. To do so we recall the following equalities relating Kelvin functions to the modified Bessel functions (Abramowitz and Stegun, 1972)
 | (A1) |
In the above Iµ( ) and Kµ( ) are the modified Bessel functions of first and second kind, respectively. Introducing Eq. [A1] into Eq. [27] leads to
 | (A2) |
For the special case of
, the modified Bessel functions of Eq. [A2] can be written as (Abramowitz and Stegun, 1972)
 | (A3) |
Then by introducing Eq. [A3] into Eq. [A2] after setting
and with some mathematical manipulation we get
|
| (A4) |
where a and b are unknowns to be determined and
 | (A5) |
Moreover, the heat flux associated with Eq. [A4] is
 | (A6) |
|
REFERENCES
|
|---|
- Abramowitz, M., and I.A. Stegun, 1972, Handbook of Mathematical Functions, Dover Publ., New York.
- Bennett C.A., Jr., Patty R.R. Thermal wave interferometry: A potential application of the photoacoustic effect. Appl. Opt. 1982;21:49-54.
- Hildebrand F.B. Advanced calculus for applications. Englewood Cliffs, NJ: Prentice-Hall, 1976.
- Karam M.A. Molecular optics approach to electromagnetic wave interactions with stratified media. J. Opt. Soc. Am. A 1996;13:2208-2218.
- Karam, M.A. 1999. Inferring sub-surface soil moisture and hydraulic properties from microwave radiometer observations. p. 411–418. In Proc. Int. Conf. on Applied Geologic Remote Sensing, 13th. Vol. II. Vancouver, BC, Canada. 1–3 Mar. 1999. Erim Int., Ann Arbor, MI.
- Massman W.J. Periodic temperature variations in an inhomogeneous soil: A comparison of approximate and exact analytical expressions. Soil Science 1993;155:331-338.
- Nehrbass J.W., Jevtic J.O., Lee R. Reducing the phase error for finite-difference methods without increasing the order. IEEE Trans. Antennas Propag. 1998;46:1194-1201.
- Nerpin S.V., Chudnovskii A.F. Heat and mass transfer in the plant–soil–air system. New Delhi: Oxonian Press Pvt. Ltd, 1984.
- Novak M.D., Black T.A. Theoretical determination of the surface energy balance and thermal regimes of bare soils. Boundary Layer Meteorol. 1985;33:313-333.
- Novak D.N. Theoretical values of daily atmospheric and soil thermal admittances. Boundary Layer Meteorol. 1986;344:17-34.
- Novak M.D. Analytical solutions for two-dimensional soil heatflow with radiation surface boundary conditions. Soil Sci. Soc. Am. J. 1993;57:30-39.[Abstract/Free Full Text]
- Rosencwaig A., Opsal J., Smith W.L., Willenborg D.L. Detection of thermal waves through optical reflectance. Appl. Phys. Lett. 1985;46:1013-1015.
- Tang D.W., Araki N. The wave characteristics of thermal conduction in metallic films irradiated by ultra-short laser pulses. J. Phys. D. Appl. Phys. 1996;29:2527-2533.
- Van Wijk W.R., Derksen W.J. Sinusoidal temperature variation in a layered soil. In: Van Wijk W.R., ed. Physics of plant environment. Amsterdam: North-Holland Publ, 1963:170-373.
- Wait, J.R. 1975. Electromagnetic waves in stratified media. Chapter II. p. 10–13. Macmillan, New York.
TOP
ABSTRACT
INTRODUCTION
