UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

The thermal effects of three dimensional groundwater flow Woodbury, Allan David 1983

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

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

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

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0052537/manifest

Comment

Related Items