SSSAJ Journal of Natural Resources and Life Sciences Education
HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
QUICK SEARCH:   [advanced]


     

Published online 21 June 2006
Published in Soil Sci Soc Am J 70:1281-1294 (2006)
DOI: 10.2136/sssaj2005.0293
© 2006 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
This Article
Right arrow Abstract Freely available
Right arrow Figures Only
Right arrow Full Text (PDF) Free
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Similar articles in this journal
Right arrow Similar articles in ISI Web of Science
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Right arrow reprints & permissions
Citing Articles
Right arrow Citing Articles via ISI Web of Science (1)
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Markle, J. M.
Right arrow Articles by Molson, J. W.
Right arrow Search for Related Content
PubMed
Right arrow Articles by Markle, J. M.
Right arrow Articles by Molson, J. W.
Agricola
Right arrow Articles by Markle, J. M.
Right arrow Articles by Molson, J. W.
Related Collections
Right arrow Soil Thermal Properties
Right arrow Ground Penetrating Radar, GPR
Right arrow Heat Transport

Soil Physics

Characterizing the Two-Dimensional Thermal Conductivity Distribution in a Sand and Gravel Aquifer

Jeff M. Marklea,*, Robert A. Schincariola, John H. Sassa and John W. Molsonb

a Dep. of Earth Sciences, Univ. of Western Ontario, London, ON, Canada, N6A 5B7
b Dep. of Civil, Geological and Mining Engineering, École Polytechnique Montréal, Montréal, QC, Canada, H3C 3A7

* Corresponding author (jmmarkle{at}uwo.ca)


   ABSTRACT
TOP
 ABSTRACT
INTRODUCTION
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES
 
Both hydrologic and thermal transport properties play a significant role in the movement of heat through permeable sedimentary material; however, the thermal conductivity is rarely characterized in detailed spatial resolution. As part of a study of the movement of thermal plumes through a sand and gravel aquifer, we have constructed a two-dimensional profile of thermal conductivity. This work consisted of: (i) measuring the thermal conductivity of the soil solids, {lambda}s, for the main stratigraphic units using the steady-state divided-bar apparatus and estimating conductivity from mineral composition; (ii) measuring the volumetric water content and porosity using crosshole ground-penetrating radar; (iii) evaluating four models used to predict the apparent thermal conductivity, {lambda}, of variably saturated soils and selecting the best model using the information-theoretic approach, (iv) calculating the {lambda} field on a 0.25-m square cell grid using measured data and the selected model, and (v) simulating thermal transport within the two-dimensional domain using a finite element numerical model. The apparent thermal conductivity in the saturated aquifer ranges from 2.14 to 2.69 W m–1 K–1 with a mean of 2.42 W m–1 K–1. Numerical simulations show that the heterogeneous thermal conductivity field results in increased thermal dispersion that is most pronounced at the plume front. Our values for {lambda} and {lambda}s may be used for glacial soils with similar mineralogy and texture. Our methods may also be used at other sites to construct the thermal conductivity distribution.

Abbreviations: AIC, Akaike's information criterion • AICC, Akaike's information criterion for small sample sizes • bgs, below ground surface • GPR, ground-penetrating radar • MOP, multiple-offset profile • ZOP, zero-offset profile


   INTRODUCTION
TOP
 ABSTRACT
 INTRODUCTION
MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES
 
THE temperature of discharging ground water is an important factor in determining if temperature-sensitive aquatic animals can be supported in ground water discharge areas such as river, streams, and wetlands (e.g., Garside, 1966; Acornley, 1999; Power et al., 1999). Many temperature-sensitive species have narrow ranges of thermal tolerance and even small increases in discharging ground water temperature can degrade the habitat. For example, the optimum spawning and incubation temperature for brook trout (Salvelinus fontinalis) lies between 6 and 9°C, with 50% mortality above 11.7°C (Hokanson et al., 1973). The temperature of the soil and ground water within the upper 10 to 20 m of the subsurface is controlled by the annual variation in the amount of heat transferred at the ground surface. Any disturbance that alters energy transfer at the ground surface may alter ground water temperature and adversely affect temperature-sensitive aquatic animals present in discharge areas. This work is part of an investigation into the potential thermal disturbance to ground water that may result from aggregate extraction operations in a glaciofluvial outwash sand and gravel aquifer in southwestern Ontario, Canada. The ground water from this aquifer discharges to a wetland and stream, supporting a cool-water fishery. The excavation of aggregate material below the water table involves removal of the forest cover and soil, followed by excavation of the unsaturated and saturated porous medium. Removal of the forest cover and the unsaturated porous medium increases the amount of solar radiation reaching the water table (e.g., Deardorff, 1978; Kaufmann et al., 2003; Nitoiu and Beltrami, 2005) and eliminates the insulating effect of the unsaturated zone. This results in an energy transfer across the air–water interface of the pond that is many times larger than the energy transfer across the water table under forested conditions. As a result, the temperature of the water in the pond is different from the ground water under forested conditions. In the summer months, the water in the pond is much warmer and in the winter it is colder. Under the influence of the hydraulic gradient in the surrounding aquifer, the water in these ponds moves back into the ground water system. It then moves through the aquifer as a series of alternating warm and cool thermal plumes. If this thermally altered ground water discharges before reaching background temperature, it may adversely affect the aquatic biota in the discharge area.

The movement of thermal plumes through the subsurface is controlled by the ground water velocity and by the aquifer thermal properties. Key thermal properties are the volumetric heat capacity and thermal conductivity. Since an aquifer is a granular medium consisting of solid, liquid, and gaseous phases, the volumetric heat capacity and thermal conductivity will depend on the volumetric proportions of these components. The volumetric heat capacity of an aquifer can be calculated accurately from the heat capacities and volume fractions of these three phases (Smith, 1939, 1942; Woodside and Messmer, 1961; de Vries, 1963). The apparent thermal conductivity, {lambda}, is more complicated to calculate. It depends mainly on the mineral composition of the aquifer solids, and the porosity and degree of saturation. To a lesser extent, it depends on the bulk density of the aquifer solids, the shapes, sizes, and arrangement of the solid particles, the contact area between the particles, the interfacial contact between the solid and liquid phases, the vapor diffusion in the unsaturated pores, and the temperature and pressure conditions (Smith 1939, 1942; de Vries, 1963; Hopmans and Dane, 1986).

There are several methods available for estimating the apparent thermal conductivity of unconsolidated porous media. The most common methods include the direct measurement of conductivity using probes, and the estimation of conductivity using either empirical or mixing models. In situ transient line source probes have been used successfully in fine-textured porous media to measure thermal conductivity (Lubimova et al., 1961; Sass et al., 1981; Bristow et al., 1994); however, none of the currently available probes are durable enough for in situ measurements in coarse-textured media with cobbles and boulders. Measurement of thermal conductivity in these materials requires the use of alternate methods.

Predicting thermal transport through the subsurface is often accomplished with numerical finite difference or finite element models (e.g., Andrews and Anderson, 1979; Molson et al., 1992). As input to numerical simulations, these models require values of the apparent thermal conductivity for a variety of porous media across saturation conditions that range from nearly dry to fully saturated. If the thermal conductivity for the individual components is known, values of the apparent conductivity can be calculated using mixing models (e.g., de Vries, 1963; Gori, 1983; Campbell et al., 1994), a number of which were evaluated in this study.

Our goal was to simulate the migration of a thermal plume (emanating from a nearby aggregate pit) through the shallow aquifer using a finite element numerical model. As input to the model, we required values of thermal conductivity for the glaciofluvial outwash sand and gravel aquifer. The main objectives of this study were to characterize the two-dimensional distribution of the apparent thermal conductivity in the aquifer, to evaluate the suitability of four candidate models for calculating the thermal conductivity, and to assess the influence of heterogeneous thermal conductivity on heat transport using numerical simulations.


   MATERIALS AND METHODS
TOP
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES
 
Site Description
The study area is located in the Tricks Creek watershed of southwestern Ontario, ~180 km west of Toronto (Fig. 1). The watershed is characterized by undulating topography. Before being cleared for agriculture, the area was covered by mixed deciduous forest. Presently, 12% of the watershed is forested. Tricks Creek lies within a wetland complex that encompasses an area of ~105 ha (4% of the watershed). Tricks Creek is characterized by cool water and supports resident brook and rainbow trout (Salmo Gairdneri) populations. The creek and wetland are situated in a former glacial outwash channel in which the upper 6 m of the subsurface consist of glaciofluvial outwash deposits of sands and gravels. The outwash material overlies 30 m of silty clay till. The outwash sands and gravels were deposited in a meltwater channel at the ice margin in the last retreat of the ice sheet during the Wisconsinan glaciation, which occurred in this area approximately 13000 yr ago (Barnett, 1992). The sands and gravels are mixtures of predominantly carbonate and quartz minerals, and form an unconfined aquifer. Within the study area, the ground water flows to the southwest and discharges to the wetland and creek (Fig. 1). Several aggregate operations are active along the western edge of the wetland and are upgradient of the creek.


Figure 1
View larger version (49K):
[in this window]
[in a new window]
  		Fig. 1. Site location map.

 
Field and Laboratory Methods
Thermal Conductivity of Soil Solids
The presence of gravel, cobbles, and boulders in the outwash aquifer made the use of in situ probes impractical. For this investigation, we determined that it was more practical to recover aquifer material during drilling and complete measurements of thermal conductivity in the laboratory. Aquifer samples were collected using a truck-mounted drill rig equipped with hollow stem augers and a split-barrel sampler. The tip of the sampler preceded the augers during drilling and a PVC (polyvinyl chloride) sleeve inside the core barrel provided for the retrieval of aquifer cores, 0.126 m in diameter and 1.52 m long, with minimal disturbance.

Within the study area, the soils are of the Humo-Ferric Podzol great group (Agriculture and Agri-Food Canada, 1998). These have developed beneath mixed deciduous forest cover with undulating terrain under well-drained conditions. In the area of the boreholes, the A horizon and much of the B horizon have been excavated in preparation for aggregate extraction. Based on observations during drilling, all that remains is ~0.05 m of the B horizon at the surface. Below the soil lies 6 m of parent geologic material composed of carbonate-rich, glaciofluvial outwash sand and gravel. The outwash can be subdivided into four stratigraphic units, based on particle size distribution. These include poorly sorted gravel with sand, and well-sorted coarse, medium, and fine sand. The outwash is underlain by glacial till, which was the lowest stratigraphic unit encountered during drilling. From each unit, we selected representative samples and measured the thermal conductivity of the aquifer solids, {lambda}s, on 41 of these samples. The mass of each sample was between 2 and 4 kg. We ground the samples to a particle size of <1 mm (Sass et al., 1971) and measured the thermal conductivity (at 20°C) using the steady-state divided-bar apparatus. This method involves filling a cylindrical cell with crushed material, saturating the sample with water under vacuum, and measuring the conductivity of the cell in the same manner as a cylinder of solid rock. The conductivity of the solid component is then backed out from the geometric mean of the water and solid mixture. Additional details of the apparatus and method verification can be found in Sass et al. (1971). Values measured with this apparatus are generally accurate to within ±5%. From a subset, 27 samples were ground to a particle size of 20 µm, and the mineralogy was obtained using the XRD (x-ray diffraction) technique. The semiquantitative estimation of relative mineral abundance was based on the integrated peak areas after the removal of the background response. This method yields estimates that are within 15% for the clay minerals and 5% for the other minerals (Mitchell, 1976). Both the divided-bar and XRD methods involve crushing and mixing the aquifer solids. As a result, the derived value of thermal conductivity represents a bulk value for the solid fraction of the aquifer. These methods do not provide information on the thermal conductivity or the mineralogy of individual particles. Furthermore, they do not account for the influence of particle size and shape.

We selected the divided-bar apparatus to measure {lambda}s and the mineral composition to estimate {lambda}s, since these two methods are commonly used for consolidated rocks, and the porous medium mineral composition is often known. Furthermore, most predictive models use a bulk thermal conductivity for the porous media solids. While the model proposed by de Vries (1963) can account for a porous medium having particles of different thermal conductivity, mineral composition, and shape, the majority of the results reported in the literature use a bulk thermal conductivity, and a common mineral composition and shape for the particles. The exception to this is the work by Tarnawski et al. (2000), in which the porous medium is modeled as a mixture of 20 unique types of particles having unique thermal conductivities, mineral compositions, and shapes. For most studies, however, particle shape and mineralogy are not known.

Volumetric Water Content and Porosity
In our study, the volumetric water content, {theta}, of the aquifer was estimated from a series of crosshole GPR (ground-penetrating radar) surveys. The surveys were completed across six boreholes, which span a 12.3- by 7.6-m portion of the aquifer, using a Sensors and Software (Mississauga, ON, Canada) pulseEKKO 100 GPR system equipped with borehole antennas (Fig. 1). We conducted three ZOP (zero-offset profile) surveys using antennas with a center frequency of 100 MHz, and five MOP (multiple-offset profile) surveys using 200 MHz antennas. We began all the ZOP surveys at the ground surface and proceeded to the bottom of the borehole with a step size of 0.125 m. For all the MOP surveys, we began just below the water table and proceeded to the bottom of the borehole with a step size of 0.25 m. For each GPR trace, we picked the arrival time of the direct wave from which the electromagnetic wave velocity was extracted. The value of {theta} was calculated from the velocities using the BHS (Bruggeman–Hanai–Sen) mixing formula (Sen et al., 1981; Feng and Sen, 1985).

The interwell velocity structure, measured with ZOPs and MOPs, is different as a result of the different geometrical configuration used for each of these surveys. The ZOPs have only horizontal ray paths collected between each station in the boreholes. Since the travel time measured for each station reflects the average of the electromagnetic wave velocity between the boreholes, only vertical variations in these average values are measured by ZOPs. In contrast, MOPs have ray paths at many different angles between the boreholes. Inversion of travel times for these ray paths yields both vertical and horizontal variations in interwell velocity. Since only ZOPs were completed in the unsaturated zone, only the vertical variation in the velocity and water content could be measured. In the saturated zone, we completed ZOPs and MOPs. The travel times from these were inverted to reconstruct the horizontal and vertical interwell velocity structure using the tomographic inversion code Pronto (Aldridge and Oldenburg, 1993). The inversions were performed with the domain divided into 0.25-m square cells. In the saturated zone, the water content is equal to the porosity, {phi}, and the BHS equation provides a direct estimate of the soil porosity. In the unsaturated zone, the moisture content is given by {theta} = {phi}Sr, where Sr is the degree of saturation, which ranges between 0 and 1. In this zone, the BHS equation provides an estimate of the water content only; therefore, we measured porosity directly on cores recovered from the unsaturated zone. The porosity was estimated from the difference in mass between saturated and oven-dried samples.

Predictive Models
We evaluated four predictive models: one empirical model (Johansen, 1975) and three mixing models (de Vries, 1963; Gori, 1983; Campbell et al., 1994). Each of these models may be used to predict the apparent thermal conductivity in saturated and unsaturated porous media under variable temperature conditions. The empirical model by Johansen (1975) uses a form of interpolation between the apparent thermal conductivity of dry and saturated sediments. The mixing models by de Vries (1963) and Campbell et al. (1994) are based on the analog to the Maxwell model for the electrical conductivity of a mixture of spheres dispersed in a continuous fluid. The mixing model by Gori (1983) models the porous medium as a cubic space with a cubic centered solid grain surrounded by a mixture of air and water. Only the basic equations for the four models are presented below.

The apparent thermal conductivity of an unsaturated porous medium is given by Johansen (1975) as

Formula 1[1]
where {lambda}sat is the thermal conductivity of the saturated porous medium, {lambda}dry is the thermal conductivity of the dry porous medium, and Ke is the Kersten number. The thermal conductivity of the saturated porous medium is

Formula 2[2]
where {lambda}s is the thermal conductivity of the solids and {lambda}w is the thermal conductivity of water. The thermal conductivity of a dry, coarse porous medium is given by

Formula 3[3]
where {rho}d is the dry bulk density (kg m–3) and {rho}s is the density of the solids (kg m–3). The form of Ke given by Johansen (1975) applies only when Sr > 0.05. Below this level, it underestimates the value of the thermal conductivity and alternate models must be used (Farouki, 1981, 1982); however, the use of different models produces discontinuities at the transition points and is cumbersome to implement. To overcome these problems, we implemented the following form of the Kersten number (Ewen and Thomas, 1987):

Formula 4[4]
where ß is a fitting parameter. While this form of Ke provides a continuous equation that applies across the full range of saturation, Eq. [1] and [2] do not yield the same {lambda} at full saturation when ß > –4.5. We eliminated this discrepancy by modifying Eq. [1] as follows

Formula 5[5]
Here {lambda}Sr=1 is evaluated using the unmodified form of Eq. [1].

Of much greater complexity than the Johansen empirical model are the mixing models. In the mixing model by de Vries (1963), the apparent thermal conductivity is calculated using

Formula 6[6]
where xi is the volume fraction of each constituent (air, water, or soil particle or mineral fraction), {lambda}i is the thermal conductivity of each constituent, and n is the number of soil constituents. The weighting factor ki is the ratio of the average temperature gradient in the ith component in the soil to the temperature gradient in the continuous medium and is related to the shape and conductivity of the component. All components with the same shape and conductivity are considered as one type and have common {lambda}is and kis. The subscript zero applies to the continuous medium surrounding the soil particles, which for dry soils is air and for moist to saturated soils is water. For the continuous medium k0 = 1, and the remaining kis are given by

Formula 7[7]
where gj are the shape factors for the ith component and {lambda}0 is the thermal conductivity of the continuous phase. The quantities gj depend on the ratio of the major axes of the ellipsoid for the soil component, and g1 + g2 + g3 = 1. Most soil particles are spheroids having g1 = g2 = mg3, where m varies from 0.1 to 100 (de Vries, 1963). Thus only one shape factor must be estimated for each component.

In unsaturated porous media, the temperature gradients cause moisture movement across the air-filled pores, which redistribute the heat across the pores. This can be described by an apparent thermal conductivity of the air-filled pores due to heat transport by conduction through dry air {lambda}a, and by the movement of vapor {lambda}vs in the pores containing moist air at a relative humidity h. Thus the apparent thermal conductivity is

Formula 8[8]
There are several different expressions for h and {lambda}vs (e.g., de Vries, 1963; Hopmans and Dane, 1986; Campbell et al., 1994; Tarnawski et al., 2000).

The most complicated aspect of implementing the de Vries model is the evaluation of gj used in Eq. [7] for the air pore-shape factors in unsaturated porous media. These shape factors are dependent on the water content and a transition occurs at the field capacity of the soil. De Vries (1963) and Hopmans and Dane (1986) provide detailed descriptions of the procedure required to evaluate gj, and we followed these procedures in our implementation of the de Vries model. Less complex implementations are available (de Vries, 1963; Farouki, 1982), but these are derived for the quartz sand considered in the work by de Vries (1963) and they may not be applicable to other soils.

To reduce the complexity of the de Vries model, Campbell et al. (1994) introduced a continuous function for the kis, which applies across the full range of water contents, and then used gj as an empirical fitting factor. While Campbell's modified form of the de Vries model is easier to implement, it introduces two new parameters, qo and xwo. The parameter xwo is the cutoff water content for liquid recirculation and gives the water content at which water starts to affect thermal conductivity. It can be calculated using the relationship for xwo (m3 m–3) given by Campbell et al. (1994) as

Formula 9[9]
where dg is the geometric mean particle diameter (µm) (Shiozawa and Campbell, 1991). The parameter qo relates to the rapidity of the transition from air- to water-dominated conductivity and is treated as a fitting parameter (Campbell et al., 1994). Additional details of the Campbell model and the continuous function for the kis are given in Campbell et al. (1994).

Gori (1983) developed a model based on a cubic grain inside a cubic space for unsaturated frozen porous media. This model has been adapted to consider latent heat transfer in unfrozen soils (Tarnawski et al., 2000). The Gori model was shown to provide good agreement with measured values of thermal conductivity for unsaturated soils at temperatures of 30 and 50°C (Tarnawski et al., 2000). The equations for this model are quite complex and are not presented here but they can be found in Tarnawski et al. (2000).

Model Evaluation and Selection
Our objectives for model evaluation and selection were to compare the apparent thermal conductivities predicted with the models to existing data, and to select the model that best represents a balance between bias (underfitting data with models having few parameters) and variance (overfitting data with models having many parameters). This was achieved by compiling applicable datasets of measured thermal conductivity from the literature, and using AIC (Akaike's information criterion; Akaike, 1973). The AIC is an information-theoretic procedure, based on Kullback–Leibler information theory, and it provides a method for objective model selection. When n/K < 40, where n is the sample size and K is the number of estimated parameters in the model, the AIC for small sample size, AICC (Table 1), is used (Burnham and Anderson, 2002). For all datasets in the literature, n/K was <40. The first term in the equation for AICC is a measure of the lack of fit. This term gets smaller as more parameters are added to the model to improve the fit to the data. As more parameters are introduced into the model, the remaining terms in the AICC equation get larger (a penalty for adding more parameters) and parsimony is enforced. As the sample size n increases these terms get smaller. Thus AIC provides a rigorous way to achieve a model of appropriate complexity for a dataset with a given sample size. Burnham and Anderson (2002) described the theoretical foundations of AIC and its application for model ranking, selection, and inference.


View this table:
[in this window]
[in a new window]
  		Table 1. Criteria for evaluating the fit of the apparent thermal conductivity, estimated with the candidate models, to the measured thermal conductivity.

 
Two statistics were used to measure the goodness-of-fit of predicted thermal conductivity against measured conductivity. These included the correlation coefficient r and the AICC (Table 1). In practice, one computes the criterion AICC for each model and selects the model with the smallest value. Two additional parameters, {Delta}i and wi, may be calculated from AICC values. The {Delta}i allow an easy ranking of the models from best to worst ({Delta} = 0 for the best model). In general, models having {Delta}i ≤ 2 are very good models, those for which 4 ≤ {Delta}i ≤ 7 have less support, and models having {Delta}i ≥ 10 can be eliminated as candidates (Burnham and Anderson, 2002; Anderson, 2003). The wi, called Akaike weights, are considered as the likelihood or weight of evidence in favor of model i being the best model for the situation being considered. These likelihoods are normalized and can be treated as probabilities. In addition, the ratio wi/wj gives the relative likelihood of model i vs. model j and is termed the evidence ratio. The evidence ratio allows us to state that there is wi/wj times more support for model i than model j (Poeter and Anderson, 2005).


   RESULTS AND DISCUSSION
TOP
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
CONCLUSIONS
 REFERENCES
 
Thermal Conductivity of Aquifer Solids
Divided-Bar Method
The thermal conductivity values measured using the divided-bar method are normally distributed and range between 3.38 and 4.81 W m–1 K–1. The mean value and 95% confidence interval are 4.08 ± 0.09 W m–1 K–1 and the standard deviation is 0.29 W m–1 K–1 (Table 2). A box-whisker plot of measured thermal conductivity values suggests that {lambda}s may be correlated with the stratigraphic units (Fig. 2). To test this hypothesis, we conducted one-way unbalanced ANOVA. We excluded the data for the soil (S) from the analysis since there are only two samples. The Brown–Forsythe modification of Levene's test for the homogeneity of variance (Brown and Forsythe, 1974) indicated that the variances could be assumed equal [F(4,34) = 2.035, p ≤ 0.112].


View this table:
[in this window]
[in a new window]
  		Table 2. Thermal conductivity of the porous media solids ({lambda}s) measured by the divided-bar apparatus and estimated from the mineral composition measured by x-ray diffraction.

 

Figure 2
View larger version (27K):
[in this window]
[in a new window]
  		Fig. 2. Box-whisker plot of the measured thermal conductivity for the solid fraction of porous media grouped by stratigraphic unit. The caps at the end of each box indicate the minimum and maximum values, the box is defined by the lower and upper quartiles (25th and 75th percentiles), and the line in the center of the box is the median. No outliers were present in the data.

 
A significant difference between the mean thermal conductivity of the units was observed (Table 3). Post hoc MCT (multiple comparison tests) using the methods of Tukey, Scheffé, and Bonferroni (StatSoft, 2001), with an overall error rate of {alpha} = 0.05, suggested that the mean thermal conductivities of till (T), gravel and sand (G), and coarse sand (CS) were less than those of medium sand (MS) and fine sand (FS), and that T ~ G ~ CS < MS ~ FS. A nonparametric MCT (Conover, 1999) suggested T ~ G < CS ~ MS ~ FS, based on the median thermal conductivities. Given the difficulty of distinguishing between coarse sand and medium sand in the field, it is more practical to group all sand units together. For till, only four measurements were made and the mean has a large standard deviation. We placed till in a separate group since additional measurements would probably decrease the standard deviation, and the MCTs would indeed identify till as a separate group. The basic statistical parameters for each of the grouped units are summarized in Table 4.


View this table:
[in this window]
[in a new window]
  		Table 3. The ANOVA table (one-way unbalanced analysis) shows the between-groups and within-groups sources of variance for the thermal conductivity data. The null hypothesis, H0, was rejected at p < 0.05.

 

View this table:
[in this window]
[in a new window]
  		Table 4. Mean value ± 95% confidence interval, standard deviation, and upper and lower quartiles of the measured thermal conductivity of the porous media solids for individual stratigraphic units and for the grouped sand units.

 
Mineral Composition Method
The mineral compositions for 27 outwash samples were primarily calcite, dolomite, quartz, and plagioclase feldspar (Table 5). Some samples contained hornblende, illite, montmorillonite, and chlorite or possibly kaolinite in minor quantities (<5% total). From the known mineral compositions, we estimated the bulk thermal conductivity for the aquifer solids {lambda}s (Table 2), having n mineral components with a volume fraction xi and conductivity {lambda}i, using the geometric mean equation given by

Formula 10[10]
A wide range of values are cited in the literature for several rock-forming minerals in their monomineralic and rock form (Table 6); we used the underlined values in our calculations. In general, the thermal conductivities given for minerals in their monomineralic form are higher than those reported for rocks. This difference may be due to the presence of other minerals in the rocks tested, intragranular porosity in the rocks that reduces the thermal conductivity, or other factors. For calcite and dolomite, we selected the values of thermal conductivity reported for rocks, as these will account for the intragranular porosity common in these rocks.


View this table:
[in this window]
[in a new window]
  		Table 5. Mineral composition of the porous media measured by x-ray diffraction and calculated using the integrated peak area method.

 

View this table:
[in this window]
[in a new window]
  		Table 6. Thermal conductivity of selected minerals and rocks. The underlined values are considered to be representative of the minerals in the glaciofluvial outwash sands and gravels.

 
Comparison of these calculated conductivities to those obtained using the divided-bar apparatus shows reasonable agreement, with the exception of the values for till. For the gravels and sands, the predicted and measured values are generally within ±10%, and have a correlation coefficient r of 0.473 that is significant at p < 0.05. For the till samples, the thermal conductivities estimated from the mineral composition are ~20% larger than the measured values. This difference may be due to incomplete saturation of the till samples when measured using the divided-bar method. Air entrapped in the pore space will decrease the apparent thermal conductivity of the sample. Therefore, with the exception of the till, our results suggest that {lambda}s for the glacial outwash can be measured using the divided-bar apparatus or estimated using the geometric mean equation and the mineral composition of the aquifer solids.

Volumetric Water Content and Porosity
The two-dimensional distribution of volumetric water content in the saturated zone, measured with the GPR MOPs, shows three distinct layers with different water contents (Fig. 3a). Between ~1.5 and 4 m bgs (below ground surface), {theta} varies from 0.23 to 0.30 m3 m–3. Between 4 and 5.6 m bgs, {theta} increases to 0.35 m3 m–3, and below 5.6 m bgs, {theta} is 0.32 m3 m–3. Comparison of the variation in the water content (Fig. 3a) to the geological cross-section (Fig. 4a) suggests that {theta} is not directly correlated to the stratigraphic units. We speculate that these differences in {theta} may be related to variations in the depositional environment during deposition of the outwash sediments. Thus the GPR provides an important direct measurement of the aquifer water content.


Figure 3
View larger version (90K):
[in this window]
[in a new window]
  		Fig. 3. (a) Two-dimensional volumetric water content tomogram for the saturated zone, calculated from the interwell velocity tomogram and the Burggeman–Hanai–Sen mixing formula (Sen et al., 1981; Feng and Sen, 1985). The white Xs indicate the transmitter and receiver station locations. (b) Measured water content ({theta}) variation vs. depth between the two boreholes on the right side of the section. In the saturated zone, the water content is equal to the porosity ({phi}). Above the water table, the porosity was measured in the laboratory on cores collected during drilling. The locations and measured values of the porosity are indicated by filled circles.

 

Figure 4
View larger version (45K):
[in this window]
[in a new window]
  		Fig. 4. (a) Geological cross-section of the area over which the GPR (ground-penetrating radar) survey was completed. The three major stratigraphic units shown are gravel and sand, sand, and till. The vertical lines are the locations of the boreholes where core samples were collected during drilling and where the crosshole GPR survey was conducted. The unsaturated zone is ~1.5 m thick. (b) Envelope of the annual temperature variation (minimum observed temperature on the left and maximum on the right) and the temperature profile on 1 July for this section of the aquifer.

 
The volumetric water content in the saturated and unsaturated zones, measured with the ZOPs, shows the large contrast between these zones (Fig. 3b). In the unsaturated zone, the water content decreases rapidly from 0.25 m3 m–3 at the water table to ~0.07 m3 m–3 above the capillary fringe. In the saturated zone, the variations in the water content that are evident in the ZOP, such as the higher water content layer between 4 and 5.6 m bgs, span the width of the tomogram (Fig. 3a). Porosity values in the unsaturated zone are sparse laterally and provide information on only the vertical variation. Porosity increases from 0.25 m3 m–3 at the water table to 0.4 m3 m–3 near the ground surface.

Evaluation of Candidate Models
We chose six datasets that were representative of the types of sediments found at our site, and that had all the input data required by each of the candidate models (Table 7). We did not evaluate the models for gravel, as a dataset with sufficient information was not available. For most datasets, the thermal conductivity measurements are reported at a single temperature that is generally between 20 and 26°C. For three soils we considered, measurements at more than one temperature are available: the quartz sand (de Vries, 1963); the loamy sand (Sepaskhah and Boersma, 1979); and the L-Soil (Campbell et al., 1994). While Campbell et al. (1994) completed measurement on 10 soils, all but the L-Soil were finer grained than those at our site and were not considered. Hopmans and Dane (1986) measured the thermal conductivity of a Norfolk sandy loam at four different temperatures; however, there were too few measurements in this dataset for use in our analysis.


View this table:
[in this window]
[in a new window]
  		Table 7. Parameter values input into the candidate models used to estimate the apparent thermal conductivity of the soil.

 
To evaluate the models, we compared the predicted apparent thermal conductivity with the measured conductivity for each soil across moisture conditions ranging from dry to saturated soil. In our analysis, we treated ß (Johansen's model), g1, qo, and xwo (Campbell's model), g1, fc and xc (de Vries' model), {theta}aw and xc (Gori's model), and {lambda}s (all models) as fitting parameters. The best fit for each model to the soil data (Table 7) was obtained using the parameter estimation techniques in PEST, Version 9.01 (Doherty, 2005).

The correlation coefficients were generally >0.95 (Table 8), indicating reasonable fits to the data with all the models. The AICC values indicate that Campbell's model is the best approximating model for five of the datasets (quartz sand at 40°C, loamy sand at 25 and 45°C, Sandfly Creek sand, and Tottori dune sand), Johansen's model is the best model for two datasets (quartz sand at 20°C and Leighton Buzzard sand), and de Vries' model and Gori's model are each the best for one dataset (L-soil at 30°C and L-soil at 50°C, respectively). For many of the datasets, there is no competitor to the top-ranked model ({Delta}i > 10). The exceptions, for which there is a competitor (4 < {Delta}i < 7), are the quartz sand at 40°C, the L-soil at 30°C, and the Tottori dune sand. For the quartz sand at 40°C, the weight for the top-ranked model (Campbell's model) is 0.99, the weight for the second-ranked model (de Vries' model) is 0.072, and the evidence ratio is 14 (there is 14 times more support for Campbell's model). For the Tottori dune sand, the weight for the top-ranked model (Campbell's model) is 0.84, the weight for the second-ranked model (Johansen's model) is 0.16, and the evidence ratio is 5.1. In these two cases, there is strong support for Campbell's model. For the L-soil at 30°C, the Campbell model has similar support to the de Vries model ({Delta}i = 0.8). Both models have similar Akaike weights and the evidence ratio for the de Vries model vs. the Campbell model is only 1.5. This is not strong evidence that the de Vries model is the best model. Since Campbell's model is the AICC-selected model for five of the nine datasets, and a strong competitor for another, we chose it as the "best-approximating model."


View this table:
[in this window]
[in a new window]
  		Table 8. Summary of Akaike Information Criterion (AICC) measures between the apparent thermal conductivity predicted by the candidate models and the measured values reported in the literature.

 
Model Predicted Apparent Thermal Conductivities
Using Campbell's model, we calculated the thermal conductivity for the glacial outwash aquifer. The three main units considered were sand, gravel and sand, and till (Fig. 4a). The unsaturated zone is approximately 1.5 m thick. For the simulations, g1 and qo were set equal to 0.100 and 4.0 for all soils. These are approximately the average values reported in Table 7. We investigated the implications of selecting these values through a sensitivity analysis with the sensitivity coefficient SC given by

Formula 11[11]
where Ob is the model outcome for the parameter base value, O is the model outcome for an alternate parameter value, pb is the parameter base value, and p is an alternate parameter value. Across the range of fitted values of g1 for the Campbell model (Table 7), the sensitivity coefficient ranged from –0.5 to –4.7 (model-predicted {lambda} decreased as g1 increased). If g1 for the outwash actually lies near the extremes of this range, then our predicted values of {lambda} may vary by ±20% for dry porous media and ±3% for saturated porous media. The model is much less sensitive to qo, with SC ranging from 0 to –0.18. For qo, {lambda} varied by ±5% in the unsaturated aquifer and less than ±0.5% in the saturated aquifer across the range of fitted values.

Using Eq. [9] and grain size analyses from 60 sediment samples, we determined xwo to be equal to 0.04 m3 m–3 for gravel, 0.05 m3 m–3 for sand, and 0.25 m3 m–3 for till. The thermal conductivity values assigned to the aquifer solids were taken from Table 4. In the saturated zone, the measured values of {phi} and {theta} (Fig. 3a) were used with the exception of the bottom portion of the section occupied by till. Since the GPR did not extend to this depth, we set {phi} and {theta} equal to 0.4 m3 m–3, which is typical of tills in this area. In the unsaturated zone, we used the laboratory-measured values of the porosity and average values of {theta} from the ZOPs. Uniform values of {phi} and {theta} were assigned laterally across the unsaturated portion of the aquifer, and we used the temperature profile measured on 1 July.

In the datasets we used for our model selection, conductivity values were all measured at temperatures between 20 and 50°C. Throughout most of the subsurface profile at our site, the annual variation in temperature is below this range (Fig. 4b). Suitable datasets, with measurements <20°C, are not available and we were unable to evaluate the candidate models at these lower temperatures; however, until data are available with which to evaluate these models at temperatures <20°C, we have no reason to reject the Campbell model as the "best-approximating model" for our analysis.

The two-dimensional distribution of the apparent thermal conductivity of the glacial outwash sand and gravel for 1 July was calculated using the above approach (Fig. 5). The thermal conductivity of the saturated gravel and sand aquifer ranges from 2.14 to 2.69 W m–1 K–1, while {lambda} in the till is 1.90 W m–1 K–1. In the saturated zone, two distinct layers of different conductivity are evident. The upper layer extends from the water table down to 4 m bgs. In this layer, the mean and standard deviation of the apparent thermal conductivity, {lambda} and s, are 2.53 and 0.08 W m–1 K–1 for the gravel and sand. The lower layer extends from 4 m bgs to the base of the section. In this layer, {lambda} and s are 2.30 and 0.06 W m–1 K–1, for the gravel and sand, and 2.41 and 0.07 W m–1 K–1 for the sand. The 10% difference in {lambda} for the gravel and sand unit between the upper and lower layers is due to the difference in porosity. The influence of {lambda}s on the apparent thermal conductivity is visible in the bottom portion of the saturated zone, where Formula 11 for the sand layers is 5% higher than for the gravel and sand layers. These values are in reasonable agreement with those obtained for aquifers of similar texture by Andrews and Anderson (1979) of 2.13 W m–1 K–1 for very fine to fine sand and 1.88 W m–1 K–1 for medium to coarse sand with gravel; by Palmer et al. (1992) of 2.1 ± 0.3 W m–1K–1 for Borden sand; and by Parr et al. (1983) of 2.29 ± 0.19 W m–1K–1 for sand and gravel. Given the relatively narrow range of {lambda} that we found, it may be sufficient to use literature-cited values of thermal conductivity and porosity for many investigations. As will be shown in the simulations below, however, for applications where small temperature differences are important, even the 5 to 10% difference in {lambda} measured in this aquifer may influence heat transport significantly and require detailed measurements of the aquifer properties.


Figure 5
View larger version (114K):
[in this window]
[in a new window]
  		Fig. 5. The two-dimensional distribution of the apparent thermal conductivity, {lambda} (W m–1 K–1) for the glaciofluvial outwash sand and gravel aquifer as calculated using the Campbell et al. (1994) model.

 
In the unsaturated zone, {lambda} drops from 2.6 W m–1K–1 at the water table to 1.4 W m–1K–1 at the ground surface. These values of {lambda} are 40 to 50% lower than in the saturated zone due to the lower water content above the capillary fringe. It must be recognized that, in areas where the subsurface temperatures are higher than those found at our site, the influence of heat transport by vapor diffusion will increase the apparent thermal conductivity in the unsaturated zone thereby decreasing the contrast between the saturated and unsaturated zones.

Numerical Simulations of Heat Transport
Using the measured thermal conductivity field, we investigated the influence of the heterogeneous thermal conductivity on heat transport within the two-dimensional study section. Simulations were completed using the finite element numerical model Heatflow (Molson et al., 1992) after modifications were made to include the Campbell model for apparent thermal conductivity. The Heatflow model accounts for density-dependent groundwater flow, thermal advection, conduction through the porous medium, thermal buoyancy, and thermal retardation. For all simulations, we assumed a uniform hydraulic conductivity of 5.4 x 10–4 m s–1 (at 10°C) so that heat transport would be influenced by only heterogeneity of the aquifer thermal properties. We expect that heterogeneities in the aquifer hydraulic conductivity will increase dispersion. A hydraulic gradient of 0.01 m m–1 was applied across the section using constant heads at the lateral boundaries. The average measured porosity in the saturated zone (Fig. 3a) was 0.30, giving a mean groundwater velocity of 1.8 x 10–5 m s–1 (1.55 m d–1) at 10°C. The top and bottom boundaries were assumed impermeable to flow. For the thermal transport simulations, a uniform fixed temperature of 10°C was assigned along the left boundary, except between 4.5 and 5.5 m below ground surface, where the temperature was set to 30°C. All other boundaries were assigned a temperature gradient of zero. Throughout the domain, the initial temperature was set to 10°C and hydrodynamic dispersivities were set to zero. Simulations were completed using a grid of 221 by 133 nodes in the x and z directions, respectively, and using time steps of 0.02 d. The temperature convergence criterion was 0.001°C.

To investigate the influence of the heterogeneous thermal conductivity field on heat transport, we completed three sets of simulations. Each set comprised one simulation using a heterogeneous thermal conductivity field and one simulation using a homogeneous or uniform field. The differences in the temperatures between these two simulations were calculated at a time of 10 d. For the homogeneous fields, thermal conductivities in the saturated and unsaturated zones were set equal to the geometric mean of the conductivities in the corresponding portions of the heterogeneous field. For the first set of simulations, the apparent thermal conductivities in the observed heterogeneous field (Fig. 5) were decreased by 10%, and the thermal conductivities in the saturated and unsaturated portion of the homogeneous field were 2.16 and 1.57 W m–1 K–1, respectively. For the second set of simulations, the observed {lambda} field was used, and conductivities in the saturated and unsaturated portion of the homogeneous field were 2.40 and 1.73 W m–1 K–1, respectively. Finally, for the third set, thermal conductivities from the observed field were increased by 10%, and conductivities in the saturated and unsaturated portion of the homogeneous field were 2.64 and 1.89 W m–1 K–1, respectively.

Temperatures predicted using the heterogeneous {lambda} fields (Fig. 6a, 6b, and 6c) were compared with the homogeneous {lambda} fields and differences were calculated (Fig. 6d, 6e, and 6f). For the heterogeneous thermal conductivity field observed at our site, the plume front, defined by the 11°C contour, had migrated approximately 9.5 m across the model domain after 10 d (Fig. 6b). The maximum temperature difference between the plumes in the heterogeneous and homogeneous fields of –0.36°C was centered at 8 m, near the front of the plume (Fig. 6e). In this area, temperatures in the heterogeneous field plume were lower than in the uniform field plume due to increased thermal dispersion through the heterogeneous thermal conductivity field. While reducing the core temperatures, this increased thermal dispersion also increased the temperature, in the area directly above the plume, relative to the uniform field plume, thereby producing a positive temperature difference "halo." This halo was focused above the plume, probably due to thermal buoyancy, which causes the simulated plume to rise.


Figure 6
View larger version (44K):
[in this window]
[in a new window]
  		Fig. 6. Simulated thermal plumes (left column) and corresponding temperature differences (right column) at 10 d using various thermal conductivity distributions. The thermal plumes on the left are shown using (a) a 10% decrease in the observed field, (b) the observed field, and (c) a 10% increase in the observed field. The corresponding temperature differences on the right were obtained by subtracting the plume temperatures simulated in the heterogeneous fields (a, b, and c) from plume temperatures simulated using equivalent mean thermal conductivity fields in the saturated and unsaturated zones of (d) 2.16 and 1.57 W m–1 K–1, (e) 2.40 and 1.73 W m–1 K–1, (f) and 2.64 and 1.89 W m–1 K–1, respectively.

 
When thermal conductivities in the heterogeneous field were decreased by 10% and the thermal plume (Fig. 6a) was compared to a uniform field case with a mean thermal conductivity in the saturated zone of 2.16 W m–1 K–1, the maximum temperature difference increased slightly to –0.42°C (Fig. 6d). The positive temperature difference "halo" concurrently became smaller. The thermal plume in the uniform field dispersed less due to lower rates of heat conduction, thereby maintaining higher temperatures in the plume core. For lower thermal conductivities, the influence of the heterogeneities on thermal dispersion was therefore enhanced. In contrast, when thermal conductivity was increased by 10%, dispersion of the thermal plume in the uniform field increased due to higher heat conduction rates. The influence of the heterogeneities on the thermal dispersion was decreased relative to the uniform case. The maximum temperature difference in this case was only –0.33°C, and the positive "halo" increased in area (Fig. 6f).

From other simulations, not shown here, we found that, as the plume transport distance increased, temperature differences increased in the frontal region of the plume. This is analogous to solute transport, where the size of the dispersion or mixing zone increases as the advective front moves farther from the source. While we considered only a two-dimensional system, we would expect increased thermal dispersion in a three-dimensional thermal conductivity field. Furthermore, we did not investigate the influence of a heterogeneous hydraulic conductivity field, which would further increase the dispersion.


   CONCLUSIONS
TOP
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
REFERENCES
 
We developed a method for constructing the two-dimensional thermal conductivity field for a section of a glaciofluvial outwash deposit. The method involved a combination of field and laboratory measurements to determine the bulk thermal conductivity of the aquifer solids, the volumetric water content, and the porosity of the aquifer, as well as a model selection procedure using the information-theoretic approach. Using the AICC, the Campbell model was selected as the best-approximating model for predicting the apparent thermal conductivity of variably saturated sands and gravels.

Thermal conductivities of aquifer solids were determined using two laboratory methods. Conductivity values measured directly with the divided-bar apparatus and estimated from the mineral composition were correlated, indicating that, where direct measurements are not available, estimating thermal conductivity from the mineral composition is a reasonable alternative. For this glacial outwash deposit, the thermal conductivities of the porous medium solids can be divided into three groups, which included fine to coarse sand having a mean thermal conductivity of 4.22 ± 0.10 W m–1 K–1, gravel and sand having a mean conductivity of 3.94 ± 0.12 W m–1 K–1, and till having a mean conductivity of 3.72 ± 0.59 W m–1 K–1.

By combining measured thermal conductivities and site stratigraphy with the measured porosity, we were able to define a two-dimensional apparent thermal conductivity field (Fig. 5) for the glacial outwash deposit as input to a numerical model for simulating heat transport. In the saturated zone, the mean value and standard deviation of apparent thermal conductivity were 2.42 and 0.13 W m–1 K–1, respectively. For the moisture and temperature conditions present, the apparent thermal conductivities in the unsaturated zone were between 40 and 50% lower than the apparent thermal conductivities in the saturated zone. Porosity strongly influenced the predicted two-dimensional conductivity field, indicating that this parameter must be defined carefully.

The numerical simulations showed that, for short transport distances, using a mean thermal conductivity in place of a fully heterogeneous field would yield temperature differences of <1°C relative to the fully heterogeneous field. For the homogeneous cases, predicted temperatures were higher in the plume core and lower along the plume fringes, indicative of reduced thermal dispersion; however, predicted temperature differences may increase with transport distance, plume scale, and in fully three-dimensional systems with heterogeneous aquifer thermal and hydraulic properties. Where small temperature differences are important, such as for temperature-sensitive aquatic environments, consideration of the heterogeneities in thermal conductivity may be necessary. These issues will be explored in future work.


   ACKNOWLEDGMENTS
 
We wish to thank Fred Grubb, USGS, Sacramento, CA, for measuring the thermal conductivity of the aquifer solids with the divided-bar apparatus; Ron Griffiths for his guidance on the appropriate use of the ANOVA analysis and multiple comparison tests; David Aldridge for providing a copy of Pronto; and Hansreudi Maurer for providing an independent inversion of the GPR velocity data. We also wish to thank Yasushi Mori for providing the thermal conductivity value for the dry Tottori dune sand. This work was funded by the Department of Fisheries and Oceans Canada, the Ontario Aggregate Resources Corporation, the Natural Science and Engineering Research Council of Canada through grants awarded to R.A. Schincariol, and student grants from the Geological Society of America, the American Association of Petroleum Geologists, and the Province of Ontario.

Received for publication September 7, 2005.


   REFERENCES
TOP
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES
 





This Article
Right arrow Abstract Freely available
Right arrow Figures Only
Right arrow Full Text (PDF) Free
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Similar articles in this journal
Right arrow Similar articles in ISI Web of Science
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Right arrow reprints & permissions
Citing Articles
Right arrow Citing Articles via ISI Web of Science (1)
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Markle, J. M.
Right arrow Articles by Molson, J. W.
Right arrow Search for Related Content
PubMed
Right arrow Articles by Markle, J. M.
Right arrow Articles by Molson, J. W.
Agricola
Right arrow Articles by Markle, J. M.
Right arrow Articles by Molson, J. W.
Related Collections
Right arrow Soil Thermal Properties
Right arrow Ground Penetrating Radar, GPR
Right arrow Heat Transport

HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
The SCI Journals Agronomy Journal Crop Science
Journal of Natural Resources
and Life Sciences Education		 Vadose Zone Journal
Journal of Plant Registrations Journal of
Environmental Quality		 The Plant Genome