THE THERMAL EFFECTS OF THREE DIMENSIONAL GROUNDWATER FLOW By ALLAN DAVID WOODBURY B.Sc. The University of B r i t i s h Columbia, 1973 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE n THE FACULTY OF GRADUATE STUDIES (The Department of Geological Sciences) We accept this Thesis as conforming to the required standard The University of B r i t i s h Columbia September 1983 Â© A l l a n David Woodbury, 1983 V In presenting requirements of British it freely agree for this thesis i n partial f o r an advanced degree Columbia, I agree that available that f o r reference permission scholarly and study. I be g r a n t e d that copying or publication shall n o t be of /ffe/o<;/C<zÂ£ The U n i v e r s i t y o f B r i t i s h 1956 Main M a l l Vancouver, Canada V6T 1Y3 (3/81) of this SC/tfA/CfS Columbia thesis of my I t i s thesis a l l o w e d w i t h o u t my permission. Department further by t h e head understood gain make copying of this o r by h i s o r h e r r e p r e s e n t a t i v e s . financial University shall department for DE-6 at the the Library f o r extensive p u r p o s e s may fulfilment of the written TO My Family iii ABSTRACT The proper measurement and understanding of heat flow values i s essential in studies of tectonics (regional heat flow), o i l maturation, low and high temperature and earthquake geothermal environments, research. Numerical solutions of the equations of f l u i d flow and heat transport in porous media are used to quantify the effects of three-dimensional regional groundwater flow on the thermal regime. The Galerkin F i n i t e Element method is used to solve the coupled equations governing f l u i d and heat transfer. Simplex tetrahedral elements are used to subdivide the region. The resulting f i n i t e element equations are non-linear with temperature and a "leap frog" i t e r a t i o n technique i s employed to solve the coupled equations. A series of computer simulations are c a r r i e d out to investigate how t y p i c a l three dimensional flow systems can influence heat flow measurements. Generic modeling i s carried out on hypothetical basins with a 40 km separation between the regional topographic highs and lows. Emphasis i s placed on understanding the conditions under which groundwater flows severely perturb the thermal f i e l d . The results of the s e n s i t i v i t y study indicate that the t r a n s i t i o n from a conduction-dominated to an advective-dominated thermal regime i s sharp, and depends on: the length, depth, and width of the basin, water table configuration, permeability, existance of extensive or discontinuous aquifers, and hydraulic anisotropy. In addition, the advective threshold i s a function of location within the basin. iv Simulations show that surface heat flow r e f l e c t s the s p a t i a l d i s t r i b u t i o n of recharging and discharging f l u i d associated with the three dimensional groundwater flow system. For the water table configurations studied, there i s a high s e n s i t i v i t y of surface heat flow to large scale variations of water table topography. A s i g n i f i c a n t difference i s noted between two and three dimensional simulations of the same water table configuration. The results underscore the need to recognize the effects of large scale r e l i e f in the d i r e c t i o n orthogonal to any section l i n e . The l a t e r a l variation in surface heat flow i s considered to be a poor indicator of active groundwater flow systems in an advectively-disturbed three dimensional basin. Suggestions are made for a two step strategy to be employed for locating boreholes for heat flow measurements. A number of c l o s e l y boreholes in several directions may spaced be needed to c o r r e c t l y interpret heat flow measurements. Near-correct basal heat flux values can be measured in areas of the basin where there i s a t r a n s i t i o n from regional recharge to discharge. At these locations heat transfer i s found to be l a t e r a l l y advective and v e r t i c a l l y conductive. V TABLE OF CONTENTS Page LIST OF TABLES ' LIST OF ILLUSTRATIONS ACKNOWLEDGEMENTS 1. INTRODUCTION 1.1 1.2 1.3 1 .4 Uses of Heat Flow Measurements Heat flow Determinations Previous Work Present Work 2. GOVERNING EQUATIONS 2.1 2.2 2.3 2.4 2.5 2.6 2.7 Conservation Equations Mass Energy Momentum Potential Constituative Relationships Summary of Equations and Assumptions 3. NUMERICAL SOLUTION 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9 Development of the Numerical Model Numerical Solution Techniques Galerkin F i n i t e Elements Basis Functions F i n i t e Element Formulation - F l u i d Flow F i n i t e Element Formulation - Heat Transfer Solution Techniques Programming Considerations Verification 4. RESULTS 4.1 4.2 4.3 4.4 vii viii xi 1 1 2 5 11 13 13 14 17 21 22 24 25 29 29 31 32 35 38 42 43 45 46 53 Introduction 53 Boundary Conditions 55 Presentation of Model Results 56 Simulations 59 Series 1. - Simple Water Table 59 Series 2. - Small Scale Variations in Upper Surface. 71 Series 3. - E f f e c t s of Permeable Units at Depth .... 81 Series 4. - Additional Topographic Examples 95 Series 5. - Anisotropic Case 109 Series 6. - Effect of Depth of Basin 111 4.5 Discussion 113 vi TABLE OF CONTENTS...Continued Page 5. SUMMARY AND CONCLUSIONS 123 REFERENCES 131 ) vii LIST OF TABLES Table Page 1. Review of available studies on thermal e f f e c t s of groundwater flow (after Smith and Chapman, 1983) 7 2. Summary of equations 27 3. Outline of s e n s i t i v i t y analysis 54 4. Parameter values for f l u i d and thermal properties held fixed for a l l simulations 5. Groundwater recharge rates to hydrologic system with root square deviations from basal heat flux 57 114 viii LIST OF ILLUSTRATIONS Figure 1. Concepts of heat and Sedimentary Basin Page f l u i d transfer in a 6 2. Density and v i s c o s i t y of water as a function of temperature 26 3. (a) Basic tetrahedral element, (adapted from Segerlind, 1976) (b) Assemblage of elements (adapted from Zienkiewicz, 1977) 33 4. Flow chart showing "leap-frog" algorithm 44 5. Typical f i n i t e element mesh used in simulations. Only element prism faces are shown 47 6. Case 4 v e r i f i c a t i o n s . 7. Case 5 v e r i f i c a t i o n s 8. (a) Conceptual i l l u s t r a t i o n - simple water table configuration with 40 m r e l i e f in the Y d i r e c t i o n , (b) Conceptual i l l u s t r a t i o n - simple water table configuration with 80 m r e l i e f in the Y d i r e c t i o n 9. Simple water table - 40 m r e l i e f in the Y direction 50 51 60 63 10. Surface heat flow - simple water table 40 m r e l i e f in the Y d i r e c t i o n 64 11. Section plots through basin for the Y equals 40m r e l i e f case, k = 5 x 10~ m 66 12. Temperature/gradient - depth plots for coordinate position (a) 67 13. Temperature/gradient - depth plots for coordinate position (b) 68 14. Root mean square deviation versus permeability for the simple water table configuration. 40 m r e l i e f in the Y d i r e c t i o n 70 15. Simple water table - 80 m r e l i e f direction 72 16 2 in the Y 16. Surface heat flow - simple water table 80 m relief 73 ix LIST OF ILLUSTRATIONS...Continued Figure Page 17. Root mean square deviation versus permeablity for the simple water table configuration - 80 m r e l i e f in the Y direction 74 18. Conceptual topography 75 illustration - " f l a t - top" 19. Surface heat flow - " f l a t - top" topography 77 20. Comparison of two and three dimensional model results on a section l i n e at Y equals 10 km - " f l a t - top " topography 78 21. Conceptual i l l u s t r a t i o n water table. 80 m r e l i e f 79 - Local rise on the 22. Surface heat flow - Local r i s e on water table 80 23. (a) Extensive aquifer - 40 m r e l i e f on the simple water table, (b) P a r t i a l aquifer case. 40 m r e l i e f on the simple water table. Aquifer terminates along ridge-line of basin, (c) Stratigraphic pinchout. 40 m r e l i e f on the simple water table. Aquifer terminates along a section at X equals 20 km 83 24. Surface heat flow for various aquifer cases 86 25. Section l i n e s through heat flow for the extensive aquifer case - 40 m r e l i e f in Y 88 26. Surface heat flow along three section l i n e s . 40 m r e l i e f , simple water table 89 27. Comparison between surface heat flow for a homogeneous case and a p a r t i a l aquifer case 91 28. Surface heat flow - Stratigraphic pinchout. 80 m r e l i e f in the Y d i r e c t i o n 93 29. Surface heat flow along three section l i n e s for the pinchout number 2 case 94 30. Conceptual i l l u s t r a t i o n - Concave water table configuration. 40 m r e l i e f in the Y d i r e c t i o n 97 X LIST OF ILLUSTRATIONS...Continued Figure Page 31. Concave water table topography 98 32. Surface heat flow - Concave water table case 99 33. Root mean square deviation versus permeabily for the concave water table 100 34. Conceptual i l l u s t r a t i o n s of the Segmented water table configuration, (a) 80 m r e l i e f (b) 40 m r e l i e f (c) 80 m r e l i e f upland - 8 m r e l i e f lowland 102 35. Segmented water table topography 105 36. Surface heat flow for the segmented water table 106 37. Surface heat flow along three sections in the 40 m r e l i e f segmented water table, k = 5 x 10" m 108 38. Surface heat flow for an anisotropic permeability case. Simple water table configuration. 110 39. Surface heat flow for a reduced depth of basin (new depth: 3.2 km ) 112 40. Plot of basin recharge versus root mean square deviation for three water table configurations 116 41. (a) Surface heat flow along three sections in the X d i r e c t i o n for the simple water table, (b) Surface heat flow along two sections in the X d i r e c t i o n for the concave water table 120 16 2 xi ACKNOWLEGDEMENTS It i s a great pleasure for me to express my appreciation to both Allan Freeze and L e s l i e Smith for their encouragement, invaluable advise and support. I am p a r t i c u l a r l y indebted to L e s l i e Smith for agreeing to be my thesis supervisor, for suggesting the research topic and for the many interactions that we have had. Many hours of discussion were>also shared with Keith Loague, Craig Forster, Grant Garven, and Charles Mase, who I would p a r t i c u l a r l y l i k e to thank for sharing many of h i s ideas with me. Tom Nicol was very helpful in reorganizing one of my computer programs to u t i l i z e some of the U.B.C. matrix solution routines. Gord Hodge produced the i l l u s t r a t i o n s in this thesis and prepared slides for various presentations. I would also l i k e to thank B.C. Hydro for providing summer employment during 1981 and 1982. This research was supported by the Natural Sciences and Engineering Research Council of Canada. 1 CHAPTER 1 INTRODUCTION 1.1 Uses of Heat Flow Measurements The proper measurement and understanding of t e r r e s t r i a l heat flow values i s essential in many Geo-Science studies. These measurements are used in such studies as tectonics (regional heat flow), investigation of low and high temperature geothermal environments, o i l maturation, and earthquake research. From previous work (Smith and Chapman,1983) there i s ample evidence to suggest that the thermal regime of a basin, and thus the surface heat flow, can be highly influenced by c i r c u l a t i n g groundwater. This thesis w i l l examine quantitatively the e f f e c t s of three dimensional groundwater flow systems on surface heat flow and heat flow measurements. Tectonic studies include global tectonics, crust movements, regional u p l i f t , and cooling-contracting studies of the earth's crust. Differences in heat production of d i f f e r e n t rock types can be s i g n i f i c a n t , thus d e t a i l s of crustal thicknesses can be inferred. An example of the applications of heat flow measurements i s found in Hyndman(1976), although there are many such examples to be found in the l i t e r a t u r e . 2 In many areas (Mase et a l . , 1978,1982) heat flow measurements are used to obtain information on geophysical anomalies, c i r c u l a t i n g hot water and magmatic a c t i v i t y . Many areas of the United States and Canada are being explored for potential geothermal exploitation (Lewis et a l . , 1964, Lachenbruch and Sass, 1977, Muffer, 1979). The transformation process of petroleum-like organic matter into petroleum requires energy, and heat combined with pressure has been advanced as an o i l maturation mechanism. A thermal "window" of 65 - 150 Â°C (Gretener, 1981) exists in which o i l maturation occurs, although the process i s sensitive to time, heat and pressure ( i . e . , lower temperatures require longer time periods of b u r i a l ) . Thus, exploration for petroleum using heat flow measurements could be possible (Gretener, 1981). Heat flow/groundwater interactions are part of a larger group of heat and mass transfer studies such as reservoir engineering (thermal recovery), nuclear waste repository studies, and hot water injection into aquifers. Because of the importance of the problems that have been outlined, an accurate understanding of the thermal regime of a basin i s e s s e n t i a l . The primary technique for determining the thermal regime i s from surface heat flow data. Unfortunately, in hydrologically perturbed environments this determination may be accurate. 1.2 Heat Flow Determinations Heat flow values are determined from temperature measurements in near surface boreholes. After a suitable time not 3 period in which downhole temperatures have s t a b i l i z e d following a drilling operation, temperature values are determined, usually with thermocoules or thermistors, at various points down a borehole. Temperature values must be corrected for various e f f e c t s . Topographic corrections are made to account for warping of isotherms near surface and corrections to account for the presence of surface bodies of water can also be made (Lewis and Beck, 1977). The geothermal gradient i s computed as the temperature difference between two points divided by the distance between those points. The thermal conductivityT\â€ž must be determined over the same interval at which the geothermal gradient i s computed and both laboratory and " i n - s i t u " techniques are used. In the laboratory the divided bar method i s used (Beck, 1957). Thin discs cut from core samples are inserted into a divided bar of brass or copper and with measured thermal gradients established across the sample discs a value of ^ can w be calculated. In the f i e l d " i n - s i t u " techniques may consist of inserting down-hole probes with heaters. Analytic formulas (see Garland, 1971) are available which relate the temperature of the probe and the quantity of heat produced to 1\m . The heat flow i s then the geothermal gradient multiplied by the thermal conductivity of the rock/water units (Levorsen,1967). The p r i n c i p a l uncertainty in c a l c u l a t i n g heat flow values is the determination of the geothermal gradient, e s p e c i a l l y in the of presence of an active groundwater flow system. Determination the gradient can be a r b i t r a r y , however, four d i f f e r e n t approaches are normally used: (1) In certain cases where groundwater flows are suspected in 4 portions of the hole, gradients from these intervals are rejected. (2) If portions of the temperature/ depth plots are curved only portions of the plots at the bottom of the hole are used. This approach can lead to results showing that heat flow depends on the depth of the hole. Smith and Chapman (1983) showed that using progressively more data down the hole in this manner can lead to innaccurate results. (3) A least squares l i n e i s f i t t e d through a l l the points on the temperature depth plots and one representative gradient value for the entire hole i s determined. Smith and Chapman (1983) warn that t h i s approach can also lead to unsatisfactory r e s u l t s . (4) One could use methods suggested by Mansure and Reiter (1979) that take into account an idealized condition of one dimensional v e r t i c a l flow between two fixed end temperatures. In a c t u a l i t y , groundwater flows are seldom simple and a complete treatment should account for the coupled nature of heat and f l u i d flow. Freeze and Cherry (1979) present a good introduction to the basics of this coupled process and the following remarks are similar to their description. The flow of water i s controlled by hydraulic gradients and hydraulic conductivity with an additional flow induced by thermal effects (buoyancy). Energy or heat i s transferred through a system by conduction and convection. Conductive transport occurs in the rock/water mass and i s controlled by thermal conductivities of the water/geologic formations and thermal gradients, while convective transport occurs under the influence of moving groundwater. Two types of convective transport are i d e n t i f i e d as 5 convection and advection (usually referred to in the l i t e r a t u r e as free and forced convection). The former i s related to the c i r c u l a t i o n of water due to forces caused by buoyancy, and the l a t t e r , the flow of heat and groundwater under some topographically driven system. It i s now recognized (Smith and Chapman,1983) that neither conduction nor advection t o t a l l y dominates in a l l geologic environments. For example, in t e c t o n i c a l l y stable areas of low r e l i e f and in areas of c r y s t a l l i n e basement rocks at surface, temperature f i e l d s would be almost t o t a l l y conductive. Contrasting with these conductive cases would be flow in permeable sedimentary units overlying basement rocks, and in t e c t o n i c a l l y active areas, where a large potion of the thermal f i e l d could be advectively dominated. Thus, measurements taken in certain areas of a basin with an active groundwater flow system could give investigators an apparent high (or low) heat flow anomaly more related to the hydrologic flow system than an indication of tectonic a c t i v i t y at depth. Figure 1 i s an i l l u s t r a t i o n of the concepts presented in this section. 1.3 Previous Work Smith and Chapman (1983) provide a review of previous work in this f i e l d . This section w i l l summerize some of the more important papers. A descriptive and comparitive study i s given in Table 1. Early heat flow investigators (Van Orstrand, 1934, 1939, Bullard, Birch, 1947) a l l found that pure heat conduction could not Surface heat flux Groundwater flow q =K + + + + + + + +â€¢+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + A + + + + + + + + + + + + +A+ + V \\\V V \\T + + + + + + + + + + + + + + + + â€¢+ F i g u r e 1. Concepts of heat and f l u i d + + + + + + + + + Basal heat flux t r a n s f e r i n a Sedimentary Basin. + + A + + -bflW- + X% I + -- ' 1 Table 1. Review of available studies on thermal effect (Adapted from Smith and Chapman, 1983) â€¢a 01 01 N AUTHOR DATE o c QJ u 01 DIMENSION CH UH o ur-. â€¢rl ID 2D 3D 1 Bredehoeft & 13 u c Ol 6 01 01 4-1 rH â€¢H â€¢H CO c c c T H â€¢H CO CH Â«4H P a p a d o p u l o s 1965 X X 01 S 4-1 â€¢H rH 0) T) 01 of groundwater flow. to 01 4J CO 4-t 4J w G 01 -rt T> CO c01o 4J c CO 1-1 CO 4-1 X 3 o 01 3 o 01 c CO CJ â€¢rl cj a â€¢H o a 4HJ o o GO I-I u o o 01 u TCtO B u01 . o CO c o J= J C â€¢rl CO 01 c 00 01 X X 3 4-> CO rH APPLICATION 3 CJ CO rH CQ CO CO CJ >. 4>J, 4J â€¢H â€¢ri o a o o r-C rH Ol 01 > > 4J â€¢H CO 4 J -H o CJ CO CO c â€¢H 01 ^ > a u B co cu oi 4J J2 X X vertical groundwater flow 2 Stallman 1967 X X X X X X X modify 3 Cartwright 1970 X X X X X X X apply r e f . 1 t o e s t i m a t e groundwater Sorey 1971 X X X X X X X apply r e f . 1 to d r i l l h o l e X X X X X X X X modify X X X X X X X temperature X X X 5 Mansure & R e i t e r 1979 6 Domenico & P a l c i a u s k a s 1973 7 Morgan e t a l . 1981 8 Bodvarsson 1972 X X X X X X X model temperature o f water X X X X X X vert. X X X X X X temperature X X X X X X 1973 1971 11 Lewis & Beck 1977 12 K e y s & Brown 1978 13 Kilty 1980 X 14 Ziagos & Blackwell 1981 X 15 Brott. eta l . 1981 X 16 Parsons 1970 X 17 Sorey 1976 18 Betcher 1977 X X 19 Andrews 1978 X X 20 Andrews & A n d e r s o n 1979 X X X 21 Smith 1983 X X X X X X 22 Woodbury X X X X X X & Chapman 1983 X X X X X X X X X X X X X X X X X r e f . 6 t o R i o Grande r i f t X X Cartwright & Chapman i n two-dimensional X Bodvarsson X fields X 9 X X X X X X X X X X temp, i n v e r s i o n s model h e a t injected of horz. flow aquifer surface water f l o w v a r i a t i o n s by assumed waters i n aquifer flow i n a q u i f e r s X X X X X transient X X X model s u b s u r f a c e t e m p e r a t u r e s X fluid f l o w and t e m p e r a t u r e s X fluid flow, d i s c h a r g e , temperature X g e n e r i c modeling o f temperature X X temperature X X X X X X X g e n e r i c modeling o f thermal regime of b a s i n s X X X g e n e r i c modeling o f thermal regime of basins X X X X X X X X X X X X X X X X X threshold aquifer effects l a y e r above ref . 8 & 9 to trace perm, into by t r a n s i e n t i n confining in a drillhole basin basins, injected m o d e l down-dip p e r c o l a t i o n o f c o l d apply profiles r e f . 1 to study heat flow v a r i a t i o n s X 10 discharge temperature X X X apply r e f . 1 f o r low v e l o c i t i e s temp, i n c o n f i n i n g b e d s s u r r o u n d i n g h o r z . i n aquifer thermal a l t e r a t i o n i n Snake R i v e r i n a glacial aquife plain complex i n Long Valley caldera field during o p e r a t i o n o f heat o f groundwater, s u r f a c e water pump seepage 8 explain observed heat flow patterns in the f i e l d . Groundwater movement was c i t e d as the disturbing factor on an otherwise conductive regime and subsequent investigators have confirmed these early observations. A number of groundwater hydrologists have recognized the importance of the overall coupling of heat and f l u i d transfer with respect to regional heat and groundwater flow. Donaldson (1962) was one of the f i r s t investigators to solve the combined h e a t / f l u i d transfer equations in a two dimensional system. He studied the growth and d i s t r i b u t i o n of convective thermal c e l l s with application to thermal regimes in New Zealand. This particular model assumed a constant flow f i e l d and constant f l u i d properties. A widely used formulation describing the thermal aspects of groundwater flow i s the one dimensional solution of Bredehoeft and Papadopulos (1965). A set of type curves are presented so that, with a match between these curves and observed temperature data, the s p e c i f i c discharge of the f l u i d can be computed. The boundary value problem solved was the steady flow of heat and f l u i d between two points with constant temperatures. Modifications to this method have been made by Stallman (1967) and Mansure and Reiter (1979). Cartwright (1971) applied the method of Bredehoeft and Papadopulos in tracing shallow groundwater flow systems (determining recharge and discharge areas ) using s o i l Parsons temperatures. (1970) used a f i n i t e difference technique to solve the heat and mass transfer equations. His model assumed two dimensional steady flow with constant f l u i d properties. Results 9 showed that groundwater flows exerted a c o n t r o l l i n g influence on the shallow temperature f i e l d s in a g l a c i a l complex in northern Ontario, Canada. Another important contribution to this f i e l d i s the presentation by Domenico and Pakiausicas (1973), who solved the h e a t / f l u i d equations a n a l y t i c a l l y for steady flow with constant f l u i d properties, in a manner r e f l e c t i v e of Toth's (1962,1963) work in regional groundwater studies. The importance of this paper i s in i t s contribution to the understanding of how basin surface topography, length and depth of flow f i e l d , and hydraulic and thermal conductivities governed the subsurface temperature d i s t r i b u t i o n . The l i m i t a t i o n s of their approach include constant f l u i d properties and the limited range of v a l i d i t y of their l i n e a r i z e d solutions. There have also been many studies by other authors (Lewis and Beck, 1977, K i l t y and Chapman, 1980, Brott et al.,1981 have attempted to deal with two dimensional problems, ) who but have assumed either a known v e l o c i t y f i e l d or temperature f i e l d in advance. Toth (1962,1963), and Freeze and Witherspoon (1967) showed that even in homogeneous, isotropic sytems, simple variations in the regional topography can result in complicated flow systems. For these reasons the approaches suggested by these other authors have l i m i t a t i o n s . In the more recent l i t e r a t u r e , a general c l a s s i f i c a t i o n of three d i f f e r e n t types of studies i s possible. One class of studies i s the high temperature, multiphase type of geothermal reservoir simulations as reported by Mercer et a l . (1975), Faust and Mercer (1979) and Huyakorn and Pinder (1977). Pinder (1979) 10 summarizes many of the existing models. In another class are the studies concerned with thermal energy storage in aquifers; such as Andrews (1978), Tsang et a l . (1981) and Sauty et a l . (1982a,b). In a broad t h i r d class are studies that deal with low temperature geothermal environments with respect to forced or topographically controlled regional systems. The next four papers, including Parsons (1970) which was previously described under early contributors, belong in t h i s group. Betcher (1977) carried out a s e n s i t i v i t y analysis using a f i n i t e element technique. His study included simulations to confirm dependancies of the thermal regime on hydraulic and thermal conductivities, porosity, and geothermal flux. However, most of the simulations assumed constant f l u i d properties. Sorey (1978) developed a numerical solution to the coupled equations using an integrated f i n i t e difference scheme. This p a r t i c u l a r paper dealt mostly with free convection. Garven and Freeze (1982) and Garven (1982), from an ore genesis perspective solved the coupled equations of f l u i d , heat, and mass (chemicaly reactive) transfer in a heterogenous media. F l u i d properties were coupled to temperatures in these studies. Smith and Chapman (1983) provided the most complete treatment to date on the thermal e f f e c t s of regional groundwater flow. The authors used a f i n i t e element technique of the same type as Frind (1980) and Garven (1982) to solve the coupled equations. Generic simulations were c a r r i e d out on a basin with varying topography t o t a l l i n g 1 km elevation change over 40 km , with a maximum depth of 5 km. Model output consisted of the 11 temperature f i e l d with surface heat flows plotted across the basin. Emphasis was placed on understanding under what conditions an advective threshold was reached (an advective threshold i s defined as when the root mean square deviation of the surface heat flow with the basal heat flow exceeds an a r b i t r a r y l i m i t ) . Results showed that a large range (0 to almost 100%) of the basin can have surface heat flow values s i g n i f i c a n t l y different from the true basal heat flux depending on the groundwater flow system. Their results confirmed the dependence of the advective threshold on permeability, basin depth, anisotropic ratios between the horizontal and v e r t i c a l permeability, and the existence and locations of aquifers. As a result of their work, i t i s clear that s i g n i f i c a n t groundwater disturbance may mask or d i s t o r t the basal heat flux, when estimated from measurements in near surface bore holes. The authors conclude: "The results of our simulations suggest that a knowledge of the complete environment of a s i t e , including the water table configuration and subsurface flow system, combined with more c l o s e l y spaced heat flow measurements may be necessary to unravel the true background heat flux in active flow regions". 1.4 Present Work The research presented in t h i s thesis consists of the development of a f u l l y three dimensional model for the purpose of: (1) Providing a numerical model which can make best use of heat flow data sets, because data are seldom c o l l e c t e d on cross 12 sections and flow systems are seldom two dimensional. (2) Performing a generic study of hypothetical three dimensional basins to determine the s e n s i t i v i t y of various physical parameters such as water table topography, permeability, etc., on the thermal regime of a basin. (3) Show how the thermal disturbances from the above simulations effect surface heat flow measurements. (4) Discuss the implications of the studies to heat flow investigations and suggest areas of future research. A numerical program w i l l be developed and checked against known one dimensional dimensional solutions, then compared with two simulations using Smith and Chapman's (1983) computer program. This v e r i f i c a t i o n stage w i l l be followed by simulations of hypothetical basins examining the three dimensional e f f e c t s of: varying permeability and water table topography, extensive and discontinuous aquifers, reduced depth of basin and anisotropic permeability. Emphasis w i l l be placed on understanding when the advective threshold has been reached and what portions of the basin are within a "minimumly disturbed area" with respect to heat flow determinations. F i e l d for heat flow deteminations in hydrologically perturbed environments are introduced, along with suggestions research. techniques for future 13 CHAPTER 2 GOVERNING EQUATIONS 2.1 Conservation Equations As described in Chapter 1, the flow of heat and groundwater is a coupled process. Therefore, any quantitative mathematical description of combined heat and f l u i d transfer must be based on the fundamental equations describing conservation of mass, energy, and momentum. The continuum approach w i l l be adopted in the development and solution of the d i f f e r e n t i a l equations presented in t h i s study. This approach avoids the d i f f i c u l t i e s of working at the pore scale or at the size of individual fractures and involves the use of averaged medium parameters over a representative elementary volume "REV" (Bear,1972) . The REV averaging approach transforms the fundamental conservation equations from a microscopic pore l e v e l to a macroscopic continuum l e v e l . The continuum approach i s v a l i d for unconsolidated sediments, sedimentary rocks, highly fractured igneous or metamorphic rocks, and in broken landslide debris, talus, etc. In media such as karstic limestone or fractured media where the scale of fracture spacings i s large, the approach may not be v a l i d (Long 14 et a l . , 1982). Many presentations and derivations of these equations can be found in the l i t e r a t u r e . Stallman (i960) developed the basic equations for a saturated porous media. Brownell et a l . (1977) develop the balance laws and constituative relations for convective hydrothermal (multi-phase) geothermal reservoirs. Rock porosity and permeability dependance on stress i s considered in their development. Faust and Mercer (1979) present mathematical models suitable for numerical simulation of geothermal reservoirs. The equations where written in terms of pressure and enthalpy. L i (1980) derives the balance equations for use in simulating transport of single-phase hot water containing dissolved s i l i c a . Bear (1972) and Bear and Corapcioglu (1981) present rigorous developments of f l u i d , momentum, and energy transfer in thermo-elastic media. The governing equations are now presented and discussed i n d i v i d u a l l y . Resonable s i m p l i f i c a t i o n s are made to each equation, pertinent to the problem being solved. A l l of the following equations are v a l i d for a single-phase Newtonian f l u i d in a saturated porous medium. Solute concentrations are assumed to be n e g l i g i b l e and cross-coupling phenomena such as the Soret and Dufour e f f e c t s (Bear, 1972, Freeze and Cherry, 1979) are also assumed to be n e g l i g i b l e . 2.2 Mass The equation of f l u i d flow follows d i r e c t l y from the theory that matter can neither be created nor destroyed. The fluid (macroscopic) mass conservation equation for saturated flow with 15 v e r t i c a l deformations and with temperature dependant properties is (Bear and Corapcioglu, 1981, p725): (2-0 v - < j r < + C*PM4)H Where: {^-^)Â¥t <|â€ž = fy- />[V S Cj^y. = the r e l a t i v e s p e c i f i c f l u i d discharge {L/t} = f l u i d s p e c i f i c discharge {L/t} c7Cp = compressibily of s o l i d at a constant temperature {Lt /M} 2 tf[ = porosity {dimensionless} V5 = average velocity of the s o l i d {L/t} y^p = fluid's coefficient of compressibility at a constant temperature {Lt /M} 2 P = f l u i d pressure {M/Lt } 2 OC r = thermal expansivity of the s o l i d at a constant pressure {T~ } 1 Pr = thermal expansivity of the f l u i d at a constant pressure {T~ } 1 ~]~ V . = temperature {T} - operate, ( J A + j ^ ,12) The f i r s t term of equation (2-1) represents the amount of f l u i d entering or leaving a control volume. The second term represents the time rate of change of f l u i d mass due to pressure changes (contributions due to expansion of water and the porous medium). The t h i r d term represents the time rate of change of 16 f l u i d temperature due to thermal expansivity of water and the s o l i d matrix (this term could also be viewed as a f l u i d source/sink term due to changes in temperature). The last term is an additional f l u i d source due to thermally induced f l u i d density changes. Only steady state pressure and temperature f i e l d s w i l l be considered. This i s an appropriate approximation where yearly fluctuations in the water table are small in comparison to the depth of the flow system (Freeze and Witherspoon, 1966) and where seasonal temperature variations effect only the near surface temperature f i e l d s . Because only regional and not l o c a l thermal f i e l d s are being considered a steady state thermal f i e l d is considered appropriate for t h i s study. This approach has been adopted by Betcher (1977), Garven (1982), and Smith and Chapman (1983). Thus, equation (2-1) reduces to: o Noting that where ( ) = density of f l u i d {M/L } and substituting (2-2A) into (2-2) leads to: 3 17 (2-3) v- (e<f) = 0 This equation describes steady-state, temperature dependant flow in a saturated porous medium with no sources or sinks. 2.3 Energy The conservation equation for thermal energy results from the f i r s t law of thermodynamics: " If a system i s c a r r i e d through a cycle, the t o t a l heat added to the system from i t s surroundings i s proportional to the work done by the system on i t s surroundings" (Welty et a l . , 1976,p.73). Bear and Corapcioglu (1981) develop an expression for thermal conservation in a saturated, deforming porous media: h(T%-p) +(i-*Tr) ]'g=o , 18 with ( ? 0 â€ž = Â«?< *< c ('-"OpÂ» s 4 c where: = heat capacity per unit volume of the porous media {M/Lt T} 2 ^ -i = s p e c i f i c heat capacity of the f l u i d v r {L /Tt } 2 2 = s p e c i f i c heat capacity of the s o l i d {L /Tt } 2 2 = density of the s o l i d {M/L } 3 = coefficient of thermal conductivity of the porous media {ML/t T} 3 = thermal conductivity of the f l u i d {ML/t T} 3 u a r "l^ s = thermal conductivity of the s o l i d {ML/t T} Â£ 3 = d i l a t a t i o n of the porous media {L/L}'= stress in the s o l i d {M/Lt } 3 0 = 2 S ^"^s This equation neglects thermal micro-dispersion, sources or sinks, and heat produced by radioactive decay or chemical reactions. I n e r t i a l terms are assumed to be smaller than time rate of changes ( i . e . Vj* * V"" 7 ^ ). porous medium and f l u i d temperatures are assumed to be in l o c a l thermal equilibrium. The f i r s t term of equation (2-4) represents the time rate of change of energy, and the second term i s the flux of thermal energy by conduction. The t h i r d term represents heat transfer by advection and the last two terms express the source of heat due to viscous dissapation and by compression (or expansion) of the porous medium. (Viscous dissapation represents the t o t a l dissapation of energy through the action of v i s c o s i t y on the 19 porous media.) At steady state (2-4) reduces to: (*-5) V-l\ -VT- ^0,^.77 m - O neglecting viscous dissapation with no sources or sinks (Garven ( 2 - 5Â°) ,1982). Pure heat conduction i s described by: V - â€¢ V T V T = O Various authors (Pickens and Grisak, 1979, Garven, Smith and Chapman, 1983) have redefined "Piâ„¢ 1982, as a combined tensor composed of two components; a purely conductive term and a v e l o c i t y dependant macro-dispersion term "ftâ„¢ = D;; + E;j. The concepts of thermal dispersion have been inferred from mass transport theory (Bear, 1972). Macro-dispersion results from large scale heterogeneities in permeability present in a flow system. The usual approach in dealing with dispersion i s to model i t as a d i f f u s i v e process. For heat transport the dispersion tensor i s (Bear, 1979) o k - 0 F; o o O O cc r = longitudinal dispersivity{L} CC- = transverse dispersivity{L} 20 This equation i s written in i t s simplest form (groundwater is flowing p a r a l l e l to the x a x i s ) . Part of the conceptual problem of modeling dispersion as a d i f f u s i v e process i s that in f i e l d situations the <J\_ 's seem to depend on the length of the flow system (Robertson, 1974). Thus, as van der Kamp (1982) points out, the predictive power of (2-6) is questionable. Its value l i e s in predicting thermal plumes at s i t e s where values of cC L r oC" 7 can be calibrated (for example see Andrews and Anderson, 1979). There are also fundamental differences between contaminant and thermal transport. Chemical species move only within void spaces while heat travels in both the s o l i d grains and with the moving water. Theoretically, because equation (2-5) assumes that the temperatures in the s o l i d matrix and water are the same, dispersion should not occur. However, one cause of observed thermal macro-dispersion probably results from this very assumption of thermal equilibrium not being reached within large heterogeneitites. Dispersion would be more important in transient problems. Other differences e x i s t . Thermal terms E,j and D,j are probably within the same magnitude, whereas E ; J for solute transport may be orders of magnitude higher than d i f f u s i v e transport. Thus, macro-dispersion could also be explained by lack of knowledge of variations in thermal conductivity at a f i e l d s i t e . Van der Kamp (1983) provides a good discussion of macro-dispersion in heat transport. F i e l d observations of this phenomena are described in Sauty et a l . (1982a,b). , Smith and Chapman (1983) modeled two dimensional systems including high d i s p e r s i v i t y values ( up to 100 m in a basin 21 40 km long) and concluded that model results were not sensitive to these values. Hence, after Sorey (1978) macro-scale thermal dispersion w i l l not be considered in the model that w i l l be developed. 2.4 Momentum In groundwater flow problems Darcy's Law s a t i s f y s the momentum requirements. This law i s b a s i c a l l y a linear flow law similar to Fourier's Law and Ohm's Law and can be developed from the Navier-Stokes momentum equations for a "creeping" Newtonian fluid (Bear, 1972, Sorey, 1978, Whitaker, 1966). "Creeping" refers to low Reynolds number flows of an upper l i m i t of between 1 and 10 which i s applicable to most groundwater problems. (1972) expresses the generalized form of Darcy's law as: r where: s p e c i f i c f l u i d discharge vector {L/t} f l u i d dynamic v i s c o s i t y {M/Lt} p f l u i d pressure {M/Lt } 2 acceleration of gravity {L/t } 2 2 elevation of reference point above a datum plane {L} k permeability tensor {L } 2 f l u i d potential function {L} Bear 22 In this form, l o c a l i n e r t i a l forces are neglected with respect to viscous forces. 2.5 Potential To u t i l i z e equation (2-7) the f l u i d potential (j) must be described. Hubbert (1940) derived an energy potential 'o for isothermal flow where: is a compressible f l u i d density{M/L } 3 The usual form of equation (2-8) i s : p 1 3 If p ( i - i Â« ' ) ( d f ) , f i s constant, then: h = * + ^ However, because we are dealing with nonisothermal flows, hydraulic head as expressed in equation (2-9) cannot be used as ^ i n equation (2-7) because a unique potential would not e x i s t . One would normally have to work with the pressure form of equation (2-7) as i s common in many engineering studies and in petroleum applications. However, pressures are not usually used 23 by groundwater hydrologists. Bear (1972) and Frind (1980) proposed the use of a pseudo*-potential h (2-10) where p c = _ L h 0 such that: 2 + i s some reference density at a constant o temperature. Frind (1980) defines h as the 0 "equivalent freshwater head" and states that numerical accuracies improved through the use of h Q are instead of pressure. However, t h i s is not the elevation that water would rise to in an open standpipe. Substituting (2-10) into where { TÂ° r (2-7) ~ r e l a t i v e density {dimensionless} The equation describing f l u i d motion i s now separated into two parts; the f i r s t r e s u l t i n g from head differences (forced convection) and the second an additional v e r t i c a l component caused by buoyancy forces of a f l u i d p a r t i c l e f l u i d of density into (2-3) p o ^ imbedded in a (free convection). By substituting (2-11) the f l u i d conservation equation: 24 we a r r i v e at the f i n a l f l u i d continuity equation. 2.6 Constitutive Relationships F l u i d density and v i s c o s i t y vary substantually with temperature. importance Strauss and Schubert (1977) have shown the of incorporating temperature dependance of density and v i s c o s i t y in the equations describing heat and f l u i d flow. Furthermore, Strauss and Schubert (1977), Sorey (1976), and Garven (1982) present evidence that the pressure dependence of these two f l u i d properties may resonably be neglected for the type of problems being addressed in t h i s study. A l l other f l u i d and material properties such as thermal conductivity of the f l u i d , porosity, etc., are assumed to be neither pressure nor temperature dependant (after Garven, 1982). Therefore only two functional relationships of the form: are required. 25 Two of the more simple relationships available describing this dependence are: (2-13) f, J2 C, / \ (a-'*; y = p. [ I - = 3 . 17 x io-" = x 2.56 = 2.4 = C, x 10" 10- r (_ / ? ' ( c6 T - T Â» ) y'(r-zy] * 1 C" 2 5 /O 2<tr.i7 / (~r + J The v i s c o s i t y and density equations are summarized in Bear and Corapcioglu (1981). These two relationships are plotted on Figure 2 . 2 . 7 Summary of Equations and Assumptions The basic equations describing non-homogeneous, anisotropic heat and steady-state f l u i d transport in a porous media have been presented (Table 2 ) . Before proceeding onto Chapter 3 i t i s useful to summarize b r i e f l y the assumptions in the governing equat ions: (1) The porous medium i s assumed to be a continuum, f u l l y saturated with a single phase Newtonian f l u i d (water) with negligible solute concentrations. F l u i d flow rates are small enough to consider the flow laminar. DENSITY 9Z XIO 3 (kg m"3 ) 27 Table 2. Summary of equations. GOVERNING EQUATIONS Fluid Continuity V* (p q) = 0 f Thermal Energy V-A -VT -p C m f v f q â€¢ VT = 0 Momentum q* = â€” T - ( V h 2 0 + p vZ) Potential Equations of State p = p(T) ^=/u(T) r 28 (2) Cross-coupling phenomena such as f l u i d pressures adding an additional driving force to heat flow or vise versa are assumed to be neglibly small. It i s recognized that density and v i s c o s i t y are non-linear functions of temperature but not, pressure. (3) The temperature f i e l d and hydraulic heads are assumed to be at steady state. Hence, consolidation and/or thermal expansion of the medium are not considered. (4) Micro- and Macroscopic thermal dispersion are not considered. Permeability and thermal conductivity can be heterogeneous and anisotropic. Viscous dissapation i s neglected. (5) Temperatures of the f l u i d and s o l i d media are assumed to be in l o c a l thermal equilibrium. (6) No sources or sinks are present in either the f l u i d or temperature fields. (7) Radioactive decay, chemical reactions, etc., are not considered as seperate thermal sources internal to the porous medium. Heat flux sources w i l l be input as boundary conditions. 29 CHAPTER 3 NUMERICAL SOLUTION 3.1 Development of the Numerical Model The formulation and solution of a Boundary V a l u e / I n i t i a l Value problem in groundwater hydrology should proceed in a four step manner as described by Freeze (1978): (1) I d e n t i f i c a t i o n of the physical problem and formulation of a conceptual model; (2) Determination of the governing equations and formulation of a mathematical model; (3) Solution of the mathematical model, given i n t i a l and boundary conditions. (4) Interpretations of the model results in l i g h t of geologic data and model c a l i b r a t i o n . Steps (1) and (2) have been presented in the l a s t two chapters and this chapter w i l l deal with the development of a suitable technique to solve the governing equations. Five methods are available for solving Boundary Value Problems in groundwater hydrology. These are: (a) Flow nets; (b) Physical models; 30 (c) E l e c t r i c analogs; (d) Analytic solutions; (e) Numerical solutions. Flow nets can be constructed for simple groundwater flow geometries (Freeze and Cherry, 1979) but are not well suited for anything but homogeneous isotropic systems with boundary conditions of the f i r s t kind ( D i r i c h l e t ) . A physical model of a f i e l d situation could be constucted in a laboratory with pertubations and response to the system measured d i r e c t l y . These models provide an excellent check on numerical solutions but are time consuming to construct and are not well suited to the varying of material properties. E l e c t r i c analog models have been constructed to solve the groundwater equations. These models are based on the anology between groundwater and current flows (see Freeze and Cherry, 1979) and have been used successfully to model a number of groundwater/contaminant problems (Bourne,1976). Prickett and Lonnquist (1968) review analog and d i g i t a l (computer) techniques in aquifer simuations. In cases such as the problems that w i l l be solved in t h i s study an analog model cannot be used. Analytic solutions to groundwater problems are usually in a closed form or in a s e r i e s / i n t e r g r a l form that has to be numerically evaluated. Analytic solutions are only available where the region of flow i s simple ( i . e . one dimensional flows or simple aquifer problems) and when the equations can be linearized. Numerical methods offer the only way in which we can solve the coupled non-linear equations presented in Chapter 2. 31 Furthermore, these models have the advantage of allowing the user to vary the flow geometry and material properties in an easy way; however, some l i m i t a t i o n s to these models are evident and these w i l l be noted later in this chapter. 3.2 Numerical Solution Techniques Some of the more widely used numerical techniques in groundwater hydrology include f i n t i t e difference, finite element, integrated f i n i t e d i f f e r e c e , boundary intergral equation method, random walk p a r t i c l e models, and the method of c h a r a c t e r i s t i c s . Of these six techniques integrated finite difference, f i n i t e difference, and f i n i t e element are the most suitable to solve the equations presented in Chapter 2. The method of c h a r a c t e r i s t i c s has been used to solve groundwater contaminant problems (Konikow and Bredeheoft, 1978), however i t s use i s more applicable to hyperbolic equations usually found in open channel or stream flow. Grove (1977) describes this method as "cumbersome, expensive, and lacking in mathematical rigor". Random walk p a r t i c l e models are generally used for transient contaminant transport problems (Garven,1982), although i t i s t h e o r e t i c a l l y possible to solve the heat transport equation by t h i s method. Boundary intergral equation methods are currently one of the newer numerical techniques, and have been applied to transient free surface flows (Liggett and L i u , 1983). However only homogeneous systems can be considered and l i t t l e work has been done in three dimensional geothermal modeling with t h i s technique. Integrated f i n i t e difference methods have been developed 32 and applied to groundwater/geothermal problems in one, two, and three dimensions (Narasimhan and Witherspoon, 1976). This method employs the simple f i n i t e difference approxiation and cannot e a s i l y incorporate general tensorial quantities. F i n i t e difference techniques have been used for many years and can handle a wide varity of problems (Remson et a l . , 1971). Freeze (1971) has applied f i n i t e difference to non-linear unsaturated flow in three dimensions. This method i s suitable for many problems, but i t i s d i f f i c u l t to d i s c r e t i z e an i r r e g u l a r l y shaped region and flow v e l o c i t i e s are not easily calculated. The f i n i t e element technique i s the most popular method for most groundwater problems and i s the method that w i l l be used to solve the governing equations in this study. 3.3 Galerkin F i n i t e Elements The Galerkin-based f i n i t e element method i s one of a number of f i n i t e element techniques that have been applied in groundwater hydrology. The interested reader i s referred to Pinder and Gray (1977) for d e t a i l s on other methods. The Galerkin method has received considerable attention in the l i t e r a t u r e . Examples of the application of the technique to advective heat transport problems include Mercer et. a l . (1975), Andrews and Anderson (1979), L i (1980), Garven (1982) and Smith and Chapman (1983). The F i n i t e Element approach involves subdividing the region of flow into a grid or mesh containing small blocks or c e l l s c a l l e d elements. Element corner points are c a l l e d nodes. Properties of the porous media such as porosity, permeability, 33 F i g u r e 3. B a s i c t e t r a h e c l r a l element. (Adapted froir. S e g e r l i n d , 1976) (a) and (b) Assemblage of elements. (Adapted from Z i e n k i e w i c z , 1977) 34 thermal conductivity etc., are assigned to each element and the solution to the dependant variables (head, temperature) i s obtained at each of the nodal points. The dependant variables are assumed to vary in some functional way quantities may between nodes so that be easily calculated at any point within an element. The simplest variation i s a linear one, although higher order variations such as cubic or quadratic polynomials are quite common (Pinder and Gray, 1977). The simplest shape of an element in three dimensions i s a tetrahedron with straight sides (Figure 3). This element i s c a l l e d a simplex three dimensional element when the nodal variation of the dependant variable i s linear (Segerlind, 1976,p.23). There have been only a few f u l l y three dimensional F i n i t e Element groundwater models developed. Both Verge (1975) and Segol (1976) use isoparametric quadrilateral elements (see Pinder and Gray, 1977). These models allow for a higher order polynomial variation between nodes, as well as allowing for curved elements. The advantage of using higher order elements i s that fewer nodes and hence less computer core storage i s needed; however, their disadvantage i s that increased computer time i s needed (see sect. 3.7). Current computer capacities are larger today than when Verge (1975) or Segol (1976) developed programs; consequently i t i s now their more economic to work with more simplex elements and use additional storage rather than increase computer time by using complex elements. Grove (1977) compared linear and hermite polynomial elements for contaminant problems and concluded that linear elements were computationally more e f f i c i e n t for two dimensional problems. Therefore, the simplex 35 tetfahedral shaped element w i l l be used for solving the boundary value problems presented 3.4 Basis in t h i s study. Functions A basis function i s a mathematical expression which has a unit value at one nodal point and zero value at a l l others. It is through these basis functions that the l i n e a r i t y of the solution between nodal points i s introduced. We can define an approximate solution to dependant variables h 0 and T as: A A/; = N- (x, where: 1 = 1,2,3 n nodes A and h ^ 2) A a and T represent approximate solutions at each nodal point. Since we have adopted a linear variation of and T between nodes: A hi = /\x f = . /\x By + t By + O + C2 t D * D where A,B,C,D are c o e f f i c i e n t s to be determined. In the next few pages the basis functions for hydraulic head h 0 36 w i l l be developed. The basis functions for temperature are exactly the same. Hence, at each node (since four nodes form a tetrahedral element): = + / U ; -- A =A i By, â€¢ %i ( t t a 4 B y , * B y Â« + D C z j + C * z D + D D or in matrix form: A X; h; 6 V, c h. 2. LD J /J The values of the linear c o e f f i c i e n t s A,B,C,D can be found by inversion of the above system, leading to: A = ' A v (Cll.h; +. 3 ^ 4 61, h K t QK L where V = Volume of the tetrahedral element and: ) (w) a; J a L a -((z *y K = +((z *y K = -((z *y- - y * ) - y * ( z = + A L b K ; z K c X K K K " Zj ) +z; * ( y L - zj ) â€¢ z;*(y + K Z - y )) L 3 - K y j + zj * U * < c - K> j - y )) L Z )) - L x * z ) - x; * ( z - z ) + z; * (x = L K L +(z *xj - Z j * X ) L - L k = 3 K +(x *(y K - = -(x; *(yj - y = +(x;*( = x y j - - Yj * ( L X < x ) + (y % L xj) x - x ) + (y * Â« - , % ) ) L x t x L (Xj k W - y. * ( x * z 2 xj) - i*y<>> 4 - x ) + y, * ( X J - K j * y < * L .â€¢ ( - L y > K K " Y.- *(xj - x ) + (y. *xj - x * } L x )) - L K < y ) K * ( z - Zj ) - z; * ( x - - y ) - yj *(x t x) - L L K -(x *(y K x; * ( z - z- ) + z; * (x - ( z * x j - Zj *x ) +x] ; V Z j *(y - z ) + zi*(y, - y ) ) L - yÂ«* 3) - yj * ( z y j K K ; ; + (z,*x - z * x ) - = c 2 L L ((z * Z j - y * *> - y * ( z L bl = bj - y, *z<) - K *y< - x K Y ( / * y j * )) * L - x *z ) + z j * ( y * c - x *z ) + z * ( y * x L K L X < <> = x; *(y *z K â€¢- u " y; * ( \ * z y * K\ z L L K ; L < x,*y ) k d K = x ; * ( z * y j - y * z ) - y; Mx. *z L L 3 L - x * z j ) + z *(y _*xj L ; ( x *yj) L d L = x;*(z *yj - y * K v K = 7^ Z j ) - y *(x *z i [ <d J ^ 4 K - x *zj> + z'My^xj K d K - o t L ] ' 38 To identify the basis functions for the linear element, we substitute equations (3-6) into (3-3) and regroup: A (l-X) h = 'Av [ ( flix + b y ; f Q j f C ; 2 +cfj) h; f + j y t Cj x b 2 i dj) hj t A or h + A/j hj = A/; h ; i A / * /)Â« V V , e.g. N = 1 at node i , 0 at j , k , l With the basis functions defined, we can proceed to the next step; the development of the f l u i d flow equation. 3.5 F i n i t e Element Formulation - F l u i d Flow Equation We define an operator L(h ) based on our f i n a l f l u i d flow 0 equation (2-12) as: We define k as = P/ fÂ°3 k At t h i s point we note that i f the f i n i t e element grid i s set up so that the p r i n c i p a l directions of anisotropy in K are a l l i g n e d with the coordinate directions, off diagonal terms in K are reduced to zero (Bear,1972). Then (2-12) reduces to: 3-<0 L (K) - 1 (K.&) ,1^4) 4fo I * K^r) - o 39 The problem i s to find an approximate solution to this equation, To accomplish this we multiply (3-9) by a weighting function and integrate over the flow domain: ^ The L (K) ^ oi2 = o reader i s referred to Segerlind (1976) for a complete treatment of t h i s method. The Galerkin requirement (Pinder Gray, 1977) for this system i s when W; and is set to the basis functions N; : A / ; dx Substituting equations (3-8), (3-9) dl = 0 into (3-10) and applying Green's f i r s t i d e n t i t y , the r e s u l t i n g set of matrix equations can be written as (Garven,1982): (*-") W;, A; k + B. = Qâ€ž 40 where A m w i s an n x n matrix and B* and Q vectors of size n. The matrix k mn n are i s usually referred to as the global s t i f f n e s s matrix and i s composed of the terms: where V i s a volume i n t e r g r a l . This matrix i s a sparse, symmetric, diagonally is banded matrix that i s dominent and positive d e f i n i t e . The bouyancy term B : .2 Mass f l u i d flux boundaries can be imposed through the ( S i s a suface intergral) vector h 41 where Jj 71 i s the area of the element on which flux i s imposed on. Elements of Q h are set to zero i f no boundaries in the f i n i t e element mesh are of the specified flux type or i f a node i s on an impermeable boundary. S p e c i f i c discharges can be calculated substituting for each element by (3-8) into (2-11). The result i s : These three quantities w i l l be needed for each element in order to solve equation (2-5) . Because the solution based on a linear interpolation for h i s 0 between nodes, v e l o c i t i e s w i l l be constant within the elements but discontinous between adjacent elements. Because the numerical solution converges to the correct solution as the element size decreases the d i s c o n t i n u i t i e s i n v e l o c i t i e s between adjacent elements gradually vanish as the element size i s reduced. Higher order basis functions, such as hermite polynomials, allow for a continuous v e l o c i t y d i s t r i b u t i o n between elements, however substantial more computational time i s required because the i n t e r g r a l s expressed in (2-12,13,14) would have to be evaluated numerically. o 42 3.6 F i n i t e Element Formulation - Heat Transfer An operator L(T) based on equation (2-5) can be defined as: L (l ) = Vâ€¢ â€¢ VT - ^ 'VT^O If we a l i g n the coordinate axes with the p r i n c i p a l directions of anisotropy in then: In an i d e n t i c a l manner to that of the f l u i d flow equation, (3-16) i s substituted into the Galerkin requirement (3-10) and Green's f i r s t identity i s applied. The f i n a l result i s another set of matrix equations: (*-'7) where i s an n x n s t i f f n e s s matrix and F^ i s a vector of applied fluxes of size n. The s t i f f n e s s matrix S be expressed as: w/V can 43 This matrix i s also banded and p o s i t i v e - d e f i n i t e but i s nonsymmetric. As terms q ,q , q , increase in magnitude the matrix S , may no longer be diagonally dominant. This behaviour l7 h can lead to innaccuracies in the solution when normal Gaussian Elimination techniques are applied to solve (3-17) (Gambolati,1980). The conductive flux vector F n where S i s expressed as: denotes a surface intergral of an element with a face on a flux surface. The normal flux i s and 9i i s the area of the element face. Like the f l u i d case vanishes on ; insulated boundaries or i f no flux i s prescribed on a surface. 3.7 Solution Techniques The two systems of equations (3-11) and (3-17) with (3-15) and the equations of state (2-13) and (2-14), define a coupled heat and mass transfer problem. The proceedure used to solve t h i s non-linear problem i s c a l l e d a "leap-frog" i t e r a t i o n technique and i s shown on Figure 4. The following algorithm i s similar to that described in Smith and Chapman (1983). The f i r s t step in the proceedure after the generation of a f i n i t e elment mesh i s to solve this system of equations assuming a conductive heat flow regime. Appropriate densities and v i s c o s i t e s are assigned to each element using equations (2-13) and (2-14). The f l u i d flow equation equivalent freshwater (3-11) i s solved to give values of heads at each nodal point, following which equations (3-15) are employed to c a l c u l a t e s p e c i f i c discharges 44 ( START ) PRINT DATA ECHO PRINT R E A D IN DATA NODAL GENERATE MESHES COORDINATES I COMPUTE TEMPERATURES D E T E R M I N E FLUID PROPERTIES PRINT H E A D S / F O R N O D E S /**"" COMPUTE HYDRAULIC HEADS PRINT E L E M E N T S DARCY VELOCITY CALCULATE r DARCY VELOCITIES \ PRINT T E M P E R A T U R E S FOR NODES COMPUTE TEMPERATURES 1 NO CHECK FOR CONVERGENCE IF N E E D E D YES ( STOP Figure 4. Flow c h a r t PLOT S U R F A C E HEAT FLOW, TOPOGRAPHY PLOT C R O S S S E C T I O N DATA THERMAL PROFILES ) showing " l e a p - f r o g " algorithm. 45 q t ,q 5 ,q ? within each element. Next, the heat transfer equation (3-17) i s solved to determine nodal temperature values. For each element, densities and v i s c o s i t i e s are updated using this new estimate of the temperature f i e l d . An i t e r a t i v e procedure i s then employed; the f l u i d flow and heat transport equations are solved repeatedly with densities and v i s c o s i t i e s updated u n t i l the maximum temperature change between i t e r a t i o n s is less than a specified tolerance. For the three dimensional simulations that w i l l be reported in Chapter 4, t h i s tolerance is 0.1Â° C. 3.8 Programming Considerations The computer code developed was written in Fortran IV and compiled under Fortran H, level 2, on the University of B r i t i s h Columbia's Amdahl 470 V/8. The key components of the computer code l i e in the system of data preparation, and the method of solution for equations (3-11) and (3-17). A'seperate mesh generation code was developed to rapidly construct a large, simple three dimensional flow system in which element prisms consisting of eight nodes are f i r s t formed and then subdivided into f i v e tetrahedral elements per prism (Figure 3b). These elements must be assembled in the correct order to guarantee element compatibiliy (Zienckiewicz, 1971). Many solution techniques are available to solve a system of equations t y p i c a l of the form (3-11) or (3-17). The interested reader i s referred to Remson et. a l . (1971), or Pinder and Gray (1977) for d e t a i l s on many of these methods. Verge (1975) recommends the use of Cholesky elimination for banded, symmetric 46 systems (flow equation) and Gaussian elimination for banded non- symmetric systems (heat t r a n s f e r ) . We have used a modified version of Verge (1975) Cholesky routine in the program (based on Neumann et. a l . ,1974) to solve (3-11) and LINPACK (Dongarra et a l . , 1979) routines to solve (3-17). However, both of these routines require the f u l l or one half bands of the system to be stored. Even with e f f i c i e n t nodal numbering, large band widths result from the choice of elements used, l i m i t i n g the size and mesh density of the problems considered. Because of these problems another version of the same program was employed U.B.C SPARSPAK (Nicol, 1982) developed that routines to reorder and compress matrices. This algorithm allowed for a substantial increase in storage and hence a larger mesh, although execution time increased. Economics proved to be the factor l i m i t i n g the simulations reported in t h i s study. Typical meshes consisted of 19 nodal spacings in the X d i r e c t i o n , 10 in Z and 8 in Y for a t o t a l of 1980 nodes and 7600 elements ( A t y p i c a l section through a f i n i t e element mesh i s shown on Figure 5. Here only prism faces are shown). This entire mesh required 6.8 of storage megabytes (computer l i m i t : 7.0 megabytes) for the LINPAK version of the program. For advectively disturbed cases, of which most of the simulations are, n o n - l i n e a r i t i e s in density and v i s c o s i t y required the program to iterate about 11 times to converge, requiring about 700 seconds of CPU time using the LINPAK routines. 3.9 Verification In this section, a summary of the approach used in testing F i g u r e 5. T y p i c a l f i n i t e element mesh used i n s i m u l a t i o n s . Only element p r i s m f a c e s a r e shown. 48 the numerical code i s presented. The numerical solution was compared to known one dimensional analytic solutions and a two dimensional numerical solution. The following comparisons were made: Cased) - Test of the f l u i d (2-3) and heat conduction equations (2-5a). In t h i s case a one dimensional grid with homogeneous, constant f l u i d properties (constant head and temperature boundaries at either end) was compared with an anayltic solution. Exact matches between the analytic and numerical solutions were noted. Case(2) - Test of the f l u i d (2-3) and heat conduction equations (2-5a). In t h i s case a one dimensional grid with homogeneous medium parameters and constant f l u i d properties with constant f l u i d properties (constant temperatures at one end and heat flux conditions at the other ,with constant heads at either ends) was compared with an analytic solution .Exact matches were noted. Case(3) - Test of the f l u i d (2-3) and heat conduction equations (2-5a). A simulation of a one dimensional heterogeneous (the system was divided into two regions with two d i f f e r e n t permeabilities and thermal conductivities) with constant head and temperature conditions at either end of the mesh was compared with an analytic solution. Again, exact matches were noted. Case(4) - Test of the combined f l u i d (2-3) and heat transfer (2-5) equations. In a series of computer runs, a one dimensional system with constant f l u i d properties and homogeneous media parameters (constant head and temperature 49 boundaries at either solutions. Permeabilities (Figure 6) t o c h e c k agreement was ends) were compared four d i f f e r e n t transfer dimensional other) system were heat transfer 14). Two c o m p a r i s o n s for both modeled. with Both case p r o p e r t i e s and system which and largest ( F i g u r e 7 ) . Good (6) - Test solutions of the f u l l y equations w e r e made b e t w e e n was m o d e l e d dimension. condition at the coupled case difference (2-13), (2- dimensional program. cases. and a homogeneous, f o r the isotropic innaccuracies i n solution. negligible allowed being case difference w a s 1.3% f o r the o f 0.004%. i n the discharge groundwater made anisotropic o f 0.14%. Maximum d i f f e r e n c e of highest A two Comparisons were t h r e s h o l d f o r t h e system i s the area (2-3), p r o p e r t i e s were was 0.05% w i t h a n a v e r a g e roundoff fluid of s t a t e the three and f l u i d i n these isotropic programs were c o n s i d e r e d numerical ends, to analytic o f t h e maximum d i f f e r e n c e s o c c u r r e d the a one (constant head a t both The maximum d i f f e r e n c e anisotopic (2-3) and flux at the advective an a v e r a g e fluid a t one end a n d a h e a t temperature a homogeneous, case; Good of runs a n d Chapman's (1983) c o m p u t e r system with fluid (2-5), and coupled model and Smith vary cases. noted. Case to velocity (2-5). In a series parameters were compared dimensional i n four simulations of the combined with constant temperature agreements - Test equations homogeneous media constant analytic noted. Case(5) heat were v a r i e d t o an end of velocities D i f f e r e n c e s i n t h e two and probably e r r o r s i n t r o d u c e d by i n c l u d i n g due t o a third 50 TEMPERATURE F i g u r e 6. Case 4 v e r i f i c a t i o n s . (Â°C) 51 TEMPERATURE (Â°C) Figure 7. Case 5 v e r i f i c a t i o n s . 52 With the six v e r i f i c a t i o n cases completed the program was considered to be f u l l y operational and i t was decided to proceed with the production runs that w i l l be presented in the next chapter. 53 CHAPTER 4 RESULTS 4.1 Introduction As described in Chapter 1, one of the goals of this work i s to assess how t y p i c a l three dimensional flow systems can influence heat flow measurements. The method used to accomplish this goal w i l l be to carry out a s e n s i t i v i t y study using the mathematical model developed in Chapter 3. This study w i l l be similar to the one c a r r i e d out by Smith and Chapman (1983) and w i l l consist of generic modeling on hypothetical basins. Such models are designed to be representative of real field situations, but are kept simple in order to more easily understand the phenomena being studied. A range of r e a l i s t i c input data are chosen, and the model results are assessed for s e n s i t i v i t y to the variation in a parameter while others are held fixed. The parameters that w i l l be investigated in t h i s study are: variations in permeability and water table configuration, intoducing more permeable units at depth, p a r t i a l and discontinuous aquifer cases, anisotropic formations, and reduced depths of active flow (Table 3). 54 Table 3. Outline of s e n s i t i v i t y a n a l y s i s . Simulation Parameter Figure Variations in basin permeability 9,13,15,16 31,32,35 Variations in water table topography 15,18,21,31,35 Permeable units at depth 23,26,27,28 Anisotropic permeability 27 Reduced depth of basin 38 55 4.2 Boundary conditions The models considered in this study consist of simple sloping topographies in three dimensions. Various topographic surfaces w i l l be modeled, but the boundary conditions w i l l be the same for a l l cases, and are similar to those adopted by Smith and Chapman (1983). The upper boundary of the flow domain is the water table, which i s assumed to be a subdued r e p l i c a of the surface topography. At any point on this surface the value of equivalent freshwater head i s equal to the elevation of that point. An isothermal temperature d i s t r i b u t i o n i s also imposed along the water table. A conductive heat flux i s prescribed on the basal boundary, which i s assumed to be impervious to f l u i d flow. At the basin divides along the sides of the model, a condition of zero l a t e r a l heat and f l u i d flux i s set; a result of the symmetry of the basin (Freeze and Cherry, 1979). The reasons for specifying a basal heat flux rather than a constant temperature at the base of the model are that temperature d i s t r i b u t i o n s are then more r e f l e c t i v e of the upper topography (Smith and Chapman, 1983) and that a fixed temperature on the lower boundary does not adequately represent the thermal regime in an advectively disturbed environment. By specifying a.fixed temperature on the water table we are assuming that some average climatic factor i s the c o n t r o l l i n g influence on the near surface thermal regime. For the simulations examined, a temperature of 10Â°C was set on the upper surface of the model. It i s recognized that t h i s temperature boundary condition i s more complicated and.represents the interplay of heat flow in the unsaturated zone, vegetation 56 types, snow l o a d i n g , g l a c i a l transient effects rebound, r a d i a t i o n , and other (Sekioka and Uhara, 1978), i n c l u d i n g e x i s t e n c e of thermal which these e f f e c t s the s p r i n g s . However, because the depths to infuence the thermal regime are relatively small with respect to the o v e r a l l depth of the b a s i n , the assumption of an isothermal temperature Values of m a t e r i a l and boundary i s v a l i d . f l u i d p r o p e r t i e s used s i m u l a t i o n s are summarized i n Table 4. Basal heat i n the flux values imposed on the model are set to 60 mWnr , a value near the 2 c o n t i n e n t a l average otherwise as r e p o r t e d by Jessop e t . a l . (1976). s t a t e d , p o r o s i t y i s assumed to decrease a piecewise manner from 0.20 of the model. The Unless with depth i n at the s u r f a c e to 0.06 at the base thermal c o n d u c t i v i t y of the s o l i d i s constant in the s i m u l a t i o n s , but the thermal c o n d u c t i v i t y of the s o l i d / w a t e r mass does vary with p o r o s i t y , and hence, 4.3 depth. P r e s e n t a t i o n of Model R e s u l t s R e s u l t s from s i m u l a t i o n s w i l l be presented i n the f o l l o w i n g manner. Plan view contour maps of the topographic be presented along with three dimensional surface w i l l schematic r e p r e s e n t a t i o n s of the s u r f a c e . Output of the computer model c o n s i s t s of contour maps of s u r f a c e heat showing s u r f a c e heat flow, c r o s s - s e c t i o n s flow versus d i s t a n c e , and contoured isotherms on s e c t i o n s . In a d d i t i o n , temperature-depth g r a d i e n t - d e p t h p l o t s w i l l be presented at s e l e c t e d and locations within'the basin. For c o n s i s t e n c y i n r e s u l t s , s u r f a c e heat flow values are computed by m u l t i p l y i n g the n e a r - s u r f a c e thermal c o n d u c t i v i t y by TABLE 4. Parameter Values for f l u i d and thermal properties held fixed for a l l simulations Property Reference Temperature Reference Density of f l u i d Specific heat of f l u i d Thermal conductivity of f l u i d Thermal conductivity of rock Basal heat flux Value 10Â°C 997.97 Kgrn' 4186.0 JKg- C" 0. 586 Wm" C 2.51 Wm-'C" 60.0 mWnr 3 1 1 1 1 2 1 58 the temperature gradient computed over a 100 m interval immediately below the water table. Ideally, these model results should be compared with f i e l d measurements over the same i n t e r v a l . Discrepancies could exist between model results and heat flows determined from deeper boreholes, however model results w i l l represent the "worst case", as heat flow results from deeper holes w i l l tend to be closer to the basal heat flux. Surface heat flow i s contoured by U.B.C.'s SURFACE contouring package. In some cases v i s u a l smoothing was necessary to f i l t e r out d i s t o r t i o n s caused by the gridding approximations. Inflow f l u i d quantities to the basin are calculated, along with root mean square deviations from the basal heat flux along each section, and for the entire upper surface. The root mean square deviation serves as a useful a i d in assessing the effects of various topographic and geologic pertubations on the thermal regime. Before proceeding to the r e s u l t s , i t i s necessary to define some terms that w i l l be used throughout the next few sections: Minimumly Disturbed Area This region i s defined as the surface area of the basin in which surface heat flow ranges between 50 - 70 mWnr , or Â± 16% 2 of the basal heat flux. The selection of these values i s a r b i t r a r y , but r e f l e c t s that portion of the basin which could be considered as f a l l i n g within a reasonable range of accuracy for heat flow values in advectively disturbed environments. The shape of the minimumly disturbed area has some physical significance but i t s primary use i s to highlight differences in surface heat flow between simulations as d i f f e r e n t hydrologic 59 disturbances are imposed on the model. Conductive Case In t h i s case permeability, and hence f l u i d v e l o c i t i e s , are so low that almost a l l the heat transfer i s by conduction. Surface heat flow then varies across the basin only because of topographic e f f e c t s . In these cases 100% of the basin i s within 57-63 mWm~ or Â± 5% of the background flux for the cases that 2 w i l l be considered. Advectively Disturbed Case In t h i s instance the modes of heat transfer are by both conduction and advection to the extent that the root mean square deviation from the basal heat flux i s greater than 10% of the imposed value or, Â± 6 mWm -2 (Smith and Chapman, 1983). Heat r e d i s t r i b u t i o n by moving f l u i d dominates the character of the thermal regime. 4.4 Simulations Series j_ - Simple Water Table The three dimensional topographic effect w i l l be examined by carrying out a s e n s i t i v i t y analysis on a very simple water table configuration (Figure 8). With a l l other parameters held fixed, the homogeneous permeability of the basin w i l l be increased in three l e v e l s : from 1 x 10" 1 6 f i n a l l y , 5 x 10" 16 to 3 x 10" 1 6 m . Heat flow values on different 2 m and 2 section l i n e s w i l l be compared with one another to assess the three dimensional e f f e c t . In a f i n a l suite of simulations (for t h i s topographic case) the t o t a l reief in the Y direction w i l l be increased to 80 m from 40 m for these permeability ranges. The 60 Figure 8. (a) Conceptual i l l u s t r a t i o n - simple water table configuration with 40 m r e l i e f i n the Y d i r e c t i o n . 61 F i g u r e 8. (Continued) (b) Conceptual i l l u s t r a t i o n - s i m p l e water t a b l e c o n f i g u r a t i o n w i t h 80 m r e l i e f i n the Y direction. 62 increase in r e l i e f w i l l demonstrate the s e n s i t i v i t y of the basin heat flows to increases in groundwater velocity in the Y direction. Figures 8a and b show the two d i f f e r e n t r e l i e f s of 40 m and 80 m over 10km in the Y d i r e c t i o n , and 500 m in 40 km in the X direction that are imposed. Elevation contours for the 40 m case are shown on Figure 9, along with the location of two f i c t i c o u s boreholes a and b. Gradient/temperature - depth plots w i l l be examined at each of these locations at a later stage. Surface heat flow for this topography 10" 16 m 2 for a mildly advective (k = 1 x ) case i s shown on Figure 10a. Heat flows determined anywhere within this basin would indicate values close to the basal heat flux. S l i g h t l y warped i s o l i n e s indicate the e f f e c t s of groundwater flows in both the X-Z and Y-Z d i r e c t i o n s . These e f f e c t s are more pronounced on Figures: 10b (k = 3 x 10" and 10c, a severely disturbed case (k = 5 x 10" 16 m 2 16 m 2 ) ). In t h i s case only 55% of the basin i s within the minimumly disturbed area, a reduction of 45% from the k = 1 x 10" 16 m 2 case in just one half an order of magnitude increase in permeablility. Because permeability can vary by thirteen orders of magnitude, the t r a n s i t i o n from a conductive to an advective-dominated system i s sharp, confirming the results of two dimensional models (Smith and Chapman, 1983). Figure 10c also indicates the three dimensional nature of the groundwater flow. The maximum heat flows are concentrated at the (40km, 0km) and (40km, 40km) positions, which are the main discharge areas for the deep, regional flow system. A near-surface heat flow value determined in this location would y i e l d substantially higher values than SIMPLE W A T E R TABLE 40 M e t r e Relief 64 SURFACE HEAT FLOW Simple Water Table 40 Metre Relief k = 1 x 10" m 16 2 X-DISTANCE (km) k - 5 x 10~ m 16 X-DISTANCE (km) 2 Contours in mWm" 2 Figure 10. Surface heat flow - simple water table - 40m r e l i e f i n the Y d i r e c t i o n . 65 the imposed b a s a l heat f l u x , i n d i c a t i v e of groundwater d i s t u r b a n c e and not a l o c a l l y higher b a s a l heat f l u x at depth. Note a l s o the lack of c o r r e l a t i o n between the minimumly d i s t u r b e d area and the contours of the water t a b l e . If one t r a c e d the 4600 m contour l i n e around the s i d e of the h i l l would encounter heat flows v a r y i n g from 55 to 75 mWnr 2 10c). Thus, a d r i l l area i n one s e c t i o n one (Figure hole p o s i t i o n w i t h i n the minimumly d i s u r b e d of the basin cannot be r e l o c a t e d w i t h i n the minimumly d i s t u r b e d area by simply f o l l o w i n g a h i l l s l o p e contour. This is a noticeable effect for t h i s the Y d i r e c t i o n and as w i l l be shown l a t e r , more pronounced when the t o t a l relief slight slope i n the e f f e c t in Y is i s much increased. We w i l l now examine s e l e c t e d s e c t i o n s of s u r f a c e heat (Figure ) is 10c). A section at the edge of the basin shown on F i g u r e s 11a, downward groundwater flow b, and c , (Y = 0 or 20 km and shows the e f f e c t s 12) reflects this temperature/depth p l o t , geothermal g r a d i e n t upward flow i n the hole a convex i n d i c a t i v e of upward groundwater For c o m p a r i t i v e purposes, of i n the recharge area and upward flows and steeper g r a d i e n t s at the d i s c h a r g e a r e a . D r i l l (Figure flow flow. a l s o p l o t t e d as the dashed l i n e is the for the c o n d u c t i v e c a s e . As can be seen, only g r a d i e n t s and thermal c o n d u c t i v i t i e s i n depth would y i e l d heat flow v a l u e s measured over 2.5 km c l o s e to the b a s a l heat flux. Figure 13 i l l u s t r a t e s l o c a t i o n of the s e c t i o n deviations the interdependance through the b a s i n , between the root mean square from the b a s a l heat f l u x and p e r m e a b i l i t y . F i g u r e 13 a l s o i n t r o d u c e s an important d i f f e r e n c e between two and three 66 SIMPLE WATER TABLE 4 0 M e t r e Relief k - 5 x 10" m 1 6 2 Section at Y = 10 km Section at Y = 5 km Section at Y = 0 km 20 25 DISTANCE ( km ) w 30 35 40 Isolines in "Centigrade 200 i DISTANCE (km) Figure 11. Section plots through basin f o r the Y equals 40 m r e l i e f case. 67 O 1 0 1 I 10 20 | 30 | 1 40 50 GEOTHERMAL GRADIENT (Â°Ckm- ) 1 60 1 Coordinate Position ( a) X = 37.5 km Y = 2.5 km k = 5 x 10~ m 16 F i g u r e 12. Temperature/gradient p o s i t i o n (.a) . 2 - depth p l o t s f o r c o o r d i n a t e 68 TEMPERATURE ( Â° C ) 10 o 1 0 30 50 70 90 110 130 1 1 1 1 1 1 20 30 10 40 50 150 60 GEOTHERMAL GRADIENT (Â°Ckrrr ) 1 Coordinate Position (b) X - 20 km Y - 10 km k - 5 x 10' m 16 F i g u r e 13. Temperature/gradient position (.b) . 2 - depth p l o t s f o r c o o r d i n a t e 69 dimensional only The results; namely t h a t a f u n c t i o n of p e r m e a b i l t y advective value the advective threshold but l o c a t i o n w i t h i n the b a s i n . t h r e s h o l d h a s been s h i f t e d t o a lower i n a s e c t i o n n e a r t h e edge o f t h e b a s i n , a s e c t i o n near the r i d g e line. i s not Sharp permeabilty i n comparison t o increases i n the root mean square d e v i a t i o n a r e noted f o r these r i d g e l i n e s , which a r e located for i n areas t h a t a r e g e n e r a l l y a c t i n g as d i s c h a r g e the system. A heat flow d e t e r m i n a t i o n 5.0 o r 15.0 km ( F i g u r e values the 11d) b e t w e e n on a s e c t i o n a t Y = 16 a n d 32km w o u l d of t h e c o r r e c t b a s a l heat f l u x . At t h i s groundwater Although t h i s heat f l u x , recharge flow areas yield s e c t i o n most o f i s i n the Y - Z d i r e c t i o n ( h o r i z o n t a l ) . s e c t i o n has v a l u e s significant c l o s e s t to the correct basal p o r t i o n s of the s e c t i o n c o n t r i b u t e t o (0 - 16km) a n d d i s c h a r g e (32 - 40km). O v e r a l l , a s e c t i o n a t Y = 7.5 o r 15.0 km h a s t h e l o w e s t r o o t mean deviations in deviations the (Figure 14) b e c a u s e o f r e d u c t i o n s b a s a l heat f l u x a t the discharge hydrologic s y s t e m . The l a t e r a l square from and r e c h a r g e zones of t h e v a r i a t i o n i n surface heat flow i s o f t e n u s e d a s an i n d i c a t o r o f an a c t i v e g r o u n d w a t e r flow system. H o w e v e r , t h e r i d g e - l i n e o f t h e s y s t e m a t 10km, w h i c h i s a of symmetry, does n o t r e p r e s e n t even though there islittle a good l o c a t i o n f o r m e a s u r e m e n t s lateral v a r i a t i o n i n heat b e t w e e n 10 a n d 35 km. M o s t o f t h e g r o u n d w a t e r plot the that t h i s e f f e c t . Here, a n e a r - l i n e a r i s i n d i c a t e d . S l i g h t downward f l o w w o u l d t e m p e r a t u r e p l o t o f a deep d r i l l t h i s t r e n d w o u l d be n o t i c e d flow f l o w s a r e downward t h r o u g h o u t most o f t h e s e c t i o n . H o l e b ( F i g u r e illustrates line 13) a t X= 20 km, temperature/depth be n o t i c e d from hole, but i t i s d o u b t f u l from a s h a l l o w bore hole. 70 60 H SIMPLE WATER TABLE 40 Metre Relief 50 H 40 H CM E E 30 cr 20 H 7.5 k m 10 k m 10 H â€” 10% q r17 10 ~i 1 1â€”iâ€”iâ€”r ~i r16 10 PERMEABILITY (m ) 1 1 h 1 1â€”r- .-15 10 2 Figure 14. Root mean square deviation versus permeability for the simple water table configuration. 40 m r e l i e f i n the Y d i r e c t i o n . o 71 Near surface heat flow values of 50 mWm" would be indicated, 2 d i f f e r i n g by 17% from the imposed value. The effects of increasing the slope in the Y direction are examined next; a surface with an 80 m drop in the Y d i r e c t i o n is shown on Figure 15. As in the 40 m case, permeabilities were varied in three levels (Figures: 16a, b, c) k = 1 x '10" 1 6 10" 16 and 5 x 10" 1 6 m . Comparing Figures 16c and 2 ,3 x 10c we see that with just a 40km increase in r e l i e f in the Y d i r e c t i o n , the size of the minimumly disturbed area has been reduced to 31% from 55% of the t o t a l area, and the root mean square for the entire surface has increased from Â±30.6 to Â±35.6 mWm". Root 2 mean square deviations for each section l i n e are shown on Figure 17. Again, the 7.5 or 12.5 km section l i n e has the lowest root mean square values; maximum values on the 0 or 20 km l i n e are Â±66 mWm" compared to Â±43 mWm" 2 2 for the 40 m r e l i e f case (Figure 14). The s h i f t i n g of the advective threshold on each l i n e i s noted as well. Thus, model results demonstrate a high s e n s i t i v i t y of surface heat flow to variations in r e l i e f in the Y direction. Series 2 ^ Small Scale Variations in Water Table In a series of two simulations the effects of introducing small scale variations in the water table, while maintaining same overall slopes, are examined. Another topographic the surface is shown on Figure 18; a " f l a t - t o p " topography. An overall r e l i e f of 40 m in 10km 5km from Y = 0 to 5km 5 to 15km. i s maintained, but t h i s now and 15 to 20 km, occurs over with a f l a t surface from The temperature f i e l d and surface heat flow for an SIMPLE W A T E R TABLE 80 M e t r e Relief r 30 X-DISTANCE (km) Figure 15. Simple water table - 80 m r e l i e f i n the Y d i r e c t i o n . 35 40 Contours in metres 73 SURFACE HEAT FLOW Simple Water Table 80 Metre Relief 1 x 10" m 16 15 20 25 X-DISTANCE (km) -r 30 35 40 3 x 10" m 16 10 15 20 - i â€” 25 30 2 35 2 40 X-DISTANCE (km) ^2 k - 5 x 10"1 6m ,o T 15 20 25 30 35 40 Contours in mWm X-DISTANCE (km) Figure 16. Surface heat flow - simple water table 8 0 m r e l i e f . -2 70 PERMEABILITY Figure 17. Root mean square deviation versus permeability - 80 m r e l i e f i n the Y d i r e c t i o n . (m ) 2 for the simple water table configuration Figure 18. Conceptual i l l u s t r a t i o n - " f l a t - top" topography. 76 advectively-disturbed case was computed with k = 5 x 10" 16 m 2 and the results are presented on Figure 19. Comparing these results to those for the simple water table case (Figure 10c), we see that while the shape of the minimumly disturbed region has changed s l i g h t l y , the area i s 51% of the t o t a l basin with an overall root mean square deviation of Â±32.3 mWm", compared to 2 55% and Â± 30.6 mWm" for the simple water table case. These 2 results show l i t t l e difference for the two water table configurations, indicating a lack of s e n s i t i v i t y of the surface heat flow to small scale variations in the water table. This case also points out some of the important differences between two and three dimensional models and between small and large scale r e l i e f . Figure 20 shows the surface heat flow at Y = 10 km for both two and three dimensional simulations of the " f l a t - t o p " topography. A substantial difference i s noted. One might have considered that with a horizontal surface over 10 km in width and a simple slope along i t s length, a two dimensional representation of the thermal regime would be appropriate. This is not the case. The results underscore the need to recognize the e f f e c t s of large scale r e l i e f in the d i r e c t i o n orthogonal to any section l i n e . Small scale e f f e c t s were examined further by introducing a small l o c a l r i s e or "hummock" onto the simple water table topography (Figure 21). The hummock i s 25 m high with a diameter of about 1 km. In this case an 80 m r e l i e f in the Y - d i r e c t i o n was imposed to further enhance any potential three dimensional e f f e c t s . The surface heat flow for t h i s case i s shown on Figure 22. Some change to the 40 mWm" contour i s evident, but 2 SURFACE HEAT F L O W 20 ~ E "Flat Top" Water Table 40 Metre Relief k = 5 x 10- 1 6 m 2 15 UJ o 10H C0 Q I >- 5H ~1 X-DISTANCE (km) F i g u r e l'J. S u r f a c e heat f l o w - " f l a t - t o p " topography. 30 r 35 40 Contours in mWm" 2 SURFACE HEAT FLOW Tlat-Top" Water Table 40 Metre Relief Section at Y = 1 0 km k - 15 CO 20 25 DISTANCE (km) 30 5 x 1(f m 1 6 35 2 40 F i g u r e 20. Comparison of two and t h r e e d i m e n s i o n a l model r e s u l t s on a s e c t i o n l i n e a t Y e q u a l s 10 km â€” " f l a t - t o p " topography. cc 79 Figure 21. Conceptual i l l u s t r a t i o n - Local r i s e on water table. 80 m r e l i e f . SURFACE H E A T F L O W Local Rise in Water Table 80 Metre Relief k = 5 x 10" m 16 2 r 30 X-DISTANCE (km) F i g u r e 22. S u r f a c e heat f l o w - L o c a l r i s e on water t a b l e . 35 40 Contours In mWm" 2 CO c 81 comparing the results to the simple water table case (Figure 16c), we see that with the hummocky topography 29.1% of the basin i s within the minimumly disturbed area, compared to 31.0% for the simple water table. Root mean square deviations for the entire basin were Â± 35.9 and 35.6 mWnr , respectively. These 2 results indicate l i t t l e change between the two cases. It i s expected that the introduction of the l o c a l r i s e in other locations would not produce any s i g n i f i c a n t changes to the surface heat flow. The o v e r a l l r e l i e f in both the X and Y directions overshadows any p o t e n t i a l e f f e c t s . A series of humps over the entire surface would introduce i n t e r f e r i n g l o c a l groundwater flows that may cause small scale pertubations on the surface heat flow. However, these l o c a l groundwater flows would not penetrate deeply enough to r e d i s t r i b u t e warmer waters, which is the c o n t r o l l i n g influence on pertubations to the heat flow regime. Hence, the thermal regime of a basin and therefore the surface heat flow, i s not sensitive to small scale v a r i a t i o n s in the water table. Series 3 z. E f f e c t s of Permeable Units At Depth So far only basins with homogeneous permeability have been modeled. However, i t i s much more common for sedimentary basins to exhibit strong layering e f f e c t s , and one or more units may behave as a p r e f e r e n t i a l hydraulic pathway. These aquifers are often t i l t e d or offset and in some instances may be 82 discontinuous (see Freeze and Cherry, 1979). Smith and Chapman (1983) examined the quantitative aspects of disturbances to the heat flow regime by aquifers and concluded that aquifers have a s i g n i f i c a n t e f f e c t , but only i f the surrounding permeability i s high enough to allow for adequate recharge to and discharge from the aquifer. They also determined that thermal pertubations are proportional to aquifer permeability, thickness and depth of burial. What i s not well known at t h i s point are how aquifers in a three dimensional setting affect the thermal regime. For instance, i s the three dimensional effect as noted in the f i r s t two series of simulations enhanced or damped out? In this section we w i l l examine effects of extensive and discontinuous horizontal aquifers. Unless otherwise noted, the aquifers modeled are horizontal, 500 m thick units at an elevation of 3.5 km with permeabilities of k= 1 x 10" 14 m. 2 Figure 23a shows the surface heat flow for an extensive aquifer in a basin with a simple water table topography and 40 m relief in the Y d i r e c t i o n . In this simulation the permeability of the unit surrounding the aquifer i s 1 x 10" 1 6 m . For t h i s 2 water table configuration, a homogeneous permeability of 1 x 10" m produces a mild distubance. However, with the 16 2 intoduction of the aquifer (Figure 24a) the root mean square deviation for the entire basin increases to Â± 12.8 mWm", which 2 is almost three times greater than the homogeneous case ( Â± 4.6 mWm" ) and about the same as the s l i g h t l y disturbed homogeneous 2 case ( k = 3 x 10" 1 6 m 2 : Â± 14.5 mWm" ). Heat flow and water 2 table contour patterns are s i m i l a r . Two other effects are noted. Â«3 Figure 23. (a) Extensive aquifer - 40 m r e l i e f on the simple water table. F i g u r e 23. (Continued) (b) P a r t i a l a q u i f e r case. 40 m r e l i e f on simple water t a b l e . the Figure 23. (Continued) (..c) Statigraphic pinchout. 40 m r e l i e f on the simple water table. Aquifer terminates along a section at Y equals 20 km. 86 SURFACE HEAT FLOW Aquifer at Depth Simple Water Table k., - 1 x 1 0 " m 1 6 2 k - 1 x 1<r m 14 a 2 Extensive Aquifer 40 m relief 30 X-DISTANCE Extensive Aquifer X-DISTANCE X-DISTANCE 40 80 m relief (km) Partial Aquifer 20 35 (km) 25 (km) 80 m relief 30 35 40 Contours in mWm Figure 24. Surface heat flow f o r various aquifer cases. -2 87 Section l i n e s taken through the basin (Figure 25) and compared to the k = 3 x 10" m case (Figure 26) show that the three 16 2 dimensional effect has been damped out. Root mean square deviations for individual section l i n e s are similar, Â± 14.0, 12.5, and 11.5 mWnr for the Y =(0, 20 km) and ( 5, 15 km) and 2 10km l i n e s , respectively, compared with Â± 17.9, 13.3, and 11.6 mWnr for the homogeneous case ( k = 3 x 1 0 2 - 1 6 m ). Furthermore, 2 in the homogeneous basin, portions of the heat flow p r o f i l e s on some section l i n e s are close to the basal heat flux value. With the aquifer introduced, about 50% of the heat flow values are above and 50% are below the basal heat flux on each l i n e . The intersection of the heat flow curve with the 60 mW/m i s o l i n e i s 2 shifted s l i g h t l y on each section l i n e . These observations can be explained by examining the flow system. Because of their higher permeability, aquifers have the effect of substantially a l t e r i n g the hydrologic system so that most of the flow i s channeled to the discharge end of the system. The preceeding simulation was repeated with 80 m r e l i e f in the Y d i r e c t i o n so that the e f f e c t s of aquifers could be examined for s e n s i t i v i t y to changes in slope in the Y d i r e c t i o n . The r e s u l t i n g surface heat flow (Figure 24b) has a root mean square deviation of Â± 13.0 mWm' , only s l i g h t l y higher than the 2 40 m case ( Â± 12.8 mWm -2 while enhancing ). This further confirms that aquifers, the thermal disturbance, tend to damp out any three dimensional e f f e c t . In some instances facies changes or faulting may cause the development of a p a r t i a l aquifer (Figure 23b) in which a portion of the basin has an aquifer at depth with the remaining portion EFFECTS OF AQUIFER AT DEPTH 15 20 DISTANCE 25 (km) Figure 25. Section l i n e s through heat flow for the extensive aquifer case - 40 m r e l i e f i n Y. cc CO SIMPLE WATER TABLE LL H â€” i 0 5 1 1 1 1 " ' ' 10 15 20 25 30 35 40 DISTANCE (km ) Figure 26. Surface heat flow along three section l i n e s . 40 m r e l i e f , simple water table. cc 90 homogeneous. This case was investigated by modeling a simple water table case having a drop of 80 m in the Y direction with an aquifer on one half of the basin terminating along the ridgel i n e . The resulting heat flow i s shown on Figure 24c. In comparing this result with the f u l l basin aquifer (Figure 24b) and the homogeneous case ( k=1 x 10" 1 6 m 2 - Figure 16a), i t can be seen that while the half section of the basin with the aquifer i s i d e n t i c a l to the f u l l basin aquifer, the remaining half does not resemble a homogeneous case. In this half (Y = 10 to 20 km) the three dimensional effect has been damped out by the aquifer in the other half of the basin (Y = 0 to 10 km). Figure 27 i l l u s t r a t e s t h i s observation by comparing a section l i n e at 12.5km for the homogeneous case (root mean square deviation Â± 4.27 mWm") with a 12.5km line for the p a r t i a l 2 aquifer case (root mean square deviation Â± 6.81 mWm"). This 2 effect can be explained by examining the flow system in the Y d i r e c t i o n , both with and without the aquifer on one half of the flow system. In many cases groundwater flow divides develop in regional flow systems along axes of symmetry (Freeze and Cherry, 1979). In the simple water table case, for the homogeneous basin, the divide would be a v e r t i c a l plane through Y = 10 km. With the introduction of the p a r t i a l aquifer this effect i s eliminated and the groundwater divide s h i f t s towards higher Y coordinate values and may no longer be v e r t i c a l . A portion of the groundwater flow from the homogeneous half of the basin would then become part of the flow system of the aquifer side. A common occurence that i s of interest to heat flow investigators i s the case of discontinuous aquifers SIMPLE WATER TABLE 80 Metre Relief Section at Y - 12.5 km 1 Homogeneous k = 1x10 m 1 6 2 Partial aquifer k = 1 x 10~ m 16 120n o H < LU I LU O < Li. I 14 a 80 2 A 1 2 40 DC CO 2 k = 1 x 10~ m _J LL CM 2 T 5 10 15 â€”r~ 20 25 DISTANCE (km) 30 35 40 F i g u r e 27. Comparison between s u r f a c e heat f l o w f o r a homogeneous case and a p a r t i a l a q u i f e r c a s e . 92 (statigraphic pinchouts). Figure 23c shows this case conceptually. This problem i s of interest because i t i s known that the pinching out of aquifers tends to d i r e c t flows upward at the terminus of an aquifer, i f located in the upper half of the basin (for example, see Freeze and Witherspoon, 1967). To what extent does this tend to a l t e r the heat flow regime at these locations? To examine t h i s question two simulations were carried out on the simple water table configuration with a discontinuous aquifer introduced at depth (Figure 23c). Background permeabilities were set to 1 x 10" aquifer to 1 x 10" 14 m 2 16 m and the (Case #1) and to 5 x 10" 2 14 m 2 (Case #2). These permeabilities were chosen to be compatible with the extensive aquifer case. The simple water table configuration has a r e l i e f of 80 m in the Y d i r e c t i o n . Results are shown on Figure 28a and b. L i t t l e difference i s noted for these two simulations. Root mean square deviations are Â± 7.5 and Â± 9.9 mW/m. Figure 29 2 shows surface heat flow on three section l i n e s for Case #2: (0, 20), (5, 15) and 10 km. A c h a r a c t e r i s t i c curve i s developed that appears to be a combination of the aquifer and homogeneous cases. It appears from the simulations that the presence of the discontinuous aquifer increases disturbance to the system by causing upward groundwater flows at the aquifer terminus (root mean square deviation i s Â± 5.2 mW/m in comparison for the 2 homogeneous case - k = 1 x 10" 16 m 2 ). The e f f e c t s of increasing the background permeability were also examined. In t h i s case the background permeability was increased to 5 x 10" 16 m and the 2 result i s shown on Figure 28c. At t h i s permeability the thermal e f f e c t s of the stratigraphic pinchout are masked by the o v e r a l l 93 SURFACE HEAT FLOW Simple Water Table Stratigraphic Pinchout 80 Metre Relief k -Aquifer 15 20 1 x 10" m 2 1 x 10" 1 6 2 1 x 10 _ 1 4 14 a m 25 X-DISTANCE (km) m 2 5 x 10" m 2 1 6 15 20 25 X-DISTANCE (km) 30 35 40 Contours in mWm" 2 Figure 28. Surface heat flow - Statigraphic pinchout. 80 m r e l i e f i n the Y d i r e c t i o n . SIMPLE WATER TABLE Stratigraphic Pinchout k - 1 x 10- m 16 5o 120 k =5 x 10 a _J _ 1 4 m 2 2 0 km 5 km 10 km 80 LU = X LU E O < LL 5 w CC CO 40 ~r 5 I 10 i 15 1 20 r~ 25 30 35 40 DISTANCE (km) Figure 29. Surface heat flow along three section l i n e s for the pinchout number 2 case. 4^ 95 disturbance caused by higher groundwater v e l o c i t i e s . A root mean square deviation of Â± 38.5 mW/m was calculated for this case 2 which i s similar to the homogeneous case (k = 5 x 10" 16 m and Â± 2 35.6 mW/m). 2 Based on model results i t can be concluded that there i s only a small "window" of contrasting permeability values in which the pinchout effect would be noticed (for a given aquifer depth and thickness). At background permeabilty than 1 x 10" 16 m, 2 values of less the system would behave as a conductive At higher background values of k = 5 x 10" 16 m, 2 case. the system behaves as an advectively disturbed homogeneous case. An increase in aquifer permeability (for a given background permeability) w i l l increase the disturbance to some maximum l i m i t . Only a further increase in background permeability will cause an increase in disturbance by allowing for more warmer water to be diverted upwards. Series 4 - Additional Topographic Examples In the preceeding series of simulations, we have demonstrated the s e n s i t i v i t y of surface heat flow to changes in relief in the Y d i r e c t i o n and also to small scale pertubations of the water table surface. In t h i s section we w i l l examine the s e n s i t i v i t y of surface heat flow to large scale changes in shape and r e l i e f of the water table in the X d i r e c t i o n . To do t h i s , two a d d i t i o n a l topographic surfaces w i l l be simulated; the 96 f i r s t , a concave surface suggested from Blackwell et. a l . (1980), and the second, a Basin and Range topography similar to one studied by Smith and Chapman (1983). The concave surface consists of a segmented water table (Figure 30 and 31) which has a drop of 200 m in the f i r s t 10km, 150 m in the second, 100 m in the t h i r d , and 50 m in the last 10km. The r e l i e f in the Y d i r e c t i o n i s 40 m. Surface heat flow maps are presented (Figure 32) at three permeability l e v e l s : 1, 3, and 5 x 10" 1 6 m . Some interesting e f f e c t s are noted when one 2 compares t h i s case to the simple water table case (Figure 10c) at k= 5 x 10" 16 m . F i r s t , because of the broad, low sloping 2 discharge region the higher heat flow values decrease in magnitude but are spread out over a larger area. The t o t a l conductive energy remains the same; 48.0 MW for both cases. The maximum values of heat flow for the simple water table case exceeded 140 mWm", whereas the maximum values for the concave 2 water table configuration are about 90 mWm". The root mean 2 square deviation for the basin i s reduced to Â± 19.9 mWm" 2 compared to Â± 30.6 mWm" for the prior case. Second, in contrast 2 to the simple water table case, the section at Y = 7.5 or 12.5 km instead of Y = 5.0 or 15.0 km has the largest portion of i t s length closest to the basal heat flux. These e f f e c t s are a direct result of changes in slope of the water table. With the concave slope, groundwater v e l o c i t i e s at depth are lower at the discharge end of the system and thus, the o v e r a l l disturbance i s reduced. Figure 33, a plot of root mean square deviation versus permeability on each section l i n e , shows how the advective threshold has shifted to a higher permeability value with the 97 Figure 30. Conceptual i l l u s t r a t i o n â€” Concave water table configuration. 40 m r e l i e f i n the Y d i r e c t i o n . C O N C A V E WATER TABLE X-DISTANCE (km) Contours in metres F i g u r e 31. Concave water t a b l e topography. CO 99 SURFACE HEAT FLOW Concave Water Table X - D I S T A N C E (km) k - 3 x 10" m 2 k = 5 x 10" m 2 16 X - D I S T A N C E (km) X - D I S T A N C E (km) 16 Contours in mWm Figure 32. Surface heat flow - Concave water table case. C O N C A V E WATER TABLE 40 M e t r e Relief PERMEABILITY F i g u r e 33. Root mean square d e v i a t i o n v e r s u s p e r m e a b i l i t y f o r t h e concave water t a b l e . 101 change in the topographic configuration of the water table. Therefore, knowledge of only"the o v e r a l l change in elevation of the slope may be inadequate to determine the l e v e l of disturbance in heat flow. The precise location of the water table may be required to properly understand the thermal regime of a basin in an advectively disturbed case. Another flow domain considered i s a water table r e p l i c a t i n g subdued Basin and Range topography, which can be represented by three linear segments of upland, slope and lowland regions. The upland and lowland regions are characterized by slopes in the X d i r e c t i o n of 17.2 m/km. The gradient of the slope connecting two segments i s 51.6 m/km. the The t o t a l r e l i e f i s 1000 m in 40 km in the X d i r e c t i o n . Smith and Chapman (1983) analyzed this configuration in d e t a i l . In three simulations, permeabilities were fixed at 5 x 10" 16 m and r e l i e f in Y was 2 m, and a f i n a l case of 80 m in the upland/slope 8 m in the lowland set at 40 m, 80 regions and only (Figures 34 a, b, c and 35). The f i r s t two simulations were c a r r i e d out to test the s e n s i t i v i t y of surface heat flow to changes in slope for t h i s type of topography. The f i n a l simulation was c a r r i e d out to understand the three dimensional thermal e f f e c t s in a more r e a l i s t i c s e t t i n g . In t y p i c a l Basin and Range topography the slope segments are usually a result of fault u p l i f t i n g , with side r e l i e f developed by erosion and stream action. The lowland regions are often dominated by large a l l u v i a l fans with much gentler slopes to the side. The resulting heat flow (Figure 36) ranges from lows of 30 mWm" to highs of over 140 mWm". These plots are characterized 2 2 by the large gradients of heat flow at the steepest part of the 102 Figure 34. Conceptual i l l u s t r a t i o n s of the segmented water table configuration, (.a) 80 m r e l i e f . 103 F i g u r e 34. (Continued) (.b) 40. m relief. 104 Figure 34. (Continued) (c) 80. m r e l i e f upland - 8 m relief lowland. 105 SEGMENTED WATER TABLE 40 Metre Relief 20 o m a 5H 4400 4200 4600 4700 / 4500 / 4300 15 20 25 X - D I S T A N C E (km) 10 15 20 25 4100 â€”i30 4000 35 40 80 Metre Relief 30 35 40 X - D I S T A N C E (km) 80 Metre Relief Upland/8 Metre Relief Lowland 15 20 T 25 X - D I S T A N C E (km) F i g u r e 35. Segmented water t a b l e topography. 30 35 40 C o n t o u r s in metres SURFACE HEAT FLOW Segmented Water Table k - 5 x 10" m 16 2 40 Metre Relief X-DISTANCE (km) Contours in mWm F i g u r e 3 b . S u r f a c e heat f l o w f o r the segmented water t a b l e . 1 07 topography. Figure 37 i l l u s t r a t e s this observation by presenting surface heat flow on three section lines for the 40 m r e l i e f case. Figure 36a also shows the effects of l o c a l recharge in the lowland region combined with heat flow values greater than the imposed basal heat flux. In Figure 36b, the minimumly disturbed area extends well into the lowland region. Root mean square deviations for these two cases are Â± 32.5 and Â± 36.1 mWm" , 2 indicating a s e n s i t i v i t y to changes in slope in the Y d i r e c t i o n similar to simple water table cases. Figure 36c shows the surface heat flow resulting from the t h i r d topographic case. Now, a heat flow low does not develop in the lowland as the decreased slope in the Y d i r e c t i o n damps out the three dimensional e f f e c t . In t h i s case the minimumly disturbed area consists of a thin band located at the maximum slope segment and follows the contours of the water table. A l o c a l high of 90 mWm" develops at the Y = 0 or 20 km section l i n e in the lowland 2 region ( X = 26 km p o s i t i o n ) . Root mean square values on each section l i n e are: Â± 32.3, 30.0, and 30.2 mWm for the Y equals 2 (0, 20), (5, 15) and 10 km l i n e s . The two additional topographic surfaces modeled in this section show that for a given water table configuration, the shape of the water table has a major influence on the surface heat flow. The Basin and Range simulations also showed that the l e v e l of disturbance in t h i s type of basin i s more sensitive to changes in Y r e l i e f in the lowland section. In discharge regions, higher temperatures exist because of upward groundwater flows. Coupling between temperature and f l u i d properties has the effect of further increasing groundwater v e l o c i t i e s . Hence, any SEGMENTED W A T E R T A B L E 0 5 10 15 20 25 30 35 40 DISTANCE (km) o cc Figure 37. Surface heat flow along three section l i n e s i n the 40 m r e l i e f segmented water table. 109 changes to driving forces (increased topographic slopes) in these locations has a major impact on the surface heat flow. Series 5 ^ Anisotropic Case Layered heterogeneties in permeability can be represented by an equivalent media in which permeability i s homogeneous but anisotropic (Freeze and Cherry, 1979). This layering i s common in sedimentary sequences. Smith and Chapman (1983) investigated the thermal e f f e c t s of anisotropic formations and found t h i s had the effect of s h i f t i n g the advective threshold to a higher permeability value by reducing v e r t i c a l v e l o c i t i e s . In order to investigate the s e n s i t i v i t y of anisotropy in permeability on flow in a basin, the simple water table configuration was simulated with values of k = k = 5 x 10' x y 16 m and k = 5 x 1 0 2 - 1 7 z m . An 80 m r e l i e f in the Y d i r e c t i o n was chosen to accentuate 2 any three dimensional e f f e c t . The result i s shown on Figure 38. The o v e r a l l l e v e l of disturbance i s reduced from the isotropic case( root mean square of Â± 9.0 mWm" compared to Â± 35.6 mWm ) 2 -2 and i s similar to the aquifer case (Â± 13.0 mWm~ ). At t h i s 2 l e v e l of disturbance the three dimensional effect has been damped out as well. The surface heat flow looks similar to the aquifer case with heat flow contours having a similar shape to the water table contours. Root mean square deviations for sections at Y equals (0, 20), (5, 15) and 10 km are Â± 10.4, 8.7 and 8.1 mWnr 2 , respectively. The explanation of t h i s effect i s SURFACE H E A T F L O W Anisotropic Permeability Simple Water Table X-DISTANCE (km) Figure 38. Surface heat f l o w f o r an a n i s o t r o p i c p e r m e a b i l i t y 80 Metre Relief Contours in mWm case. Simple water t a b l e configuration. 111 that the anisotropic case i s a more general form of the aquifer cases simulated. The aquifer'case could have been modeled as an equivalent anisotropic medium with about the same ratios of to k as the present case being considered. Like the aquifer, z flows tend to be channelled more in the X rather than the Y d i r e c t i o n . Therefore, anisotropic cases w i l l cause advective thresholds to be shifted to higher permeabilty values and three dimensional effects w i l l be subdued. Series 6 - Effect of Depth of the Basin The s e n s i t i v i t y of the surface heat flow on the reduction of depth of the hydrologic system was investigated by simulating the simple water table configuration with a reduced flow depth of 3.2 km. This was accomplished by setting several rows at the base of the f i n i t e element mesh to permeabilities of 1 x 10" m. In the rest of the system permeabilities were set to 5 x 2 10" 21 16 m 2 , and the results are shown on Figure 39. The o v e r a l l l e v e l of disturbance was reduced from the 5 km depth case (root mean square deviation of Â± 9.0 mWm" as opposed to Â± 35.6 mWm" 2 ). This s h i f t to higher permeability values i s caused by the elimination of flow paths to the deeper, warmer areas of the basin. Consequently, permeabilities and therefore v e l o c i t i e s would have to be increased to achieve the same l e v e l of disturbance. It i s noted however, that the three dimensional effect has not been damped. 2 SURFACE H E A T F L O W X-DISTANCE (km) Contours in mWm Figure 39. Surface heat flow for a reduced depth of basin (new depth: 3.2 km). 113 4.5 Discussion The purpose of this section w i l l be to discuss some of the implications of the s e n s i t i v i t y study. Suggestions for future research in coupled groundwater/heat flow w i l l be made, along with stratagies for locating boreholes for thermal measurements. It may be possible to estimate (within broad ranges) thermal disturbance in a basin prior to extensive f i e l d exploration. Van der Ramp (1983) and van der Kamp and Betcher (1983) have attemped to relate groundwater recharge rates, dimensionless numbers based on c h a r a c t e r i s t i c lengths and depths of flow systems, and levels of heat flow disturbance. However, the results of Sections 4.4 suggest that the c a l c u l a t i o n of a single number to describe the disturbance may not be possible. A complete description of the water table configuration may required before thermal disturbances may be be estimated. To examine the idea of computing a single number to represent heat flow disturbance more c l o s e l y , we w i l l examine recharge rates and levels of disturbance for d i f f e r e n t water table configurations. The basin recharge as computed by the model, i s defined as the sum of a l l f l u i d fluxes into of the i n d i v i d u a l surface elements multiplied by individual element areas. It therefore represents the t o t a l quantity of flow into the e n t i r e basin that i s required to maintain the elevation of the water table. Table 5 i s a summary of basin recharge rates and l e v e l s of disturbance (root mean square deviations) for a l l the simulations reported here. Comparing the concave water table case (root mean square Â± 19.9 mWm" ; inflow 2.3 x 10~ 2 2 m /sec) 3 114 TABLE 5. Groundwater recharge rates to hydrologic root mean square deviations from basal heat flux. N Simulation Recharge (m /sec) 3 1 S.W.T. 40m 2 S.W.T. 40m 3 S.W.T. 40m 4 S.W.T. 40m 5 S.W.T. 40m 6 S.W.T. 40m 7 S.W.T. 80m 8 S.W.T. 80m 9 S.W.T. 80m 10 S.W.T. 80m 1 1 C.W.T. 40m 12 C.W.T. 40m 13 C.W.T. 40m 1 4C.W.T. 40m 15 B & R 80m 16 B & R 40m 17 B & R 80m 18 B & R 80/8 19 ANIST 80m 20 FLAT 40m 21 REDUC 80m 22 HUMM 80m 23 AQUIF 40m 24 AQUIF 80m 25 PARTA 80m 26 PINC1 80m 27 PINC2 80m 28 PINC3 80m k=1x10- 2 1 k=1x10" 1 7 k = 3x10" 1 7 k=1x10" 1 6 k=3x10" 1 6 k=5x10" 1 6 k=3x10" 1 7 k=1x10- 1 6 k=3x10" 1 6 k=5x10- 1 6 k=3x10" 1 7 k=1x10" 1 6 k=3x10- 1 e k=5x10- 1 6 k=1x10" 2 1 k=5x10" 1 6 k=5x10" 1 6 k=5x10" 1 6 k=5x10" 1 6 k=5x10- 1 6 k=5x10" 1 6 k=5x10- 1 6 k=5x10" 1 6 k=5x10" 1 6 k=1x10" 1 6 k=1x10" 1 6 k=1x10- 1 6 k=5x10" 1 6 4.3X10" RMS deviation (mWm") 8 3.8x10-" 1 .2X10" 3 4.4x10" 3 1.3X10- 2 1.6X10" 3 2.4x10" 2 5.3X10" 3 1.7x10-2 3.1X10- 2 1.3x10- 3 4.3X10" 3 1.3x10-2 2.3x10- 2 9.4X10" 8 4.3X104.7X10" 2 2 4.4X10-2 1.0X10- 2.6xl0' 1.9X10" 2 2 2 3.4x10-2 1.8x10-2 1.9x10-2 1.3x10- 2 1.1x10-2 1.7x10-2 4.1x10-2 Where: S.W.T. = Simple water table C.W.T. = Concave water table B & R = Basin and Range (segmented water table) ANIST = Anisotropic FLAT = "Flat-top" water table REDUC = Reduced depth of flow f i e l d to 3.2 km HUMM = "Hummocky" - l o c a l r i s e in water table AQUIF = Extensive aquifer PARTA = P a r t i a l aquifer PINCn = Discontinuous aquifer case n=1,2 or 3 RMS = Root mean square system with 2 0.8 1 .1 1 .8 4.6 14.5 30.6 2.2 5.2 16.3 35.6 1 .9 4.6 12.1 19.9 1 .3 32.5 36. 1 30.5 9.0 32.3 10.0 35.9 12.8 13.0 9.9 7.5 9.9 38.5 115 with the simple water table case (root mean square Â± 30.6 mWnr 2 ;inflow 2.4xl0~ 2 m /sec), we see almost identical recharge rates 3 to the system, but levels of disturbance that d i f f e r by 35%. The simulations of the simple and concave water table slopes showed that for a given r e l i e f in the X and Y directions the shape of the slope has a major influence on the level of disturbance to the surface heat flow. However, for a given water table configuration, the recharge rates and root mean square deviation values appear to be l o g r i t h m i c a l l y related. Figures 40a and b show this r e l a t i o n s h i p for two different cases. The significance of this observation topographic i s that i t may be possible to determine i f a basin i s hydrologically disturbed (within broad ranges) by estimating recharge and empirically r e l a t i n g this to root mean square deviation (If the water table was known or could be estimated). The only other way to make an estimation of the thermal disturbance would be either to estimate permeabilties from e x i s t i n g hydrologic data or conduct f i e l d permeability t e s t s . Permeabilities could then be used in a mathematical model to estimate the thermal regime. However, i t is d i f f i c u l t to determine permeabilities to an order of magnitude from f i e l d tests (Freeze and Cherry, 1979). Because the advective threshold i s sharp, any inaccuracies in permeability could result in an estimation of thermal disturbance when none existed. Estimations of recharge rates and subsequent calculations of hydraulic conductivity and head within a system have been c a r r i e d out sucessfully (for example, Freeze and Witherspoon, 1967, Jamieson and Freeze, Therefore, 1983). the estimation of recharge rates provides an RECHARGE (x 10" ) m /sec 2 QTT 3 RECHARGE (x 10 ) m /sec -2 o in 3 118 alternate method for determining basin thermal disturbances. An empirical relation between different water table configurations in basins of d i f f e r i n g dimensions, root mean square deviations and recharge rates might be established through the use of type curves. These curves could be developed by solving dimensionless forms of the heat and mass transfer equations (for example, Hoist and Aziz, 1972). For a given dimensionless basin and water table configuration, dimensionless recharge parameters could be calculated, along with root mean square deviations. A heat flow investigator might then be able to estimate heat flow disturbance for a basin (given an estimation of the water table) using these dimensionless type curves. A second suggestion for future research may be to model hydrologic systems with a free-surface approach. That i s , for various recharge rates to the system, the position of the water table in the lower part of the basin would be determined by an i t e r a t i v e approach. This type of modeling would allow for computing the thermal regime in basins where the precise location of the water table was not well defined. F i n a l l y , what i s evident from the results i s that a two step strategy should be employed in any f i e l d program of heat flow measurements. The f i r s t step (as was discussed in the preceeding paragraphs) i s to estimate whether or not a thermal regime may be hydrologically disturbed. If i t i s , the second stage would be to employ f i e l d methodology to deal with the problem. Results from the s e n s i t i v i t y study suggest a general approach for locating boreholes for thermal measurements in 119 hydrologically-disturbed environments. As we have seen, a minimumly disturbed area i s present somewhere within a basin. Continued modelling of d i f f e r e n t topographic shapes may y i e l d information about where minimumly disturbed areas may be located within general water table types ( i . e . concave, convex, simple, e t c . ) . This c l a s s i f i c a t i o n may help investigators to quickly determine where highly disturbed areas are located within a basin. A number of shortly boreholes may spaced then be needed in several directions to c o r r e c t l y interpret heat flow measurements. For instance, a series of closely spaced holes in the X d i r e c t i o n for the simple water table case (Figure 16c) would show l i t t l e l a t e r a l variation in heat flow from 15 to 30 km. A series of closely spaced holes in the Y d i r e c t i o n would indicate a suface heat flow.pattern l i k e Figure 41a, and shows how complicated the heat flow pattern r e a l l y i s because of three dimensional e f f e c t s . These figures show surface heat flow along various sections in the Y d i r e c t i o n for the concave and simple water table cases at k = 5 x 10" m. 2 Interestingly, Figure 41a and b show how 1 6 the correct basal heat flux appears to be located at the main i n f l e x t i o n points on these curves. At these locations the basal heat flux i s unaffected by moving groundwater. These locations also correspond to areas that acts as neither recharge nor discharge to the hydrologic system, and are probably areas of near horizontal groundwater flow. Here, the thermal regime i s v e r t i c a l l y conductive and surface heat flow values measured at these points would give values close to the basal heat flux. A suggested area of research i s to therefore integrate the f i e l d s SURFACE HEAT FLOW Simple Water Table 80 Metre Relief k Â« 5 x 10- m 16 Section at X * 10 km 100 -i X-N Section at X s 5 20 km 2 Section at X - 30 km 80 Basal Heat Flux Heat Flow Approx. Straight Line Segment 10 0 Y-COORDINATE 10 0 10 (km) Figure 41. (a) Surface heat flow along three sections i n the X d i r e c t i o n for the simple water SURFACE HEAT FLOW Concave Water Table 40 Metre Relief k - 5 x 10" m 1 6 100 Section at X - 30 km 2 Section at X - 21.8 km CM E <: o U_ 40 '< UJ I Approx. Straight Line Segment Y-COORDINATE (km) F i g u r e 41. (.Continued) (b) S u r f a c e heat f l o w a l o n g two s e c t i o n s i n t h e X d i r e c t i o n f o r t h e concave water t a b l e . 122 of s u r f a c e w a t e r hydrology, t o i n v e s t i g a t e t h e use approach and heat of d e t e r m i n i n g h o r i z o n t a l i f h y d r o l o g i c d a t a was g r a d i e n t s of heat flow i n the s p a t i a l (X,Y) used t o determine the c o r r e c t b a s a l heat moving lacking, c o o r d i n a t e s might a s e r i e s of b a s i n type c u r v e s shape, d i m e n s i o n l e s s l e n g t h , depth r e c h a r g e , and r o o t mean s q u a r e and be flux. I n summary, f u t u r e r e s e a r c h s h o u l d be c a r r i e d o u t develop This flow h e n c e , p o t e n t i a l a r e a s u n a f f e c t e d by groundwater. Conversely, (1) flow of mapping r e c h a r g e / d i s c h a r g e a r e a s . c o u l d be a way l o c a t i o n s , and groundwater hydrology relating width to d e v i a t i o n from to: topographic dimensionless the b a s a l heat flux; (2) develop field s t a t a g i e s f o r l o c a t i n g b o r e h o l e s and measurements i n h y d r o l o g i c a l l y d i s t u r b e d e n v i r o n m e n t s ; (3) i n t e g r a t e the f i e l d s of s u r f a c e water h y d r o l o g y , and heat and, groundwater f l o w t o d e v e l o p methods f o r mapping r e c h a r g e / d i s c h a r g e a r e a s and measurements. hydrology, thermal relating these t o heat flow regional 123 CHAPTER 5 SUMMARY AND CONCLUSIONS "Any solution i s possible in the absence of fact." - W. A. Meneley, 1974 Heat flow measurements in near surface boreholes have long been used for estimating thermal regimes of the subsurface. It has been recognized that groundwater movement can act as a disturbing factor on otherwise conductive environments. However, what has not been well understood are the quantitative e f f e c t s of groundwater flow on the thermal regime. Groundwater hydrologists have just begun to understand these interactions with the use of two dimensional mathematical models. This study focused on three dimensional effects and included: (1) Providing a numerical model which can make best use of heat flow data sets, as data are seldom c o l l e c t e d on cross sections and flow systems are seldom two dimensional. (2) Performing a generic study of hypothetical three dimensional basins to determine the s e n s i t i v i t y of various physical parameters such as water table topography, permeability, etc., on the thermal regime of a basin. (3) Showing how the thermal disturbances from the above 124 simulations effect surface heat flow measurements. (4) Discussing the implications of the studies to heat flow investigations and suggest areas of future research. The f i r s t goal of this study was to develop a three dimensional mathematical model that would solve the coupled f l u i d and heat transfer equations and allow for temperaturedependant f l u i d properties. The fundamental equations were presented in Chapter 2, and are applicable to steady-state, single-phase f l u i d and heat transfer in a saturated, nondeforming porous media. F l u i d and porous media matrix materials are assumed to be in local thermal equilibrium. The governing equations were solved using a f i n i t e element technique, which was described in Chapter 3. The resulting numerical model was v e r i f i e d against known analytic solutions and an existing numerical model. A second goal of this research was to assess how typical three dimensional flow systems can infuence the thermal regime, and hence, heat flow measurements. The method used to accomplish t h i s task was to carry out a s e n s i t i v i t y study by generic modeling of hypothetical basins. The properties that were investigated are: variations in permeability and water table configuration, introducing more permeable units at depth, p a r t i a l and discontinuous aquifer cases, anisotropic formations, and reduced depths of active flow. The results of the s e n s i t i v i t y study are summerized below: 125 Variations in Basin 1. Permeability. A simple sloping topography in three dimensions was modeled at three levels of permeability to produce conductive to advectively disturbed thermal regimes. The t r a n s i t i o n from a conductive to an advective-dominated system was sharp, confirming the results of two dimensional found to be models. Also, surface heat flow r e f l e c t s the s p a t i a l d i s t r i b u t i o n of recharging and discharging f l u i d associated with the three dimensional 2. groundwater flow system. A near-surface heat flow value in a hydrologically disturbed basin could be substantially higher than the imposed basal heat flux, indicative of groundwater disturbance and not a l o c a l l y higher basal heat flux at depth. There was also a lack of correlation between contours of water table topography and surface heat flow. A d r i l l hole position within the minimumly disurbed area in one section of the basin cannot be relocated within the minimumly disturbed area by simply following a h i l l s l o p e contour. In areas of the basin that were advectively disturbed, only gradients and thermal conductivities measured over 2.5 km in depth would y i e l d heat flow values close to the basal heat flux. 3. Important differences exist between two and dimensional three r e s u l t s ; namely that the advective threshold i s not only a function of permeabilty but location within the basin. The advective threshold has been shifted to a lower permeabilty value in a section near the edge of the basin, in comparison to a section near the r i d g e - l i n e . 4. The l a t e r a l v a r i a t i o n in surface heat flow i s often used as 126 an indicator of an active groundwater flow system. However, the ridge-line of the simple water table system, which i s a line of symmetry, did not represent a good location for measurements even though there is l i t t l e l a t e r a l v a r i a t i o n in heat flow. Variations in Water Table Topography 1. For the simple water table configuration, there i s a high s e n s i t i v i t y of surface heat flow to variations in r e l i e f in the Y d i r e c t i o n . However, the introduction of small scale variations (for a given r e l i e f in X and Y), such as a l o c a l r i s e or a " f l a t - t o p " type of surface, did not produce s i g n i f i c a n t l y d i f f e r e n t surface heat flow than the model without the v a r i a t i o n . Therefore, the surface heat flow i s insensitive to these small scale v a r i a t i o n s . 2. Important differences exist between two and three dimensional models and between small and large scale e f f e c t s . A s i g n i f i c a n t difference was dimensional simulations configuration. The noted between two and three of the " f l a t - t o p " water table results underscored the need to recognize the e f f e c t s of large scale r e l i e f in the d i r e c t i o n orthogonal to any section 3. line. Modeling a concave shaped water table with the same r e l i e f as the simple water table cases showed that knowledge of only the o v e r a l l change in elevation of the slope i s inadequate to determine the level of disturbance geometry of the water table may in heat flow. The be required to properly understand the thermal regime of a basin in an disturbed case. precise advectively 127 4. Another flow domain that was considered was a segmented water table r e p l i c a t i n g subdued Basin and Range topography. Three cases were considered; two to investigate s e n s i t i v i t y of the thermal regime to varying r e l i e f in the X and Y d i r e c t i o n s , and a f i n a l simulation was carried out to understand the three dimensional thermal e f f e c t s in a more r e a l i s t i c setting. It i s possible under certain topographic conditions to have l o c a l recharge in the lowland region combined with heat flow values greater than the imposed basal heat flux. Surface heat flow in segmented basins has a similar s e n s i t i v i t y to changes in slope in the Y d i r e c t i o n to the simple water table configuration. 5. The segmented water table simulations also showed that the l e v e l of disturbance in t h i s type of basin i s more sensitive to changes in Y r e l i e f in the lowland section. In discharge regions, higher temperatures exist because of upward groundwater flows. Coupling between temperature and f l u i d properties has the effect of further increasing groundwater v e l o c i t i e s . Hence, any changes to driving forces (increased topographic slopes) that effect these locations has a major impact on the surface heat flow. E f f e c t s of Introducing Permeable Units at Depth 1. Because of their higher permeability, extensive aquifers have the effect of substantially a l t e r i n g the hydrologic system so that most of the flow i s channeled to the discharge end of the system. Extensive aquifers, while enhancing the thermal disturbance, tend to damp out any three dimensional 2. effect. For a p a r t i a l aquifer terminating along the ridge-line of a 128 basin, the half section of the basin with the aquifer had i d e n t i c a l heat flow patterns to the f u l l basin aquifer, but the remaining half did not resemble a homogeneous case. 3. For a discontinuous aquifer terminating at depth along a section at X equals 20 km, a c h a r a c t e r i s t i c surface heat flow curve i s developed that appears to be a combination of the aquifer and homogeneous cases. It i s concluded from the simulations that the presence of the discontinuous aquifer increases disturbance to the system by causing upward groundwater flows at the aquifer terminus. However, there i s only a small "window" of contrasting permeability values in which the pinchout depth and effect would be noticed for a given aquifer thickness. Anisotropic Permeabilities 1. A basin with a simple water table configuration was modeled with the X and Y d i r e c t i o n permeabilities ten times greater than in the Z d i r e c t i o n . This anisotropic permeability caused the advective threshold to be shifted to a higher permeabilty and three dimensional value e f f e c t s to be subdued. Reduced Depth of Groundwater Flow 1. The s e n s i t i v i t y of the surface heat flow on the reduction of depth of the hydrologic system was investigated by simulating the simple water table configuration with a reduced flow depth of 3.2 km. The advective threshold shifted to a higher permeability value because of the elimination of flow paths to the deeper, warmer areas of the basin. It was noted however, that three The dimensional final implications These 1. goal of Simulations dimensionless disturbance water of the f i n d i n g s are research discussed number not the estimated. prior extensive to configuration, values appear recharge to be observation basin the be empirically r a t e s and logrithmically 2. What should be relating be from employed i n any m e a s u r e m e n t s . The thermal second first r e g i m e may stage problem. general recharge i s evident would approach the boreholes correctly interpret be related. possible The root needed heat flow water table i s that program of field a heat two boreholes i f a recharge step whether study or not a i t i s , the to deal suggest for thermal A strategy flow methodology sensitivity of mean s q u a r e d e v i a t i o n . to be basin estimating basin be the in a significance hydrologically disturbed. If employ to possible to determine to the thermal a given i s to estimate for locating may thermal mean s q u a r e d e v i a t i o n h y d r o l o g i c a l l y - d i s t u r b e d environments. spaced single disturbance root results field step from the measurements. be Results some o f d e s c r i p t i o n of i t may thermal i t may of a required before i s h y d r o l o g i c a l l y d i s t u r b e d by and level e x p l o r a t i o n . For i s that flow of complete However, ranges) field to heat calculation possible. A t a b l e c o n f i g u r a t i o n may are to discuss below: that be subdued. was study to describe (within broad in were not sensitivity estimate the this showed may disturbances this effects with a measurements number o f closely in several d i r e c t i o n s to m e a s u r e m e n t s . The correct basal 130 heat flux i s located at one of the main i n f l e x t i o n points of the surface heat flow map. At these locations the thermal regime i s v e r t i c a l l y conductive, because of horizontal groundwater flow. Surface heat flow values measured at these points would give values close to the basal heat flux. Gradients of heat flow in the s p a t i a l (X,Y) coordinates might be used to determine the correct basal heat flux. Smith and Chapman (1983) concluded that an understanding of the hydrogeology of a s i t e combined with more closely spaced d r i l l holes may be required to correctly interpret near surface heat flow values. It i s apparent that Smith and Chapman's conclusions are even more relevant when three dimensional effects are considered. It i s hoped that the work begun in t h i s study w i l l encourage more interaction between hydrogeologists and geophysists i n thermal investigations. 131 REFERENCES Andrews, C.B., The impact of"the use of heat pumps on ground water temperatures, Ground Water, 16(6), 437-443, 1978. Andrews, C.B., and M.P. Anderson, Thermal a l t e r a t i o n of groundwater caused by seepage from a cooling lake, Water Resour. Res., 15(3), 595-602, 1979. Bear, J . , Dynamics of Fluids in Porous Media , E l s e v i e r , 1972. Bear, J . , Hydraulics of Groundwater McGraw-Hill, New York, 1979. Bear, J . and M.Y. Corapcioglu, A mathematical model for consolidation in a thermoelastic aquifer due to hot water injection or pumping, Water Resour. Res., 17, 723-736, 1981.. Beck, A.E., Steady-state method for the rapid measurement of thermal conductivity of rocks, J . S c i . Instrum., 34, 186-189, 1957. Betcher, R.N., Temperature d i s t r i b u t i o n in deep groundwater flow systems-a f i n i t e element model, M.Sc. thesis, Univ. of Waterloo, 1977. Birch, F., Temperature and heat flow in a well near Colorado Springs, Am. J . S c i . , 245, 733-753, 1947. Blackwell, D.D., J.L. Steele and C A . Brott, The t e r r a i n effect on t e r r e s t r i a l heat flow, J . Geophys. Res., 85(B9), p.4757-4772, 1982. Bourne, D.R., An e l e c t r i c analog simulation of groundwater flow patterns at a potash waste disposal pond located near Esterhazy, Saskatchewan, Unpubl. Msc. Thesis, Univ. B r i t i s h Columbia, 1976. Brott, C.A., D.D. Blackwell, and J.P. Ziagos, Thermal and tectonic implications of heat flow in the Eastern Snake River P l a i n , J . Geophys. Res., 86, 11709-11734, 1981. Brownell, D.A., S.K. Garg and J.W. P r i t c h e t t , Governing equations of geothermal reservoirs, Water Resour. Res., 13, 929935 , 1977. Bullard, E . C , Heat flow in South A f r i c a , Proc. Royal Soc. London, A173(995), 474-502, 1939. Bredehoeft, J.D., and I.S. Papadopulos, Rates of v e r t i c a l groundwater movement estimated from the earth's thermal p r o f i l e , Water Resour. Res., 1(2), 325-328, 1965. Cartwright, K., Redistribution of geothermal heat by a shallow aquifer, Geol. Soc. Am. B u l l . , 82, 3197-3200, 1971. Domenico, P.A., and V.V. Palciauskas, Theoretical analysis of forced convective heat transfer in regional ground-water flow, 1 32 Geol. Soc. Amer. B u l l . , 84, 3803-3814, 1973. Donaldson, J.G., Temperature gradients in the upper layers of the earth's crust due to convective water flows, J . Geophys. Res., 67, 3449-3459, 1962. Donnaga, J . J . , C B . Moler, J.R. Bunch, and G.W. Stewart, LINPAK User's Guide, Soc. Int. Appl. Math., Philadephia, 1979. Faust, C.R. and J.W. Mercer, Geothermal reservoir simulation: 1. Models for liquid-dominated and vapour-dominated hydrothermal systems, Water Resour. Res., 15, p.23-30, 1979. Freeze, R.A., Three-dimensional, transient, saturatedunsaturated flow in a groundwater basin, Water Resour. Res., 7, 347-366, 1971. Freeze, R.A., Mathematical models of h i l l s l o p e hydrology, in H i l l s l o p e Hydrology , (Chapter 6), edited by M.J. Kirkby, 177225, Wiley-Interscience, New York, 1978. Freeze, R.A., and J.A. Cherry, Groundwater, Prenice-Hall, 1979. Freeze, R.A., and P.A. Witherspoon, Theoretical analysis of regional groundwater flow: 1. A n a l y t i c a l and numerical solutions to the mathematical model, Water Resour. Res., 2, 641-656, 1966. Freeze, R.A., and P.A. Witherspoon, Theoretical analysis of regional groundwater flow: 2. Effect of water table configuration and subsurface permeability variation, Water Resour. Res., 3, 623-634, 1967. Frind, E.O., Seawater intrusion in continuous coastal aquiferaquitard systems, Proc. Third Inter. Conf. on F i n i t e Elements in Water Resources, Univ. Miss.,Oxford, May, 1980. Gambolati, G., Perspective on the modified conjugate method for the f i n i t e element solution of linear subsurface equations, Proc. Third Inter. Conf. on F i n i t e Elements in Water Resources, Univ. Miss., Oxford, May, 1980. Garland, G.D., Introduction to Geophysics: Mantle Core and Crust , W.B. Saunders, London, 1971. Garven, G., The role of groundwater flow in the genesis of stratabound ore deposits: A quantitative analysis, Unpulbl. Phd Thesis, University of B r i t i s h Columbia, 1982. Garven, G. and R.A. Freeze, The role of regional groundwater flow in the formation of ore deposits in sedimentary basins: A quantitative analysis, paper presented to the 2nd. Nat. Hydrogeol. Conf., Winnipeg, Manitoba, Canada, 1982. Gretener, P.E., Geothermics: Using temperatures in hydrocarbon exploration, Education Course Notes Series #17, AAPG Annual Meeting, San Fransico, 1981. 133 Grove, D.B., The use of galerkin finite-element methods to solve mass transport equations, USGS Water-Resour. Invest. 77-49, 1977. Hoist, P.H., and K. Aziz, Transient three-dimensional natural convection in confined porous media, Int. J . Heat Mass Transfer, 15, 73-90, 1972. Hubbert, M.K., The theory of groundwater motion, J . Geol., 48, 785-944, 1940. Huyakorn, P.S., and G.F. Pinder, A pressure-enthalpy f i n i t e element model for simulating hydrothermal reservoirs, Advances in Computer Methods for P a r t i a l D i f f e r e n t i a l Equations, IMACS, 284-293, 1977. Hyndman, R.D., Heat flow measurments in the i n l e t s of southwestern B r i t i s h Columbia, J . Geophys. Res., 81, 337-349, 1976. Jamieson, G.R. and R.A. Freeze, Determining hydraulic conductivity d i s t r i b u t i o n s in a mountainous area using mathematical modeling, Groundwater, 21(2), 168-177, 1983. Jessop, A.M., M.A. Hobart, and J.G. Sclater, The world heat flow data c o l l e c t i o n - 1975, Geothermal Service of Canada, Geotherm. Ser. 5, 1976. K i l t y , K., and D.S. Chapman, Convective heat transfer in selected geologic situations, Ground Water, 18(4), 386-394, 1980. Konikow, L.F., and J.D. Bredehoeft, Computer model of twodimensional solute transport and dispersion in ground water, USGS, Tech. of Water Res. Inv., Ch. C2, 1978. Lachenbruch, A.H., and J.H. Sass, Heat flow in the.United States, in The Earth's Crust, Geophysical Monograph 20, edited by J.G. Heacock, 626-675, AGU, Washington, D.C, 1977. Levorsen, A.I., Geology of Petroleum San Fransisco, 1967. , 2nd. edition , Freeman, Lewis, T.J., and A.E. Beck, Analysis of heat flow data detailed observations in many holes in a small area, Tectonophysics, 41, 41-59, 1977. Lewis, T.J., A.S. Judge and J.G. Souther, Possible geothermal resources in the coast plutonic complex of southern B r i t i s h Columbia, Canada, Pageoph., 117, 172-179, 1979. Li,T.M.C, Axisymmetric numerical simulation of hydrothermal systems including changes in porosity and permeability due to the quartz-water reaction, Ph.D. thesis, The Pennsylvania State University, 240p., 1980. 134 Liggett, J.A. and P.L-F. L i u , The Boundary Intergral Equation Method for Porous Media Flow , Allen & Unwin, New York, 1983. Long, J.C.S., J.S. Remer, CR. Wilson and P.A. Witherspoon, Porous media equivalents for networks of discontinuous fractures, Water Resour. Res., 18, 645-658, 1982. Mansure, A.J., and M. Reiter, A v e r t i c a l groundwater movement correction for heat flow, J . Geophy. Res., 84(B7), 3490-3496, 1979. Mase, C.W., D.S. Chapman and S.H. Ward, Geophysical study of the Monroe-Red H i l l geothermal system, Topical Rep. EY-76-S-07-1601, Dept. of Geol. and Geophys., Univ. of Utah, Salt Lake City, Utah, U.S.A., 1978. Mase, C.W., J.H. Sass, A.H. Lachenbruch and R.J. Munroe, Preliminary heat-flow investigations of the C a l i f o r n i a Cascades, USGS open f i l e report 82-150, 1982. Mercer, J.W., Pinder, G.F., and I.G. Donaldson, A Galerkin f i n i t e element analysis of the hydrothermal system at Wairakei New Zealand, J . Geophy. Res., 80(17), 2608-2621, 1975. Muffer, L.J.P., editor, Assessment of geothermal resources of the United States - 1978, USGS c i r c u l a r 790, 1979. Narasimhan, T.N. and P.A. Witherspoon, An integrated f i n i t e difference method for analyzing f l u i d flow in porous media, Water Resour. Res., 12(1), 57-64, 1976. Neumann, S.P., R.A. Feddes and E. Bresler, Development of methods, tools, and solutions for unsaturated flow, Third Annual Report, Project No. A10-SWC-77, Hydraulic Eng. Lab., Technion, Haifa, I s r e a l , 1974. N i c o l , T., UBC SPARSPAK: Solving sparse systems of linear equations, Computiong Centre, Univ. of B r i t i s h Columbia, 1982. Parsons, M.L., Groundwater thermal regime in a g l a c i a l Water Resour. Res., 6(6), 1701-1720, 1970. complex, Pickens, J.F. and G.E. Grisak, Finite-element analysis of l i q u i d flow, heat transport and solute transport in a groundwater flow system: 1. Governing equations and model formulation, Tech. Rep. TR-81, AECL, 1979. Pinder, G.F., State-of-the-art review of geothermal reservoir engineering, Lawrence Berkeley Lab., LBL-9093, 1979. Pinder, G.F., and W. Gray, F i n i t e element simulation in surface and subsurface hydrology, Academic Press, 1977. Prickett, T.A. and C.G. Lonnquist, Comparison between analog and d i g i t a l simulation techniques for aquifer evaluation, Proc. Symp. Use Analog D i g i t a l Computers Hydrol., Inter. Assoc. S c i . 135 Hydrology, Publ. 81, 625-634, 1968. Remson, I., G.M. Hornberger and F.J. Molz, Numerical Methods in Subsurface Hydrology , Wiley-Interscience, New York, 1971. Robertson, J.B., D i g i t a l modeling of radioactive and chemical waste transport in the Snake River Plain Aquifer at the National Reactor Testing Station, Idaho, USGS open f i l e report IDO-22054, 1974. Sauty, J.P., Gringarten, A.C., Menjoz, A., and P.A. Landel, Sensible energy storage in aquifers 1. Theoretical study, Water Resour. Res., 18(2), 245-252, 1982a. Sauty, J.P., Gringarten, A.C., Fabris, H., Thiery, D., Menjoz, A., and P.A. Landel, Sensible energy storage in aquifers 2. F i e l d experiments and comparisons with theoretical results, Water Resour. Res., 18(2), 253-265, 1982b. Segerlind, L.J., Applied F i n i t e Element Analysis , Wiley, New York, 1976. Segol, G., A three-dimensional galerkin f i n i t e element model for the analysis of contaminent transport in variably saturated porous media, User's Guide, Dept. of Earth Sciences, Univ. of Waterloo, 1976. Sekioka, M. and K. Yuhara, Heat flux estimation in geothermal areas based on the heat balance of the ground surface, J . Geophys. Res., 79, 2053-2058, 1978. Smith, L. and D.S. Chapman, On the thermal effects of groundwate flow: 1. Regional scale systems, J . Geophys. Res., 88, 593-608, 1983. Sorey, M.L., A model of the hydrothermal system of Long Valley caldera, C a l i f o r n i a , Summaries Second Workshop Geothermal Reservoir Engineering, Stanford University, Stanford, C a l i f o r n i a , 324-338, 1976. Sorey, M.L., Numerical modeling of l i q u i d geothermal systems, U.S. Geol. Surv. Professional Paper 1044-D, 1978. Stallman, R.W., Notes on the use of temperature data for computing ground water v e l o c i t y , Societe Hydrotechnique de France, Nancy, France, 6th Assembly on Hydraulics, Rapport 3 question 1 , 1-7, 1960. Stallman, R.W., Flow in the zone of aeration, in Advances in Hydroscience, V o l . 4, edited by Ven Te Chow, pp.151-195, Academic, New York, 1967. Straus, J.M., and G. Schubert, Thermal convection of water in a porous medium: effects of temperature and pressure dependent thermodynamic and transport properties, J . Geophy. Res., 82(2), 325-333, 1977. 136 Toth, J . , A theory of groundwater motion in small drainage basins in central Alberta, J, Geophys. Res., 67, 4375-4387, 1962. Toth, J . , A theoretical analysis of groundwater flow in small drainage basins, J . Geophys. Res., 68(16), 4795-4812, 1963. Tsang, C.F., Buscheck, T., and C. Doughty, Aquifer thermal energy storage: A numerical simulation of Auburn University f i e l d experiments, Water Resour. Res., 17(3), 647-658, 1981. van der Kamp, G., Interactions between heat flow and groundwater flow - A review, Waterloo Research Institute, University of Waterloo, Project No. 109-17, 1982. van der Kamp, G. and R.N. Betcher, Dimensionless numbers to characterize interactions between heat flow and groundwater flow, paper presented to the Joint Annual Meeting of the GSC, MAC, and CGU, V i c t o r i a , B.C., Canada, May, 1983. van Orstrand, C.E., Temperature gradients, in Problems of Petroleum Geology, 989-1021, Am. Assoc. P e t r o l . Geol., Tulsa, Okla., 1934. Verge, M-J., A three-dimensional saturated-unsaturated groundwater flow model for p r a c t i c a l purposes, Unpubl. Phd. Thesis, Univ. of Waterloo, 1975. Welty, J.R., C.W. Wicks and R.E. Wilson, Fundamentals of Momentum, Heat, and Mass Transfer , 2nd. edition, Wiley, New York, 1976. Whitaker, S., The equation of motion in porous media, Chem. Eng. S c i . , 21, 291-300, 1966. Zienkiewicz, O.C., The F i n i t e Element Method , 3rd. ed., McGrawH i l l , New York, 1977.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- The thermal effects of three dimensional groundwater...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
The thermal effects of three dimensional groundwater flow Woodbury, Allan David 1983
pdf
Page Metadata
Item Metadata
Title | The thermal effects of three dimensional groundwater flow |
Creator |
Woodbury, Allan David |
Date Issued | 1983 |
Description | The proper measurement and understanding of heat flow values is essential in studies of tectonics (regional heat flow), oil maturation, low and high temperature geothermal environments, and earthquake research. Numerical solutions of the equations of fluid flow and heat transport in porous media are used to quantify the effects of three-dimensional regional groundwater flow on the thermal regime. The Galerkin Finite Element method is used to solve the coupled equations governing fluid and heat transfer. Simplex tetrahedral elements are used to subdivide the region. The resulting finite element equations are non-linear with temperature and a "leap frog" iteration technique is employed to solve the coupled equations. A series of computer simulations are carried out to investigate how typical three dimensional flow systems can influence heat flow measurements. Generic modeling is carried out on hypothetical basins with a 40 km separation between the regional topographic highs and lows. Emphasis is placed on understanding the conditions under which groundwater flows severely perturb the thermal field. The results of the sensitivity study indicate that the transition from a conduction-dominated to an advective-dominated thermal regime is sharp, and depends on: the length, depth, and width of the basin, water table configuration, permeability, existence of extensive or discontinuous aquifers, and hydraulic anisotropy. In addition, the advective threshold is a function of location within the basin. Simulations show that surface heat flow reflects the spatial distribution of recharging and discharging fluid associated with the three dimensional groundwater flow system. For the water table configurations studied, there is a high sensitivity of surface heat flow to large scale variations of water table topography. A significant difference is noted between two and three dimensional simulations of the same water table configuration. The results underscore the need to recognize the effects of large scale relief in the direction orthogonal to any section line. The lateral variation in surface heat flow is considered to be a poor indicator of active groundwater flow systems in an advectively-disturbed three dimensional basin. Suggestions are made for a two step strategy to be employed for locating boreholes for heat flow measurements. A number of closely spaced boreholes in several directions may be needed to correctly interpret heat flow measurements. Near-correct basal heat flux values can be measured in areas of the basin where there is a transition from regional recharge to discharge. At these locations heat transfer is found to be laterally advective and vertically conductive. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-04-22 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0052537 |
URI | http://hdl.handle.net/2429/24058 |
Degree |
Master of Science - MSc |
Program |
Geological Sciences |
Affiliation |
Science, Faculty of Earth, Ocean and Atmospheric Sciences, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-UBC_1983_A6_7 W66.pdf [ 7.2MB ]
- Metadata
- JSON: 831-1.0052537.json
- JSON-LD: 831-1.0052537-ld.json
- RDF/XML (Pretty): 831-1.0052537-rdf.xml
- RDF/JSON: 831-1.0052537-rdf.json
- Turtle: 831-1.0052537-turtle.txt
- N-Triples: 831-1.0052537-rdf-ntriples.txt
- Original Record: 831-1.0052537-source.json
- Full Text
- 831-1.0052537-fulltext.txt
- Citation
- 831-1.0052537.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0052537/manifest