THE THERMAL EFFECTS OF THREE DIMENSIONAL GROUNDWATER FLOW By ALLAN DAVID WOODBURY B.Sc. The University of British Columbia, 1973 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE THE FACULTY OF (The Department of n GRADUATE STUDIES Geological Sciences) We accept this Thesis as conforming to the required standard The University of British Columbia September 1983 Â©Allan David Woodbury, 1983 V I n p r e s e n t i n g t h i s t h e s i s i n p a r t i a l f u l f i l m e n t o f t h e r e q u i r e m e n t s f o r an a d v a n c e d d e g r e e a t t h e U n i v e r s i t y o f B r i t i s h C o l u m b i a , I a g r e e t h a t t h e L i b r a r y s h a l l make i t f r e e l y a v a i l a b l e f o r r e f e r e n c e a n d s t u d y . I f u r t h e r a g r e e t h a t p e r m i s s i o n f o r e x t e n s i v e c o p y i n g o f t h i s t h e s i s f o r s c h o l a r l y p u r p o s e s may be g r a n t e d by t h e h e a d o f my d e p a r t m e n t 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 . I t i s u n d e r s t o o d t h a t c o p y i n g o r p u b l i c a t i o n o f t h i s t h e s i s f o r f i n a n c i a l g a i n s h a l l n o t be a l l o w e d w i t h o u t my w r i t t e n p e r m i s s i o n . D e p a r t m e n t o f /ffe/o<;/C<zÂ£ SC/tfA/CfS The U n i v e r s i t y o f B r i t i s h C o l u m b i a 1956 Main M a l l V a n c o u v e r , C a n a d a V6T 1Y3 DE-6 (3/81) TO My Family i i i ABSTRACT The proper measurement and understanding of heat flow values is essential in studies of tectonics (regional heat flow), o i l maturation, low and high temperature geothermal environments, and earthquake research. Numerical solutions of the equations of flu 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 Finite 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" 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 f i e l d . 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, existance of extensive or discontinuous aquifers, and hydraulic anisotropy. In addition, the advective threshold is a function of location within the basin. i v 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. V TABLE OF CONTENTS Page LIST OF TABLES ' v i i LIST OF ILLUSTRATIONS v i i i ACKNOWLEDGEMENTS x i 1. INTRODUCTION 1 1.1 Uses of Heat Flow Measurements 1 1.2 Heat flow Determinations 2 1.3 Previous Work 5 1 .4 Present Work 11 2. GOVERNING EQUATIONS 13 2.1 Conservation Equations 13 2.2 Mass 14 2.3 Energy 17 2.4 Momentum 21 2.5 Potential 22 2.6 Constituative Relationships 24 2.7 Summary of Equations and Assumptions 25 3. NUMERICAL SOLUTION 29 3.1 Development of the Numerical Model 29 3.2 Numerical Solution Techniques 31 3.3 Galerkin Finite Elements 32 3.4 Basis Functions 35 3.5 Finite Element Formulation - Fluid Flow 38 3.6 Finite Element Formulation - Heat Transfer 42 3.7 Solution Techniques 43 3.8 Programming Considerations 45 3.9 Verification 46 4. RESULTS 53 4.1 Introduction 53 4.2 Boundary Conditions 55 4.3 Presentation of Model Results 56 4.4 Simulations 59 Series 1. - Simple Water Table 59 Series 2. - Small Scale Variations in Upper Surface. 71 Series 3. - Effects 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 v i TABLE OF CONTENTS...Continued Page 5. SUMMARY AND CONCLUSIONS 123 REFERENCES 131 ) v i i LIST OF TABLES Table Page 1. Review of available studies on thermal effects of groundwater flow (after Smith and Chapman, 1983) 7 2. Summary of equations 27 3. Outline of sensitivity analysis 54 4. Parameter values for fluid and thermal properties held fixed for a l l simulations 57 5. Groundwater recharge rates to hydrologic system with root square deviations from basal heat flux 114 v i i i LIST OF ILLUSTRATIONS Figure Page 1. Concepts of heat and fluid transfer in a Sedimentary Basin 6 2. Density and viscosity 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 finite element mesh used in simulations. Only element prism faces are shown 47 6. Case 4 verifications . 50 7. Case 5 verifications 51 8. (a) Conceptual illustration - simple water table configuration with 40 m relief in the Y direction, (b) Conceptual illustration - simple water table configuration with 80 m relief in the Y direction 60 9. Simple water table - 40 m relief in the Y direction 63 10. Surface heat flow - simple water table 40 m relief in the Y direction 64 11. Section plots through basin for the Y equals 40m relief case, k = 5 x 10~ 1 6 m2 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 relief in the Y direction 70 15. Simple water table - 80 m relief in the Y direction 72 16. Surface heat flow - simple water table 80 m relief 73 Figure LIST OF ILLUSTRATIONS...Continued ix Page 17. Root mean square deviation versus permeablity for the simple water table configuration - 80 m relief in the Y direction 74 18. Conceptual illustration - " f l a t - top" topography 75 19. Surface heat flow - "f l a t - top" topography 77 20. Comparison of two and three dimensional model results on a section line at Y equals 10 km - " flat - top " topography 78 21. Conceptual illustration - Local rise on the water table. 80 m relief 79 22. Surface heat flow - Local rise on water table 80 23. (a) Extensive aquifer - 40 m relief on the simple water table, (b) Partial aquifer case. 40 m relief on the simple water table. Aquifer terminates along ridge-line of basin, (c) Stratigraphic pinchout. 40 m relief 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 lines through heat flow for the extensive aquifer case - 40 m relief in Y 88 26. Surface heat flow along three section lines. 40 m r e l i e f , simple water table 89 27. Comparison between surface heat flow for a homogeneous case and a partial aquifer case 91 28. Surface heat flow - Stratigraphic pinchout. 80 m relief in the Y direction 93 29. Surface heat flow along three section lines for the pinchout number 2 case 94 30. Conceptual illustration - Concave water table configuration. 40 m relief in the Y direction 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 illustrations of the Segmented water table configuration, (a) 80 m relief (b) 40 m relief (c) 80 m relief upland - 8 m relief 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 relief segmented water table, k = 5 x 10" 1 6 m2 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 direction for the simple water table, (b) Surface heat flow along two sections in the X direction for the concave water table 120 x i ACKNOWLEGDEMENTS It is a great pleasure for me to express my appreciation to both Allan Freeze and Leslie Smith for their encouragement, invaluable advise and support. I am particularly indebted to Leslie 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 particularly like to thank for sharing many of his 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 illustrations in this thesis and prepared slides for various presentations. I would also like 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 terrestrial heat flow values is 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 is ample evidence to suggest that the thermal regime of a basin, and thus the surface heat flow, can be highly influenced by circulating groundwater. This thesis w i l l examine quantitatively the effects 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 different rock types can be significant, thus details of crustal thicknesses can be inferred. An example of the applications of heat flow measurements is found in Hyndman(1976), although there are many such examples to be found in the literature. 2 In many areas (Mase et a l . , 1978,1982) heat flow measurements are used to obtain information on geophysical anomalies, circulating hot water and magmatic activity. 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 is sensitive to time, heat and pressure (i.e., lower temperatures require longer time periods of burial). 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 is essential. The primary technique for determining the thermal regime is from surface heat flow data. Unfortunately, in hydrologically perturbed environments this determination may not be accurate. 1.2 Heat Flow Determinations Heat flow values are determined from temperature measurements in near surface boreholes. After a suitable time 3 period in which downhole temperatures have stabilized following a d r i l l i n g operation, temperature values are determined, usually with thermocoules or thermistors, at various points down a borehole. Temperature values must be corrected for various effects. 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 is 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 is computed and both laboratory and "in-situ" techniques are used. In the laboratory the divided bar method is 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 ^ w can be calculated. In the f i e l d "in-situ" 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 is then the geothermal gradient multiplied by the thermal conductivity of the rock/water units (Levorsen,1967). The principal uncertainty in calculating heat flow values is the determination of the geothermal gradient, especially in the presence of an active groundwater flow system. Determination of the gradient can be arbitrary, however, four different 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 line is fitted through a l l the points on the temperature depth plots and one representative gradient value for the entire hole is determined. Smith and Chapman (1983) warn that this approach can also lead to unsatisfactory results. (4) One could use methods suggested by Mansure and Reiter (1979) that take into account an idealized condition of one dimensional vertical flow between two fixed end temperatures. In actuality, groundwater flows are seldom simple and a complete treatment should account for the coupled nature of heat and fl 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 is controlled by hydraulic gradients and hydraulic conductivity with an additional flow induced by thermal effects (buoyancy). Energy or heat is transferred through a system by conduction and convection. Conductive transport occurs in the rock/water mass and is 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 identified as 5 convection and advection (usually referred to in the literature as free and forced convection). The former is related to the circulation of water due to forces caused by buoyancy, and the latter, the flow of heat and groundwater under some topographically driven system. It is now recognized (Smith and Chapman,1983) that neither conduction nor advection totally dominates in a l l geologic environments. For example, in tectonically stable areas of low relief and in areas of crystalline basement rocks at surface, temperature fields would be almost totally conductive. Contrasting with these conductive cases would be flow in permeable sedimentary units overlying basement rocks, and in tectonically 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 activity at depth. Figure 1 is an ill 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 will summerize some of the more important papers. A descriptive and comparitive study is given in Table 1. Early heat flow investigators (Van Orstrand, 1934, Bullard, 1939, Birch, 1947) a l l found that pure heat conduction could not q = K Surface heat flux Groundwater flow + + + +â€¢+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + A + + + + + + + + + + + A + + + + + + + + + + + + + + A+ -bflW- + + V + \ \ \ V + V + \ \ T + + Basal heat flux + X % + + + + + + + + + + + + â€¢ + I -1- ' Figure 1. Concepts of heat and f l u i d transfer i n a Sedimentary Basin. Table 1. Review of available studies on thermal effect (Adapted from Smith and Chapman, 1983) N AUTHOR DATE DIMENSION ID 2D 3D 1 Bredehoeft & Papadopulos 1965 X X X X X X X 2 Stallman 1967 X X X X X X X 3 C a r t w r i g h t 1970 X X X X X X X Sorey 1971 X X X X X X X 5 Mansure & R e i t e r 1979 X X X X X X X X 6 Domenico & P a l c i a u s k a s 1973 X X X X X X X X 7 Morgan e t a l . 1981 X X X X X X X X 8 Bodvarsson 1972 X X X X X X X 9 Bodvarsson 1973 X X X X X X X 10 C a r t w r i g h t 1971 X X X X X X X 11 Lewis & Beck 1977 X X X X X X X X 12 Keys & Brown 1978 X X X X X X X 13 K i l t y & Chapman 1980 X X X X X X X X 14 Ziagos & B l a c k w e l l 1981 X X X X X X X 15 B r o t t . e t a l . 1981 X X X X X X 16 Parsons 1970 X X X X X X X 17 Sorey 1976 X X X X X X X X X 18 Betcher 1977 X X X X X X X 19 Andrews 1978 X X X X X X X X X 20 Andrews & Anderson 1979 X X X X X X X X X 21 Smith & Chapman 1983 X X X X X X X X X X 22 Woodbury 1983 X X X X X X X X X X â€¢a 01 01 4-> o 13 CO c 01 rH QJ 4-1 S 3 u c to 3 CJ 01 Ol 01 3 CO rH CH 6 4 J CO o CJ CQ CO UH 01 CO 3 01 â€¢rl CO CJ â€¢H rH 4-t 4 J o c cj a o T) 0) w G 01 01 â€¢H o >. >, 4 J â€¢rl 01 01 c 0 0 a H 4 J 4J â€¢H ^ u 01 -rt 01 o o 4 J â€¢H â€¢ri CO 4 J r-. u 4-1 T> CO GO u I-I o o a o -H rH â€¢H â€¢H co c o 01 u CO o o CJ CO CO c c 01 CO B u . o Tt r-C rH CO c c TH â€¢H 4 J 1-1 o 01 CO c Ol 01 â€¢H 01 CO CH Â«4H CO 4-1 J= J C â€¢rl CO > > > a u B co cu oi 4 J J2 of groundwater flow. APPLICATION v e r t i c a l groundwater flow modify r e f . 1 f o r low v e l o c i t i e s apply r e f . 1 to estimate groundwater d i s c h a r g e apply r e f . 1 to d r i l l h o l e temperature p r o f i l e s modify r e f . 1 to study heat flow v a r i a t i o n s i n a d r i l l h o l e temperature f i e l d s i n two-dimensional b a s i n apply r e f . 6 to Rio Grande r i f t b a s i n s , perm, t h r e s h o l d model temperature of water i n j e c t e d i n t o a q u i f e r v e r t . temp, i n v e r s i o n s by t r a n s i e n t e f f e c t s of horz. flow temperature i n c o n f i n i n g l a y e r above a q u i f e r model down-dip p e r c o l a t i o n of c o l d s u r f a c e waters apply r e f . 8 & 9 to tr a c e i n j e c t e d water i n a q u i f e r model heat flow v a r i a t i o n s by assumed flow i n a q u i f e r s t r a n s i e n t temp, i n c o n f i n i n g beds surrounding h o r z . a q u i f e model subsurface temperatures i n Snake R i v e r p l a i n f l u i d flow and temperatures i n a g l a c i a l complex f l u i d flow, d i s c h a r g e , temperature i n Long V a l l e y c a l d e r a g e n e r i c modeling of temperature f i e l d temperature i n a q u i f e r d u r i n g o p e r a t i o n of heat pump thermal a l t e r a t i o n of groundwater, s u r f a c e water seepage g e n e r i c modeling of thermal regime of b a s i n s g e n e r i c modeling o f thermal regime of b a s i n s 8 explain observed heat flow patterns in the f i e l d . Groundwater movement was cited 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 fluid 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 heat/fluid transfer equations in a two dimensional system. He studied the growth and distribution of convective thermal cells with application to thermal regimes in New Zealand. This particular model assumed a constant flow f i e l d and constant fluid properties. A widely used formulation describing the thermal aspects of groundwater flow is 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 specific discharge of the fluid can be computed. The boundary value problem solved was the steady flow of heat and fluid 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 so i l temperatures. Parsons (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 fluid properties. Results 9 showed that groundwater flows exerted a controlling influence on the shallow temperature fields in a glacial complex in northern Ontario, Canada. Another important contribution to this f i e l d is the presentation by Domenico and Pakiausicas (1973), who solved the heat/fluid equations analytically for steady flow with constant fluid properties, in a manner reflective of Toth's (1962,1963) work in regional groundwater studies. The importance of this paper is 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 distribution. The limitations of their approach include constant fluid properties and the limited range of validity of their linearized solutions. There have also been many studies by other authors (Lewis and Beck, 1977, Kilty and Chapman, 1980, Brott et al.,1981 ) who have attempted to deal with two dimensional problems, but have assumed either a known velocity 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 limitations. In the more recent literature, a general classification of three different types of studies is possible. One class of studies is 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 third 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 this group. Betcher (1977) carried out a sensitivity 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 fl u i d properties. Sorey (1978) developed a numerical solution to the coupled equations using an integrated fi n i t e difference scheme. This particular 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. Fluid properties were coupled to temperatures in these studies. Smith and Chapman (1983) provided the most complete treatment to date on the thermal effects 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 carried out on a basin with varying topography totalling 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 is defined as when the root mean square deviation of the surface heat flow with the basal heat flow exceeds an arbitrary l i m i t ) . Results showed that a large range (0 to almost 100%) of the basin can have surface heat flow values significantly 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 vertical permeability, and the existence and locations of aquifers. As a result of their work, i t is clear that significant groundwater disturbance may mask or distort 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 site, including the water table configuration and subsurface flow system, combined with more closely 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 this thesis consists of the development of a fully 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 collected on cross 1 2 sections and flow systems are seldom two dimensional. (2) Performing a generic study of hypothetical three dimensional basins to determine the sensitivity 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 will be developed and checked against known one dimensional solutions, then compared with two dimensional simulations using Smith and Chapman's (1983) computer program. This verification stage will be followed by simulations of hypothetical basins examining the three dimensional effects of: varying permeability and water table topography, extensive and discontinuous aquifers, reduced depth of basin and anisotropic permeability. Emphasis wi 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. Field techniques for heat flow deteminations in hydrologically perturbed environments are introduced, along with suggestions for future research. 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 fluid transfer must be based on the fundamental equations describing conservation of mass, energy, and momentum. The continuum approach wi l l be adopted in the development and solution of the differential equations presented in this 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 level to a macroscopic continuum level. The continuum approach is valid 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 is large, the approach may not be valid (Long 14 et a l . , 1982). Many presentations and derivations of these equations can be found in the literature. 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 is 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. Li (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 individually. Resonable simplifications are made to each equation, pertinent to the problem being solved. A l l of the following equations are valid for a single-phase Newtonian fluid in a saturated porous medium. Solute concentrations are assumed to be negligible and cross-coupling phenomena such as the Soret and Dufour effects (Bear, 1972, Freeze and Cherry, 1979) are also assumed to be negligible. 2.2 Mass The equation of fl u i d flow follows directly from the theory that matter can neither be created nor destroyed. The fl u i d (macroscopic) mass conservation equation for saturated flow with 1 5 vertical deformations and with temperature dependant properties is (Bear and Corapcioglu, 1981, p725): (2-0 v - < j r < C*PM4)H + {^-^)Â¥t Where: <|â€ž = fy- />[VS Cj^y. = the relative specific fluid discharge {L/t} = flu i d specific discharge {L/t} c7Cp = compressibily of solid at a constant temperature {Lt2/M} tf[ = porosity {dimensionless} V 5 = average velocity of the solid {L/t} y^ p = fluid's coefficient of compressibility at a constant temperature {Lt2/M} P = fluid pressure {M/Lt2} OCr = thermal expansivity of the solid at a constant pressure {T~1} Pr = thermal expansivity of the flu i d at a constant pressure {T~1} ~]~ = temperature {T} V . - operate, ( J A + j ^ ,12) The f i r s t term of equation (2-1) represents the amount of fluid entering or leaving a control volume. The second term represents the time rate of change of fl u i d mass due to pressure changes (contributions due to expansion of water and the porous medium). The third term represents the time rate of change of 16 fluid temperature due to thermal expansivity of water and the solid matrix (this term could also be viewed as a fluid source/sink term due to changes in temperature). The last term is an additional f l u i d source due to thermally induced fl u i d density changes. Only steady state pressure and temperature fields will be considered. This is 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 fields. Because only regional and not local thermal fields are being considered a steady state thermal f i e l d is considered appropriate for this 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 fl u i d {M/L3} and substituting (2-2A) into (2-2) leads to: 1 7 (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 is carried through a cycle, the total heat added to the system from i t s surroundings is proportional to the work done by the system on its 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*< 4 ('-"OpÂ» cs where: = heat capacity per unit volume of the porous media {M/Lt2T} ^v-i = specific heat capacity of the flu i d {L 2/Tt 2} r = specific heat capacity of the solid {L 2/Tt 2} = density of the solid {M/L3} = coefficient of thermal conductivity of the porous media {ML/t3T} = thermal conductivity of the fl u i d {ML/t3T} " l ^ s = thermal conductivity of the solid {ML/t3T} Â£ = dilatation of the porous media {L/L}'- = ^ "^s u a r 3 0S = stress in the solid {M/Lt2} This equation neglects thermal micro-dispersion, sources or sinks, and heat produced by radioactive decay or chemical reactions. Inertial terms are assumed to be smaller than time rate of changes ( i.e. Vj* * V"7" ^ ). porous medium and fl u i d temperatures are assumed to be in local 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 is the flux of thermal energy by conduction. The third 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 total dissapation of energy through the action of viscosity on the 19 porous media.) At steady state (2-4) reduces to: ( *-5) V-l\m-VT- ^0,^.77 - O neglecting viscous dissapation with no sources or sinks (Garven ,1982). Pure heat conduction is described by: ( 2 - 5Â°) V - T V â€¢ V T = O Various authors (Pickens and Grisak, 1979, Garven, 1982, Smith and Chapman, 1983) have redefined "Piâ„¢ as a combined tensor composed of two components; a purely conductive term and a velocity 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 is to model i t as a diffusive process. For heat transport the dispersion tensor is (Bear, 1979) o O k - 0 F ; O o o cc r CC-= longitudinal dispersivity{L} = transverse dispersivity{L} 20 This equation is written in its simplest form (groundwater is flowing parallel to the x axis). Part of the conceptual problem of modeling dispersion as a diffusive process is 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 sites where values of cCL 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 solid grains and with the moving water. Theoretically, because equation (2-5) assumes that the temperatures in the solid 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 exist. 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 diffusive transport. Thus, macro-dispersion could also be explained by lack of knowledge of variations in thermal conductivity at a f i e l d site. Van der Kamp (1983) provides a good discussion of macro-dispersion in heat transport. Field observations of this phenomena are described in Sauty et a l . (1982a,b). , Smith and Chapman (1983) modeled two dimensional systems including high dispersivity 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 will be 2.4 Momentum In groundwater flow problems Darcy's Law satisfys the momentum requirements. This law is basically 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 limit of between 1 and 10 which is applicable to most groundwater problems. Bear (1972) expresses the generalized form of Darcy's law as: developed. r where: specific f l u i d discharge vector {L/t} p f l u i d dynamic viscosity {M/Lt} flu i d pressure {M/Lt2} acceleration of gravity {L/t2} 2 elevation of reference point above a datum plane {L} k permeability tensor {L2} flui d potential function {L} 22 In this form, local in 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 flu i d potential (j) must be described. Hubbert (1940) derived an energy potential 'o for isothermal flow where: is a compressible fl u i d density{M/L3} The usual form of equation (2-8) is : p 1 ( d f 3 ) , f If p is constant, then: ( i - i Â« ' ) h = * + ^ However, because we are dealing with nonisothermal flows, hydraulic head as expressed in equation (2-9) cannot be used as ^ in equation (2-7) because a unique potential would not exist. One would normally have to work with the pressure form of equation (2-7) as is 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 0 such that: (2-10) h c = _ L + 2 where p o is some reference density at a constant temperature. Frind (1980) defines h 0 as the "equivalent freshwater head" and states that numerical accuracies are improved through the use of h Q instead of pressure. However, this is not the elevation that water would rise to in an open standpipe. Substituting (2-10) into (2-7) where {r TÂ° ~ r e l a t i v e density {dimensionless} The equation describing fluid motion is now separated into two parts; the f i r s t resulting from head differences (forced convection) and the second an additional vertical component caused by buoyancy forces of a f l u i d particle ^ imbedded in a fluid of density po (free convection). By substituting (2-11) into (2-3) the fluid conservation equation: 24 we arrive at the final fluid continuity equation. 2.6 Constitutive Relationships Fluid density and viscosity vary substantually with temperature. Strauss and Schubert (1977) have shown the importance of incorporating temperature dependance of density and viscosity in the equations describing heat and fluid flow. Furthermore, Strauss and Schubert (1977), Sorey (1976), and Garven (1982) present evidence that the pressure dependence of these two fl u i d properties may resonably be neglected for the type of problems being addressed in this study. A l l other fl u i d and material properties such as thermal conductivity of the fl 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. 2 5 Two of the more simple relationships available describing this dependence are: (2-13) f , = p . [ I - / ? ' ( T - T Â» ) * y'(r-zy] J2 = 3 . 17 x i o - " c - 1 = 2 . 5 6 x 1 0 " 6 C"2 C, = 2 . 4 x 1 0 - 5 / \ r 2<tr.i7 / (~r + ( a - ' * ; y = C, (_ /O J The viscosity 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 fl u i d transport in a porous media have been presented (Table 2 ) . Before proceeding onto Chapter 3 i t is useful to summarize briefly the assumptions in the governing equat ions: ( 1 ) The porous medium is assumed to be a continuum, fully saturated with a single phase Newtonian fluid (water) with negligible solute concentrations. Fluid flow rates are small enough to consider the flow laminar. D E N S I T Y X I O 3 (kg m"3 ) 9Z 27 Table 2. Summary of equations. GOVERNING EQUATIONS Fluid Continuity V * (p f q) = 0 Thermal Energy V - A m - V T - p f C v f q â€¢ V T = 0 Momentum q* = â€” T 2 - ( V h 0 + p r v Z ) Potential Equations of State p = p(T) ^=/u(T) 28 (2) Cross-coupling phenomena such as fluid pressures adding an additional driving force to heat flow or vise versa are assumed to be neglibly small. It is recognized that density and viscosity 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 is neglected. (5) Temperatures of the fl u i d and solid media are assumed to be in local thermal equilibrium. (6) No sources or sinks are present in either the flu 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 wi 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 Value/Initial Value problem in groundwater hydrology should proceed in a four step manner as described by Freeze (1978): (1) Identification 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 light of geologic data and model calibration. Steps (1) and (2) have been presented in the last two chapters and this chapter will 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) Electric 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 (Dirichlet). 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 directly. 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. Electric 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 di g i t a l (computer) techniques in aquifer simuations. In cases such as the problems that will be solved in this study an analog model cannot be used. Analytic solutions to groundwater problems are usually in a closed form or in a series/intergral form that has to be numerically evaluated. Analytic solutions are only available where the region of flow is 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 limitations to these models are evident and these will 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, fi n i t e element, integrated finite differece, boundary intergral equation method, random walk particle models, and the method of characteristics. Of these six techniques integrated f i n i t e difference, fi n i t e difference, and finite element are the most suitable to solve the equations presented in Chapter 2. The method of characteristics has been used to solve groundwater contaminant problems (Konikow and Bredeheoft, 1978), however it s use is 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 particle models are generally used for transient contaminant transport problems (Garven,1982), although i t is theoretically possible to solve the heat transport equation by this method. Boundary intergral equation methods are currently one of the newer numerical techniques, and have been applied to transient free surface flows (Liggett and Liu, 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 this technique. Integrated finite 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 fin i t e difference approxiation and cannot easily incorporate general tensorial quantities. Finite 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 finite difference to non-linear unsaturated flow in three dimensions. This method is suitable for many problems, but i t is d i f f i c u l t to discretize an irregularly shaped region and flow velocities are not easily calculated. The fin i t e element technique is the most popular method for most groundwater problems and is the method that will be used to solve the governing equations in this study. 3.3 Galerkin Finite Elements The Galerkin-based finite element method is one of a number of fi n i t e element techniques that have been applied in groundwater hydrology. The interested reader is referred to Pinder and Gray (1977) for details on other methods. The Galerkin method has received considerable attention in the literature. Examples of the application of the technique to advective heat transport problems include Mercer et. a l . (1975), Andrews and Anderson (1979), Li (1980), Garven (1982) and Smith and Chapman (1983). The Finite Element approach involves subdividing the region of flow into a grid or mesh containing small blocks or cells called elements. Element corner points are called nodes. Properties of the porous media such as porosity, permeability, 33 Figure 3. Basic tetraheclral element. (Adapted froir. Segerlind, 1976) (a) and (b) Assemblage of elements. (Adapted from Zienkiewicz, 1977) 34 thermal conductivity etc., are assigned to each element and the solution to the dependant variables (head, temperature) is obtained at each of the nodal points. The dependant variables are assumed to vary in some functional way between nodes so that quantities may be easily calculated at any point within an element. The simplest variation is 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 is a tetrahedron with straight sides (Figure 3). This element is called a simplex three dimensional element when the nodal variation of the dependant variable is linear (Segerlind, 1976,p.23). There have been only a few fully three dimensional Finite 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 is that fewer nodes and hence less computer core storage is needed; however, their disadvantage is that increased computer time is needed (see sect. 3.7). Current computer capacities are larger today than when Verge (1975) or Segol (1976) developed their programs; consequently it is now 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 efficient for two dimensional problems. Therefore, the simplex 35 tetfahedral shaped element will be used for solving the boundary value problems presented in this study. 3.4 Basis Functions A basis function is 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 linearity of the solution between nodal points is introduced. We can define an approximate solution to dependant variables h 0 and T as: A where: A/; = N- (x, ^ 2 ) 1 = 1,2,3 n nodes A A and ha and T represent approximate solutions at each nodal point. Since we have adopted a linear variation of h0 and T between nodes: A hi = / \ x + By + O t D f = . /\x t By + C2 * D where A,B,C,D are coefficients to be determined. In the next few pages the basis functions for hydraulic head 36 will be developed. The basis functions for temperature are exactly the same. Hence, at each node (since four nodes form a tetrahedral element): = / U ; + B y , t a + D -- A % i â€¢ B y , 4 C z j + D = A i ( t B y Â« * C z * + D or in matrix form: h ; h . X ; V, 2 . / J D A 6 c L D J The values of the linear coefficients A,B,C,D can be found by inversion of the above system, leading to: A = ' A v ( C l l . h ; +. 3 ^ 4 61, h K t QLK ) where V = Volume of the tetrahedral element and: ( w ) a; - ( ( z A * y K - y, *z<) - Z j * ( y L - y K)) a J = +((z L*y K - yL*2*> - y ; * ( z L - z K) + zi*(y, - y K)) = -((z L*y- - y L* Z j) - y ; * ( z L " Zj ) + z; *(y L - y 3)) aL = + ( ( z K * y j - yÂ«* z3) - yj * ( z K - zj ) â€¢ + z;*(y K - y j )) bl = + (z,*x K - z K*x c) - X j * <Zc - ZK> + zj * UL - xK) b j = x L*z K) - x; * (z L - z K) + z; * (xL - xK)) b K +(z L*xj - Z j * X L ) - x; * ( z L - z- ) + z; * (xL - xj) = -(z k*xj - Zj *xK ) + x] * ( z K - Zj ) - z; * (x K - xj) c ; = -(x 3*(y K - y t ) - yj *(x< - x L) + (y4 % - xi*y<>> +(x ;*(y K - y L ) - Yj * ( X < - x L) + (yt * xÂ« - x , % ) ) = -(x; *(yj - y L } " Y.- *(xj - x L) + (y. *xj - x L* Y ( /)) = + (x;*( y j - y K > - y , * ( X J - x k) + (X j *y< - x K * y j * * = xj* (y<* 2L . â€¢ W - y. *(x K*z L - xL *z K) + z j * ( y L * X < V <> = x; *(yK*zu â€¢ - yL*zK\ " y; * ( \ * z c - xL *z K) + z ; *(y L*x < x,*yk) d K = x;*(z L*yj - y L*z 3 ) - y; Mx. *zL - x L*zj) + z ;*(y (_*xj -x L*yj) d L = x;*(z K*yj - y K * Z j ) - y i * ( x J * z K - x K*zj> + z'My^xj ' v = 7^ [ d<- ^ 4 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 C; 2 + c f j ) h ; f f Q j x + b j y t Cj 2 i dj) h j t A or h = A/; h ; + A/j hj 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 fl u i d flow equation. 3.5 Finite Element Formulation - Fluid Flow Equation We define an operator L(h 0) based on our final f l u i d flow equation (2-12) as: We define k as = P/ k fÂ°3 At this point we note that i f the f i n i t e element grid is set up so that the principal directions of anisotropy in K are alligned 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 is to find an approximate solution to this equation, To accomplish this we multiply (3-9) by a weighting function W;, and integrate over the flow domain: ^ L (K) ^ oi2 = o The reader is referred to Segerlind (1976) for a complete treatment of this method. The Galerkin requirement (Pinder and Gray, 1977) for this system is when W; is set to the basis functions N; : A / ; dx dl = 0 Substituting equations (3-8), (3-9) into (3-10) and applying Green's f i r s t identity, the resulting set of matrix equations can be written as (Garven,1982): (*-") A; k + B. = Qâ€ž 40 where A m w is an n x n matrix and B* and Q n are vectors of size n. The matrix kmn is usually referred to as the global stiffness matrix and is composed of the terms: where V is a volume intergral. This matrix is a sparse, symmetric, banded matrix that is diagonally dominent and positive definite. The bouyancy term B h i s : .2 Mass fl u i d flux boundaries can be imposed through the vector ( S is a suface intergral) 41 where 71 is the area of the element on which flux Jj is imposed on. Elements of Q h are set to zero i f no boundaries in the fini t e element mesh are of the specified flux type or i f a node is on an impermeable boundary. Specific discharges can be calculated for each element by substituting (3-8) into (2-11). The result i s : These three quantities will be needed for each element in order to solve equation (2-5) . Because the solution for h 0 is based on a linear interpolation between nodes, velocities 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 discontinuities in velocities between adjacent elements gradually vanish as the element size is reduced. Higher order basis functions, such as hermite polynomials, allow for a continuous velocity distribution between elements, however substantial more computational time is required because the intergrals expressed in (2-12,13,14) would have to be evaluated numerically. o 42 3.6 Finite 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 align the coordinate axes with the principal directions of anisotropy in then: In an identical manner to that of the fluid flow equation, (3-16) is substituted into the Galerkin requirement (3-10) and Green's f i r s t identity is applied. The final result is another set of matrix equations: ( * - ' 7 ) where is an n x n stiffness matrix and F^ is a vector of applied fluxes of size n. The stiffness matrix S w / V can be expressed as: 43 This matrix is also banded and positive-definite but is nonsymmetric. As terms q , q , q, increase in magnitude the matrix Sl7,h may no longer be diagonally dominant. This behaviour can lead to innaccuracies in the solution when normal Gaussian Elimination techniques are applied to solve (3-17) (Gambolati,1980). The conductive flux vector Fn is expressed as: the area of the element face. Like the fl u i d case; vanishes on insulated boundaries or i f no flux is 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 this non-linear problem is called a "leap-frog" iteration technique and is shown on Figure 4. The following algorithm is similar to that described in Smith and Chapman (1983). The f i r s t step in the proceedure after the generation of a fi n i t e elment mesh is to solve this system of equations assuming a conductive heat flow regime. Appropriate densities and viscosites are assigned to each element using equations (2-13) and (2-14). The flui d flow equation (3-11) is solved to give values of equivalent freshwater heads at each nodal point, following which equations (3-15) are employed to calculate specific discharges where S denotes a surface intergral of an element with a face on a flux surface. The normal flux is and 9i is 44 PRINT DATA E C H O PRINT N O D A L C O O R D I N A T E S ( START ) R E A D IN DATA G E N E R A T E M E S H E S I C O M P U T E T E M P E R A T U R E S DETERMINE FLUID P R O P E R T I E S PRINT H E A D S / FOR N O D E S /**"" PRINT E L E M E N T S D A R C Y V E L O C I T Y PRINT T E M P E R A T U R E S F O R N O D E S r C O M P U T E HYDRAULIC H E A D S C A L C U L A T E D A R C Y V E L O C I T I E S \ C O M P U T E T E M P E R A T U R E S 1 NO C H E C K FOR C O N V E R G E N C E IF N E E D E D P L O T S U R F A C E HEAT F L O W , T O P O G R A P H Y Y E S P L O T C R O S S S E C T I O N DATA T H E R M A L P R O F I L E S ( S T O P ) Figure 4 . Flow chart showing "leap-frog" algorithm. 45 q t , q 5 , q ? within each element. Next, the heat transfer equation (3-17) is solved to determine nodal temperature values. For each element, densities and viscosities are updated using this new estimate of the temperature f i e l d . An iterative procedure is then employed; the f l u i d flow and heat transport equations are solved repeatedly with densities and viscosities updated until the maximum temperature change between iterations is less than a specified tolerance. For the three dimensional simulations that will be reported in Chapter 4, this 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 British 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 five 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 typical of the form (3-11) or (3-17). The interested reader is referred to Remson et. a l . (1971), or Pinder and Gray (1977) for details 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 transfer). 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 efficient nodal numbering, large band widths result from the choice of elements used, limiting the size and mesh density of the problems considered. Because of these problems another version of the same program was developed that employed U.B.C SPARSPAK (Nicol, 1982) 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 limiting the simulations reported in this study. Typical meshes consisted of 19 nodal spacings in the X direction, 10 in Z and 8 in Y for a total of 1980 nodes and 7600 elements ( A typical section through a f i n i t e element mesh is shown on Figure 5. Here only prism faces are shown). This entire mesh required 6.8 megabytes of storage (computer limit: 7.0 megabytes) for the LINPAK version of the program. For advectively disturbed cases, of which most of the simulations are, non-linearities in density and viscosity 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 Figure 5. Typical f i n i t e element mesh used i n simulations. Only element prism faces are shown. 48 the numerical code is 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 fluid (2-3) and heat conduction equations (2-5a). In this case a one dimensional grid with homogeneous, constant fluid 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 fl u i d (2-3) and heat conduction equations (2-5a). In this case a one dimensional grid with homogeneous medium parameters and constant f l u i d properties with constant fl 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 fluid (2-3) and heat conduction equations (2-5a). A simulation of a one dimensional heterogeneous (the system was divided into two regions with two different 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 flu i d (2-3) and heat transfer (2-5) equations. In a series of computer runs, a one dimensional system with constant fl u i d properties and homogeneous media parameters (constant head and temperature 49 b o u n d a r i e s a t e i t h e r e n d s ) were c o m p a r e d t o an a n a l y t i c s o l u t i o n s . P e r m e a b i l i t i e s were v a r i e d i n f o u r s i m u l a t i o n s ( F i g u r e 6) t o c h e c k f o u r d i f f e r e n t v e l o c i t y c a s e s . Good a g r e e m e n t was n o t e d . C a s e ( 5 ) - T e s t o f t h e c o m b i n e d f l u i d ( 2 - 3 ) and h e a t t r a n s f e r e q u a t i o n s ( 2 - 5 ) . I n a s e r i e s o f r u n s a one d i m e n s i o n a l s y s t e m w i t h c o n s t a n t f l u i d p r o p e r t i e s a n d homogeneous m e d i a p a r a m e t e r s ( c o n s t a n t h e a d a t b o t h e n d s , c o n s t a n t t e m p e r a t u r e a t one end a n d a h e a t f l u x c o n d i t i o n a t t h e o t h e r ) were c o m p a r e d t o a n a l y t i c s o l u t i o n s ( F i g u r e 7 ) . Good a g r e e m e n t s were n o t e d . C a s e ( 6 ) - T e s t o f t h e f u l l y c o u p l e d f l u i d ( 2 - 3 ) , h e a t t r a n s f e r ( 2 - 5 ) , a n d c o u p l e d e q u a t i o n s o f s t a t e ( 2 - 1 3 ) , ( 2 -1 4 ) . Two c o m p a r i s o n s were made b e t w e e n t h e t h r e e d i m e n s i o n a l m o d e l a n d S m i t h a n d Chapman's ( 1 9 8 3 ) c o m p u t e r p r o g r a m . A two d i m e n s i o n a l s y s t e m was m o d e l e d a n d f l u i d p r o p e r t i e s were a l l o w e d t o v a r y w i t h t e m p e r a t u r e i n t h e s e c a s e s . C o m p a r i s o n s w ere made f o r a homogeneous, i s o t r o p i c c a s e a n d a homogeneous, a n i s o t r o p i c c a s e ; b o t h a t t h e a d v e c t i v e t h r e s h o l d f o r t h e s y s t e m b e i n g m o d e l e d . The maximum d i f f e r e n c e f o r t h e i s o t r o p i c c a s e was 1.3% w i t h an a v e r a g e d i f f e r e n c e o f 0.14%. Maximum d i f f e r e n c e f o r t h e a n i s o t o p i c c a s e was 0.05% w i t h an a v e r a g e d i f f e r e n c e o f 0.004%. B o t h 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 i n t h e d i s c h a r g e e nd o f t h e s y s t e m w h i c h i s t h e a r e a o f h i g h e s t g r o u n d w a t e r v e l o c i t i e s a n d l a r g e s t i n n a c c u r a c i e s i n s o l u t i o n . D i f f e r e n c e s i n t h e two p r o g r a m s were c o n s i d e r e d n e g l i g i b l e a n d p r o b a b l y due t o n u m e r i c a l r o u n d o f f 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 a t h i r d d i m e n s i o n . 50 T E M P E R A T U R E ( Â° C ) Figure 6. Case 4 v e r i f i c a t i o n s . 51 T E M P E R A T U R E ( Â° C ) Figure 7. Case 5 verifications. 52 With the six verification cases completed the program was considered to be fully operational and i t was decided to proceed with the production runs that will 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 is to assess how typical three dimensional flow systems can influence heat flow measurements. The method used to accomplish this goal will be to carry out a sensitivity study using the mathematical model developed in Chapter 3. This study w i l l be similar to the one carried out by Smith and Chapman (1983) and wil l consist of generic modeling on hypothetical basins. Such models are designed to be representative of real f i e l d situations, but are kept simple in order to more easily understand the phenomena being studied. A range of rea l i s t i c input data are chosen, and the model results are assessed for sensitivity to the variation in a parameter while others are held fixed. The parameters that w i l l be investigated in this study are: variations in permeability and water table configuration, intoducing more permeable units at depth, partial and discontinuous aquifer cases, anisotropic formations, and reduced depths of active flow (Table 3). 54 Table 3. Outline of sensitivity analysis. 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 will be modeled, but the boundary conditions will 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 is assumed to be a subdued replica of the surface topography. At any point on this surface the value of equivalent freshwater head is equal to the elevation of that point. An isothermal temperature distribution is also imposed along the water table. A conductive heat flux is prescribed on the basal boundary, which is assumed to be impervious to fl u i d flow. At the basin divides along the sides of the model, a condition of zero lateral heat and fluid flux is 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 distributions are then more reflective 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 is the controlling 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 is recognized that this temperature boundary condition is more complicated and.represents the interplay of heat flow in the unsaturated zone, vegetation 56 types, snow loading, g l a c i a l rebound, radiation, and other transient e f f e c t s (Sekioka and Uhara, 1978), including the existence of thermal springs. However, because the depths to which these effects infuence the thermal regime are r e l a t i v e l y small with respect to the overa l l depth of the basin, the assumption of an isothermal temperature boundary i s v a l i d . Values of material and f l u i d properties used in the simulations are summarized in Table 4. Basal heat flux values imposed on the model are set to 60 mWnr2, a value near the continental average as reported by Jessop et. a l . (1976). Unless otherwise stated, porosity is assumed to decrease with depth in a piecewise manner from 0.20 at the surface to 0.06 at the base of the model. The thermal conductivity of the s o l i d i s constant in the simulations, but the thermal conductivity of the solid/water mass does vary with porosity, and hence, depth. 4.3 Presentation of Model Results Results from simulations w i l l be presented in the following manner. Plan view contour maps of the topographic surface w i l l be presented along with three dimensional schematic representations of the surface. Output of the computer model consists of contour maps of surface heat flow, cross-sections showing surface heat flow versus distance, and contoured isotherms on sections. In addition, temperature-depth and gradient-depth plots w i l l be presented at selected locations within'the basin. For consistency in re s u l t s , surface heat flow values are computed by multiplying the near-surface thermal conductivity by TABLE 4. Parameter Values for flu i d and thermal properties held fixed for a l l simulations Property Reference Temperature Reference Density of fluid Specific heat of fluid Thermal conductivity of flu i d Thermal conductivity of rock Basal heat flux Value 10Â°C 997.97 Kgrn'3 4186.0 JKg-1C"1 0. 586 Wm" 1 C 1 2.51 Wm-'C"1 60.0 mWnr2 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 interval. Discrepancies could exist between model results and heat flows determined from deeper boreholes, however model results will represent the "worst case", as heat flow results from deeper holes will tend to be closer to the basal heat flux. Surface heat flow is contoured by U.B.C.'s SURFACE contouring package. In some cases visual smoothing was necessary to f i l t e r out distortions caused by the gridding approximations. Inflow fl 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 aid in assessing the effects of various topographic and geologic pertubations on the thermal regime. Before proceeding to the results, i t is necessary to define some terms that w i l l be used throughout the next few sections: Minimumly Disturbed Area This region is defined as the surface area of the basin in which surface heat flow ranges between 50 - 70 mWnr2, or Â± 16% of the basal heat flux. The selection of these values is arbitrary, but reflects that portion of the basin which could be considered as fal 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 its primary use is to highlight differences in surface heat flow between simulations as different hydrologic 59 disturbances are imposed on the model. Conductive Case In this case permeability, and hence fluid velocities, are so low that almost a l l the heat transfer is by conduction. Surface heat flow then varies across the basin only because of topographic effects. In these cases 100% of the basin is within 57-63 mWm~2 or Â± 5% of the background flux for the cases that will be considered. Advectively Disturbed Case In this 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 is greater than 10% of the imposed value or, Â± 6 mWm-2 (Smith and Chapman, 1983). Heat redistribution 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 sensitivity 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 levels: from 1 x 10" 1 6 to 3 x 10" 1 6 m2 and fin a l l y , 5 x 10" 1 6 m2. Heat flow values on different section lines w i l l be compared with one another to assess the three dimensional effect. In a fi n a l suite of simulations (for this topographic case) the total 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 in the Y direction. 61 Figure 8. (Continued) (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 i n the Y d i r e c t i o n . 62 increase in relief will demonstrate the sensitivity of the basin heat flows to increases in groundwater velocity in the Y direction. Figures 8a and b show the two different rel i e f s of 40 m and 80 m over 10km in the Y direction, 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 ficticous boreholes a and b. Gradient/temperature - depth plots will be examined at each of these locations at a later stage. Surface heat flow for this topography for a mildly advective (k = 1 x 10" 1 6 m2 ) case is shown on Figure 10a. Heat flows determined anywhere within this basin would indicate values close to the basal heat flux. Slightly warped isolines indicate the effects of groundwater flows in both the X-Z and Y-Z directions. These effects are more pronounced on Figures: 10b (k = 3 x 10" 1 6 m2 ) and 10c, a severely disturbed case (k = 5 x 10" 1 6 m2 ). In this case only 55% of the basin is within the minimumly disturbed area, a reduction of 45% from the k = 1 x 10" 1 6 m2 case in just one half an order of magnitude increase in permeablility. Because permeability can vary by thirteen orders of magnitude, the transition from a conductive to an advective-dominated system is 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 yield substantially higher values than SIMPLE WATER TABLE 40 Metre Relief 64 SURFACE HEAT FLOW Simple Water Table 40 Metre Relief k = 1 x 10" 16m2 X-DISTANCE (km) k - 5 x 10~16m2 X-DISTANCE (km) Contours in mWm"2 Figure 10. Surface heat flow - simple water table - 40m r e l i e f in the Y direction. 65 the imposed basal heat f lux , i n d i c a t i v e of groundwater disturbance and not a l o c a l l y higher basal heat flux at depth. Note a l so the lack of c o r r e l a t i o n between the minimumly dis turbed area and the contours of the water t a b l e . If one traced the 4600 m contour l i n e around the side of the h i l l one would encounter heat flows varying from 55 to 75 mWnr2 (Figure 10c). Thus, a d r i l l hole pos i t ion within the minimumly disurbed area in one section of the basin cannot be relocated within the minimumly disturbed area by simply fol lowing a h i l l s l o p e contour. This i s a not iceable e f fect for th i s s l i g h t slope in the Y d i r e c t i o n and as w i l l be shown l a t e r , the ef fect i s much more pronounced when the t o t a l r e l i e f in Y is increased. We w i l l now examine selected sections of surface heat flow (Figure 10c). A sect ion at the edge of the basin (Y = 0 or 20 km ) i s shown on Figures 11a, b, and c , and shows the e f fec ts of downward groundwater flow in the recharge area and upward flows and steeper gradients at the discharge area. D r i l l hole a (Figure 12) r e f l e c t s th i s upward flow in the convex temperature/depth p l o t , i n d i c a t i v e of upward groundwater flow. For comparitive purposes, a lso p lo t ted as the dashed l i n e i s the geothermal gradient for the conductive case. As can be seen, only gradients and thermal c o n d u c t i v i t i e s measured over 2.5 km in depth would y i e l d heat flow values close to the basal heat f lux . Figure 13 i l l u s t r a t e s the interdependance between the locat ion of the sect ion through the bas in , root mean square deviat ions from the basal heat f lux and permeabi l i ty . Figure 13 a lso introduces an important d i f ference between two and three 66 S I M P L E W A T E R T A B L E 4 0 M e t r e R e l i e f k - 5 x 10 " 1 6 m 2 Section at Y = 10 km Section at Y = 5 km Section at Y = 0 km 20 25 30 35 40 DISTANCE ( km ) Isolines in "Centigrade w 200 i DISTANCE (km) Figure 11. Section plots through basin for the Y equals 40 m r e l i e f case. 67 O 1 I 1 | | 1 1 0 1 0 2 0 3 0 4 0 5 0 6 0 GEOTHERMAL GRADIENT ( Â°Ckm- 1 ) Coordinate Position ( a) X = 37.5 km Y = 2.5 km k = 5 x 10~ 1 6m 2 Figure 12. Temperature/gradient - depth p l o t s f o r coordinate p o s i t i o n (.a) . 68 TEMPERATURE ( Â° C ) 1 0 3 0 5 0 7 0 9 0 1 1 0 1 3 0 1 5 0 o 1 1 1 1 1 1 1 0 1 0 2 0 3 0 4 0 5 0 6 0 GEOTHERMAL GRADIENT ( Â°Ckrrr 1 ) Coordinate Position (b) X - 20 km Y - 10 km k - 5 x 10' 1 6 m 2 Figure 13. Temperature/gradient - depth p l o t s f o r coordinate p o s i t i o n (.b) . 69 d i m e n s i o n a l r e s u l t s ; namely t h a t the a d v e c t i v e t h r e s h o l d i s not o n l y a f u n c t i o n of p e r m e a b i l t y but l o c a t i o n w i t h i n the b a s i n . The a d v e c t i v e t h r e s h o l d has been s h i f t e d t o a lower p e r m e a b i l t y v a l u e i n a s e c t i o n near the edge of the b a s i n , i n comparison t o a s e c t i o n near the r i d g e l i n e . Sharp i n c r e a s e s i n the r o o t mean square d e v i a t i o n a r e n o t e d f o r t h e s e r i d g e l i n e s , which a r e l o c a t e d i n a r e a s 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 a r e a s f o r the system. A heat f l o w d e t e r m i n a t i o n on a s e c t i o n a t Y = 5.0 or 15.0 km ( F i g u r e 11d) between 16 and 32km would y i e l d v a l u e s of the c o r r e c t b a s a l heat f l u x . At t h i s s e c t i o n most of the groundwater f l o w 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 ) . A l t h o u g h t h i s s e c t i o n has v a l u e s c l o s e s t t o the c o r r e c t b a s a l heat f l u x , s i g n i f i c a n t 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 re c h a r g e (0 - 16km) and 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 at Y = 7.5 or 15.0 km has the lowe s t r o o t mean square d e v i a t i o n s ( F i g u r e 14) because of r e d u c t i o n s i n d e v i a t i o n s from the b a s a l heat f l u x a t the d i s c h a r g e and r e c h a r g e zones of the h y d r o l o g i c system. The l a t e r a l v a r i a t i o n i n s u r f a c e heat f l o w i s o f t e n used as an i n d i c a t o r of an a c t i v e groundwater f l o w system. However, the r i d g e - l i n e of the system a t 10km, which i s a l i n e of symmetry, does not r e p r e s e n t a good l o c a t i o n f o r measurements even though t h e r e i s l i t t l e l a t e r a l v a r i a t i o n i n heat f l o w between 10 and 35 km. Most of the groundwater f l o w s a r e downward throughout most of the s e c t i o n . Hole b ( F i g u r e 13) at X= 20 km, i l l u s t r a t e s t h i s e f f e c t . Here, a n e a r - l i n e a r t e m p e r a t u r e / d e p t h p l o t i s i n d i c a t e d . S l i g h t downward f l o w would be n o t i c e d from the t emperature p l o t of a deep d r i l l h o l e , but i t i s d o u b t f u l t h a t t h i s t r e n d would be n o t i c e d from a s h a l l o w bore h o l e . 7 0 CM E 6 0 H 5 0 H 4 0 H E 3 0 cr 20 H 10 H 10 r17 SIMPLE WATER TABLE 40 Metre Relief 7 . 5 k m 1 0 k m â€” 1 0 % q h ~i 1 1â€”iâ€”iâ€”r ~i 1 1 1 1â€”r-10 r16 10 .-15 P E R M E A B I L I T Y ( m 2 ) Figure 14. Root mean square deviation versus permeability for the simple water table configuration. 40 m relief in the Y direction. o 71 Near surface heat flow values of 50 mWm"2 would be indicated, differing 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 direction 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 ,3 x 10" 1 6 and 5 x 10" 1 6 m2. Comparing Figures 16c and 10c we see that with just a 40km increase in relief in the Y direction, the size of the minimumly disturbed area has been reduced to 31% from 55% of the total area, and the root mean square for the entire surface has increased from Â±30.6 to Â±35.6 mWm"2. Root mean square deviations for each section line are shown on Figure 17. Again, the 7.5 or 12.5 km section line has the lowest root mean square values; maximum values on the 0 or 20 km line are Â±66 mWm"2 compared to Â±43 mWm" 2 for the 40 m relief case (Figure 14). The shifting of the advective threshold on each line is noted as well. Thus, model results demonstrate a high sensitivity of surface heat flow to variations in relief 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 the same overall slopes, are examined. Another topographic surface is shown on Figure 18; a "flat-top" topography. An overall relief of 40 m in 10km is maintained, but this now occurs over 5km from Y = 0 to 5km and 15 to 20 km, with a flat surface from 5 to 15km. The temperature f i e l d and surface heat flow for an SIMPLE WATER TABLE 80 Metre Relief X-DISTANCE (km) r 30 35 40 Contours in metres Figure 15. Simple water table - 80 m r e l i e f in the Y direction. 73 SURFACE HEAT FLOW Simple Water Table 80 Metre Relief 1 x 10"16m2 15 20 X-DISTANCE (km) 25 -r 30 35 40 3 x 10"16m2 10 15 20 X-DISTANCE (km) - i â€” 25 30 35 40 k - 5 x 10",om 16^2 Figure 16. Surface heat flow -T 15 20 25 X-DISTANCE (km) simple water table 80 m r e l i e f . 30 35 40 Contours in mWm-2 70 P E R M E A B I L I T Y ( m 2 ) Figure 17. Root mean square deviation versus permeability for the simple water table configuration - 80 m rel i e f in the Y direction. Figure 18. Conceptual ill 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" 1 6 m2 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 slightly, the area is 51% of the total basin with an overall root mean square deviation of Â±32.3 mWm"2, compared to 55% and Â± 30.6 mWm"2 for the simple water table case. These results show l i t t l e difference for the two water table configurations, indicating a lack of sensitivity 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 "flat-top" topography. A substantial difference is 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 effects of large scale relief in the direction orthogonal to any section line. Small scale effects were examined further by introducing a small local rise or "hummock" onto the simple water table topography (Figure 21). The hummock is 25 m high with a diameter of about 1 km. In this case an 80 m relief in the Y - direction was imposed to further enhance any potential three dimensional effects. The surface heat flow for this case is shown on Figure 22. Some change to the 40 mWm"2 contour is evident, but SURFACE HEAT FLOW "Flat Top" Water Table 40 Metre Relief k = 5 x 10 - 1 6 m 2 20 ~ 15 E UJ o C0 Q I >-10H 5H X-DISTANCE (km) Figure l'J. Surface heat flow - " f l a t - top" topography. ~1 r 30 35 40 Contours in mWm"2 CO SURFACE HEAT FLOW Tlat-Top" Water Table 40 Metre Relief Section at Y = 1 0 km k - 5 x 1 ( f 1 6 m 2 15 20 25 DISTANCE (km) 30 35 40 Figure 20. Comparison of two and three dimensional model r e s u l t s on a section l i n e at Y equals 10 km â€” " f l a t - top" topography. cc 79 Figure 21. Conceptual ill u s t r a t i o n - Local rise on water table. 80 m r e l i e f . SURFACE HEAT FLOW Local Rise in Water Table 80 Metre Relief k = 5 x 10" 1 6m 2 X-DISTANCE (km) r 30 35 40 Contours In mWm"2 Figure 22. Surface heat flow - Local r i s e on water table. 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 is 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 mWnr2, respectively. These results indicate l i t t l e change between the two cases. It is expected that the introduction of the local rise in other locations would not produce any significant changes to the surface heat flow. The overall relief in both the X and Y directions overshadows any potential effects. A series of humps over the entire surface would introduce interfering local groundwater flows that may cause small scale pertubations on the surface heat flow. However, these local groundwater flows would not penetrate deeply enough to redistribute warmer waters, which is the controlling influence on pertubations to the heat flow regime. Hence, the thermal regime of a basin and therefore the surface heat flow, is not sensitive to small scale variations in the water table. Series 3 z. Effects of Permeable Units At Depth So far only basins with homogeneous permeability have been modeled. However, i t is much more common for sedimentary basins to exhibit strong layering effects, and one or more units may behave as a preferential 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 significant effect, but only if the surrounding permeability is 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 is not well known at this point are how aquifers in a three dimensional setting affect the thermal regime. For instance, is 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 will 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" 1 4 m2. 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 direction. In this simulation the permeability of the unit surrounding the aquifer is 1 x 10" 1 6 m2. For this water table configuration, a homogeneous permeability of 1 x 10" 1 6 m2 produces a mild distubance. However, with the intoduction of the aquifer (Figure 24a) the root mean square deviation for the entire basin increases to Â± 12.8 mWm"2, which is almost three times greater than the homogeneous case ( Â± 4.6 mWm"2 ) and about the same as the slightly disturbed homogeneous case ( k = 3 x 10" 1 6 m2 : Â± 14.5 mWm"2 ). Heat flow and water table contour patterns are similar. Two other effects are noted. Â«3 Figure 23. (a) Extensive aquifer - 40 m r e l i e f on the simple water table. Figure 23. (Continued) (b) P a r t i a l aquifer case. 40 m r e l i e f on the simple water table. 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 10 " 1 6 m 2 k a - 1 x 1<r14m2 Extensive Aquifer 40 m relief 30 35 40 X-DISTANCE (km) Extensive Aquifer 80 m relief X-DISTANCE (km) Partial Aquifer 80 m relief 20 25 X-DISTANCE (km) 30 35 40 Contours in mWm - 2 Figure 24. Surface heat flow for various aquifer cases. 87 Section lines taken through the basin (Figure 25) and compared to the k = 3 x 10" 1 6 m2 case (Figure 26) show that the three dimensional effect has been damped out. Root mean square deviations for individual section lines are similar, Â± 14.0, 12.5, and 11.5 mWnr2 for the Y =(0, 20 km) and ( 5, 15 km) and 10km lines, respectively, compared with Â± 17.9, 13.3, and 11.6 mWnr2 for the homogeneous case ( k = 3 x 10 - 1 6 m2). Furthermore, in the homogeneous basin, portions of the heat flow profiles on some section lines 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 line. The intersection of the heat flow curve with the 60 mW/m2 isoline is shifted slightly on each section line. These observations can be explained by examining the flow system. Because of their higher permeability, aquifers have the effect of substantially altering the hydrologic system so that most of the flow is channeled to the discharge end of the system. The preceeding simulation was repeated with 80 m relief in the Y direction so that the effects of aquifers could be examined for sensitivity to changes in slope in the Y direction. The resulting surface heat flow (Figure 24b) has a root mean square deviation of Â± 13.0 mWm'2 , only slightly higher than the 40 m case ( Â± 12.8 mWm-2 ). This further confirms that aquifers, while enhancing the thermal disturbance, tend to damp out any three dimensional effect. In some instances facies changes or faulting may cause the development of a partial 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 25 D I S T A N C E (km) Figure 25. Section lines through heat flow for the extensive aquifer case - 40 m r e l i e f in Y. cc CO SIMPLE WATER TABLE LL H â€” i 1 1 1 1 " ' ' 0 5 10 15 20 25 30 35 40 DISTANCE (km ) Figure 26. Surface heat flow along three section lines. 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 ridge-line. The resulting heat flow is 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 m2 - Figure 16a), i t can be seen that while the half section of the basin with the aquifer is identical 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 illustrates this observation by comparing a section line at 12.5km for the homogeneous case (root mean square deviation Â± 4.27 mWm"2) with a 12.5km line for the partial aquifer case (root mean square deviation Â± 6.81 mWm"2). This effect can be explained by examining the flow system in the Y direction, 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 vertical plane through Y = 10 km. With the introduction of the partial aquifer this effect is eliminated and the groundwater divide shifts towards higher Y coordinate values and may no longer be vertical. 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 is of interest to heat flow investigators is the case of discontinuous aquifers o _J LL H < LU I LU O < Li. DC CO CM I 120n 80 A 40 T 5 SIMPLE WATER TABLE 80 Metre Relief Section at Y - 12.5 km 10 1 2 Homogeneous k = 1x10 1 6 m 2 Partial aquifer k = 1 x 10~16m2 ka = 1 x 10~14m2 15 â€”r~ 20 25 DISTANCE (km) 30 35 1 2 40 Figure 27. Comparison between surface heat flow for a homogeneous case and a p a r t i a l aquifer case. 92 (statigraphic pinchouts). Figure 23c shows this case conceptually. This problem is of interest because i t is known that the pinching out of aquifers tends to direct 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 alter the heat flow regime at these locations? To examine this 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" 1 6 m2 and the aquifer to 1 x 10" 1 4 m2 (Case #1) and to 5 x 10" 1 4 m2 (Case #2). These permeabilities were chosen to be compatible with the extensive aquifer case. The simple water table configuration has a relief of 80 m in the Y direction. Results are shown on Figure 28a and b. L i t t l e difference is noted for these two simulations. Root mean square deviations are Â± 7.5 and Â± 9.9 mW/m2. Figure 29 shows surface heat flow on three section lines for Case #2: (0, 20), (5, 15) and 10 km. A characteristic curve is 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 is Â± 5.2 mW/m2 in comparison for the homogeneous case - k = 1 x 10" 1 6 m2 ). The effects of increasing the background permeability were also examined. In this case the background permeability was increased to 5 x 10" 1 6 m2 and the result is shown on Figure 28c. At this permeability the thermal effects of the stratigraphic pinchout are masked by the overall 93 SURFACE HEAT FLOW Simple Water Table Stratigraphic Pinchout 80 Metre Relief -Aquifer k a - 1 x 1 0 " 1 4 m 2 1 x 10" 1 6 m 2 15 20 25 X-DISTANCE (km) 1 x 1 0 _ 1 4 m 2 5 x 1 0 " 1 6 m 2 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 in the Y direction. SIMPLE WATER TABLE Stratigraphic Pinchout 5 o _ J LU = X 5 LU E O w < LL CC CO 120 80 40 k - 1 x 10- 1 6 m 2 k a = 5 x 1 0 _ 1 4 m 2 ~r 5 I 10 i 1 r~ 15 20 25 DISTANCE (km) 30 35 0 km 5 km 10 km 40 Figure 29. Surface heat flow along three section lines for the pinchout number 2 case. 4 ^ 95 disturbance caused by higher groundwater velocities. A root mean square deviation of Â± 38.5 mW/m2 was calculated for this case which is similar to the homogeneous case (k = 5 x 10" 1 6 m2 and Â± 35.6 mW/m2). Based on model results i t can be concluded that there is 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 values of less than 1 x 10" 1 6 m2, the system would behave as a conductive case. At higher background values of k = 5 x 10" 1 6 m2, the system behaves as an advectively disturbed homogeneous case. An increase in aquifer permeability (for a given background permeability) will increase the disturbance to some maximum limit. Only a further increase in background permeability w i l l 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 sensitivity of surface heat flow to changes in relief in the Y direction and also to small scale pertubations of the water table surface. In this section we wi l l examine the sensitivity of surface heat flow to large scale changes in shape and relief of the water table in the X direction. To do this, two additional 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 third, and 50 m in the last 10km. The relief in the Y direction is 40 m. Surface heat flow maps are presented (Figure 32) at three permeability levels: 1, 3, and 5 x 10" 1 6 m2. Some interesting effects are noted when one compares this case to the simple water table case (Figure 10c) at k= 5 x 10" 1 6 m2. Firs t , because of the broad, low sloping discharge region the higher heat flow values decrease in magnitude but are spread out over a larger area. The total 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"2, whereas the maximum values for the concave water table configuration are about 90 mWm"2. The root mean square deviation for the basin is reduced to Â± 19.9 mWm"2 compared to Â± 30.6 mWm"2 for the prior case. Second, in contrast 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 effects are a direct result of changes in slope of the water table. With the concave slope, groundwater velocities at depth are lower at the discharge end of the system and thus, the overall disturbance is reduced. Figure 33, a plot of root mean square deviation versus permeability on each section line, shows how the advective threshold has shifted to a higher permeability value with the 97 Figure 30. Conceptual il 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 direction. CONCAVE WATER TABLE X-DISTANCE (km) Contours in metres Figure 31. Concave water table topography. CO 99 SURFACE HEAT FLOW Concave Water Table X-DISTANCE (km) k - 3 x 10"16m2 X-DISTANCE (km) k = 5 x 10"16m2 X-DISTANCE (km) Contours in mWm Figure 32. Surface heat flow - Concave water table case. CONCAVE WATER TABLE 40 Metre Relief PERMEABILITY Figure 33. Root mean square deviation versus permeability for the concave water table. 101 change in the topographic configuration of the water table. Therefore, knowledge of only"the overall change in elevation of the slope may be inadequate to determine the level 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 is a water table replicating 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 direction of 17.2 m/km. The gradient of the slope connecting the two segments is 51.6 m/km. The total relief is 1000 m in 40 km in the X direction. Smith and Chapman (1983) analyzed this configuration in detail. In three simulations, permeabilities were fixed at 5 x 10" 1 6 m2 and relief in Y was set at 40 m, 80 m, and a final case of 80 m in the upland/slope regions and only 8 m in the lowland (Figures 34 a, b, c and 35). The f i r s t two simulations were carried out to test the sensitivity of surface heat flow to changes in slope for this type of topography. The final simulation was carried out to understand the three dimensional thermal effects in a more rea l i s t i c setting. In typical Basin and Range topography the slope segments are usually a result of fault uplifting, with side relief 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"2 to highs of over 140 mWm"2. These plots are characterized by the large gradients of heat flow at the steepest part of the 102 Figure 34. Conceptual illustrations of the segmented water table configuration, (.a) 80 m r e l i e f . 103 Figure 34. (Continued) (.b) 40. m r e l i e f . 104 Figure 34. (Continued) (c) 80. m r e l i e f upland -8m r e l i e f lowland. 105 SEGMENTED WATER TABLE 20 40 Metre Relief o m a 5H 4600 4400 4700 / 4500 / 4300 4200 4100 15 20 25 X - D I S T A N C E (km) â€” i -30 4000 35 40 80 Metre Relief 10 15 20 25 X - D I S T A N C E (km) 30 35 40 80 Metre Relief Upland/8 Metre Relief Lowland 15 20 X - D I S T A N C E (km) T 25 30 35 40 Contours in metres Figure 35. Segmented water table topography. SURFACE HEAT FLOW Segmented Water Table k - 5 x 10" 1 6m 2 40 Metre Relief X-DISTANCE (km) Contours in mWm Figure 3 b . Surface heat flow for the segmented water table. 1 07 topography. Figure 37 illustrates this observation by presenting surface heat flow on three section lines for the 40 m relief case. Figure 36a also shows the effects of local 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 sensitivity to changes in slope in the Y direction similar to simple water table cases. Figure 36c shows the surface heat flow resulting from the third topographic case. Now, a heat flow low does not develop in the lowland as the decreased slope in the Y direction damps out the three dimensional effect. In this 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 local high of 90 mWm"2 develops at the Y = 0 or 20 km section line in the lowland region ( X = 26 km position). Root mean square values on each section line are: Â± 32.3, 30.0, and 30.2 mWm2 for the Y equals (0, 20), (5, 15) and 10 km lines. 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 level of disturbance in this type of basin is more sensitive to changes in Y relief in the lowland section. In discharge regions, higher temperatures exist because of upward groundwater flows. Coupling between temperature and fluid properties has the effect of further increasing groundwater velocities. Hence, any SEGMENTED WATER TABLE 0 5 10 1 5 2 0 2 5 3 0 3 5 4 0 DISTANCE (km) o cc Figure 37. Surface heat flow along three section lines in 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 is homogeneous but anisotropic (Freeze and Cherry, 1979). This layering is common in sedimentary sequences. Smith and Chapman (1983) investigated the thermal effects of anisotropic formations and found this had the effect of shifting the advective threshold to a higher permeability value by reducing vertical velocities. In order to investigate the sensitivity of anisotropy in permeability on flow in a basin, the simple water table configuration was simulated with values of kx = ky = 5 x 10' 1 6 m2 and kz = 5 x 10 - 1 7 m2. An 80 m relief in the Y direction was chosen to accentuate any three dimensional effect. The result is shown on Figure 38. The overall level of disturbance is reduced from the isotropic case( root mean square of Â± 9.0 mWm"2 compared to Â± 35.6 mWm-2 ) and is similar to the aquifer case (Â± 13.0 mWm~2 ). At this level 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 mWnr2 , respectively. The explanation of this effect is SURFACE HEAT FLOW Anisotropic Permeability Simple Water Table 80 Metre Relief X-DISTANCE (km) Contours in mWm F i g u r e 38. S u r f a c e 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 c a s e . S i m p l e w a t e r t a b l e c o n f i g u r a t i o n . 111 that the anisotropic case is 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 z as the present case being considered. Like the aquifer, flows tend to be channelled more in the X rather than the Y direction. Therefore, anisotropic cases will cause advective thresholds to be shifted to higher permeabilty values and three dimensional effects will be subdued. Series 6 - Effect of Depth of the Basin The sensitivity 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 fini t e element mesh to permeabilities of 1 x 10" 2 1 m2. In the rest of the system permeabilities were set to 5 x 10" 1 6 m2 , and the results are shown on Figure 39. The overall level of disturbance was reduced from the 5 km depth case (root mean square deviation of Â± 9.0 mWm"2 as opposed to Â± 35.6 mWm"2 ). This shift to higher permeability values is caused by the elimination of flow paths to the deeper, warmer areas of the basin. Consequently, permeabilities and therefore velocities would have to be increased to achieve the same level of disturbance. It is noted however, that the three dimensional effect has not been damped. SURFACE HEAT FLOW 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 will be to discuss some of the implications of the sensitivity study. Suggestions for future research in coupled groundwater/heat flow will 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 characteristic lengths and depths of flow systems, and levels of heat flow disturbance. However, the results of Sections 4.4 suggest that the calculation of a single number to describe the disturbance may not be possible. A complete description of the water table configuration may be required before thermal disturbances may be estimated. To examine the idea of computing a single number to represent heat flow disturbance more closely, we will examine recharge rates and levels of disturbance for different water table configurations. The basin recharge as computed by the model, is defined as the sum of a l l fluid fluxes into of the individual surface elements multiplied by individual element areas. It therefore represents the total quantity of flow into the entire basin that is required to maintain the elevation of the water table. Table 5 is a summary of basin recharge rates and levels 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"2 ; inflow 2.3 x 10~2 m3/sec) 114 TABLE 5. Groundwater recharge rates to hydrologic system with root mean square deviations from basal heat flux. N Simulation Recharge RMS deviation (m3/sec) (mWm"2) 1 S.W.T. 40m k=1x10- 2 1 4.3X10"8 0.8 2 S.W.T. 40m k=1x10" 1 7 3.8x10-" 1 .1 3 S.W.T. 40m k = 3x10" 1 7 1 .2X10" 3 1 .8 4 S.W.T. 40m k=1x10" 1 6 4.4x10"3 4.6 5 S.W.T. 40m k=3x10" 1 6 1.3X10-2 14.5 6 S.W.T. 40m k=5x10" 1 6 2.4x10"2 30.6 7 S.W.T. 80m k=3x10" 1 7 1.6X10"3 2.2 8 S.W.T. 80m k=1x10- 1 6 5.3X10" 3 5.2 9 S.W.T. 80m k=3x10" 1 6 1.7x10-2 16.3 10 S.W.T. 80m k=5x10- 1 6 3 . 1 X 1 0 - 2 35.6 1 1 C.W.T. 40m k=3x10" 1 7 1.3x10-3 1 .9 12 C.W.T. 40m k=1x10" 1 6 4.3X10"3 4.6 13 C.W.T. 40m k=3x10- 1 e 1.3x10-2 12.1 1 4 C.W.T. 40m k=5x10- 1 6 2.3x10-2 19.9 15 B & R 80m k=1x10" 2 1 9.4X10" 8 1 .3 16 B & R 40m k=5x10" 1 6 4.3X10-2 32.5 17 B & R 80m k=5x10" 1 6 4.7X10"2 36. 1 18 B & R 80/8 k=5x10" 1 6 4.4X10-2 30.5 19 ANIST 80m k=5x10" 1 6 1.0X10-2 9.0 20 FLAT 40m k=5x10- 1 6 2.6xl0'2 32.3 21 REDUC 80m k=5x10" 1 6 1.9X10"2 10.0 22 HUMM 80m k=5x10- 1 6 3.4x10-2 35.9 23 AQUIF 40m k=5x10" 1 6 1.8x10-2 12.8 24 AQUIF 80m k=5x10" 1 6 1.9x10-2 13.0 25 PARTA 80m k=1x10" 1 6 1.3x10-2 9.9 26 PINC1 80m k=1x10" 1 6 1.1x10-2 7.5 27 PINC2 80m k=1x10- 1 6 1.7x10-2 9.9 28 PINC3 80m k=5x10" 1 6 4.1x10-2 38.5 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" - local rise in water table AQUIF = Extensive aquifer PARTA = Partial aquifer PINCn = Discontinuous aquifer case n=1,2 or 3 RMS = Root mean square 115 with the simple water table case (root mean square Â± 30.6 mWnr2 ;inflow 2.4xl0~2 m3/sec), we see almost identical recharge rates to the system, but levels of disturbance that differ by 35%. The simulations of the simple and concave water table slopes showed that for a given relief 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 logrithmically related. Figures 40a and b show this relationship for two different topographic cases. The significance of this observation is that i t may be possible to determine i f a basin is hydrologically disturbed (within broad ranges) by estimating recharge and empirically relating 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 existing hydrologic data or conduct f i e l d permeability tests. 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 is 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 carried out sucessfully (for example, Freeze and Witherspoon, 1967, Jamieson and Freeze, 1983). Therefore, the estimation of recharge rates provides an RECHARGE (x 10"2) m3/sec QTT RECHARGE (x 10 - 2) m3/sec o in 118 alternate method for determining basin thermal disturbances. An empirical relation between different water table configurations in basins of differing 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 iterative 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. Finally, what is evident from the results is that a two step strategy should be employed in any fi e l d program of heat flow measurements. The f i r s t step (as was discussed in the preceeding paragraphs) is to estimate whether or not a thermal regime may be hydrologically disturbed. If i t i s , the second stage would be to employ fi e l d methodology to deal with the problem. Results from the sensitivity study suggest a general approach for locating boreholes for thermal measurements in 119 hydrologically-disturbed environments. As we have seen, a minimumly disturbed area is present somewhere within a basin. Continued modelling of different topographic shapes may yield information about where minimumly disturbed areas may be located within general water table types (i.e. concave, convex, simple, etc.). This classification may help investigators to quickly determine where highly disturbed areas are located within a basin. A number of shortly spaced boreholes may then be needed in several directions to correctly interpret heat flow measurements. For instance, a series of closely spaced holes in the X direction for the simple water table case (Figure 16c) would show l i t t l e lateral variation in heat flow from 15 to 30 km. A series of closely spaced holes in the Y direction would indicate a suface heat flow.pattern like Figure 41a, and shows how complicated the heat flow pattern really is because of three dimensional effects. These figures show surface heat flow along various sections in the Y direction for the concave and simple water table cases at k = 5 x 10" 1 6 m2. Interestingly, Figure 41a and b show how the correct basal heat flux appears to be located at the main inflextion points on these curves. At these locations the basal heat flux is 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 is vertically conductive and surface heat flow values measured at these points would give values close to the basal heat flux. A suggested area of research is to therefore integrate the fields SURFACE HEAT FLOW Simple Water Table 80 Metre Relief k Â« 5 x 10- 1 6m 2 Section at X * 10 km Section at X s 5 20 km Section at X - 30 km 100 -i X - N 80 Basal Heat Flux 10 0 10 0 Y - C O O R D I N A T E ( k m ) Heat Flow Approx. Straight Line Segment 10 Figure 41. (a) Surface heat flow along three sections in the X direction for the simple water 100 SURFACE HEAT FLOW Concave Water Table 40 Metre Relief k - 5 x 1 0 " 1 6 m 2 Section at X - 30 km Section at X - 21.8 km CM E <: o U_ 40 '< UJ I Approx. Straight Line Segment Y-COORDINATE (km) Figure 41. (.Continued) (b) Surface heat flow along two sections i n the X d i r e c t i o n f o r the concave water table. 122 of s u r f a c e water h y d r o l o g y , groundwater h y d r o l o g y and heat f l o w t o i n v e s t i g a t e the use of mapping r e c h a r g e / d i s c h a r g e a r e a s . T h i s approach c o u l d be a way of d e t e r m i n i n g h o r i z o n t a l f l o w l o c a t i o n s , and hence, 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 moving groundwater. C o n v e r s e l y , i f h y d r o l o g i c d a t a was l a c k i n g , g r a d i e n t s of heat f l o w i n the s p a t i a l (X,Y) c o o r d i n a t e s might be used t o d e t e r m i n e the c o r r e c t b a s a l heat f l u x . 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 out t o : (1) d e v e l o p a s e r i e s of b a s i n type c u r v e s r e l a t i n g t o p o g r a p h i c shape, d i m e n s i o n l e s s l e n g t h , depth and w i d t h t o d i m e n s i o n l e s s r e c h a r g e , and r o o t mean square d e v i a t i o n from the b a s a l heat f l u x ; (2) d e v e l o p f i e l d 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 t h e r m a l 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 ; and, (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 , groundwater h y d r o l o g y , and heat f l o w t o d e v e l o p methods f o r mapping r e g i o n a l r e c h a r g e / d i s c h a r g e a r e a s and r e l a t i n g t h e s e t o heat f l o w measurements. 123 CHAPTER 5 SUMMARY AND CONCLUSIONS "Any solution is 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 effects 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 collected on cross sections and flow systems are seldom two dimensional. (2) Performing a generic study of hypothetical three dimensional basins to determine the sensitivity 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 fl u i d and heat transfer equations and allow for temperature-dependant fluid 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, non-deforming porous media. Fluid 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 verified 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 this task was to carry out a sensitivity 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, partial and discontinuous aquifer cases, anisotropic formations, and reduced depths of active flow. The results of the sensitivity study are summerized below: 125 Variations in Basin Permeability. 1. A simple sloping topography in three dimensions was modeled at three levels of permeability to produce conductive to advectively disturbed thermal regimes. The transition from a conductive to an advective-dominated system was found to be sharp, confirming the results of two dimensional models. Also, surface heat flow reflects the spatial distribution of recharging and discharging fl u i d associated with the three dimensional groundwater flow system. 2. 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 locally 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 hillslope contour. In areas of the basin that were advectively disturbed, only gradients and thermal conductivities measured over 2.5 km in depth would yield heat flow values close to the basal heat flux. 3. Important differences exist between two and three dimensional results; namely that the advective threshold is 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 ridge-line. 4. The lateral variation in surface heat flow is often used as 126 an indicator of an active groundwater flow system. However, the ridge-line of the simple water table system, which is a line of symmetry, did not represent a good location for measurements even though there is l i t t l e lateral variation in heat flow. Variations in Water Table Topography 1. For the simple water table configuration, there is a high sensitivity of surface heat flow to variations in relief in the Y direction. However, the introduction of small scale variations (for a given relief in X and Y), such as a local rise or a "flat-top" type of surface, did not produce significantly different surface heat flow than the model without the variation. Therefore, the surface heat flow is insensitive to these small scale variations. 2. Important differences exist between two and three dimensional models and between small and large scale effects. A significant difference was noted between two and three dimensional simulations of the "flat-top" water table configuration. The results underscored the need to recognize the effects of large scale relief in the direction orthogonal to any section line. 3. Modeling a concave shaped water table with the same relief as the simple water table cases showed that knowledge of only the overall change in elevation of the slope is inadequate to determine the level of disturbance in heat flow. The precise geometry of the water table may be required to properly understand the thermal regime of a basin in an advectively disturbed case. 127 4. Another flow domain that was considered was a segmented water table replicating subdued Basin and Range topography. Three cases were considered; two to investigate sensitivity of the thermal regime to varying relief in the X and Y directions, and a final simulation was carried out to understand the three dimensional thermal effects in a more reali s t i c setting. It is possible under certain topographic conditions to have local 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 sensitivity to changes in slope in the Y direction to the simple water table configuration. 5. The segmented water table simulations also showed that the level of disturbance in this type of basin is more sensitive to changes in Y relief in the lowland section. In discharge regions, higher temperatures exist because of upward groundwater flows. Coupling between temperature and fluid properties has the effect of further increasing groundwater velocities. Hence, any changes to driving forces (increased topographic slopes) that effect these locations has a major impact on the surface heat flow. Effects of Introducing Permeable Units at Depth 1. Because of their higher permeability, extensive aquifers have the effect of substantially altering the hydrologic system so that most of the flow is channeled to the discharge end of the system. Extensive aquifers, while enhancing the thermal disturbance, tend to damp out any three dimensional effect. 2. For a partial aquifer terminating along the ridge-line of a 128 basin, the half section of the basin with the aquifer had identical 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 characteristic surface heat flow curve is developed that appears to be a combination of the aquifer and homogeneous cases. It is 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 is only a small "window" of contrasting permeability values in which the pinchout effect would be noticed for a given aquifer depth and thickness. Anisotropic Permeabilities 1. A basin with a simple water table configuration was modeled with the X and Y direction permeabilities ten times greater than in the Z direction. This anisotropic permeability caused the advective threshold to be shifted to a higher permeabilty value and three dimensional effects to be subdued. Reduced Depth of Groundwater Flow 1. The sensitivity 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, t h a t t h r e e d i m e n s i o n a l e f f e c t s were n o t s u b d u e d . The f i n a l g o a l o f t h i s r e s e a r c h was t o d i s c u s s some o f t h e i m p l i c a t i o n s o f t h e s e n s i t i v i t y s t u d y t o h e a t f l o w m e a s u r e m e n t s . T h e s e f i n d i n g s a r e d i s c u s s e d b e l o w : 1. S i m u l a t i o n s showed t h a t t h e c a l c u l a t i o n o f a s i n g l e d i m e n s i o n l e s s number t o d e s c r i b e t h e l e v e l o f t h e r m a l d i s t u r b a n c e may n o t be p o s s i b l e . A c o m p l e t e d e s c r i p t i o n o f t h e w a t e r t a b l e c o n f i g u r a t i o n may be r e q u i r e d b e f o r e t h e r m a l d i s t u r b a n c e s a r e e s t i m a t e d . H o w e v e r , i t may be p o s s i b l e t o e s t i m a t e ( w i t h i n b r o a d r a n g e s ) t h e r m a l d i s t u r b a n c e i n a b a s i n p r i o r t o e x t e n s i v e f i e l d e x p l o r a t i o n . F o r a g i v e n w a t e r t a b l e c o n f i g u r a t i o n , r e c h a r g e r a t e s a n d r o o t mean s q u a r e d e v i a t i o n v a l u e s a p p e a r t o be l o g r i t h m i c a l l y r e l a t e d . The s i g n i f i c a n c e o f t h i s o b s e r v a t i o n i s t h a t i t may be p o s s i b l e t o d e t e r m i n e i f a b a s i n 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 e s t i m a t i n g b a s i n r e c h a r g e a n d e m p i r i c a l l y r e l a t i n g r e c h a r g e t o r o o t mean s q u a r e d e v i a t i o n . 2. What i s e v i d e n t f r o m t h e r e s u l t s i s t h a t a two s t e p s t r a t e g y s h o u l d be e m p l o y e d i n any f i e l d p r o g r a m o f h e a t f l o w m e a s u r e m e n t s . The f i r s t s t e p i s t o e s t i m a t e w h e t h e r o r n o t a t h e r m a l r e g i m e may be h y d r o l o g i c a l l y d i s t u r b e d . I f i t i s , t h e s e c o n d s t a g e w o u l d be t o e m p l o y f i e l d m e t h o d o l o g y t o d e a l w i t h t h e p r o b l e m . R e s u l t s f r o m t h e s e n s i t i v i t y s t u d y s u g g e s t a g e n e r a l a p p r o a c h f o r l o c a t i n g b o r e h o l e s f o r t h e r m a l m e a s u r e m e n t s 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 . A number o f c l o s e l y s p a c e d b o r e h o l e s may be n e e d e d i n s e v e r a l d i r e c t i o n s t o c o r r e c t l y i n t e r p r e t h e a t f l o w m e a s u r e m e n t s . The c o r r e c t b a s a l 130 heat flux is located at one of the main inflextion points of the surface heat flow map. At these locations the thermal regime is vertically 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 spatial (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 site combined with more closely spaced d r i l l holes may be required to correctly interpret near surface heat flow values. It is apparent that Smith and Chapman's conclusions are even more relevant when three dimensional effects are considered. It is hoped that the work begun in this study w i l l encourage more interaction between hydrogeologists and geophysists in 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 alteration 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 , Elsevier, 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. Sci. Instrum., 34, 186-189, 1957. Betcher, R.N., Temperature distribution in deep groundwater flow systems-a fin 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. Sci., 245, 733-753, 1947. Blackwell, D.D., J.L. Steele and CA. Brott, The terrain effect on terrestrial heat flow, J. Geophys. Res., 85(B9), p.4757-4772, 1982. Bourne, D.R., An electric analog simulation of groundwater flow patterns at a potash waste disposal pond located near Esterhazy, Saskatchewan, Unpubl. Msc. Thesis, Univ. British Columbia, 1976. Brott, C.A., D.D. Blackwell, and J.P. Ziagos, Thermal and tectonic implications of heat flow in the Eastern Snake River Plain, J. Geophys. Res., 86, 11709-11734, 1981. Brownell, D.A., S.K. Garg and J.W. Pritchett, Governing equations of geothermal reservoirs, Water Resour. Res., 13, 929-935 , 1977. Bullard, E.C, Heat flow in South Africa, Proc. Royal Soc. London, A173(995), 474-502, 1939. Bredehoeft, J.D., and I.S. Papadopulos, Rates of vertical groundwater movement estimated from the earth's thermal profile, Water Resour. Res., 1(2), 325-328, 1965. Cartwright, K., Redistribution of geothermal heat by a shallow aquifer, Geol. Soc. Am. Bull., 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. Bull., 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., CB. 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, saturated-unsaturated flow in a groundwater basin, Water Resour. Res., 7, 347-366, 1971. Freeze, R.A., Mathematical models of hillslope hydrology, in Hillslope Hydrology , (Chapter 6), edited by M.J. Kirkby, 177-225, 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. Analytical 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 aquifer-aquitard systems, Proc. Third Inter. Conf. on Finite Elements in Water Resources, Univ. Miss.,Oxford, May, 1980. Gambolati, G., Perspective on the modified conjugate method for the fin i t e element solution of linear subsurface equations, Proc. Third Inter. Conf. on Finite 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 British 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 finite element model for simulating hydrothermal reservoirs, Advances in Computer Methods for Partial Differential Equations, IMACS, 284-293, 1977. Hyndman, R.D., Heat flow measurments in the inlets of southwestern British Columbia, J. Geophys. Res., 81, 337-349, 1976. Jamieson, G.R. and R.A. Freeze, Determining hydraulic conductivity distributions 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 collection - 1975, Geothermal Service of Canada, Geotherm. Ser. 5, 1976. Kilty, 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 two-dimensional 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 , 2nd. edition , Freeman, San Fransisco, 1967. 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 British 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. Liu, 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 vertical 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 California Cascades, USGS open f i l e report 82-150, 1982. Mercer, J.W., Pinder, G.F., and I.G. Donaldson, A Galerkin fin 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 circular 790, 1979. Narasimhan, T.N. and P.A. Witherspoon, An integrated fi 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, Isreal, 1974. Nicol, T., UBC SPARSPAK: Solving sparse systems of linear equations, Computiong Centre, Univ. of British Columbia, 1982. Parsons, M.L., Groundwater thermal regime in a glacial complex, Water Resour. Res., 6(6), 1701-1720, 1970. Pickens, J.F. and G.E. Grisak, Finite-element analysis of liquid 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, Finite element simulation in surface and subsurface hydrology, Academic Press, 1977. Prickett, T.A. and C.G. Lonnquist, Comparison between analog and di g i t a l simulation techniques for aquifer evaluation, Proc. Symp. Use Analog Digital Computers Hydrol., Inter. Assoc. Sci. 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., Digital 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. Field experiments and comparisons with theoretical results, Water Resour. Res., 18(2), 253-265, 1982b. Segerlind, L.J., Applied Finite Element Analysis , Wiley, New York, 1976. Segol, G., A three-dimensional galerkin finite 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, California, Summaries Second Workshop Geothermal Reservoir Engineering, Stanford University, Stanford, California, 324-338, 1976. Sorey, M.L., Numerical modeling of liquid 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 velocity, 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, Vol. 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, Victoria, B.C., Canada, May, 1983. van Orstrand, C.E., Temperature gradients, in Problems of Petroleum Geology, 989-1021, Am. Assoc. Petrol. Geol., Tulsa, Okla., 1934. Verge, M-J., A three-dimensional saturated-unsaturated groundwater flow model for practical 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. Sci., 21, 291-300, 1966. Zienkiewicz, O.C., The Finite Element Method , 3rd. ed., McGraw-H 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
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | The thermal effects of three dimensional groundwater flow |
Creator |
Woodbury, Allan David |
Publisher | University of British Columbia |
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 |
AggregatedSourceRepository | 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}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0052537/manifest