UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Interaction of groundwater flow systems and thermal regimes in mountainous terrain : a numerical study Forster, Craig Burton 1987

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

Item Metadata

Download

Media
831-UBC_1987_A1 F67.pdf [ 10.31MB ]
Metadata
JSON: 831-1.0052356.json
JSON-LD: 831-1.0052356-ld.json
RDF/XML (Pretty): 831-1.0052356-rdf.xml
RDF/JSON: 831-1.0052356-rdf.json
Turtle: 831-1.0052356-turtle.txt
N-Triples: 831-1.0052356-rdf-ntriples.txt
Original Record: 831-1.0052356-source.json
Full Text
831-1.0052356-fulltext.txt
Citation
831-1.0052356.ris

Full Text

INTERACTION OF GROUNDWATER FLOW SYSTEMS AND THERMAL REGIMES IN MOUNTAINOUS TERRAIN: A NUMERICAL STUDY By CRAIG BURTON FORSTER B . S c , The U n i v e r s i t y of B r i t i s h Columbia, 1975 M . S c , The U n i v e r s i t y of Waterloo, 1979 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY i n THE FACULTY OF GRADUATE STUDIES (Department of G e o l o g i c a l S c i e n c e s ) We a c c e p t t h i s t h e s i s as conforming t o the r e q u i r e d s t a n d a r d THE UNIVERSITY OF BRITISH COLUMBIA October 1987 © C r a i g Burton F o r s t e r , 1987 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of GEOLOGICAL SCIENCES The University of British Columbia 1956 Main Mall Vancouver, Canada V6T 1Y3 Date October 15, 1987 DE-6(3/81) ABSTRACT I t i s w i d e l y r e c o g n i z e d t h a t t o p o g r a p h i c a l l y - d r i v e n groundwater f low can p e r t u r b c o n d u c t i v e thermal reg imes . H i g h - r e l i e f topography a m p l i f i e s the impact of f a c t o r s c o n t r o l l i n g groundwater flow and a d v e c t i v e heat t r a n s f e r . A f i n i t e element method i s deve loped to model the i n f l u e n c e of g e o l o g y , c l i m a t e , s u r f a c e topography and r e g i o n a l heat f l u x on s teady groundwater f low and heat t r a n s f e r . Because f l u i d v i s c o s i t y (hence f l u i d f l u x ) depends upon t e m p e r a t u r e , groundwater f low i s i n f l u e n c e d by the r e g i o n a l heat f l u x . As a consequence , i s o t h e r m a l approaches to model ing deep groundwater flow i n mountains may be i n a p p r o p r i a t e . U s i n g a f r e e - s u r f a c e a p p r o a c h , the water t a b l e i s r e p r e s e n t e d as an i n t e r n a l c h a r a c t e r i s t i c of the groundwater flow system, r a t h e r than the upper boundary f o r f l u i d f l o w . T h i c k u n s a t u r a t e d zones are expec ted i n h i g h - p e r m e a b i l i t y t e r r a i n ( g r e a t e r than 10"15 m 2 ) w i t h a r i d c l i m a t e , or where groundwater recharge i s r e s t r i c t e d by e x t e n s i v e a l p i n e g l a c i e r s . Only v e r t i c a l f l u i d flow i s assumed to occur i n the u n s a t u r a t e d zone, t h e r e f o r e , heat t r a n s f e r above the water t a b l e i s r e p r e s e n t e d by o n e - d i m e n s i o n a l a d v e c t i o n and t w o - d i m e n s i o n a l c o n d u c t i o n . S i m u l a t i o n r e s u l t s i n d i c a t e t h a t water t a b l e e l e v a t i o n s are h i g h l y s e n s i t i v e to changes i n the c o n t r o l l i n g f a c t o r s , but have l i t t l e impact on the thermal reg ime . C o n d u c t i v e thermal regimes are p r e d i c t e d i n l o w - p e r m e a b i l i t y t e r r a i n ( l e s s than 1 0 " 1 8 m 2 ) or in h i g h - p e r m e a b i l i t y t e r r a i n w i t h a r i d c l i m a t e (recharge r a t e s i i l e s s than 1 0 " 1 1 m / s e c ) . S t r o n g a d v e c t i v e heat t r a n s f e r masks the r e g i o n a l heat f l u x when p e r m e a b i l i t y exceeds 1 0 " 1 6 m 2 i n t e r r a i n w i t h r e l i e f of 2 km over a h o r i z o n t a l d i s t a n c e of 6 km. Less than one percent of t y p i c a l mean annual p r e c i p i t a t i o n i s t r a n s m i t t e d through deep groundwater flow systems under these c o n d i t i o n s . Asymmetric s u r f a c e topography c o m p l i c a t e s e f f o r t s to i n t e r p r e t c h e m i c a l and t h e r m a l da ta c o l l e c t e d near the v a l l e y f l o o r . F r a c t u r e zones o u t c r o p p i n g a t the v a l l e y f l o o r can c a p t u r e a l a r g e percentage of groundwater f l o w i n g through the system and a s i g n i f i c a n t percentage of the b a s a l heat f l u x . Maximum s p r i n g temperatures are i n d i c a t e d when bulk p e r m e a b i l i t y i s between 1 0 " 1 7 m 2 and 1 0 " 1 5 m 2 . O u t s i d e t h i s range , s p r i n g temperatures approach ambient a i r t e m p e r a t u r e . T o p o g r a p h i c a l l y - d r i v e n groundwater flow can d i s t o r t and o b l i t e r a t e f r e e - c o n v e c t i o n c e l l s tha t might o therwise d e v e l o p w i t h i n a mountain m a s s i f . TABLE OF CONTENTS Page A b s t r a c t i i T a b l e of Contents i v L i s t of T a b l e s v i i L i s t of F i g u r e s v i i i Acknowledgements x i CHAPTER 1 - INTRODUCTION 1 References 9 F i g u r e s 13 CHAPTER 2 - DEVELOPMENT OF THE NUMERICAL METHOD 14 2.1 I n t r o d u c t i o n 14 2.2 Conceptua l Model 15 2.3 Mathemat i ca l Model 20 2.4 N u m e r i c a l Method 30 2.4 .1 F r e e - S u r f a c e Models 30 2 .4 .2 S o l u t i o n Procedure 32 2.5 Summary of Assumptions 40 2.6 Implementat ion of the Model 42 2.7 Summary 46 References 50 T a b l e s . . . . . . 53 F i g u r e s 55 i v TABLE OF CONTENTS - Cont'd Page CHAPTER 3 - GROUNDWATER FLOW IN MOUNTAINOUS TERRAIN 58 3.1 I n t r o d u c t i o n 58 3.2 Method of P r e s e n t a t i o n 60 3.3 F a c t o r s C o n t r o l l i n g Groundwater Flow i n Mountains . . . . 6 6 3 .3 .1 S lope P r o f i l e 66 3 . 3 . 2 I n f i l t r a t i o n R a t i o 67 3 . 3 . 3 T o p o g r a p h i c R e l i e f and B a s i n Depth 73 3 . 3 . 4 T o p o g r a p h i c Symmetry 75 3 . 3 . 5 Permeable H o r i z o n s and F r a c t u r e Zones 78 3 . 3 . 6 G l a c i e r s 81 3 . 3 . 7 Non-Uni form I n f i l t r a t i o n 84 3 . 3 . 8 B a s a l Heat Flow 86 3 . 3 . 9 Thermal C o n d u c t i v i t y • • • • 87 3 .3 .10 S u r f a c e Temperature 87 3. 4 D i s c u s s i o n 88 3.5 C o n c l u s i o n s 96 Re ferences 99 T a b l e s 101 F i g u r e s 110 v TABLE OF CONTENTS - Cont'd Page CHAPTER 4 - THERMAL REGIMES IN MOUNTAINOUS TERRAIN 122 4.1 I n t r o d u c t i o n 122 4.2 S i m u l a t i o n Parameters 124 4.3 R e s u l t s of N u m e r i c a l S i m u l a t i o n s 126 4 .3 .1 P a t t e r n of A d v e c t i v e Thermal D i s t u r b a n c e 126 4 . 3 . 2 Magnitude of A d v e c t i v e Thermal D i s t u r b a n c e . . . . 1 3 0 4 . 3 . 3 I n f l u e n c e of Water T a b l e C o n f i g u r a t i o n 134 4 . 3 . 4 Onset of A d v e c t i v e Thermal D i s t u r b a n c e 137 4 . 3 . 5 F r e e - C o n v e c t i o n i n Mountainous T e r r a i n 140 4 . 3 . 6 I n f l u e n c e of Mountain Topography . . . . . 1 4 4 4 . 3 . 7 Permeable F r a c t u r e Zones and Thermal S p r i n g s . . 1 4 8 4.4 I m p l i c a t i o n s for Heat Flow S t u d i e s 154 4.5 C o n c l u s i o n s 161 References 166 T a b l e s 169 F i g u r e s 171 CHAPTER 5 - SUMMARY OF CONCLUSIONS 182 Cumula t ive Reference L i s t 187 APPENDIX I - Nomenclature 195 APPENDIX II - D i s t i n c t i o n Between P o i n t of detachment and Hinge P o i n t on a Seepage Face 198 v i LIST OF TABLES Table Page 2.1 T y p i c a l S i m u l a t i o n Parameters - Chapter 2 53 2.2 I n f l u e n c e of P e r m e a b i l i t y ku and I n f i l t r a t i o n / on Maximum Temperatures in Flow System T x 54 3.1 S i m u l a t i o n Summary and Guide to I l l u s t r a t i o n s - Chapter 3 101 3.2 T y p i c a l S i m u l a t i o n Parameters - Chapter 3 102 3.3 S i m u l a t i o n R e s u l t s 103 4.1 S i m u l a t i o n Summary and Guide to I l l u s t r a t i o n s - Chapter 4 169 4.2 T y p i c a l S i m u l a t i o n Parameters - Chapter 4 170 v i i LIST OF FIGURES Figure Page 1 .1 H y p o t h e t i c a l groundwater flow systems for homogeneous p e r m e a b i l i t y 13 2.1 C o n c e p t u a l model for groundwater flow and heat t r a n s f e r i n mountainous t e r r a i n 55 2.2 T y p i c a l f i n i t e element meshes 56 2.3 Temperature f i e l d s and groundwater flow p a t t e r n s as a f u n c t i o n of upper zone p e r m e a b i l i t y ku f o r f i x e d water t a b l e c o n f i g u r a t i o n s i n concave and convex s l o p e p r o f i l e s 57 3.1 Groundwater flow p a t t e r n s , water t a b l e c o n f i g u r a t i o n s and r e c h a r g e - d i s c h a r g e p r o f i l e s as a f u n c t i o n of s lope p r o f i l e 110 * 3.2 I n f l u e n c e of i n f i l t r a t i o n r a t i o / on water t a b l e c o n f i g u r a t i o n s , flow p a t t e r n s and r e c h a r g e - d i s c h a r g e p r o f i l e s 111 3.3 Water t a b l e e l e v a t i o n WT„„V as a f u n c t i o n of ma x i n f i l t r a t i o n r a t i o I and d i m e n s i o n l e s s t o t a l flow Q* 112 3.4 M o d i f i e d s u r f a c e t o p o g r a p h i e s 113 3 .5 I n f l u e n c e of r a d i a l symmetry on water t a b l e c o n f i g u r a t i o n s , flow p a t t e r n s and r e c h a r g e - d i s c h a r g e p r o f i l e s 114 3 .6 I n f l u e n c e of r i d g e asymmetry on water t a b l e c o n f i g u r a t i o n , flow p a t t e r n s and r e c h a r g e - d i s c h a r g e p r o f i l e 115 v i i i LIST OF FIGURES - Cont'd Figure Page 3.7 I n f l u e n c e of v a l l e y asymmetry on water t a b l e c o n f i g u r a t i o n , f low p a t t e r n s and r e c h a r g e - d i s c h a r g e p r o f i l e 116 3.8 I n f l u e n c e of a permeable h o r i z o n 100 m t h i c k wi th p e r m e a b i l i t y 10 t imes ku l o c a t e d 100 m below the v a l l e y f l o o r 117 3 .9 I n f l u e n c e of a s t e e p l y d i p p i n g f r a c t u r e zone tha t o u t c r o p s at the v a l l e y f l o o r 118 3.10 I n f l u e n c e of a g l a c i e r e x t e n d i n g from an e l e v a t i o n of 2000 m to 1500 m above the v a l l e y f l o o r on water t a b l e c o n f i g u r a t i o n , flow p a t t e r n and r e c h a r g e - d i s c h a r g e p r o f i l e 119 3.11 I n f l u e n c e of b a s a l heat flow on water t a b l e c o n f i g u r a t i o n , f low p a t t e r n and thermal regime 120 3.12 Summary of the i n f l u e n c e of f a c t o r s c o n t r o l l i n g water t a b l e e l e v a t i o n (WTmax) and t o t a l f low (Q) . . . . .121 4.1 P a t t e r n s of groundwater flow and heat t r a n s f e r i n convex and concave topography 171 4.2 Contour p l o t s of temperature r e s i d u a l 172 4.3 M a t c h i n g water t a b l e c o n f i g u r a t i o n s i n convex topography 173 4.4 I n f l u e n c e of i n f i l t r a t i o n Iz and p e r m e a b i l i t y ku on water t a b l e e l e v a t i o n s and a d v e c t i v e t h r e s h o l d s 174 i x LIST OF FIGURES - Cont'd Figure Page 4.5 Thermal regimes and p a t t e r n s of groundwater flow in mixed f r e e - and f o r c e d - c o n v e c t i o n s c e n a r i o s 175 4.6 Thermal regimes beneath mountain v a l l e y s s i m u l a t e d wi th r e f e r e n c e c o n d i t i o n s of T a b l e 4.1 176 4.7 Comparison of thermal regimes i n r a d i a l and p l a n a r symmetry 1 77 4.8 I n f l u e n c e of a s t e e p l y d i p p i n g f r a c t u r e zone w i t h u n i f o r m kj--b on t h e r m a l regimes w i t h i n the asymmetric topography of F i g u r e 4.5a 178 4.9 S p r i n g temperature as a f u n c i o n of upper zone p e r m e a b i l i t y ku, t r a n s m i s s i v i t y of the f r a c t u r e zone kj'b, and b a s a l heat f low H^ 179 4.10 Contour p l o t s of t emperature g r a d i e n t r a t i o 180 4.11 I n f l u e n c e of a g l a c i e r m a n t l i n g a convex mounta in , e x t e n d i n g from 1500 m to 2000 m above the v a l l e y f l o o r 181 I I . 1 Mapping water t a b l e c o n f i g u r a t i o n s and flow p a t t e r n s from p h y s i c a l to hodograph p l a n e s - APPENDIX II 202 x ACKNOWLEDGEMENTS I t i s w i t h grea t p l e a s u r e t h a t I acknowledge the generous s u p p l y of i n s i g h t , encouragement and f i n a n c i a l support p r o v i d e d by L e s l i e Smith throughout the c o u r s e of t h i s p r o j e c t . T h i s r e s e a r c h has b e n e f i t t e d from the many hours spent d i s c u s s i n g the thermal i m p l i c a t i o n s w i t h D a v i d Chapman. I would a l s o l i k e to thank A l l a n F r e e z e and Gedeon Dagan f o r s h a r i n g t h e i r views on the development of seepage faces i n mountainous t e r r a i n and the use of hodograph methods. B r i a n F a i r b a n k and the s t a f f of Nevin S a d l i e r - B r o w n Goodbrand C o . L t d . (NSBG) are to be thanked for t a k i n g an a c t i v e i n t e r e s t i n t h i s r e s e a r c h and p r o v i d i n g acces s to t h e i r o f f i c e f a c i l i t i e s and t e c h n i c a l l i b r a r y . The e f f o r t s of my s u p e r v i s o r y committee Tom Brown, A l l a n Freeze and John Ross are g r a t e f u l l y acknowledged. Thanks are extended to C h a r l e s Mase for p r o v i d i n g the sof tware to compute water p r o p e r t i e s and to Gord Hodge for d r a f t i n g the f i g u r e s . T h i s work was funded by g r a n t s from the N a t u r a l S c i e n c e s and E n g i n e e r i n g Research C o u n c i l of Canada (NSERC) and a GREAT award p r o v i d e d by the B r i t i s h Columbia S c i e n c e C o u n c i l ( i n c o o p e r a t i o n wi th NSBG). Computat ions were c a r r i e d out on an FPS 164/MAX A r r a y P r o c e s s o r supported by an NSERC Major I n s t a l l a t i o n Grant to the U n i v e r s i t y of B r i t i s h C o l u m b i a . x i 1 CHAPTER 1 INTRODUCTION Recent i n t e r e s t i n the thermal regimes of mountains has deve loped i n response to c o n t i n u i n g s t u d i e s of the e a r t h ' s thermal s t a t e and the d e s i r e to e x p l o i t geothermal energy r e s o u r c e s . A number of thermal data se t s c o l l e c t e d i n mountainous t e r r a i n i n d i c a t e a s t r o n g a d v e c t i v e d i s t u r b a n c e of c o n d u c t i v e thermal regimes by groundwater flow (Black et al. , 1983; Mase et al. , 1982; Reader and F a i r b a n k , 1983). T h i s d i s t u r b a n c e c o m p l i c a t e s the i n t e r p r e t i o n of b o r e h o l e temperature logs used i n e s t i m a t i n g r e g i o n a l heat f low and may mask the thermal s i g n a t u r e of u n d e r l y i n g geothermal systems. When i n t e r p r e t i n g a thermal da ta s e t , or d e v e l o p i n g c o n c e p t u a l models of geothermal systems i n mountainous t e r r a i n , i t i s e s s e n t i a l t h a t the p o t e n t i a l impact of groundwater flow be a d d r e s s e d . C o n d u c t i v e thermal regimes are m o d i f i e d when heat i s r e d i s t r i b u t e d by f l u i d f low. F l u i d flow may be d r i v e n by t o p o g r a p h i c r e l i e f i n r e g i o n a l groundwater f low systems ( a d v e c t i o n or f o r c e d - c o n v e c t i o n ) , by t e m p e r a t u r e - i n d u c e d d e n s i t y c o n t r a s t s i n r e g i o n s of h i g h heat flow ( thermal b u o y a n c y - d r i v e n or f r e e - c o n v e c t i o n ) or by a c o m b i n a t i o n of both mechanisms. S t u d i e s of b u o y a n c y - d r i v e n flow systems in 2 f l a t topography are rev iewed by Cheng (1978) . These s t u d i e s , common i n the geothermal l i t e r a t u r e , emphasize f l u i d c i r c u l a t i o n w i t h i n a l o c a l i z e d r e g i o n of the f low system (or w i t h i n a geothermal r e s e r v o i r ) and n e g l e c t the p o s s i b l e i n f l u e n c e of a d v e c t i v e heat t r a n s f e r by r e g i o n a l groundwater f low. E a r l y e f f o r t s to examine the i n f l u e n c e of r e g i o n a l flow are d e s c r i b e d by P r a t s (1966) and E l d e r (1967) . U s i n g n u m e r i c a l methods, they show how f r e e - c o n v e c t i o n c e l l s can be m o d i f i e d and o b l i t e r a t e d when a r b i t r a r i l y - d e f i n e d e x t e r n a l f lows are imposed on s imple b u o y a n c y - d r i v e n flow systems. Model s t u d i e s tha t attempt to l i n k , i n a q u a n t i t a t i v e manner, a thermal d i s t u r b a n c e wi th the r e g i o n a l groundwater flow system are p r e s e n t e d by Faus t et al. , (1984), B l a c k w e l l (1985) , and Sorey (1985a) , among o t h e r s . I n t e r e s t in the thermal regimes of sed imentary b a s i n s has prompted f i e l d a n a l y s e s , and model s t u d i e s , tha t r e c o g n i z e the i n f l u e n c e of t o p o g r a p h i c a l l y - d r i v e n groundwater flow on heat t r a n s f e r i n l o w - r e l i e f t e r r a i n . T y p i c a l f i e l d s t u d i e s are d e s c r i b e d by Majorowicz et al. (1985) and G o s n o l d (1985) . Model s t u d i e s are p r e s e n t e d by Domenico and P a l c i a u s k a s (1973) , Smith and Chapman (1983, 1985), Garven and Freeze (1984) , and Woodbury and Smith (1985). I t i s commonly r e c o g n i z e d that t o p o g r a p h i c a l l y - d r i v e n groundwater flow causes a d v e c t i v e heat t r a n s f e r in . 3 mountainous t e r r a i n . The d e t a i l s of groundwater flow w i t h i n the mountain m a s s i f , however, a r e p o o r l y u n d e r s t o o d . Upper r e g i o n s of flow have been e x p l o r e d , to a l i m i t e d e x t e n t , i n f i e l d s t u d i e s tha t emphasize the i n t e r f a c e between s u r f a c e h y d r o l o g y and sha l low groundwater flow ( H a l s t e a d , 1969; S k l a s h and F a r v o l d e n , 1979; B o r t o l a m i et al. , 1979; M a r t i n e c et al. , 1982; Smart , 1985). A l t h o u g h these s t u d i e s p r o v i d e i n s i g h t i n t o the h y d r o l o g y of a l p i n e watersheds and the r e l a t i o n s h i p s between water t a b l e f l u c t u a t i o n s and s e a s o n a l snowmelt, they y i e l d l i t t l e i n f o r m a t i o n on the deep flow systems that c o n t r o l a d v e c t i v e heat t r a n s f e r i n mountainous t e r r a i n . I n f o r m a t i o n on the c h a r a c t e r and d i s t r i b u t i o n of permeable zones w i t h i n a mountain mass i f i s a v a i l a b l e i n r e p o r t s d e s c r i b i n g i n f l o w s to a l p i n e t u n n e l s (Schadt , 1905; Fox , 1907; Hennings , 1910; K e a y s , 1928; Mears , 1932). U n f o r t u n a t e l y , measurements of f l u i d p r e s s u r e that c o u l d a i d in d e f i n i n g the na ture of mountain flow systems are g e n e r a l l y l a c k i n g . Jamier (1975) a s se s sed the h y d r a u l i c c h a r a c t e r i s t i c s of f r a c t u r e d c r y s t a l l i n e rock deep w i t h i n Mount B lanc (France) on the b a s i s of geochemica l and h y d r a u l i c data o b t a i n e d d u r i n g c o n s t r u c t i o n of a highway t u n n e l , ye t an i n t e g r a t e d d e s c r i p t i o n of the flow system w i t h i n the mountain mass i f was not a t t empted . Summit water t a b l e and h y d r a u l i c head data are r a r e because most w e l l s and b o r e h o l e s are l o c a t e d on the lower f l a n k s of mountain s l o p e s . Two summit water l e v e l measurements are noted in the 4 l i t e r a t u r e ; at a depth of 30 m i n f r a c t u r e d c r y s t a l l i n e rock a t M t . Kobau, B . C . ( H a l s t e a d , 1969) and at a depth of 488 m in the b a s a l t s of M t . K i l a u e a , Hawai i ( Z a b l o c k i et al. , 1974). N u m e r i c a l s t u d i e s of m o u n t a i n - s c a l e f low systems have been p r e s e n t e d by Jamieson and F r e e z e (1983) and I n g e b r i t s e n and Sorey (1985). Jamieson and F r e e z e (1983) used a f r e e - s u r f a c e model and a water budget approach to e s t i m a t e the range of h y d r a u l i c c o n d u c t i v i t y that might be found w i t h i n Meager M o u n t a i n , B r i t i s h C o l u m b i a . I n g e b r i t s e n and Sorey (1985) i n c o r p o r a t e the i n f l u e n c e of h i g h - r e l i e f topography i n t h e i r n u m e r i c a l model of a d v e c t i o n - d o m i n a t e d heat t r a n s f e r at M t . L a s s e n , C a l i f o r n i a . In u s i n g a b a s a l source of heated f l u i d to r e p r e s e n t the c i r c u l a t i o n and h e a t i n g of groundwater recharged on the mountain f l a n k s , the r e l a t i o n s h i p between the magnitude of the f l u i d source and r a t e s of groundwater recharge at the ground s u r f a c e i s p o o r l y d e f i n e d . In t h i s t h e s i s , a model ing approach i s adopted tha t s p e c i f i e s recharge to the flow system as i n f i l t r a t i o n at the ground s u r f a c e and s o l v e s for the temperature of groundwater c i r c u l a t i n g through a mountain m a s s i f . The nature of deep groundwater flow i s of p a r t i c u l a r i n t e r e s t in s t u d i e s of geothermal systems in mountainous t e r r a i n . Groundwater samples o b t a i n e d from s p r i n g s and b o r e h o l e s d u r i n g geothermal e x p l o r a t i o n o f t e n p r o v i d e 5 geochemica l i n d i c a t i o n s of a r e s o u r c e at d e p t h . I d e n t i f y i n g the source of a c h e m i c a l s i g n a t u r e r e q u i r e s an u n d e r s t a n d i n g of the r a t e s and p a t t e r n s of groundwater f low. E f f o r t s to i d e n t i f y a geothermal r e s o u r c e a l s o r e l y on temperature data c o l l e c t e d i n shal low b o r e h o l e s . U n f o r t u n a t e l y , a d v e c t i v e d i s t u r b a n c e of c o n d u c t i v e thermal regimes by groundwater flow can c o m p l i c a t e the i n t e r p r e t a t i o n of b o r e h o l e temperature l ogs and may mask the thermal s i g n a t u r e of an u n d e r l y i n g r e s o u r c e . Geochemica l and thermal data have been used by s e v e r a l workers (Lahsen and T r u j i l l o , 1975; B l a c k w e l l and S t e e l e , 1983; S o r e y , 1985b) to form g e n e r a l i z a t i o n s on the n a t u r e of mountain h y d r o t h e r m a l systems. A comprehens ive , q u a n t i t a t i v e a n a l y s i s of groundwater flow systems i n mountainous t e r r a i n has ye t to be r e p o r t e d in the l i t e r a t u r e . The o b j e c t i v e of t h i s t h e s i s i s t o i n v e s t i g a t e the nature of t o p o g r a p h i c a l l y - d r i v e n groundwater flow systems in mountainous t e r r a i n and t h e i r i n f l u e n c e on c o n d u c t i v e thermal reg imes . A n u m e r i c a l model ing approach p r o v i d e s a q u a n t i t a t i v e b a s i s for examining the. i n f l u e n c e of f a c t o r s c o n t r o l l i n g the r a t e s and p a t t e r n s of groundwater flow and heat t r a n s f e r . In t h i s s t u d y , i d e a l i z e d mountainous t e r r a i n i s modeled for a range of c o n d i t i o n s r e p r e s e n t a t i v e of the Western C o r d i l l e r a in N o r t h A m e r i c a . Mountainous t e r r a i n i s d e f i n e d as rugged topography w i t h l o c a l r e l i e f i n excess of 600 m (Thompson, 1964). In the Coast Mountains of B r i t i s h Columbia and the c e n t r a l Cascades of the P a c i f i c Northwes t , 6 t o p o g r a p h i c r e l i e f of 2 km over a h o r i z o n t a l d i s t a n c e of 6 km i s t y p i c a l . In the Rocky Mounta ins of Canada and U . S . A . , a more subdued r e l i e f of 1 km over 6 km i s not uncommon. V e r t i c a l s e c t i o n s and schemat ic groundwater p a t h l i n e s r e p r e s e n t a t i v e of the Coast Mounta ins of B r i t i s h Columbia and the Rocky Mounta ins at the A l b e r t a - B r i t i s h Columbia border are shown i n F i g u r e s 1.1a and 1.1b. For c o m p a r i s o n , F i g u r e 1.1c shows flow systems i n a l o w - r e l i e f topography s i m i l a r to those d e s c r i b e d by Freeze and Witherspoon (1967). Groundwater flow systems in mountains d i f f e r from those found i n l o w - r e l i e f t e r r a i n i n two important r e s p e c t s : 1. F o r a g i v e n set of c o n d i t i o n s , w i t h g r e a t e r t o p o g r a p h i c r e l i e f a g r e a t e r range i n water t a b l e e l e v a t i o n and form i s p o s s i b l e . In l o w - r e l i e f t e r r a i n , water t a b l e c o n f i g u r a t i o n s can be d e f i n e d wi th r e a s o n a b l e a c c u r a c y u s i n g water l e v e l e l e v a t i o n s and h y d r a u l i c head d a t a o b t a i n e d from b o r e h o l e s and w e l l s l o c a t e d a c r o s s the r e g i o n of i n t e r e s t . In many i n s t a n c e s , e s t i m a t e d water t a b l e e l e v a t i o n s are used i n d e f i n i n g the upper boundary of r e g i o n a l f low sys tems . In mountainous t e r r a i n measured water t a b l e e l e v a t i o n s and h y d r a u l i c head data are sparse and , where a v a i l a b l e , u s u a l l y c o n c e n t r a t e on the lower f l a n k s of mountain s l o p e s . T h i s r e s t r i c t e d d i s t r i b u t i o n of data l eads to c o n s i d e r a b l e u n c e r t a i n t y in d e f i n i n g water t a b l e c o n f i g u r a t i o n s beneath 7 mountain summits. 2. H i g h - r e l i e f t e r r a i n enhances groundwater c i r c u l a t i o n to depths where e l e v a t e d s u b s u r f a c e temperatures may be e n c o u n t e r e d . War ing , 1965; Souther and H a l s t e a d , 1973 invoke t h i s mechanism to e x p l a i n the presence of thermal s p r i n g s in mountain v a l l e y s where the r e g i o n a l heat f l u x d i f f e r s l i t t l e from the n o r m a l l y a c c e p t e d median of 60 mW/m 2. V a r i a t i o n in temperature has a s t r o n g e f f e c t on f l u i d d e n s i t y and v i s c o s i t y t h a t , in t u r n , have an important i n f l u e n c e on the r a t e s and p a t t e r n s of groundwater f l o w . T h e r m a l l y - i n d u c e d d i f f e r e n c e s i n f l u i d d e n s i t y produce a b u o y a n c y - d r i v e n component of f l u i d flow t h a t enhances v e r t i c a l movement of groundwater . In a d d i t i o n , reduced f l u i d v i s c o s i t y i n r e g i o n s of e l e v a t e d temperature c o n t r i b u t e s to i n c r e a s e d r a t e s of groundwater f low. Chapter 2 e x p l a i n s how the c h a r a c t e r of mountainous t e r r a i n i s i n c o r p o r a t e d i n a c o n c e p t u a l model for f l u i d flow and heat t r a n s f e r . The n u m e r i c a l method used here d i f f e r s from those of p r e v i o u s work because a f r e e - s u r f a c e approach i s used i n a n o n - i s o t h e r m a l f o r m u l a t i o n to e s t i m a t e water t a b l e e l e v a t i o n s w i t h i n a mountain m a s s i f . Important a s p e c t s of the c o n c e p t u a l and n u m e r i c a l models are summarized at the c o n c l u s i o n of Chapter 2. In Chapters 3 and 4, the n u m e r i c a l model i s used to examine the i n t e r a c t i o n of groundwater flow and heat t r a n s f e r in mountainous t e r r a i n . N u m e r i c a l r e s u l t s p r e s e n t e d i n Chapter 3 i l l u s t r a t e the i n f l u e n c e of geo logy , climate, surface topography and regional heat flux on water table elevations and the pattern and magnitude of groundwater flow. Numerical results presented in Chapter 4 i l l u s t r a t e the impact of topographically-driven groundwater flow on thermal regimes and the role that advective heat transfer plays in masking the regional heat f l u x . Conditions leading to the development of free-convection c e l l s within a mountain massif are also considered. Conclusions presented at the end of Chapters 3 and 4 are summarized in Chapter 5. 9 REFERENCES Black, G.L., D.D. Blackwell, and J.L. Steele, Heat flow in the Oregon Cascades, in Geology and Geothermal Resources of the Central Oregon Cascades Range, ed. by G.R. Priest and B.F. Vogt, Oregon Dept. of Geol. and Min. Ind., Special Paper 15, 69-76, 1983. Blackwell, D.D., A transient model of the geothermal system of Long Valley Calder, C a l i f o r n i a , Jour. Geophys. Res., 90, 11229-11242, 1985. Blackwell, D.D., and J.L. Steele, A summary of heat flow studies in the Cascade Range, Geoth. Resources Council Trans., 7, 233-236, 1983. Bortolami, G.S., B. R i c c i , G.G. Susella, and G.M. Zuppi, Hydrogeochemistry of the Corsaglia Valley, Maritime Alps Piedmont It a l y , Jour, of Hydrol., 44, 57-79, 1979. Cheng, P., Heat transfer in geothermal systems, Adv. in Heat Transfer, Vol. 14, 1-105, 1978. Domenico,P.A. and V.V. Palciauskas, Theoretical analysis of forced convective heat transfer in regional ground-water flow, Geol. Soc. Amer., B u l l e t i n 84, 3303-3814, 1973. Elder, J.W., Steady free convection in a porous medium heated from below, J . F l u i d Mech., 27, 29-48, 1967. Faust, C.R., J.W. Mercer, S.D. Thomas, and W.P. Balleau, Quantitative analysis of existing conditions and production strategies for the Baca Geothermal System, New Mexico, Water Resour. Res., 20(5), 601-618, 1984. Fox, F.M., The Simplon tunnel, Minutes of the Proceedings of the Inst, of C i v i l Engineers, 168, 61-86, 1907. Freeze, R.A., and P.A. Witherspoon, Theoretical analysis of regional groundwater flow: 2. Effect of water-table configuration and subsurface permeability v a r i a t i o n , Water Resour. Res., 3, 623-634, 1967. Garven, G. and R.A. Freeze, Theoretical analysis of the role of groundwater flow in the genesis of stratabound ore -deposits: 2. Quantitative results, Am. Jour. S c i . , 284, 1125-1174, 1984. Gosnold, W.D., Heat flow and ground water flow in the great plains of the United States, Jour, of Geodynamics, 4, 247-264, 1985. Halstead, E.C., Groundwater investigation, Mount Kobau, B r i t i s h Columbia, Canada Inland Waters Branch, Dept. of Energy Mines and Resources, Tech. B u l l . 17, 1969. Hennings, F., On the question of long railway tunnels Construction, v e n t i l a t i o n and operation, B u l l , of the Internl. Railway Congress Assoc., 24(10), 943-983, 1910 Ingebritsen, S.E., and M.L. Sorey, A quantitative analysis of the Lassen hydrothermal system, North Central C a l i f o r n i a , Water Resour. Res., 21(6), 853-868, 1985. Jamier, D., Etude de l a f i s s u r a t i o n , de 1'hydrogeologie et de l a geochemie des eaux profondes des massifs de l ' A r p i l l e et du Mont Blanc, PhD. d i s s e r t a t i o n , Faculte des Sciences, Universite de Neuchatel, Switzerland (in French), 1975. 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. Keays, R.N., Construction methods in the Moffat Tunnel, Amer. Soc. C i v i l Eng., 62, 63-112, 1928. Lahsen, A., and P. T r u j i l l o , The geothermal f i e l d of E l Tatio, Chile, Proc. 2nd. U.N. Symp. on Development and Use of Geothermal Resources, San Francisco, C a l i f . , Vol 1, 157-175, 1975. Majorowicz, J.A., F.W. Jones, H.L. Lam, and A.M. Jessop, T e r r e s t r i a l heat flow and geothermal gradients in r e l a t i o n to hydrodynamics in the Alberta Basin, Canada, Jour. Geodynamics, 4, 265-283, 1985. Martinec, J., H. Oeschager, U. Schotterer, and U. Siegenthaler, Snowmelt and groundwater storage in an alpine basin, in, Hydrological Aspects of Alpine and High Mountainous Terrain, ed. J. W. Glen, IAHS Publ. 138, 169-175, 1982. 11 Mase, C.W., J.H. Sass, A.H. Lachenbruch, and J.R. Munroe, Preliminary heat-flow investigations of the C a l i f o r n i a Cascades, U.S. Geol. Survey, Open-File Report 82-150. 1982 Mears, F., The eight mile Cascade Tunnel, Great Northern Railway Part I I . Surveys, construction methods and a comparison of routes, Amer. Soc. C i v i l Eng., Trans., 96, 915-1004, 1932. Prats, M. , The e f f e c t of horizontal f l u i d flow on thermally induced convection currents in porous mediums, Jour. Geophys. Res., 71(2), 4835-4837, 1966. Reader, J.F., and Fairbank, B.D., Heat flow in the v i c i n i t y of the Meager Volcanic Complex, Southwestern B r i t i s h Columbia, Geoth. Resour. Council Trans., Vol. 7, 535-539, 1983. Schadt, M.H., Les resultats s c i e n t i f i q u e s du percement du tunnel du Simplon, B u l l e t i n Technique de l a Suisse Romande, 125-178, 1905. S.klash, M.G., and R.N. Farvolden, The role of groundwater in storm runoff, Jour, of Hydrol., 43, 45-65, 1979. Smart, C.C., The hydrology of the Castleguard Karst, Columbia I c e f i e l d s , Alberta, Canada, A r c t i c and Alpine Research, 15(4), 471-486, 1985. Smith, L., and D.S. Chapman,The influence of water table configuration on the near-surface thermal regime, Jour, of Geodynamics, 4, 183-198, 1985. Smith, L., and D.S. Chapman, On the thermal effects of groundwater flow 1. Regional scale systems, Jour. Geophys. Res., 88(B1), 593-608, 1983. Sorey, M.L., Evolution and present state of the hydrothermal system in Long Valley Caldera, Jour. Geophys. Res., 90, 11219-11228, 1985a. Sorey, M.L., Types of hydrothermal convection systems in the Cascade Range of C a l i f o r n i a and Oregon, Proc. Workshop on Geothermal Resources of the Cascade Range, ed. M. Guffanti and L.J.P. Muffler, 63-67, 1985b. Souther, J.G., and Halstead, E.C., Mineral and thermal waters of Canada, Canada Dept. of Energy, Mines and Resources, Paper 73-18, 1973. Thompson, W.T., How and why to d i s t i n g u i s h between mountains and h i l l s , Prof. Geographer, 16, 6-8, 1964. Waring, G.A., Thermal springs of the United States and other countries of the world - a summary, U.S. Geol. Survey Prof. Paper 492, 1965. Woodbury, A.D., and L. Smith, On the thermal effects of three-dimensional groundwater flow, Jour. Geophys. Res., 90, 759-767, 1985. Zablocki. C.J., R.I. T i l l i n g , D.W. Peterson, R.L. Christiansen, G.V. K e l l e r , and J.C. Murray, A deep research hole at the summit of an active volcano, Kilauea, Hawaii, Geoph. Res. Lett., 1(7), 323-326, 1974. 13 Figure 1 . 1 . Hypothetical groundwater flow systems for homogeneous permeability; a. Coast Mountains of B r i t i s h Columbia, b. Rocky Mountains of B.C.-Alberta and c. conventional low-relief terrain after Freeze and Witherspoon, 1967. CHAPTER 2 DEVELOPMENT OF THE NUMERICAL METHOD 2.1 INTRODUCTION This chapter describes the conceptual model, mathematical formulation and numerical method developed to simulate steady groundwater flow and heat transfer in mountainous t e r r a i n . The method adopted here d i f f e r s from those used by previous workers because water table elevations are found as part of a solution process that accounts for f l u i d flux and heat transfer in the unsaturated zone. This approach i s p a r t i c u l a r l y useful when studying regions where limited groundwater recharge causes deep water tables and thick unsaturated zones. In thi s approach, recharge rates are cont r o l l e d as an input variable, rather than i m p l i c i t l y calculated in the solution process. As a consequence, recharge rates consistent with the clim a t i c conditions within the region of interest can be ensured. In addition, r e s u l t s of s e n s i t i v i t y analyses performed with this method are unaffected by the changes in recharge rate that would normally result when changing simulation parameters such as permeability. 15 2.2 CONCEPTUAL MODEL A schematic of the boundary value problem adopted in th i s study i s shown, without v e r t i c a l exaggeration, in Figure 2.1. The upper boundary of the domain i s the bedrock surface. Erosional processes operating in mountainous te r r a i n often promote development of a thin cover of discontinuous s u r f i c i a l deposits (often less than 10 m thick) over upland areas of mountain slopes. In t h i s conceptual model, these deposits are thought of as a thin skin of variable thickness that i s not e x p l i c i t l y included in the model. Subsurface flow within t h i s skin, in addition to overland flow and evapotranspiration, i s lumped in a single runoff term. These processes are strongly affected by sp a t i a l variations in p r e c i p i t a t i o n , slope angle and s o i l permeability as well as by temporal variations in pr e c i p i t a t i o n events. F i e l d observations of the complex nature and interaction of these factors are lacking for individual mountain slopes. As a f i r s t approximation, therefore, a lumped steady state approach i s adopted and an available i n f i l t r a t i o n rate is defined. Available i n f i l t r a t i o n Uz) represents the maximum rate of recharge possible at the bedrock surface for specified c l i m a t i c , geologic and topographic conditions. In the absence of detailed s i t e data, the available i n f i l t r a t i o n rate is best thought of as a percentage of the mean annual p r e c i p i t a t i o n 16 rate. Although the units of / are m/sec, t h i s parameter actually expresses the volumetric flux of f l u i d across a horizontal plane (m3/sec per m 2). Recharge to the flow system r e f l e c t s the magnitude of the available i n f i l t r a t i o n rate, the capacity for f l u i d flow through the system and the nature of the thermal regime. In high permeability t e r r a i n , groundwater flow systems may accept a l l the available i n f i l t r a t i o n and produce a water table that l i e s below the bedrock surface (Figure 2.1). In lower permeability t e r r a i n , where recharge accepted by the flow system i s less than available i n f i l t r a t i o n , the water table w i l l be found close to the bedrock surface. Conventional approaches to modeling regional groundwater flow (Freeze and Witherspoon, 1967) recognize that a t r a n s i t i o n from groundwater recharge to groundwater discharge occurs at a specified point on the upper boundary. In this study, t h i s point i s termed the hinge point (HP shown in Figure 2.1). Because a s i m p l i f i e d representation of the unsaturated zone i s included in t h i s conceptual model i t is also necessary to i d e n t i f y the point where the water table meets the bedrock surface. This point i s defined as a point of detachment (POD shown in Figure 2.1). In developing numerical models to analyze seepage through earth dams, an exit point is commonly defined that has the properties of both the POD and the HP. In h i g h - r e l i e f t e r r a i n , the usual exit point cannot be i d e n t i f i e d because the POD and HP no 17 longer coincide. In the region between the point of detachment and the hinge point, the water table coincides with the bedrock surface and the saturated zone i s recharged d i r e c t l y from overlying s u r f i c i a l deposits or by p r e c i p i t a t i o n f a l l i n g on the bedrock surface. This suprising result implies that recharge can occur on what i s usually considered to be the seepage face. This behavior i s discussed in more d e t a i l in a subsequent section and in Appendix I I . Upslope of the point of detachment, the water table l i e s below the bedrock surface and i n f i l t r a t i o n is transferred to the water table by unsaturated flow. Although s i g n i f i c a n t l a t e r a l flow in the unsaturated zone may cause patterns of recharge at the water table to d i f f e r from patterns of i n f i l t r a t i o n at the bedrock surface, the nature and magnitude of l a t e r a l flow in unsaturated regions of mountain flow systems i s poorly understood. In this study, a one-dimensional model of v e r t i c a l flow from the bedrock surface to the water table i s adopted. This approach c l e a r l y oversimplifies the complex factors c o n t r o l l i n g f l u i d flow in regional-scale unsaturated zones. Before a more r e a l i s t i c approach can be adopted, however, additional research involving both field-based and t h e o r e t i c a l studies must be c a r r i e d out. Boundaries of the thermal regime coincide with those-of the f l u i d flow system. A basal conductive heat flow supplies thermal energy to the mountain flow system. The basal heat flow, that represents heat transfer from deeper levels in the earth's crust, i s a ch a r a c t e r i s t i c of the tectonic environment in which the mountain i s located. Chapman and Rybach (1985) report that representative heat flow values range from 30 to 150 mW/m2 with a median of 61 mW/m2. Following the example of Birch (1950), i t is assumed that temperatures at the ground surface r e f l e c t an a l t i t u d i n a l gradient in surface temperature (thermal lapse rate) with mean annual temperatures at the ground surface a few degrees warmer than mean annual a i r temperature. Smith (1975) and Barry (1981, pg. 74) indicate that the difference between mean annual s o i l temperature and mean annual a i r temperature i s strongly controlled by microclimate. As a consequence, t h i s temperature difference may range from less than 1 K to in excess of 15 K, on a mean annual basis. In the s e n s i t i v i t y studies performed in thi s thesis, mean annual ground temperatures are assumed to be controlled only by variations in elevation within the range 0 to 15°C. In t h i s study, fixed temperatures are not sp e c i f i e d where fracture zones outcrop to produce groundwater springs. In such cases, temperature at the ground surface i s assumed to r e f l e c t the temperature of groundwater flowing in the fracture zone, rather than the a i r temperature. In developing the f l u i d flow component of this conceptual model i t i s assumed that only a thin skin of s u r f i c i a l deposits i s l i k e l y to be found on mountain slopes. The resistance to heat transfer through the s u r f i c i a l deposits i s assumed to cause only minor changes in temperature, therefore, temperatures at the bedrock surface are presumed to match temperatures at the ground surface. Heat transfer occurs by conduction and advection both above and below the water table, therefore, the influence of f l u i d flow on the thermal regime in the unsaturated zone must be considered. In the model presented here, only the thermal e f f e c t s of moisture movement in the l i q u i d phase are considered here. Thermal e f f e c t s of condensation, evaporation and heat transfer by vapor movement within the unsaturated zone are neglected. In developing a conceptual model of deep unsaturated zones, Ross (1984) suggests that moisture transport by vapor movement becomes negli g i b l e when recharge rates exceed about 10" 1 2 m/sec. Unsaturated zones with recharge rates less than t h i s amount are l i k e l y found only in the most a r i d of climates. Because such a r i d conditions are few in mountainous t e r r a i n of the Western C o r d i l l e r a , the influence of vapor movement on moisture and heat transport in the unsaturated zone i s neglected and a minimum i n f i l t r a t i o n rate of 10" 1 2 m/sec is assumed. Although fractures are assumed to provide the primary pathways for groundwater flow through a mountain massif, an equivalent porous medium approach i s adopted. This approach is reasonable where the length and spacing of individual • fractures i s generally much smaller than the scale of the mountain massif. An equivalent porous media approximation i s often adopted in studies where limited f i e l d data preclude a more rigorous approach. In keeping with this lumped approach, groundwater discharge is represented as seepage di s t r i b u t e d across the discharge area (downslope of the HP in Figure 2.1) rather than at isolated groundwater springs. Local variations in surface topography, permeability, and thickness of s u r f i c i a l deposits c o n t r o l l i n g the d i s t r i b u t i o n of springs are assumed to have l i t t l e e f f e c t on the overall pattern and magnitude of groundwater flow. Exceptions to this assumption are major through-going fracture zones represented as discrete permeable fracture zones with a homogeneous permeability kj- and width b (Figure 2.1). An important component of thi s study involves assessing the impact of through-going fracture zones on groundwater flow systems and springs. 2.3 MATHEMATICAL MODEL The mathematical model for groundwater flow and heat transfer in mountains i s expressed by two coupled p a r t i a l d i f f e r e n t i a l equations describing f l u i d flow and heat transfer, by equations of state for the f l u i d properties, and by boundary conditions for both the f l u i d and thermal problems. Note that in the conceptual model outlined in Figure 2.1, a free-surface exists where the water table l i e s below the bedrock surface. Appendix I contains a table of nomenclature used in the following discussion. The steady state equation of f l u i d mass conservation in the absence of internal sources and sinks i s given by: y p f x i + hte/'iii = ° ( 2 - i ) where py is f l u i d density and q x , qz are horizontal and v e r t i c a l components of f l u i d flux ( s p e c i f i c discharge). F l u i d flow i s driven by gradients in f l u i d pressure and thermally-induced density contrasts. F l u i d potential h is defined in terms of an equivalent freshwater head: h = ^i + Z i 2 ' 2 ) where p0 i s a reference f l u i d density at a s p e c i f i e d temperature, g i s acceleration due to gravity, p i s f l u i d pressure and z i s the elevation where the freshwater head i s calculated. Frind (1982) advocates the use of an equivalent freshwater head in order to describe potential gradients exclusive of s t a t i c f l u i d pressures. F l u i d flux through the saturated porous matrix i s given by: */ - - V [<vHy + . ( ' / - ( 2 - 3 ) where a i s the dynamic v i s c o s i t y of the f l u i d and ki . i s the permeability tensor for the porous matrix. This equation i s s i m p l i f i e d by defining a r e l a t i v e density Pr = y- - 1 ( 2 . 4 ) to obtain: „ - ~ki ,• p ° g r o h + Pr0Z 1 Ci Z.\ qi - lJ-~[T>Xj + rJ5j] (2*5) In thin fracture zones and discrete fractures, f l u i d flux in d i r e c t i o n s p a r a l l e l to the fracture-matrix boundary i s given: -I where kj- i s the permeability of material within a discrete fracture or fracture zone. For open fractures, a p a r a l l e l plate model i s used to define the fracture permeability in terms of an e f f e c t i v e aperture b where kj- = 62/12 (Freeze and Cherry, 1979). In t h i s study, only fractures within the saturated region of flow are considered. The equation describing the d i s t r i b u t i o n of freshwater head in the porous matrix i s obtained by substituting equation (2.5) into equation (2.1): 3 r k„„pfpos bh , k„_pfpo8 bh-i 77 1 x x f — 37 + x z f — 77 ] + h t*« P A* H + (|| + Pr)] = 0 (2.7) The cross-terms for the permeability tensor ( k x z , kzx) are included to preserve the generality of the subsequent numerical formulation. The analagous equation in one-dimensional l o c a l coordinates for a fracture zone i s h ^kfpfP-£ [U + p ' § 7 » = 0 ( 2 - 8 ) Boundary conditions for the f l u i d flow system shown in Figure 2.1 are as follows: d h for -4 km < z < z j 7* = 0 at x = XQ and x = x^ (2.9a) dh = 0 f o r x0 ~ x ~ XL 7 7 at z = -4 km (2.9b) h = zx f o r x0 * x ~ XL for z = z * (2.9c) where xQ equals 0 . 0 km, L i s the horizontal length of the flow system ( 6 . 0 km), and z ^ i s the water table elevation at the s p e c i f i e d x-coordinate. In Figure 2.1, the free-surface is that portion of the water table that l i e s below the bedrock surface. Equation (2.9c) r e f l e c t s the fact that freshwater head equals the water table elevation everywhere at the upper boundary of the saturated region of flow ( f l u i d pressure i s atmospheric at the water table). On the free-surface segment an a d d i t i o n a l constraint i s required to l i m i t the i n f i l t r a t i o n a vailable for recharge to the bedrock flow system. F l u i d flux qn , directed along the unit normal ni to the free-surface, i s driven by p o t e n t i a l gradients and thermal buoyancy where; In keeping with the assumption of only v e r t i c a l flow of f l u i d through the unsaturated zone, the available i n f i l t r a t i o n rate Iz i s applied d i r e c t l y on the free-surface by computing the flux of f l u i d normal to a free-surface that slopes at an angle 6 from horizontal. This boundary condition d i f f e r s from that used by Neuman and Witherspoon (1970) in two respects. F i r s t , f l u i d flux and freshwater head conditions are s p e c i f i e d on the free-surface segment while only freshwater head i s s p e c i f i e d on the seepage face (Figure 2.1). This approach allows recharge to occur on the seepage face between the hinge point (HP) and the point of detachment (POD). Below the HP, only discharge occurs. Second, a buoyancy term stated in terms of the r e l a t i v e buoyancy pr i s included inside the square brackets of equation (2.10) because thermal e f f e c t s are incorporated in the formulation. Appendix II contains a discussion of the conditions that promote separation of the POD and HP and outlines use of a hodograph to demonstrate the theoretical basis for such a separation. The steady state balance of thermal energy in a variably-saturated porous medium with no internal heat sources or sinks i s described by: 77 l xxtt + XZTF] + JI[ ZXTZ + ZZTZ] ' PfCf(qxE + ' H S T > = 0 ( 2 . 1 1 ) where T = temperature Cy = s p e c i f i c heat capacity of f l u i d \ej = thermal conductivity tensor for solid-fluid-vapor composite The l a s t term on the l e f t side of equation ( 2 . 1 1 ) describes the advective transfer of heat by f l u i d flow both above and below the water table. In developing equation ( 2 . 1 1 ) , thermal equilibrium i s assumed to exist between f l u i d , s o l i d and vapor. Modeling heat transfer in the unsaturated zone ensures a water table temperature and thermal regime consistent with temperature conditions specified at the bedrock surface. By adopting a free-surface method, the rate of v e r t i c a l f l u i d flow in the unsaturated zone is i m p l i c i t l y defined to be the available i n f i l t r a t i o n rate. In solving the equations of energy transport, one-dimensional advective heat transfer above the water table is modeled by setting the horizontal f l u i d flux qx to zero and the v e r t i c a l f l u i d flux qz to the available i n f i l t r a t i o n rate J in equation (2.11). Heat conduction and thermal dispersion by f l u i d flow in the solid-liquid-vapor composite i s represented by the f i r s t two terms of equation (2.11) using a thermal conductivity tensor defined by: Xf. = SnDtJ + (/ - S)n\v + (1 - «)Xf y (2.12) where n = porosity of the porous matrix Xv = thermal conductivity of vapor Di . = conduction-dispersion tensor for f l u i d X? . = thermal conductivity tensor for s o l i d S = degree of saturation S = 1 below water table 0 < S < '\ above water table In t h i s formulation, the degree of saturation i s s p e c i f i e d only to calculate the thermal conductivity of the solid-fluid-vapor composite found above the water table. Saturation may vary from less than 0.2 to 1.0 within the unsaturated zone and r e f l e c t s the r e l a t i v e proportion of vapor and f l u i d within the pore space. Because the thermal conductivity of a i r (vapor) i s about 5 percent that of water, higher levels of saturation S contribute to greater thermal conductivity. Rocks most commonly found in mountainous t e r r a i n are l i k e l y to have porosity values less than about 15 percent. Therefore, the degree of saturation has l i t t l e impact on values calculated for thermal conductivity of the s o l i d - v a p o r - f l u i d composite. In the absence of information regarding the variation of saturation in mountainous t e r r a i n , the pore space i s assumed to be f i l l e d only with vapor (S = 0). Although a minimum degree of saturation i s required to permit f l u i d flow in the unsaturated zone, this assumption provides an upper l i m i t for the influence of saturation on mountain flow systems by contributing to a minimum estimate of thermal conductivity above the water table. Note that the thermal conductivity for a composite porous medium i s defined using the arithmetic mean to account for the thermal conductivity of . each material. This approach i s commonly adopted in the hydrogeology l i t e r a t u r e (Sorey, 1978; Smith and Chapman, 1983) and d i f f e r s from approaches found in the heat flow l i t e r a t u r e that r e f l e c t d i f f e r e n t views regarding the influence of porosity on thermal conductivity (Walsh and Decker, 1966; Scharli and Rybach, 1984). Simulation results presented in subsequent sections i l l u s t r a t e the influence of di f f e r e n t values of thermal conductivity on advective heat transfer. These results suggest that differences between thermal conductivity values calculated would have l i t t l e impact on groundwater flow systems and advective disturbance of mountain thermal regimes. Bear (1972) uses a f l u i d conduction-dispersion tensor D[j to account for heat conduction in the f l u i d and thermal dispersion by mechanical mixing. The nature of thermal dispersion in a porous medium i s described by Sauty et a l . (1982). Bear (1979) shows that the dispersion terms can be expanded for an isot r o p i c porous medium as follows: nDxx = pjCf{atq\/q + atq\/q) + n\f (2.13a) nDzz = PfCf(at42x/4 + al4\/D + n x / (2.13b) nDxz ' nDzx = PfCf^l «tUxqz/q (2.13c) Here, and a( are longitudinal and transverse d i s p e r s i v i t i e s , X-^  i s the thermal conductivity of the f l u i d and q i s the magnitude of the f l u i d flux. The analagous energy balance equation for a f u l l y saturated fracture zone i s written using one-dimensional l o c a l coordinates as: %j[nfDs + (1 - nf)\s]™ - P f C f q s g = 0 (2.14) where nfDs = p^Cf(atqs + n f \ f ) (2.15) In equations (2.14) and (2.15), /iy is the porosity of the fracture where 0 < nj- < 1 for fracture zones f i l l e d with porous material and nj equals 1 for open fractures. Where n i s less than 1, an isotropic thermal conductivity Xs i s assumed for the fracture f i l l i n g . Boundary conditions for the heat transfer problem are: for -4 km <. z < z5 at x = XQ and x = x^ (2.16a) X z z b T - H X° ~ X * % L 7 7 b at z = -4 km (2.16b) T = T r + G l z x f o r x0 * x S XL at z = 2* (2.16c) where i s the conductive basal heat flow, Tr i s a reference surface temperature s p e c i f i e d at the valley f l o o r , G[ i s a thermal lapse rate applied at the bedrock surface and zsx i s the surface elevation at the spe c i f i e d x coordinate p o s i t i o n . Thermal springs are modeled where fracture zones outcrop at the bedrock surface by assuming that the heat flux from the fracture zone i s equal to the advective flux of heat along the fracture zone. With only advective flux presumed to occur at the fracture outcrop, a sp e c i f i e d temperature is not required and the bedrock temperature is dictated by the temperature of the f l u i d discharging from the fracture zone (a third-type boundary condition) . The f l u i d i s assumed to be pure l i q u i d water with density and v i s c o s i t y a function only of temperature and bT Ix = 0 pressure. Water properties are evaluated using the relationships of Keenan et al. (1978) for density and of Watson et al. (1981) for v i s c o s i t y . Specific heat capacity and thermal conductivity of water i s assumed constant, as i s the thermal conductivity of vapor in the unsaturated zone. Porosity, permeability and thermal conductivity vary in space but are to be unaffected by changes in temperature and pressure. 2.4 NUMERICAL METHOD The Galerkin f i n i t e element method i s used to solve the coupled equations of steady f l u i d flow and heat transfer. Two-dimensional v e r t i c a l sections of porous media with planar or r a d i a l symmetry are represented by linear triangular elements. Thin, high-permeability fracture zones are included by embedding one-dimensional l i n e elements in the f i e l d of triangular elements using the method described by Baca et al . (1984). 2.4.1 Free-surface Models Free-surface methods have been developed to solve unconfined groundwater flow problems in which the position of the water table i s i n i t i a l l y unknown. Most free-surface methods involve successive estimation and correction of water table positions in an i t e r a t i v e procedure (for example; Neuman and Witherspoon, 1970). The majority of free-surface codes have been applied to r e l a t i v e l y small-scale problems of isothermal flow through porous media in systems with low topographic r e l i e f . Although Jamieson and Freeze (1983) applied the steady FREESURF model of Neuman and Witherspoon (1970) to the rugged topography of Meager Mountain, their solutions suggest that the seepage face i s poorly represented. Non-isothermal free-surface models are described by Home and 0'Sullivan (1978) and Bodvarsson and Pruess (1983). These authors applied their models to regions of f l a t topography near producing geothermal f i e l d s . An important difference between the numerical method used in thi s study and other methods l i e s in the approach used to solve the free-surface problem in t e r r a i n with high topographic r e l i e f . Conventional techniques assume that groundwater recharge occurs only above the point of detachment (POD), that i s i m p l i c i t l y assumed to coincide with the hinge point (HP). In developing t h i s modeling approach, i t was found that the POD need not coincide with the HP in steep te r r a i n with rocks of r e l a t i v e l y low permeability (Appendix I I ) . Therefore, recharge is allowed to occur downslope of the POD on the seepage face (Figure 2.1). Although this approach to viewing the nature of seepage faces has been revealed through study of mountainous te r r a i n , i t i s equally applicable when examining systems in low-relief t e r r a i n . Failure to adopt this approach when studying groundwater flow in high-reief t e r r a i n may produce solutions with u n r e a l i s t i c water table configurations. 2.4.2 Solution Procedure The coupled equations of f l u i d flow and heat transfer are solved using an i t e r a t i v e procedure. The freshwater head and temperature d i s t r i b u t i o n s are calculated using the following steps: 1. Individual f i n i t e element meshes are constructed for the f l u i d flow and thermal problems. 2. Mesh accuracy i s tested and meshes are reconstructed as necessary. 3. A temperature f i e l d i s computed for heat conduction alone to obtain i n i t i a l values for the temperature-dependent properties of water. 4. Iterations are performed to solve for freshwater head, estimate a new water table configuration, deform the mesh used in solving for freshwater head, reconstruct the mesh for the thermal problem, and solve for an updated temperature f i e l d . 5. Iterations terminate when successive water table elevations d i f f e r by no more than a spec i f i e d tolerance. The number of iterations required t y p i c a l l y ranges from 7 to 15. 6. Flux of f l u i d mass and thermal energy across boundaries of the domain are calculated to evaluate mass and energy balances. Mesh const met i on F i n i t e element meshes that can deform in a v e r t i c a l d i r e c t i o n are constructed for the f l u i d flow and thermal problems using an internal mesh generator. The f l u i d flow mesh i s generated only within the saturated region of flow. The upper boundary of the i n i t i a l mesh i s defined using a reasonable guess of the water table configuration. In subsequent i t e r a t i o n s , the current estimate of the position of the water table i s used. An example of a mesh used in the solution for freshwater head, obtained at the f i n a l water table configuration, i s shown in Figure 2.2a. Throughout the solution procedure, the structure of the mesh i s unchanged below the elevation of the valley f l o o r . Although the upper region of the mesh deforms to conform to successive estimates of the water table configuration, the number of nodes (1137) and elements (2139) remains unchanged during the i t e r a t i v e procedure. Coupling between heat transfer and f l u i d flow i s f a c i l i t a t e d by ensuring a one-to-one correspondence between nodes located below the water table in each mesh. Heat transfer above the water table i s modeled by extending the mesh used in solving the thermal problem to the bedrock surface (Figure 2.2b). Because the mesh for the f l u i d flow problem does not extend above the water table, f l u i d flow in the unsaturated zone i s i m p l i c i t l y defined by the available i n f i l t r a t i o n rate applied at the bedrock surface. Although the boundaries of the thermal problem remain fixed, the thermal mesh must be deformed i n t e r n a l l y (or reconstructed) to match changes in the f l u i d mesh. A simple deformation of the mesh for the thermal problem i s unable to account for l a t e r a l pinching out of the unsaturated zone as the length of the seepage face increases and the POD moves upslope on the bedrock surface during the i t e r a t i v e procedure. Therefore, the thermal mesh must be reconstructed in i t s entirety at each i t e r a t i o n to accomodate upslope or downslope movement of the POD. Each time the mesh i s reconstructed, the number of nodes and elements may change depending upon the degree of POD movement. A f i n a l thermal mesh, that corresponds to the f l u i d mesh of Figure 2.2a, is shown in Figure 2.2b and contains 1230 nodes and 2201 elements. Because the numerical formulation does not include a method for maintaining the position of geologic structure within deforming regions of the mesh, fracture zones and heterogeneities must be r e s t r i c t e d to non-deforming regions of the mesh. The accuracy of each i n i t i a l mesh i s evaluated by solving the isothermal f l u i d flow problem and computing a t o t a l balance of f l u i d mass crossing a l l system boundaries and a f l u i d mass balance across the top surface. Because v e r t i c a l g r i d lines in the deforming mesh are rarely orthogonal to the sloping top boundary, computing f l u i d flux normal to the surface using freshwater heads can be inaccurate i f the mesh i s too coarse. Flux inaccuracies are minimized in a two-step process. F i r s t , the free-surface method i s used in an isothermal mode to solve for freshwater head and to obtain the water table configuration. F l u i d fluxes are calculated at each boundary using the freshwater head solution. Second, the problem is reformulated using stream functions, with the previously-defined water table configuration forming the upper boundary of flow. This stream function solution provides an alternative method for c a l c u l a t i n g f l u i d fluxes normal to the water table. Large differences between boundary fluxes obtained using the two solution methods often indicate regions where the water table configuration may be poorly estimated and the f i n i t e element mesh must be refined. Acceptable grids are defined when mass flux balances for flow across the water table d i f f e r by less than 1 percent and t o t a l mass balances d i f f e r by less than 5 percent. In addition to obtaining a good match in mass flux balance, i t i s also important to ensure that patterns of f l u i d flux on the upper boundary are the same for the two d i f f e r e n t formulations. Therefore, an acceptable g r i d must also provide hinge point positions that match for each solution method. Contour plots of the isothermal stream function solutions are p a r t i c u l a r l y useful for v i s u a l i z i n g patterns of groundwater flow and for providing insight into the nature of mountain flow systems. Iterative procedure A coupled solution i s obtained by solving equations (2.7) and (2.11) in an i t e r a t i v e procedure controlled by the free-surface method. At the upper boundary of the f l u i d flow problem a mixed boundary condition i s defined. An available i n f i l t r a t i o n rate i s applied where the water table l i e s below the bedrock surface (that i s , on the free-surface) and freshwater head i s specified where the water table coincides with the bedrock surface. Prior to i n i t i a t i n g the i t e r a t i o n s , equation (2.11) is solved for the conductive case to obtain the i n i t i a l temperature f i e l d . These i n i t i a l temperatures are used to compute f l u i d properties required in the f i r s t free-surface i t e r a t i o n . The numerical solution proceeds by updating the f l u i d properties using the latest estimate of the temperature f i e l d , solving equation (2.7) for freshwater head, obtaining a new estimate of the water table configuration to calculate s p e c i f i c discharge, and updating the temperature f i e l d by solving equation (2.11). The solution begins by solving equation (2.7) for freshwater head and extrapolating a new water table position using, the method described by Neuman and Witherspoon (1970). The f i n i t e element mesh i s then deformed to conform to the shape of the new upper boundary. At each i t e r a t i o n , the POD i s moved along the top surface to maintain a balance between recharge and discharge. As the position of the POD changes, the length of the free-surface changes and the region where the available i n f i l t r a t i o n rate i s applied v a r i e s . In t h i s modified free-surface method, i t i s assumed that the point of detachment (POD) marks the uppermost point where freshwater heads are known to equal the elevation of the bedrock surface, rather than the uppermost point of groundwater discharge. Although recharge rates are not i n i t i a l l y known for the region between the POD and the hinge point (HP), they can be computed from the freshwater head solution. After each i t e r a t i o n , recharge rates computed using freshwater heads in the v i c i n i t y of the water table are compared to available i n f i l t r a t i o n rates. If recharge exceeds available i n f i l t r a t i o n near the POD, additional i t e r a t i o n s are performed to refine further the water table configuration. This approach d i f f e r s from that of Neuman and Witherspoon (1970) in two respects. F i r s t , f l u i d fluxes calculated at the seepage face are not e x p l i c i t l y used in solving equation (2.7). Rather, freshwater head is spec i f i e d where the water table coincides with the bedrock surface and f l u i d flux i s spec i f i e d only on the free-surface. This minimizes the dependence of the solution on calculated f l u i d fluxes and allows the system of equations to be solved only once for each i t e r a t i o n , rather than the two-step process used by Neuman and Witherspoon (1970). Second, because only v e r t i c a l mesh deformation i s allowed, n e a r - v e r t i c a l topographic slopes are poorly represented in t h i s modified formulation. In developing t h i s approach, i t i s assumed that available i n f i l t r a t i o n i s controlled by processes acting at the bedrock surface (evapotranspiration, subsurface stormflow, surface runoff) rather than by the position of the water table or by processes acting in the unsaturated zone. Therefore, when varying parameters such as permeability in a numerical s e n s i t i v i t y analysis, the available i n f i l t r a t i o n rate i s assumed to be unaffected by varying depth to water table. Where the flow system accepts the entire available i n f i l t r a t i o n , the water table w i l l be predicted to l i e below the bedrock surface and recharge to the saturated region of flow i s exactly equal to the available i n f i l t r a t i o n rate. Where only a portion of the available i n f i l t r a t i o n rate i s accepted, the predicted water table w i l l coincide with the bedrock surface and recharge to the flow system w i l l be less than the available i n f i l t r a t i o n rate. Solutions to the thermal transport equation (equation 2.11) may be poorly approximated by conventional f i n i t e element techniques when advective heat transfer dominates conduction. Such conditions are defined on the basis of a gri d Peclet number that describes, on a l o c a l scale, the r a t i o of heat transfer by advection to that by conduction: ( 2 . 1 7 ) where L£ i s a c h a r a c t e r i s t i c length for an individual element. Patankar (1980) and Huyakorn and Pinder (1983) report that r e l i a b l e solutions to both one-dimensional and two-dimensional transport problems are possible when Pe ^ 2. Although adequate mesh d i s c r e t i z a t i o n can ensure r e l i a b l e solutions, a large number of nodes may be required. In some cases however, upstream weighting procedures may reduce the number of nodes that are required. An acceptable triangular element mesh i s constructed by ensuring that the gr i d Peclet numbers are less than 2 for each element. In representing fracture zones as l i n e elements within the triangular mesh, the fracture permeability must exceed that of the adjacent rock mass by at least four orders of magnitude (Baca et al. , 1984). Such large permeability contrasts imply that large f l u i d flux contrasts are possible. In such cases, solutions with triangular element Pe's less than 2 may also have l i n e element Pe's much greater than 2. This problem is handled by defining an exponential basis function, the form of which depends upon Pe at the previous i t e r a t i o n , for incorporation in the f i n i t e element equations and for c a l c u l a t i n g average f l u i d properties within each l i n e element. The method adopted i s similar to that described for f i n i t e difference grids by Spalding (1972), Raithby and Torrance (1974) and Patankar (1980). Patankar (1980) states that this approach provides good comparisons to one-dimensional a n a l y t i c a l solutions for any value of Pe and any number of grid points. The i t e r a t i v e procedure i s terminated when water table elevations are resolved within a tolerance of 1 m and calculated recharge i s less than available i n f i l t r a t i o n near the POD. In flow systems with topographic r e l i e f of 2 km, th i s tolerance of 1 m i s 0.05 percent of the maximum possible water table r e l i e f . Once the solution i s terminated, a f i n a l check is performed to ensure that nodal temperatures obtained at f i n a l and penultimate it e r a t i o n s d i f f e r by no more than 0.1 °C. 2.5 SUMMARY OF ASSUMPTIONS In developing the numerical method described in the previous section, a number of assumptions are invoked. Those of deemed to be of greatest importance are summarized below: 1. Two-dimensional coupled f l u i d flow and heat transfer in r i g i d v e r t i c a l sections with planar or axisymmetric symmetry. 2. Steady f l u i d flow and heat transfer. Diurnal and seasonal transients are neglected and a "dynamic" steady state i s assumed. The impact of long-term climate change or ongoing igneous a c t i v i t y are not considered. 41 3. An equivalent porous medium i s assumed except where major through-going fracture zones are represented as discrete e n t i t i e s with permeability at least 10" times that of the surrounding rock mass. Permeability and thermal conductivity may be heterogeneous and anisotropic. 4. Thermal equilibrium between f l u i d , s o l i d and vapor. 5. V e r t i c a l boundaries for f l u i d flow and heat transfer are symmetry boundaries (impermeable and insulated). 6. The basal boundary i s horizontal and impermeable with a conductive heat flux applied along the boundary. 7. The upper boundary of the domain i s the bedrock surface where temperatures are s p e c i f i e d using a thermal lapse rate. 8. One-dimensional v e r t i c a l flow of f l u i d through the unsaturated zone causes one-dimensional advective heat transfer (combined with two-dimensional heat conduction) above the water table. V e r t i c a l f l u i d flux above the water table equals the available i n f i l t r a t i o n rate. Heat and moisture transport by vapor movement in the unsaturated zone are assumed n e g l i g i b l e . 9. F l u i d density and v i s c o s i t y vary as a function of temperature and pressure while thermal conductivity and s p e c i f i c heat capacity of the f l u i d are assumed constant. 2.6 IMPLEMENTATION OF THE MODEL The numerical model described in the preceding section is used in Chapter 3 to examine the factors c o n t r o l l i n g groundwater flow through a mountain massif. To demonstrate the basic form of the results, an example i s presented here that characterizes the thermal regime in mountainous t e r r a i n . Simulations are performed by assigning a set of f l u i d flow and thermal parameters (Table 2.1) within a geometry similar to that of Figure 2.1. The importance of surface topography i s examined by considering two extremes in slope p r o f i l e ; one convex and one concave (Figure 2.3). Convex p r o f i l e s are t y p i c a l of glaciated c r y s t a l l i n e t e r r a i n while concave p r o f i l e s are often found in folded mountain belts and at volcanic cones. A basal zone of low permeability occupies the lower 2 km of the system to provide a region of conduction-dominated heat transfer. The remainder of the system is occupied by a higher permeability unit where advective heat transfer may dominate. This configuration allows advective thermal disturbances in the upper unit to propagate into the basal conductive regime ensuring a reasonable t r a n s i t i o n between advection-dominated and conduction-dominated thermal regimes. This region i s s u f f i c i e n t l y thick that the majority of simulation results show isotherms near the basal boundary to be sub-parallel to the boundary. As a consequence, the v e r t i c a l conductive heat flux applied at the basal boundary i s transferred in a consistent manner to the thermal regime. Both upper and lower zones have homogeneous and isotropic permeability (ku, k^) and uniform porosity (nu, n^). Rock thermal conductivity (X s) is uniform throughout the system, however, varying porosity and saturation produce contrasts in thermal conductivity for the s o l i d - f l u i d composite ( X e ) . Because we consider the steady state problem, porosity only has an indirect influence on the flow system through i t s impact on thermal conductivity of the s o l i d - f l u i d composite. Porosity can be expected to play an important role in the transient flow of groundwater in mountainous t e r r a i n , e s p e c i a l l y in fractured c r y s t a l l i n e rock. A uniform heat flow (Hb) i s applied at the base of the flow system and surface temperature conditions are defined in terms of a reference surface temperature (Tr) and a thermal lapse rate (G[). Longitudinal and transverse thermal d i s p e r s i v i t i e s {a,, a() are uniform throughout the system and held constant for a l l simulations. Temperature f i e l d s , r e f l e c t i n g the interaction between advective and conductive heat transfer, are influenced by both f l u i d flow and thermal parameters. Figure 2.3 shows temperature f i e l d s obtained using the values of thermal parameters specified in Table 2.1 and values of / required to maintain the same water table configurations within each system as ku varies from 10" 1 8 to 10" 1 5 m2 (Table 2.2). By maintaining the same water table position, the d i s t r i b u t i o n of thermal conductivity i s unchanged and temperature f i e l d s shown in Figure 2.3 d i f f e r only in the way that advective heat transfer, as influenced by f l u i d flow parameters alone, affects the thermal regime. Pathlines represent the track of a p a r t i c l e entering the flow system at a spec i f i e d point on the bedrock surface. Pathline spacing i s inversely proportional to the flux of f l u i d ( s p e c i f i c discharge) through flowtubes bounded by each pair of pathlines. Maintaining the same water table position produces patterns of f l u i d flow that are v i r t u a l l y i d e n t i c a l for a given slope p r o f i l e , despite a wide v a r i a t i o n in the temperature f i e l d s . Despite their obvious s i m i l i a r i t y , pathlines should not be confused with the streamlines generated from contour plots of a suitably defined stream function or vel o c i t y p o t e n t i a l . These approaches d i f f e r because f l u i d flow through each streamtube (bounded by each pair of streamlines) i s fixed throughout the domain while f l u i d flow through flowtubes defined on the basis of pathlines may d i f f e r from flowtube to flowtube. Flowtubes carry the same f l u i d flow only where they intersect the uniform.fluid flux boundary (the free-surface) with equal spacing. Water table configurations are maintained by increasing / by the same increment that ku i s increased. Maintaining the water table position as permeability i s increased implies that increasingly humid climates support the increased i n f i l t r a t i o n rates. The relat i o n s h i p between i n f i l t r a t i o n rate, permeability and water table elevation i s examined in d e t a i l in Chapter 3. Figures 2.3a to 2.3c show the t r a n s i t i o n from conduction-dominated to advection-dominated thermal regimes as a function of increasing permeability. Upper zone permeabilities (ku) less than 10" 1 8 m2 produce a purely conductive thermal regime with isotherms sub-parallel to the surface topography (Runs C1 and C3 in Figure 2.3). An increase of two orders of magnitude in ku to 10" 1 6 m2 produces a weak advective disturbance as evidenced by the warping of isotherms (Runs C2 and C4 in Figure 2.3). An additional one order of magnitude increase in k u , to 10" 1 5 m2, creates a strongly disturbed thermal regime (Runs B2 and B7 in Figure 2.3). Smith and Chapman (1983) found that groundwater flux i S r s u f f i c i e n t l y large to cause a t r a n s i t i o n from conduction-dominated heat transfer to a weakly disturbed thermal regime when permeability i s about 10" 1 7 m2 in low-relief (1 km over 40 km) sedimentary basins. In regions of hi g h - r e l i e f mountainous topography, f l u i d flux i s s u f f i c i e n t l y enhanced by the greater hydraulic gradients to cause th i s t r a n s i t i o n to occur at a lower permeability of about 10" 1 8 m2. The maximum temperature attained along the base of each flow system (T ) is tabulated in Table 2.2 for each run shown in Figure 2.3. As kti increases, T„„v decreases from 3 u max about 130 °C to about 85 °C. Higher permeabilities enhance groundwater flow, cool the subsurface and reduce conductive thermal gradients in regions of strong downward f l u i d flow. Reduced thermal gradients r e f l e c t the a b i l i t y of groundwater to absorb and re d i s t r i b u t e heat that would otherwise cause increased rock temperatures. The temperature f i e l d s shown in Figure 2.3 i l l u s t r a t e the range of thermal conditions that are encountered in the s e n s i t i v i t y analyses c a r r i e d out in Chapter 3. A more detailed examination of the nature of thermal regimes in mountainous t e r r a i n i s presented in Chapter 4. 2.7 SUMMARY 1. Groundwater flow systems in mountainous t e r r a i n d i f f e r from those in low-relief t e r r a i n in two key respects: (1) for a given set of hydrogeologic conditions, a greater range in water table elevation and form i s possible; (2) hi g h - r e l i e f t e r r a i n enhances groundwater c i r c u l a t i o n to depths where s i g n i f i c a n t heating can occur, implying that thermal effects influence the patterns and rates of groundwater flow. 2. A conceptual model has been outlined to describe groundwater flow systems and thermal regimes in mountainous t e r r a i n . The problem i s viewed as an overlay of two boundary value problems, one for f l u i d flow and one for heat transfer. The uppermost boundary of the domain i s the bedrock surface where temperatures are assumed to be a few degrees warmer than mean annual a i r temperature. Where the water table coincides with the bedrock surface, freshwater head equals the bedrock elevation. Where the water table l i e s below the bedrock surface, a uniform available i n f i l t r a t i o n rate is applied that represents one-dimensional f l u i d flow through the to the water table by one-dimensional flow through the unsaturated zone. Heat supplied by a regional heat flux i s transferred through the system by advection and conduction in both saturated and unsaturated regions of flow. F l u i d flow and heat transfer through a thin cover of discontinuous s u r f i c i a l deposits on upland areas of mountain slopes i s not included in the model. F l u i d moving within these deposits i s lumped with overland flow as a runoff term. The remaining f l u i d available for recharge i s termed an available i n f i l t r a t i o n rate which, in the absence of detailed f i e l d data, i s best thought of as a percentage of the mean annual p r e c i p i t a t i o n rate. 3. Tra d i t i o n a l free-surface approaches incorrectly assume that recharge cannot occur downslope of the point where the water table meets the bedrock surface (commonly defined as the exit point marking the upper l i m i t of the seepage face). In developing the approach used in this study, i t was necessary to identi f y two separate points on the upper boundary; the point where the water table meets the bedrock .surface (point of detachment or POD) and the point that marks the t r a n s i t i o n between recharge and discharge (hinge point or HP). Between these points the water table coincides with the bedrock surface. In t h i s region, the boundary condition i s sp e c i f i e d by setting hydraulic head equal to the elevation of the upper boundary and allowing recharge to occur at a rate defined by the nature of the groundwater flow system. In steep t e r r a i n , with rocks of low permeability, recharge can occur downslope of the POD in a region that might normally be considered part of the seepage face. In regions of reduced topographic r e l i e f , the separation between HP and POD becomes smaller and these points merge at the more commonly defined exit point. 4. The f l u i d and thermal regimes are modeled using a Galerkin f i n i t e element technique in conjunction with a free-surface approach to estimate the position of the water table. Two f i n i t e element grids are required in solving t h i s problem. The mesh for f l u i d flow i s generated only within the saturated zone. The mesh for heat transfer extends from the basal boundary to the bedrock surface, incorporating both saturated and unsaturated regions of flow. Coupling between heat transfer and f l u i d flow i s f a c i l i t a t e d by ensuring a one-to-one correspondence between nodes located below the water table. 5. Inaccuracies inherent in c a l c u l a t i n g boundary fluxes while solving the free-surface problem are minimized in a two-step procedure. The isothermal problem i s f i r s t solved for freshwater head and to obtain the water table configuration. The problem i s then reformulated using stream functions. Large differences between boundary fluxes computed using the two solution methods often indicate regions where the water table configuration i s poorly represented and the g r i d needs refinement. A further check on gr i d accuracy i s supplied by establishing that matching hinge point positions are computed using each solution method. 6. The numerical method developed in this chapter provides the means to examine the factors that control groundwater flow and heat transfer in mountainous t e r r a i n ; geology, surface topography, climate and regional heat flow. Simulation results indicate that each factor influences both advective heat transfer and the position of the water table. Adopting the free-surface approach i s advantageous in performing s e n s i t i v i t y analyses because recharge to the flow system i s controlled as an input variable rather than i m p l i c i t l y calculated in the solution procedure. REFERENCES Baca, R.G., R.C. Arnett, and D.W. Langford, Modelling f l u i d flow in fractured-porous rock masses by f i n i t e element techniques, International Journal for Numerical Methods in F l u i d s , 4, 337-348, 1984. Barry, R.G., Mountain weather and climate, Methuen, 1981. Bear, J . , Dynamics of f l u i d s in porous media, American El s e v i e r , 1972. Bear, J . , Hydraulics of Groundwater, McGraw-Hill, 1979. Birch, F., Flow of heat in the Front Range, Colorado, B u l l . Geol. Soc. Amer., 61, 567-630, 1950. Bodvarsson, G.S., and K. Pruess, Modeling studies of geothermal systems with a free water surface, Stanford Ninth Annual Workshop on Geothermal Resources Engineering, 1983. Chapman, D.S., and L. Rybach, Heatflow anomalies and the i r interpretation, Jour, of Geodynamics, 4, 3-37, 1985. Freeze, R.A., and J.A. Cherry, Groundwater, Prentice-Hall Inc., 1979. Freeze, R.A., and P.A. Witherspoon, Theoretical analysis of regional groundwater flow: 2. Effect of water-table configuration and subsurface permeability v a r i a t i o n , Water Resour. Res., 3, 623-634, 1967. Frind, E.O., Simulations of long-term transient density-dependant transport in groundwater, Adv. Water Resourc, 5, 73-88, 1982. Home, R.N., and M.J. O'Sullivan, Numerical modeling of a desaturating geothermal reservoir, Numerical Heat Transfer, 1, 203-216, 1978. Huyakorn, P.S., and G.F. Pinder, Computational methods in subsurface flow, Academic Press, 1983. 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. Keenan, J.H., F.G. Keyes, P.G. H i l l , and J.G. Moore, Steam Tables, 162 pp., John Wiley, 1978. Neuman, S.P., and P.A. Witherspoon, F i n i t e element method of analyzing steady seepage with a free-surface, Water Resour. Res., 6(3), 889-897, 1970. Patankar, S.V., Numerical heat transfer and f l u i d flow, McGraw H i l l , 1980. Raithby, G.D., and K.E. Torrance, Upstream weighting differencing schemes and their application to e l l i p t i c problems involving f l u i d flow, Comput. Fluids, 2, 191-206, 1974. Ross, B., A conceptual model of deep unsaturated zones with negligible recharge, Water Resour. Res., 20(11), 1627-1629, 1984. Sauty, J.P., A.C. Gringarten, H. Fabris, D. Thiery, A. Menjoz, and P.A. Landel, Sensible energy storage in aquifers, 2, F i e l d experiments and comparisons with theoretical r e s u l t s , Water Resour. Res., 18(2), 245-252, 1982. S c h a r l i , U., and L. Rybach, On the thermal conductivity of low-porosity c r y s t a l l i n e rocks, Tectonophysics, 103, 307-313, 1984. Smith, L., and D.S. Chapman, On the thermal effects of groundwater flow 1. Regional scale systems, Jour. Geop. Res., 88(B1), 593-608, 1983. Smith, M.W., Microclimatic influences on ground temperatures and permafrost d i s t r i b u t i o n , Mackenzie Delta, Northwest T e r r i t o r i e s , Canadian Journal of Earth Science, 12, 1421-1438, 1975. Sorey, M.L. Numerical modeling of l i q u i d geothermal systems, U.S. Geol. Surv. Prof. Paper 1044-D, 1978. Spalding, D.B., A novel f i n i t e difference formulation for d i f f e r e n t i a l equations involving both f i r s t and second derivatives Int. J. Numer. Methods Eng., 4, 551-559, 1972. V e r r u i j t , A., Theory of Groundwater Flow, Macmillan Press, 1970. Walsh, J.B., and E.R., Decker, E f f e c t of pressure and saturating f l u i d on the thermal conductivity of compact rock, Jour. Geoph. Res., 71, 3053-3061, 1966. Watson, J.T.R, R.J. Basu, and J.V. Sengers, An improved representative equation for the dynamic v i s c o s i t y of water substance, Jour. Phys. Chem. Ref. Data, 9(3), 1255-1279, 1981. TABLE 2.1 TYPICAL SIMULATION PARAMETERS - CHAPTER 2 Fluid Flow Parameters kb permeability of basal unit 1 .0x10" 2 2 m2 ku permeability of upper unit 1.0x10- 1 5 m2 'z v e r t i c a l i n f i l t r a t i o n rate 2. 0x10"9 m/s Thermal Parameters Hb basal heat flow 60.0 mW/m2 Gl thermal lapse rate 5 K/km Tr reference surface temperature 10 °C porosity of basal unit 0.01 V porosity of upper unit 0.10 \5 s o l i d thermal conductivity 2.50 W/mK \f f l u i d thermal conductivity 0.58 W/mK Xv vapor thermal conductivity 0.024 W/mK Cf s p e c i f i c heat capacity of water 4186.0 J/kgK s saturation above water table 0.0 al longitudinal thermal d i s p e r s i v i t y 100.0 m at transverse thermal d i s p e r s i v i t y 10.0 m TABLE 2.2 INFLUENCE OF PERMEABILITY it, AND INFILTRATION / u z ON MAXIMUM TEMPERATURES IN FLOW SYSTEM Tmax RUN PROFILE k I 7" u z max (m2) (m/sec) (°C) C1 X 1 0 - 1 8 2x10- 1 2 138 C2 X 10-1 6 2x10" 1 0 130 B2 X 10" 1 5 2X10' 9 81 C3 V 10" 1 B 5x10" 1 2 122 C4 V 10- 1 6 5x10" 1 0 117 B7 V 1 0 - 1 5 5x10-' 87 Notes: 1 ) Slope P r o f i l e ; X = convex, V = concave 2) Rock thermal conductivity X5 = 2.5 W/mK 3) Basal heat flow Hb = 60.0 mW/m2 55 2 + < IX X 3 -I u . o " _ l U . Mean Annual Precipitation Available Infiltration 1 V Rate z g < > Hb 0b 1 t t t t I t T • Z £ Basal Heat Flow ' — Free-surface ^ Insulated and impermeable boundary 77777- Impermeable boundary [illjjlj Basal low-permeability unit Figure 2.1. Conceptual model for groundwater flow in mountainous t e r r a i n . »:nJH ' l l" i iB?iiB^iSISiSBiiS^S!ISi!il5ii3iaiSiSiS"5i l ", l "ii"i»"i: ::::,ikvi>,i>"ik i^Bisii9iSiiiiiSii9Qii9^^>S[i?i&ifisiiiieiii9'i."i>,'i:'/>v: : : : :UUMMBBB0B0BSa8*ffiaSW&SBff«WH'jJ:: - 4 0 2 4 6 0 2 4 6 DISTANCE (km) Figure 2.2. Typical f i n i t e element meshes; a. f l u i d flow problem (1137 nodes, 2139 elements), and b. heat transfer problem (1230 nodes, 2201 elements). a. z g r-< > 111 _l HI Run C1 Run C3 2 4 6 DISTANCE (km) b. 65 - 95 --35 --95 -Run C2 .125-Run C4 — 65-DISTANCE (km) Run B2 Run B7 0 2 4 6 DISTANCE (km) Groundwater pathline — 5- Isotherm (°C) Water table Figure 2.3. Temperature f i e l d s and groundwater flow patterns as a function of upper zone permeability ku for fixed water table configurations in concave and convex slope p r o f i l e s ; a. ku = 10" 1 8 m2 {Iz = 2X10' 1 2 m/sec for convex Run C1 and 5x10" 1 2 m/sec for concave Run C3), Run C2 and 5x10" 1 0 and, c. ku = 10" 1 5 m2 (lz = 2X10" 9 m/sec for convex Run B2 and 5x10 _ 9 m/sec for concave Run B7). b. ku = 10" 1 6 m2 (Iz = 2x10" 1 0 m/sec for convex m/sec for concave Run C4) CHAPTER 3 GROUNDWATER FLOW SYSTEMS IN MOUNTAINOUS TERRAIN 3.1 INTRODUCTION Conceptual and numerical models describing f l u i d flow and heat transfer in mountainous t e r r a i n are presented in Chapter 2. This approach d i f f e r s in two respects from those used previously to simulate regional groundwater flow systems. F i r s t , a non-isothermal formulation provides the means to account for heating of groundwater deep within the mountain massif. Secondly, the free-surface approach avoids making a priori assumptions on the location of the water table. This approach i s advantageous because, with sparse f i e l d data, i t i s most appropriate to estimate water table configurations within constraints imposed by i n f i l t r a t i o n rates, surface topography, rock permeabilty, and regional heat flux. Accounting for both conductive and advective heat transfer in the unsaturated zone allows the interaction of f l u i d flow and heat transfer to be examined, at a regional-scale, throughout both saturated and unsaturated regions of flow. In t h i s chapter, the factors c o n t r o l l i n g patterns and magnitudes of groundwater flow in mountainous ter r a i n are examined in d e t a i l . Using the numerical method described in Chapter 2, idealized mountain flow systems are modeled for a range of conditions representative of the Western C o r d i l l e r a in North America. Factors to be considered are surface topography, geology, climate, and regional heat flow. Elements of surface topography include slope p r o f i l e , r e l i e f and three-dimensional form. The implications of the three-dimensional topography of mountainous t e r r a i n are explored by modeling simple planar and axisymmetric mountain forms with the two-dimensional method. Two extremes in slope p r o f i l e are examined; one convex and one concave. Convex p r o f i l e s are t y p i c a l of glaciated c r y s t a l l i n e t e r r a i n while concave p r o f i l e s can be found in folded mountain belts and at volcanic cones. The geologic environment controls s p a t i a l v a r i a b i l i t y of permeability, porosity, and thermal conductivity. Climatic factors influence available i n f i l t r a t i o n , the presence and extent of alpine g l a c i e r s , and surface temperature conditions. Regional heat flow i s a c h a r a c t e r i s t i c of the tectonic environment within which the mountain i s located. Table 3.1 summarizes the simulations car r i e d out to develop a quantitative understanding of the importance of each factor and provides a guide to the figures that i l l u s t r a t e each e f f e c t . In t h i s chapter we focus primarily on the hydrogeologic regime. Readers are referred to Figure 2.3 in Chapter 2 for examples of temperature d i s t r i b u t i o n s t y p i c a l of the systems to be described in t h i s chapter. The primary feature to r e c a l l from th i s diagram i s the t r a n s i t i o n from a conduction-dominated to an advectively-disturbed thermal regime. Higher permeabilities enhance groundwater flow which, in turn, cools the subsurface and reduces conductive thermal gradients in regions of downward flow. In the example presented, maximum temperatures were reduced by as much as 45 °C from those found in a low-permeability, conductive regime. The nature of mountain thermal regimes w i l l be considered in greater d e t a i l in Chapter 4. 3.2 METHOD OF PRESENTATION Before discussing the s e n s i t i v i t y analyses, the method of presenting the simulation r e s u l t s i s outlined. A t y p i c a l set of results is shown in Figure 3.1 for both convex and concave topographic p r o f i l e s . A basal zone of low permeability occupies the lower 2 km of the domain to provide a region of conduction-dominated heat transfer. This region i s s u f f i c i e n t l y thick that the majority of simulation results show isotherms near the basal boundary to be sub-parallel to the boundary. As a consequence, the v e r t i c a l conductive heat flux applied at the basal boundary is transferred in a consistent manner to the thermal regime. The remainder of the system i s occupied by a higher permeability unit where advective heat transfer may dominate.. In most simulations, both upper and lower zones have homogeneous and isotropic permeability (ku, k^) and uniform porosity (nu, n^). Although the thermal conductivity of the s o l i d matrix (X s) is uniform throughout the system, s p a t i a l v a r i a t i o n in porosity and saturation produce contrasts in thermal conductivity of the sol i d - v a p o r - f l u i d composite ( X e ) . Because we consider the steady state problem, porosity only has an indire c t influence on the flow system through i t s impact on thermal conductivity. Porosity can be expected to play an important role in transient mountain flow systems, especially in fractured c r y s t a l l i n e rock. A uniform heat flow {H^) i s applied at the base of the flow system and surface temperature conditions are defined in terms of a reference surface temperature (Tr) and a thermal lapse rate (<J^  ). Longitudinal and transverse thermal d i s p e r s i v i t i e s (a^, a() are uniform throughout the system and held constant for a l l simulations. F l u i d and thermal properties held constant in the simulations are l i s t e d in Table 3.2. Recharge-discharge p r o f i l e s shown in Figure 3.1 depict the mass flux of f l u i d mfl in kg/m2sec, normal to the water table, projected on a horizontal surface. In h i g h - r e l i e f systems with variable water table slope, the areas under recharge and discharge portions of the p r o f i l e do not provide a dire c t estimate of net recharge and discharge for the flow system. The t o t a l flow through each planar flow system (Q) i s calculated within the numerical formulation by summing the recharge normal to a unit width of the upper boundary in units of kg/sec. Total flow (Q) for axisymmetric systems i s calculated by summing the recharge to the upper c i r c u l a r boundary of the ra d i a l system, also with units of kg/sec. The v e r t i c a l sections in Figure 3.1 i l l u s t r a t e water table configurations and patterns of groundwater flow. Pathlines represent the track of a p a r t i c l e entering the flow system at a specified point on the bedrock surface. Pathline spacing i s inversely proportional to the flux of f l u i d ( s p e c i f i c discharge) through flowtubes bounded by each pair of pathlines. Despite their obvious s i m i l i a r i t y , pathlines should not be confused with the streamlines generated from contour plots of a suitably defined stream function or v e l o c i t y p o t e n t i a l . These approaches d i f f e r because flow of f l u i d through each streamtube (bounded by each pair of streamlines) i s fixed throughout the domain while f l u i d flow through flowtubes defined on the basis of pathlines may d i f f e r from flowtube to flowtube. Flowtubes carry the same f l u i d flow only where they intersect the uniform f l u i d flux boundary (the free-surface) with equal spacing. Because a uniform available i n f i l t r a t i o n rate is applied in each case, the magnitude of / provides a convenient reference for comparing the variation of f l u i d f l u x . For example, pathline spacing in Figures 3.1a and 3.1b varies in response to the fact that f l u i d flux in the convex domain increases approximately f i v e f o l d from a value of 2x10 _ 6 kg/sec at the free-surface to about 10"5 kg/sec (on average) at the discharge area. F l u i d flux is more uniform in the concave domain (Figure 3.1b), decreasing from 2x10-* kg/sec at the free-surface to about 1.6x10~6 kg/sec. The r e l a t i v e l y narrow range of f l u i d flux within each domain shown in Figure 3.1 suggests that the available i n f i l t r a t i o n rate can be viewed as a c h a r a c t e r i s t i c s p e c i f i c discharge for the domain under consideration, when the water table l i e s below the bedrock surface. As a consequence, c h a r a c t e r i s t i c t r a n s i t times can be calculated for each flow system where t c = nu/IzLp' Here, Lp i s the length of a spec i f i e d pathline and nu i s the porosity of the upper permeable zone. Transit times for the lowermost pathlines shown in Figure 3.1 are about 8x10 5 years for domains with bulk permeability of 10" 1 5 m2. Simulation results are quantified by defining the following parameters: maximum elevation of water table, usually located below the mountain summit. absolute difference between the value of WTmax obtained in a previous simulation and that obtained in the simulation run of inte r e s t . percentage change in t o t a l flow between a designated reference Q and that obtained in the run. £WTmax and bQ provide a convenient measure of the s e n s i t i v i t y of mountain flow systems to differences in WTmax <m> mTmax <m> AQ (%) permeability, i n f i l t r a t i o n rate, surface topography and thermal regime. For a s p e c i f i e d set of thermal parameters, an i n f i l t r a t i o n r a t i o I* = I„/K„ can be defined where K„ is Z 0 0 a reference hydraulic conductivity: k p e Ko-^TJ- « 3 . 1 ) with f l u i d density pQ and v i s c o s i t y u0 s p e c i f i e d at the reference surface temperature Tr . A wide range in permeability and i n f i l t r a t i o n rates can be explored with a r e l a t i v e l y small number of simulations by recognizing the it a p p l i c a b i l i t y of the dimensionless i n f i l t r a t i o n r a t i o / . Results are presented in Table 3.3 for sixteen simulation series designated A through P. Each series i l l u s t r a t e s the influence of a s p e c i f i e d c o n t r o l l i n g factor such as slope p r o f i l e , i n f i l t r a t i o n rate or permeability. Data tabulated in Table 3.3 (convex Runs C l , C2, B2 and concave Runs C3, C4, B7) indicate that maintaining a constant I yields differences in water table elevations that do not exceed 20 m as ku varies over three orders of magnitude, despite a s i g n i f i c a n t v a r i a t i o n in temperature f i e l d s (see Figure 2.3 in Chapter 2). Moreover, the patterns of f l u i d flow are e s s e n t i a l l y unchanged as ku varies. This result r e f l e c t s an approximate linear r elationship between Iz and ku for a given set of thermal parameters. Results shown in Table 3.3 indicate that increasing permeability and i n f i l t r a t i o n by a sp e c i f i e d increment causes t o t a l flow Q to increase by the same increment (convex Runs C1, C2, B2 and concave Runs C3, C4, B7). This * suggests that a dimensionless t o t a l flow Q can be defined * in a manner similar to that used to define I ; where / i s normalized with respect to KQ (equation 3.1). Dimensionless * t o t a l flow Q is defined as follows: Q* = Q (3.2) LpnKW O 0 where L i s the length of the section and If i s a unit width. Simulation results indicate that unique values of Q and I are found for each water table position within a specified system geometry. For example, a maximum water table elevation of about 1550 m i s found in the convex domain when I equals 0.26 (Runs C1, C2 and B2 in Table 3) and Q equals 0.22 although both permeability and i n f i l t r a t i o n vary over three orders of magnitude in the simulation series. Although it dimensionless t o t a l flow Q lacks an obvious physical meaning, i t i s a useful parameter for viewing the r e l a t i v e changes in t o t a l flow calculated in subsequent numerical simulations. In developing dimensionless parameters to characterize the simulation results a dimensionless water table elevation was considered but not implemented. Such a parameter might be useful in assessing the influence of topographic r e l i e f . Examining the absolute magnitude of changes in water table elevation, however, yields greater insight into the influence of factors c o n t r o l l i n g groundwater flow in mountainous t e r r a i n . 3.3 FACTORS CONTROLLING GROUNDWATER FLOW IN MOUNTAINS 3.3.1 Slope P r o f i l e Examining simulation results obtained for simple topographic slope p r o f i l e s provides insight into flow systems found in complex mountain topographies. Simulation results for a system with convex topographic p r o f i l e (Run A1) and a concave p r o f i l e (Run A2) are shown in Figure 3.1 and tabulated in Table 3.3. The recharge-discharge p r o f i l e s of Figure 3.1 indicate the strong influence of slope p r o f i l e on the pattern and magnitude of f l u i d flux across the bedrock surface. In convex topography, the recharge area i s more extensive than the discharge area and discharge fluxes can greatly exceed recharge fluxes (Figure 3.1a). In the concave topography recharge and discharge areas are approximately equal and discharge fluxes are similar to recharge fluxes (Figure 3.1b). A key factor in determining the influence of slope p r o f i l e i s the volume of rock within the mountain massif that acts to r e s i s t f l u i d flow. For example, convex p r o f i l e s contain a greater volume of porous rock above the valley f l o o r than concave p r o f i l e s (Figure 3.1). Therefore greater mechanical energy, reflected by higher water table elevations, i s required to drive groundwater flow through convex slopes. The higher water table elevations found in the convex slope, in turn, contribute to a more extensive region of recharge. As a consequence, applying the same available i n f i l t r a t i o n rate in each case causes a t o t a l flow through the convex system twice that of the concave system (Table 3.3). Differences between water table elevations calculated for each slope p r o f i l e highlight the advantage of using a free-surface approach to estimate water table configurations that are consistent with assumed available i n f i l t r a t i o n rates, thermo-hydrologic properties of the region and surface topography. 3.3.2 I n f i l t r a t i o n Ratio Available i n f i l t r a t i o n rate (7 Z) and the permeability of the mountain massif (*u) influence water table elevations and rates of f l u i d flow through the i n f i l t r a t i o n r a t i o (/ ). Prior to describing the influence of the i n f i l t r a t i o n r a t i o on flow systems, reasonable l i m i t s are defined for i n f i l t r a t i o n rates and permeability in mountainous t e r r a i n . Mean annual p r e c i p i t a t i o n i s the basic determinant of available i n f i l t r a t i o n rates. Review of climate summaries (Bryson and Hare, 1974; Ruffner and B l a i r , 1975) indicates that p r e c i p i t a t i o n rates at mountain summits of western North America are unlikely to exceed 6 m/year, except perhaps in southwest Alaska. In the low permeability t e r r a i n and humid climate of the Alaska Panhandle, recharge to mountain flow systems i s unlikely to exceed approximately 10 percent of t h i s summit p r e c i p i t a t i o n rate. Therefore, i t seems reasonable to assume an upper l i m i t for available i n f i l t r a t i o n rates of about 0.6 m/yr (2.0 x 10"8 m/sec). In more a r i d climates, with t e r r a i n of higher permeability, i n f i l t r a t i o n rates are unlikely to exceed t h i s upper l i m i t although a greater percentage of the p r e c i p i t a t i o n rate may be available for i n f i l t r a t i o n . Recall that, in developing t h i s conceptual model (Chapter 2), a minimum i n f i l t r a t i o n rate of 10" 1 2 m/sec i s assumed. The reference lz of 2x10" 9 m/sec shown in Table 3.2 i s selected to represent the maximum rate of i n f i l t r a t i o n that might be available in a climate t r a n s i t i o n a l between semi-arid and humid. Values of Iz used in t h i s study vary from 2x10' 1 2 m/sec to 2X10 - 8 m/sec (Table 3.3) to provide a r e a l i s t i c range of i n f i l t r a t i o n rate. Freeze and Cherry (1979) note that rock permeabilities can range from less than 10 _ 2° m2 to in excess of 10 - 9 m2. The bulk permeability of mountainous t e r r a i n considered in th i s study i s assumed to be less than 10 - 1* m2. Freeze and Cherry (1979) and Neuzil (1986) suggest t h i s i s a reasonable upper l i m i t for moderately fractured c r y s t a l l i n e and argillaceous rocks. Simulation results show that a sevenfold increase in * the i n f i l t r a t i o n r a t i o I , from 0.13 to 0.92, causes an increase in water table elevation in excess of 1000 m (concave Runs A2, B5 through B8 and convex Runs A1 , B1 through B4 in Table 3.3) as / increases from 1.0xl0 - 9 to 7.0x10"9 m/s (0.032 to 0.224 m/yr). If t h i s difference in I* were to r e f l e c t only differences in mean annual p r e c i p i t a t i o n , a sevenfold change in / indicates that a substantial variation in climate i s required to e f f e c t a large change in water table elevation. A sevenfold increase in i n f i l t r a t i o n r a t i o could equally well correspond to an equivalent decrease in permeability k u . This, variation i s only a small portion of the possible range of permeability or i n f i l t r a t i o n rates that might be encountered in mountainous t e r r a i n . Therefore, i t is important to obtain reasonable estimates of each factor i f r e l i a b l e predictions of water table elevations are desired. Water table configurations and recharge-discharge p r o f i l e s are shown in Figure 3.2 for three i n f i l t r a t i o n rates that span the sevenfold range in /* discussed in the previous paragraph (convex Runs A1, B1, B4 and concave Runs A2, B5 and B8 in Table 3.3). Pathlines are shown for two sets of paired runs (A1 and B1, A2 and B5) in Figures 3.2a and 3.2b with the dashed l i n e s indicating the runs with lower i n f i l t r a t i o n rate (B1 and B5). These paired runs * r e f l e c t a threefold difference in I from 0.13 to 0.40 that produces differences in water table elevations in excess of 500 m. Flow patterns within each p r o f i l e can be compared by recognizing that individual pathlines drawn for each pair originate at the same point on the bedrock surface. Although a large difference in water table elevation i s found for the paired runs in each system, the overa l l patterns of flow are r e l a t i v e l y unaffected. Note the assumed v e r t i c a l flow above the water table contributes to the deeper pathlines in systems with lower water table elevations. If l a t e r a l f l u i d flow in the unsaturated zone were to be considered, the resulting pathlines might resemble those shown for the cases with higher water table elevation (Figure 3.2). Recharge-discharge p r o f i l e s shown in Figure 3.2 indicate that the pattern and magnitude of mass flux in both recharge and discharge areas is strongly influenced by varying the i n f i l t r a t i o n r a t i o . In each case shown in Figure 3.2, increasing I causes the HP (hinge point marking t r a n s i t i o n between groundwater recharge and discharge) to move upslope on the bedrock surface, a shorter free-surface segment and an increase in rates of groundwater discharge. Furthermore, the difference in the form of recharge-discharge p r o f i l e s obtained for convex and concave cases provides further evidence for the way that the influence of slope p r o f i l e on flow systems i s amplified in mountainous t e r r a i n . * A plot of WT x versus I (Figure 3.3a) demonstrates the v a r i a t i o n in.water table elevation as a function of available i n f i l t r a t i o n rate (or permeability) for both slope p r o f i l e s . The s o l i d c i r c l e s (Curve 1 ) represent results obtained for the concave p r o f i l e while the open c i r c l e s (Curve 2) represent results obtained for the convex p r o f i l e . For a sp e c i f i e d set of thermal parameters and system geometry, a threshold i n f i l t r a t i o n r a t i o can be defined that represents the maximum recharge rate accepted by the flow system. When the threshold r a t i o i s applied, the water table i s everywhere at the bedrock surface and the maximum flow through the system i s attained. In Figure 3.3a, threshold i n f i l t r a t i o n ratios are the values of I* required to raise the water table to i t s maximum possible elevation of 2000 m. * The threshold I is about 1.0 for the concave topography and * about 0.50 for the convex topography. If I i s increased above the threshold, the flow system i s unaffected as * greater runoff must occur. Run B4, with I = 0.92, greatly exceeds the threshold rate for the convex p r o f i l e while Run A1 with / = 0.40 i s s l i g h l t y less than the threshold. Despite the large difference in applied i n f i l t r a t i o n rate (and hence in I ) the water table configurations and recharge-discharge p r o f i l e s found for Runs B4 and A1 are * indistinguishable in Figure 3.2a because I for Run A1 approaches the threshold value. In the companion plot of WTmax v e r s u s dimensionless t o t a l flow Q* (Figure 3.3b), maximum dimensionless flows of about 0.14 and 0.26 occur for concave and convex p r o f i l e s when the threshold i n f i l t r a t i o n r a t i o i s exceeded. Figure 3.3a indicates that water table elevations vary * as a function of / in a non-linear fashion. The overall * slope of the I against WT x plot d i f f e r s for the concave (Curve 1) and convex (Curve 2) cases. The steeper slope of the convex case indicates that water table elevations in convex topography are more sensitive to rock permeability or i n f i l t r a t i o n rate than water table elevations in concave topography. Note that, regardless of slope p r o f i l e , large differences in water table elevation may correspond to r e l a t i v e l y small variations in i n f i l t r a t i o n r a t i o . For example, a 100 m increase in water table elevation from 1600 to 1700 m corresponds to an increase of less, than 20 percent in / for each slope p r o f i l e (Figure 3a). This result suggests that measured water table elevations could aid in constraining estimates of i n f i l t r a t i o n rates and bulk permeability for the mountain massif. Unfortunately, as w i l l be shown in a later section, water table configurations are strongly influenced by the thermal regime. Therefore, a knowledge of the thermal conditions i s also required i f hydraulic properties of mountain flow systems are to be constrained using measured water table elevations. Total flow Q and dimensionless Q for each run are * tabulated in Table 3.3 and Q i s plotted versus WT x in Figure 3.3b. In a concave slope p r o f i l e with high water table elevations, the t o t a l flow through the flow system * (represented by Q ), i s v i r t u a l l y independent of water table elevation and i n f i l t r a t i o n r a t i o . Therefore, modeling studies designed to obtain estimates of t o t a l flow through the bedrock flow system of concave slope p r o f i l e s need not make accurate estimates of water table configurations i f the water table i s located high in the mountain massif. As water table elevation decreases, t o t a l flow i s increasingly * sensitive to changes in W7V „ and hence to I . A similar J max it relationship exists for the convex p r o f i l e , however, Q for high water elevations i s more sensitive to the position of the water table in the mountain massif (Figure 3.3b). 3.3.3 Topographic Relief and Basin Depth The influence of reduced topographic r e l i e f can be assessed by truncating the slope p r o f i l e s described previously at an elevation of 1000 m (Figure 3.4). In t h i s manner topographic r e l i e f i s reduced by a factor of one-half, while the p r o f i l e of the bedrock surface on the lower slope i s maintained. The results of two runs, one for each truncated p r o f i l e it (Runs 02 and 03), are shown in the plots of WT x versus I it and WTmax versus Q (Figures 3.3a and 3.3b). The square symbols represent the truncated concave case (Curve 3) and the truncated convex case (Curve 4). These curves indicate that lower r e l i e f causes an increase in water table elevation and a small decrease in t o t a l flow under conditions that d i f f e r only in topographic r e l i e f from those used to obtain Curves 1 and 2. Performing t h i s numerical exercise under isothermal conditions would show that the response of the truncated topography replicates the response of the o r i g i n a l slope p r o f i l e s . The differences between Curves 1 and 3 and between Curves 2 and 4 shown in Figures 3.3a and 3.3b result from the advective transfer of heat in the unsaturated zone. This contrast in results indicates the importance of accounting for advective heat transfer in thick unsaturated zones when attempting to predict water table elevations in mountainous flow systems. Modifying the geometry of the lower slope has a s i g n i f i c a n t impact on water table configurations and t o t a l flow rates. Curves with triangular symbols (lab e l l e d Curve 5 in Figures 3.3a and 3.3b) represent the influence of modifying the truncated concave p r o f i l e to provide a f u l l y concave p r o f i l e from ridge to v a l l e y . The re s u l t i n g slope p r o f i l e i s sketched in Figure 3.4b and simulation results for Run 03 are tabulated in Table 3.3. Results plotted for the modified concave topography (Curve 5 in Figures 3a and 3b) d i f f e r s i g n i f i c a n t l y from those found for the low-relief system with a f l a t upland area (Curve 3) because the length of the free-surface segment i s reduced in the modified concave topography. As a consequence, t o t a l flow i s s i g n i f i c a n t l y reduced in the low-relief f u l l y concave case and water table elevations are less sensitive to changes in i n f i l t r a t i o n r a t i o . This result highlights the way that h i g h - r e l i e f t e r r a i n amplifies the impact of factors c o n t r o l l i n g groundwater flow. In this case, the factors are permeability and available i n f i l t r a t i o n rate. A decrease in basin depth i s simulated by increasing the thickness of the basal low-permeability zone from 2 km to 3 km while maintaining the external geometry used in the previous concave runs. The corresponding decrease in the thickness of the upper permeable zone reduces the t o t a l flow by 21 percent (Run P1, Table 3.3) from that of Run B7. Despite t h i s s i g n i f i c a n t reduction in Q, the water table elevation increases by only 94 m as basin depth i s reduced by 1 km. 3.3.4 Topographic Symmetry Axisymmetric Topography Simple three-dimensional mountain forms amenable to two-dimensional characterization are axisymmetric conical features and linear ridges. Axisymmetric topography i s representative of isolated mountain massifs elevated above a more uniform low-relief topography. The influence of topographic symmetry is examined by modeling a s p e c i f i e d slope p r o f i l e twice; assuming planar symmetry in the f i r s t case (Run 11 and 13) and r a d i a l symmetry in the second case (Run 12 and 14). An isolated concave mountain located on a f l a t p l a i n i s shown in Figure 3.5 with flow patterns for planar Run 13 and axisymmetric Run 14. The axis of symmetry i s located at the lefthand boundary beneath the ridgetop. Because the surface area available for discharge increases with increasing r a d i a l distance, rates and patterns of f l u i d flow (shown by differences in pathline patterns) within the r a d i a l system d i f f e r from those of the equivalent planar system. If axisymmetric mountain topography i s represented using planar symmetry, water table elevations are overestimated to a s i g n i f i c a n t degree, v e r t i c a l components of f l u i d flux are underestimated throughout the region of flow, and pathlines are displaced to shallower depths in the flow system (Figure 3.5). Although s i g n i f i c a n t differences in recharge rates occur on the mountain flank, only minor differences in discharge rates occur over much of the f l a t p l a i n . Similar results are obtained, but not shown, for a convex mountain adjacent to a f l a t plain (Runs 11 and 12 in Table 3.3). Ridge and Valley Asymmetry If a groundwater flow divide i s assumed to correspond to the crest of an asymmetric ridge, errors in estimating the magnitude of groundwater flow to each valley can r e s u l t . Results obtained for an asymmetric ridge (Run J1 shown in Figure 3.6) exhibit a divide displacement of 1.9 km. The corresponding flow pattern and recharge-discharge p r o f i l e indicate that about one-third of the i n f i l t r a t i o n applied on the free-surface is applied to the right of the ridgecrest yet transmitted to the lefthand v a l l e y . This represents a s i g n i f i c a n t r e d i s t r i b u t i o n of i n f i l t r a t i o n from the righthand to the lefthand watershed. Without using the free-surface approach, i t would be d i f f i c u l t to predict, a priori, the position of the flow divide and the appropriate d i s t r i b u t i o n of i n f i l t r a t i o n to each watershed. If a groundwater divide i s assumed to correspond to the fl o o r of an asymmetric v a l l e y , the inferred source of thermal or chemical properties of groundwater samples obtained in springs and boreholes may be i n c o r r e c t l y located beneath the wrong slope. Run K1, shown in Figure 3.7a, depicts an asymmetric mountain valley with mirrored v e r t i c a l boundaries at the ridgecrests. Opposing ridges have the same r e l i e f , but d i f f e r e n t slope p r o f i l e s . In the. v a l l e y discharge area, the groundwater divide i s displaced about 0.5 km to the right and deviates s l i g h t l y from the v e r t i c a l . Reducing the r e l i e f of the righthand slope to one-half that of the lefthand slope and reducing the angle of the linear slope causes an additional s h i f t of 1 km to the right (Run K2 shown in Figure 3.7b). The d i r e c t i o n of divide displacement r e f l e c t s the r e l a t i v e magnitude of flow through each slope p r o f i l e with displacement occurring away from the slope with maximum flow. It was pointed out e a r l i e r that, for a s p e c i f i e d r e l i e f , the t o t a l flow through a convex p r o f i l e exceeds that of the corresponding linear or concave p r o f i l e . Therefore, in Runs K1 and K2, the divide i s shifted away from the convex ridge towards the linear ridge. Flow patterns shown in Figure 3.7 suggest that groundwater samples c o l l e c t e d above the valley floor at springs issuing from the linear slope could r e f l e c t chemical conditions beneath the opposing convex slope. 3.3.5 Permeable Horizons and Fracture Zones Concave slope p r o f i l e s similar to that shown in Figure 3.8 are common in volcanic t e r r a i n and regions with folded sedimentary units. In such regions, thin permeable horizons may be encountered below the valley f l o o r . Figure 3.8 demonstrates how flow patterns are modified by a 100 m thick horizon, with a permeability 10 times that of the surrounding rock, at a depth of 100 m below the valley f l o o r . Total flow through the system increases by 19 percent while the water table elevation declines by 93 m (Run D4 in Table 3.3). The recharge-discharge p r o f i l e i s r e l a t i v e l y unaffected because the influence of the permeable zone i s dis t r i b u t e d evenly across the section. As permeability or thickness of the horizon is increased, i t s influence i s greater. Simulating a permeable horizon within a convex slope p r o f i l e produces a similar impact on the flow system (Run D4 in Table 3.3). Sub-vertical and dipping fracture zones, however, are more common than permeable horizons in the glaciated c r y s t a l l i n e terrain associated with convex topographic p r o f i l e s . The response of water table elevation and t o t a l flow to a v e r t i c a l fault zone outcropping at the valley floor in a convex topography i s summarized in Table 3.3 (Runs E1 and E2). A 100 m wide fault zone is modeled in Run E1 with a permeability 10 times that of the rock mass k u . In Run E2, a 0.1 m wide fault zone is simulated with a permeability 10* times that of the surrounding rock mass. In both cases, the fa u l t zone i s modeled with the transmissivity of the fracture kfb equal to ^03^ku in units of m 2«m. Note that t h i s d e f i n i t i o n of transmissivity d i f f e r s from the product of hydraulic conductivity times thickness usually adopted in isothermal approaches. Run E1 was carr i e d out using triangular elements to represent the fracture zone, while l i n e elements were used in run E2. The resulting influence on water table configurations and t o t a l flow i s similar for both runs, causing water table declines of about 270 m and changes in t o t a l flow of about 10 percent, r e l a t i v e to a reference case without a v e r t i c a l fracture zone. Figure 3.9 demonstrates the influence of a steeply dipping fracture zone which outcrops at the floor of an asymmetric valley (Run E4). Simulation results for the case without the fracture zone are shown in Figure 3.7a (Run K1). The fracture zone i s modeled with a value of kfb equal to 10 3«& U (m 2«m). The presence of the permeable zone causes a 27 percent increase in t o t a l flow. The recharge-discharge p r o f i l e shown in Figure 3.9 indicates a s i g n i f i c a n t reduction in discharge from the surrounding rock mass because the fracture zone captures about 50 percent of the groundwater flowing through the domain. It i s of interest to note that the permeable fracture zone also influences the position of the valley groundwater divide and compensates for the divide displacement caused by the asymmetric topography (Figure 3.7a). A tenfold increase in fracture zone transmissivity to 1 0 * • /: w (m 2«m) causes additional displacement of the divide with the result that the divide i s found at the location of the fracture zone. Furthermore, th i s increased permeability causes a 75 percent increase in t o t a l flow above that of the unfractured case (Run K1 in Table 3.3) and enables the fracture zone to capture a greater portion of flow through the domain (about 85 percent). This y i e l d s a s i g n i f i c a n t impact on both the patterns and magnitudes of groundwater flow. Further increasing the transmissivity of the fracture to ^05^ku (m 2«m), however, has l i t t l e additional influence on the groundwater flow system. As a consequence, below a threshold value of kfb equal to about 10* (m2-m) the fracture zone r e s t r i c t s f l u i d flow in the surrounding rock mass. At values of kfb in excess of this value, f l u i d flux in the system i s r e s t r i c t e d only by the permeability of the surrounding rock mass. Thermal aspects of spring discharges from through-going fracture zones are discussed in Chapter 4. Permeable horizons and fau l t zones influence mountain flow systems by transmitting a reduced f l u i d potential deep within the mountain massif. This reduced potential causes a re d i s t r i b u t i o n of f l u i d flow allowing the permeable zone to capture a portion of flow. The maximum influence of permeable zones occurs in convex slope p r o f i l e s where a greater portion of the flow system l i e s above the valley floor and i s more susceptible to the impact of permeable horizons. Sub-horizontal features outcropping at the valley floor have the greatest impact on mountain flow systems while increasing angle of dip and increasing distance between the valley floor and the outcrop of the permeable zone reduces the influence of the zone on the flow system. The influence of a permeable horizon on the flow system i s also reduced as depth to the horizon increases (Runs D1 and D3 in Table 3.3). In situations where fracture density i s large r e l a t i v e to the scale of a given flow system, an equivalent anisotropic porous medium may be representative of the fracture network. Run F1 (Table 3.3) models a sub-vertical fracture network using an equivalent anisotropic permeability with the v e r t i c a l permeability enhanced by a r a t i o of 5:1 above the horizontal permeability. This r e l a t i v e l y small increase in v e r t i c a l permeability causes a s i g n i f i c a n t decline in water table elevation (444 m) and a small increase in t o t a l flow (10 percent). Although anisotropic permeability i s more l i k e l y the rule rather than the exception, i t is d i f f i c u l t to measure v e r t i c a l permeability in f i e l d settings. 3.3.6 G l a c i e r s G l a c i e r s occupying upper mountain slopes w i l l modify i n f i l t r a t i o n rates and thus, groundwater flow systems. In temperate mid-latitude climates, many gl a c i e r s can be assumed to have basal temperatures at the melting point of water (Paterson 1981), thus, permafrost conditions are absent and groundwater flow can occur at the base of the g l a c i e r . The influence of a glacier on i n f i l t r a t i o n rates i s simulated by computing sub-glacial recharge during the solution procedure. In this procedure, recharge i s assumed to result only from ice melted at the bedrock surface by conductive heat flow. Using the latent heat of fusion for ice and the computed magnitude of conductive heat flux at the base of the g l a c i e r , i n f i l t r a t i o n rates are calculated as the volume of ice melted per unit time. Figure 3.10 compares the recharge-discharge p r o f i l e and water table configuration obtained with a g l a c i e r located between an elevation of 1500 and 2000 m on a concave slope (Run G2 in Table 3.3) and the results obtained for an equivalent system without the influence of the g l a c i e r (Run M6 in Table 3.3). The flow pattern for Run M6 i s similar to those shown for the concave p r o f i l e s of Figure 2.3. In each case, a uniform Iz of 5.0x10"9 m/s i s applied downslope of the g l a c i e r . In Run M6 t h i s uniform rate i s also applied in the region where the glacier is shown. In Run G2 a reduced / beneath the gl a c i e r i s computed with a minimum value of 7.0x10' 1 1 m/sec. This reduced i n f i l t r a t i o n rate d i s t o r t s the flow system, modifies the recharge-discharge p r o f i l e (Figure 3.10) and reduces the t o t a l flow (Run G2 in Table 3.3). The humped form of the water table r e f l e c t s the fact that recharge is n e g l i g i b l e beneath the g l a c i e r , therefore, a portion of the recharge applied downslope of the glacier margin moves beneath the glacier before c i r c u l a t i n g at depth to exit at the valley f l o o r . In the corresponding convex case, a glacier located at the same elevation covers a greater portion of the recharge area, and thus there i s a greater impact on t o t a l flow and water table elevation (Run G1 in Table 3.3). It should be noted that similar results would be obtained i f groundwater recharge i s r e s t r i c t e d by extensive units of low-permeability rock found in upper regions of the domain. Runs G1 and G2 have an upper zone permeability ku equal to 1 0 - 1 5 m2. At lower rock mass permeabilities, s l i g h t l y greater i n f i l t r a t i o n rates are calculated because reduced f l u i d flux allows higher rock temperatures and enhanced melting of g l a c i a l i c e . The combination of higher i n f i l t r a t i o n rate and lower permeability cause an increase it in i n f i l t r a t i o n r a t i o I . Therefore, at lower rock permeabilities, water table elevations approach the bedrock surface and flow systems are less affected by the g l a c i e r . it The calculated i n f i l t r a t i o n rate y i e l d s an I at the threshold r a t i o (and a water table configuration everywhere at the bedrock surface) when ku is less than about 1 0 - 1 7 m2. Simulation results suggest that i n f i l t r a t i o n rates and t o t a l flow may be s i g n i f i c a n t l y influenced by glaciers covering recharge areas of mountain flow systems. In more northerly and inland regions of the Western C o r d i l l e r a (above a latitude of 50 degrees), the l i k e l i h o o d of finding temperate glaciers i s reduced and permafrost conditions may be found beneath gl a c i e r s and on slopes that are normally snow- or ice-free in winter. Under such conditions, recharge to groundwater flow systems may be further reduced because the permafrost acts as an impermeable blanket that may produce l i t t l e f l u i d by melting. 3.3.7 Non-uniform I n f i l t r a t i o n I n f i l t r a t i o n patterns are controlled by a number of factors including slope angle, p r e c i p i t a t i o n patterns, evapotranspiration and s o i l permeability. The interaction of these factors is complex and few data are available for representative mountain slopes (Barry, 1981). Two factors can be evaluated using simple techniques; slope angle and orographically-controlled p r e c i p i t a t i o n . Steep slopes reduce available i n f i l t r a t i o n rates by enhancing surface runoff. An a r b i t r a r y non-uniform i n f i l t r a t i o n rate i s incorporated in Runs L1 and L2 (Table 3.3) by multiplying / times the cosine of the angle of the bedrock surface 6 from horizontal. Modifying Iz using t h i s approach causes a subtle increase in water table slope that modifies patterns of f l u i d flow to a minor degree near the water table and produces less than 3 percent change in t o t a l flow Q (Runs L1 and L2 in Table 3.3). Barry (1981) summarizes the results of Lauscher (1976) who shows that orography causes p r e c i p i t a t i o n in mountainous terr a i n to vary with a l t i t u d e in a manner that depends upon climate. P r e c i p i t a t i o n in t r o p i c a l climates increases with elevation to a maximum at about 1 to 1.5 km above the valley f l o o r , and decreases with increasing elevation. In polar climates, p r e c i p i t a t i o n rates tend to decrease s l i g h t l y with increasing elevation. In t h i s study, mid-latitude climates are of interest and, following the example of Lauscher, p r e c i p i t a t i o n is assumed to increase with increasing a l t i t u d e . Few data are available to describe p r e c i p i t a t i o n gradients on mountain slopes. Moreover, co r r e l a t i o n between pr e c i p i t a t i o n and available i n f i l t r a t i o n rates i s lacking at the mountain scale. Review of climate data (Bryson and Hare, 1974; Ruffner and B l a i r , 1975) and reported a l t i t u d i n a l gradients in p r e c i p i t a t i o n (Schermerhorn, 1967; Storr and Ferguson, 1972; Slaymaker and Zeman, 1975) suggest that a threefold increase in p r e c i p i t a t i o n rate for 2000 m of elevation gain i s not unusual. In t h i s study, a linear i n f i l t r a t i o n gradient i s assumed to r e f l e c t a threefold increase in p r e c i p i t a t i o n rate. The i n f i l t r a t i o n gradient is applied to both convex and concave p r o f i l e s using the maximum lz of 2x10 - 8 m/s at the ridgecrest and ku of I x l O - 1 " m2. Results for Runs L3 through L6, tabulated in Table 3.3, indicate that these conditions cause d e f i n i t e differences in predicted water table elevations and a minor difference in t o t a l flow. Slope effects and orographic e f f e c t s are defined with reference to variation in surface topography, therefore, th e i r maximum influence is f e l t when water table elevations are low and there is a substantial variation in surface topography upslope of the point of detachment (POD). Therefore, as ku i s increased or I i s reduced (more a r i d climate with r e l a t i v e l y high permeability rock), water table elevations are reduced and the influence of these e f f e c t s increases. 3.3.8 Basal Heat Flow Increased basal heat flow causes a warming of conduction-dominated and f l u i d upflow regions, an increase in t o t a l flow and a decline in water table elevation (Figure 3.11). Patterns of f l u i d flow, however, are affected only to a minor degree. Simulation results for equal to 30, 60 and 120 mW/m2 are tabulated in Table 3.3 for the convex p r o f i l e (Runs H1, B2 and H2) and for the concave p r o f i l e (Runs H3, B7 and H4). Flow patterns and temperature f i e l d s for basal heat flows of 60 and 120 mW/m2 are shown for convex Runs B2 and H2 in Figure 11. F l u i d flow patterns are similar despite changes in water table elevation and H^. A small downward s h i f t in pathline position i s found in the greater heat flow case because f l u i d flux i s enhanced in regions of warming (Figure 3.11). 87 3.3.9 Thermal Conductivity Changes in thermal conductivity of the s o l i d matrix Xs influence water table elevations and flow rates in a manner similar to that described for variations in basal heat flow, but to a lesser degree. Reduced conductive heat transfer occurs at lower \ 5 to cause heating of the thermal regime in regions of f l u i d upflow and conduction dominated heat transfer. For example, Runs N3 and N4 (Table 3.3) indicate that increasing Xs from 2.5 to 3.5 W/mK causes increased water table elevations and a minor increase in t o t a l flow. Therefore, changes in thermal conductivity have a d e f i n i t e , but r e l a t i v e l y small, impact on mountain flow systems. 3.3.10 Surface Temperature Review of climate data summaries indicates a range in mean annual valley temperatures from about 5 °C to 15 °C in western North America. Incorporating t h i s range of valley temperatures in a series of simulations y i e l d s s i g n i f i c a n t differences in water table elevations and minor differences in t o t a l flow. As T i s increased, upper regions of flow are warmed (Runs M5 through M8; Table 3.3) causing lower water table elevations and increased f l u i d flow. Barry (1981) suggests that atmospheric thermal lapse rates may vary from less than 2 K/km to greater than 8 °/km. In heat flow studies, a linear thermal lapse rate G[ of 5 K/km i s often assumed. In t h i s study, non-linear thermal lapse rates are used in cases where lin e a r thermal lapse rates predict sub-zero temperatures. Two extremes in thermal lapse rate G^ are simulated; 2 K/km and 8 K/km. The greater lapse rate y i e l d s lower temperatures at higher elevations. This, in turn, causes higher water table elevations and s l i g h t l y reduced t o t a l flow. (Runs M1 through M4; Table 3.3). 3.4 DISCUSSION Numerical results described in the preceding section provide insight into the influence of topography, geology, climate and regional heat flow on the patterns and magnitude of groundwater flow in mountainous t e r r a i n . These results are summarized in Figure 3.12 by bar graphs of simulated response of water table elevation AWT x and t o t a l flow bQ to each factor that has been considered. This information provides a means for assessing the r e l a t i v e importance of each factor. Bar graphs shown in Figure 3.12 indicate that bulk permeability has the greatest impact on water table elevations and f l u i d flow while i n f i l t r a t i o n rate, g l a c i e r s , slope p r o f i l e and basal heat flow play important, but lesser, r o l e s . In most cases, each factor i s varied over a large portion of the possible range of uncertainty associated with estimates of parameter values. Results shown in Figure 3.12 suggest upper bounds for errors in predicting water table elevations and t o t a l flow given an inherent uncertainty in the parameter considered. Exceptions to t h i s generalization include the influence of bulk permeability, permeable fault zones and permeable horizons where only a small portion of the possible range of variation is tested. For example, only a f i v e - f o l d variation in bulk permeability i s shown in Figure 3.12 while uncertainties in permeability estimates are l i k e l y to exceed one or two orders of magnitude. Such large variations contribute to uncertainties in estimated water table elevations that might e a s i l y exceed 1000 m. Furthermore, c h a r a c t e r i s t i c s for permeable fault zones and horizons are selected to indicate the minimum thickness and permeability contrast required to influence mountain flow systems. Clearly, thicker zones with higher permeability that intersect other permeable zones w i l l have a much greater influence than that shown in Figure 3.12. Results shown in Figure 3.12 indicate the possible error associated with each individual factor, the t o t a l error associated with a p a r t i c u l a r simulation would r e f l e c t the combined influence of each individual error. Geology influences steady state mountain flow systems through the s p a t i a l v a r i a b i l i t y of rock permeability and, to a lesser degree, thermal conductivity. Permeability has a greater impact, in part, because thermal conductivity has a much narrower range of v a r i a b i l i t y (bulk permeability may range from less than 10" 2 0 m2 to about 10" 1 4 m2 while thermal conductivity may range from about 0.5 to 5 W/mK). In addition, numerical results indicate a smaller flow system response occurs when thermal conductivity is varied by the , same factor as bulk permeability. The strong influence of rock permeablility is also noted in the influence of thin permeable.horizons (100 m) and fracture zones (0.1 m). Because sparsely distributed features contribute l i t t l e to the bulk permeability of a mountain massif, estimates of to t a l flow through the mountain massif are only s l i g h t l y influenced by the absence of information on thin permeable zones. Water table configurations and patterns of groundwater flow, however, may be strongly influenced by unidentified permeable horizons and fracture zones. Numerical results indicate that surface topography has a pervasive influence on underlying mountain flow systems. Variations in slope p r o f i l e , topographic r e l i e f and three-dimensional form control the ov e r a l l patterns of groundwater flow in addition to influencing water table elevations and t o t a l flow. Of p a r t i c u l a r importance i s the a b i l i t y of asymmetric topography to displace groundwater flow divides from the corresponding topographic divides. In using model studies to assess the source of chemical signatures in water samples from springs or in estimating inter-basin groundwater flow, the region modeled should, extend s u f f i c i e n t distance to include ridge and valley topography adjacent to the mountain slope of inte r e s t . Although surface topography can e a s i l y be resolved using readily available topographic information, incorporating complex three-dimensional topography into an analysis of mountain flow systems i s less e a s i l y accomplished. Thermal energy introduced into mountain flow systems by regional heat flow plays an important role in defining subsurface temperatures that, in turn, influence the nature of groundwater flow in mountains. Thermal regimes computed for flow systems with homogeneous isotropic permeability and normal regional heat flow (60 mW/m2) suggest that temperatures 2 km below the valley floor may range from 90°C to 40°C as the bulk permeability of the mountain massif i s varied from less than 1 0 - 1 8 m2 to 10~ 1 5 m2. Doubling the basal heat flow causes a doubling of temperatures in conduction-dominated and fluid-upflow zones while temperatures in other regions of the flow system are e s s e n t i a l l y unaffected. Because f l u i d density and v i s c o s i t y depend upon temperature, increased basal heat flow leads to higher subsurface temperatures, increased rates of groundwater flow and reduced water table elevations. Although basal heat flow i s normally assumed to range from about 30 to 120 mW/m2, advective heat transfer by groundwater flow in mountainous t e r r a i n often i n h i b i t s accurate estimates of heat flow. This uncertainty in basal heat flow, therefore, contributes to uncertainty in estimating water table elevation and t o t a l flow. This dependence of the groundwater flow system on regional heat flow suggests that isothermal approaches to modeling groundwater flow in mountainous t e r r a i n may be inappropriate. Climate influences the nature of mountain flow systems through variations in available i n f i l t r a t i o n rates, variations in surface temperature and the presence or absence of alpine g l a c i e r s . Numerical results indicate that large variations in mean annual surface temperature (5 to 15 °C) and thermal lapse rate (2 to 8 K/km) exert only a minor influence on mountain flow systems. Therefore, although surface temperature measurements can be made with reasonable accuracy, only rough estimates may be needed when assessing water table configurations and groundwater flow rates in mountainous t e r r a i n . Although available i n f i l t r a t i o n rates are expected to vary over a much narrower range than rock permeability, each parameter can exert a similar influence on water table elevations in mountain flow systems. Unfortunately, i t is d i f f i c u l t to make measurements of available i n f i l t r a t i o n rates that are meaningful at a mountain scale. Numerical results indicate that, i f a reasonable estimate of i n f i l t r a t i o n rate can be made (perhaps as a percentage of mean annual p r e c i p i t a t i o n ) , the s p a t i a l v a r i a t i o n of i n f i l t r a t i o n on the upper surface contributes l i t t l e to variations in the overall nature of the mountain flow system. Therefore, in most cases a uniform i n f i l t r a t i o n rate can be assumed. Climatic conditions that contribute to development of alpine glaciers can have an important impact on flow systems in permeable mountainous t e r r a i n where glaciers act as barriers to groundwater recharge. The character and s p a t i a l d i s t r i b u t i o n of alpine g l a c i e r s is r e l a t i v e l y e a s i l y defined by f i e l d mapping, however, the d e t a i l s of sub-glacial groundwater recharge are d i f f i c u l t to assess without a complete water budget analysis within the basin of inte r e s t . Calculated t o t a l flow through each hypothetical flow system represents a deep subsurface component of baseflow that is often neglected in water budget computations. The magnitude of the t o t a l flow i s dominated by bulk permeability (represented here by the upper zone permeability ku) and i n f i l t r a t i o n rate I z . In thi s study, an upper l i m i t for i n f i l t r a t i o n to deep flow systems i s estimated to be 2x10~8 m/sec (0.6 m/yr). This large i n f i l t r a t i o n can only be transmitted through high permeability t e r r a i n (in excess of 10' 1 5 m2) in a mountain slope with r e l i e f of 2 km over 6 km and a normal regional heat flow (60 mW/m2). As bulk permeability is reduced, the maximum i n f i l t r a t i o n that can be transmitted by the flow system i s reduced. Consider an i n f i l t r a t i o n rate of 2x10" 1 0 m/sec (0.006 m/yr). This i n f i l t r a t i o n rate i s a small portion of t y p i c a l p r e c i p i t a t i o n rates (in excess of 1.0 m/yr in humid coastal climates and as low as 0.4 m/yr in semi-arid climates) yet can only be transmitted through te r r a i n with permeability in excess of 10" 1 6 m2. Therefore, in t e r r a i n with permeability less than 10" 1 6 m2 and with a humid climate, only a small portion of the mean annual p r e c i p i t a t i o n i s transmitted through deep flow systems. Recharge to deep flow systems may be approximated by measuring v e r t i c a l hydraulic gradients and. permeability below the bedrock surface. Such measurements are costly and usually represent conditions only in the immediate v i c i n i t y of the measurement point. Furthermore, extrapolating measured values to other parts of the flow region i s fraught with uncertainty. In many instances, thin s u r f i c i a l deposits and open fractures near the bedrock surface may transmit a substantial portion of groundwater that appears as baseflow in water budget calculations yet i s not transmitted through deep flow systems. Therefore, estimates of recharge to deep flow systems may be poorly approximated by a water balance approach that assumes a l l baseflow travels through the deep flow system. The water table elevation i s a sensitive indicator of the influence of factors c o n t r o l l i n g groundwater flow systems. This high degree of s e n s i t i v i t y suggests that predicting water table configurations in mountainous terrain w i l l be d i f f i c u l t due to uncertainties inherent in resolving the magnitude of many of the c o n t r o l l i n g factors. Mapping water table elevations high in the mountain massif through d r i l l i n g or geophysical surveys, however, may a s s i s t in constraining parameters such as rock permeability and available i n f i l t r a t i o n rate. It seems reasonable to assume that water table elevations obtained from such a c t i v i t i e s w i l l r e f l e c t an averaging of permeability and i n f i l t r a t i o n rate that would be d i f f i c u l t to obtain from a series of point measurements of each parameter. Unfortunately, uncertainties in basal heat flow (in addition to d i f f i c u l t i e s inherent in defining the character and d i s t r i b u t i o n of thin permeable horizons and f a u l t zones) i n h i b i t defining a unique set of c h a r a c t e r i s t i c s that might cause the observed water table configuration. Throughout th i s Chapter, steady state hydrologic and thermal conditions have been assumed. Diurnal and seasonal variations in climate cause periodic fluctuations in surface temperature and i n f i l t r a t i o n rates. Buntebarth (1984) shows that seasonal changes in surface temperature are less than 10 - 3 °C at depths in the order of 30 m in conduction-dominated thermal regimes. In mountain-scale systems, such shallow penetration of seasonal temperature changes w i l l not induce s i g n i f i c a n t variations to o v e r a l l patterns and magnitudes of groundwater flow. Seasonal changes in i n f i l t r a t i o n rate are reported to cause seasonal fluctuations in mountaintop water table elevations of about 10 to 20 m (Halstead, 1969 and Province of B r i t i s h Columbia, 1974). Numerical results described in the previous section indicate that minimal changes in flow patterns and t o t a l flow rates are associated with such small changes in water table elevation. Therefore, periodic diurnal and seasonal transients (both thermal and hydraulic) can be neglected when studying deep groundwater flow in mountainous t e r r a i n . Longer-term variations in thermal and hydrologic conditions may result from ongoing processes of igneous a c t i v i t y , orogeny, erosion and climate change. In assessing the geochemical history and t r a n s i t times of groundwater samples, the influence of these processes may become important. Recall that a c h a r a c t e r i s t i c t r a n s i t time along the lowermost pathlines in Figure 3.1 can be estimated by assuming the available i n f i l t r a t i o n rate Iz i s a reasonable estimate of volumetric s p e c i f i c discharge within the domain (Section 3.2). For a system with permeability 10~ 1 5 m2, tra n s i t times are about 8x105 years. G l a c i a l events and plutonic a c t i v i t y (Norton and Knight, 1977) can cause s i g n i f i c a n t perturbations to the hydrologic and thermal regimes during t h i s time period. Major g l a c i a l advances and retreats have also occurred over the same time period. F l u i d t r a n s i t times in the order of 105 years, and longer, suggest that f l u i d entering the system as recharge at a par t i c u l a r point in time may bear the imprint of hydrologic, thermal and chemical conditions that no longer e x i s t . 3.5 CONCLUSIONS 1. Simulation results indicate that slope p r o f i l e , rock permeability, i n f i l t r a t i o n rate, g l a c i e r s and basal heat flow exercise a strong influence on water table configurations and the rates of groundwater flow. For example, a three-fold increase in permeability (or a three-fold decrease in available i n f i l t r a t i o n rate) can cause at least 500 m decline in water table elevation and a 45 percent change in t o t a l flow through a system with r e l i e f of 2 km over 6 km. Variations of similar magnitude are found to result when; a) two extremes of slope p r o f i l e are compared (concave and convex), b) alpine gl a c i e r s act to r e s t r i c t groundwater recharge, c) thin permeable fracture zones and horizons are incorporated in the flow system and d) two extremes of basal heat flow are compared (30 to 120 mW/m2). It should be noted that the dependence of water table elevation and f l u i d flux on regional heat flow suggests that model studies conducted using an isothermal approach may be inappropriate, even in regions with normal basal heat flow (60 mW/m2). 2. The high topographic r e l i e f of mountainous t e r r a i n amplifies the impact of slope p r o f i l e and thin permeable zones on patterns of groundwater flow. For example, asymmetry in ridge topography alone can cause s i g n i f i c a n t displacement of upland groundwater divides that, in turn, may enhance inter-basin groundwater flow. Furthermore, lowland groundwater divides displaced by asymmetry in valley topography can lead to uncertainties in defining the source of chemica.l signatures found in groundwater samples obtained from springs and shallow boreholes. Thin permeable fault zones and horizons (thickness of 0.1 m and permeability 10" times the surrounding rock or thickness of 100 m and permeability 10 times the surrounding rock) exert a strong influence on groundwater flow patterns. This influence further enhances the p o s s i b i l i t y of inter-basin groundwater flow and complicates accurate interpretation of chemical sampling r e s u l t s . 3. In modeling groundwater flow in mountainous t e r r a i n , i t is useful to f i r s t simulate a f u l l y saturated system and compare computed recharge rates with an estimated minimum value for the threshold i n f i l t r a t i o n r a t i o . If computed recharge rates in the saturated case approach or exceed the threshold a free-surface approach should be considered. 4 . Simulation results suggest that water table positions are l i k e l y to be found at r e l a t i v e l y high elevations in most mountainous t e r r a i n . Exceptions to t h i s generalization should be expected in terrain with r e l a t i v e l y high permeability (bulk permeability in excess of 1 0 " 1 5 m2) and ar i d climate ( i n f i l t r a t i o n less than 1 0 " 1 1 m/sec), or where alpine g l a c i e r s r e s t r i c t groundwater recharge. In lower-permeability t e r r a i n , less than one percent of t y p i c a l mean annual p r e c i p i t a t i o n rates may be transmitted through deep regions of groundwater flow. REFERENCES Barry, R.G., Mountain weather and climate, Methuen, 1981. Bryson, R.A., and F.K. Hare, Climates of North America, E l s e v i e r , 1974. Buntebarth, B., Geothermics, Springer-Verlag, 1984. Freeze, R.A., and J.A. Cherry, Groundwater, Prentice-Hall Inc., 1979. Halstead, E.C., Groundwater investigation, Mount Kobau, B r i t i s h Columbia, Canada Inland Waters Branch, Dept. of Energy Mines and Resources, Tech. B u l l . 17, 1969. Lauscher, F., Wettweite typen der hohenabgangigkeit des Niederschlags, Wetter U. Leben., 28, 80-90, 1976. Neuzil, C.E., Groundwater flow in low-permeability environments, Water Resour. Res., 22(8), 1163-1196, 1986. Norton, D., and J. Knight, Transport phenomena in hydrothermal systems: Cooling plutons, Amer. Jour, of S c i . , 277, 937-981, 1977. Paterson, W.S.B., The physics of g l a c i e r s , Pergamon, 1981. Province of B r i t i s h Columbia, Groundwater observation wells of B r i t i s h Columbia, B r i t i s h Columbia Water Resources Service, Water Investigations Branch, 1974. Ruffner, J.A., and F.E. B l a i r , Climates of the States with current tables of normals 1941-1970 and means and extremes to 1975, National Ocean and Atmospheric Adimin., Vols. I and I I , 1975. Schermerhorn, V.P., Relations bewteen topography and annual p r e c i p i t a t i o n in Western Oregon and Washington, Water Resour. Res., 3(3), 707-711, 1967. Slaymaker, H.O., and L.J. Zeman, Influence of a l t i t u d e and con t i n e n t a l i t y on watershed hydrology in the Coast Mountains of B r i t i s h Columbia, Proc. Canadian Hydrology Symposium, National Research Council, 1975. Storr, D., and H.L. Ferguson, The d i s t r i b u t i o n of p r e c i p i t a t i o n in some mountainous Canadian watersheds, i n , D i s t r i b u t i o n of Pr e c i p i t a t i o n in Mountainous Areas, Geilo Symposium, Norway, Proc. World Meteor. Assoc., WMO/OMM No. 326, Vol. II, 243-263, 1972. TABLE 3.1 SIMULATION SUMMARY AND GUIDE TO ILLUSTRATIONS - CHAPTER FACTOR SIMULATION FIGURE SERIES * TOPOGRAPHY Slope P r o f i l e A 3.1 Radial Symmetry I 3.5 Ridge Asymmetry J 3.6 Valley Asymmetry K 3.7 Topographic Relief 0 NS GEOLOGY Bulk Permeability C 3.2 Permeable Horizons D 3.8 Permeable Fracture Zones E 3.9 Basin Depth P NS Thermal Conductivity N NS Anisotropic Permeability F NS CLIMATE Available I n f i l t r a t i o n B 3.2 Alpine Glaciers G 3.10 Non-Uniform I n f i l t r a t i o n L NS Surface Temperature M NS THERMAL REGIME Basal Heat Flow H 3.11 Notes: * Simulation results are presented in Table 3.3 NS Not Shown 102 TABLE 3.2 TYPICAL SIMULATION PARAMETERS - CHAPTER 3 Fluid Flow Parameters kb permeability of basal unit 1 .0X10- 2 2 m2 ku * permeability of upper unit 1 .0x10- 1 5 m2 *z * v e r t i c a l i n f i l t r a t i o n rate 2.0x10"9 m/s< Thermal Parameters »b * basal heat flow 60.0 mW/m2 G l * thermal lapse rate 5 K/km T r * reference surface temperature 10 °C nb porosity of basal unit 0.01 nu * porosity of upper unit 0.10 \s * s o l i d thermal conductivity 2.50 W/mK \f f l u i d thermal conductivity 0.58 W/mK X? vapor thermal conductivity 0.024 W/mK Cf s p e c i f i c heat capacity of water 4186.0 J/kgK S saturation above water table 0.0 a l longitudinal thermal d i s p e r s i v i t y 100.0 m at transverse thermal d i s p e r s i v i t y 10.0 m Note: * Denotes parameters changed in simulation series. RUN PROFILE *.« '* * I # TYPE (m2) (m/sec) A. SLOPE PROFILE A1 X 1 0 - 1 5 3X10" 9 0.40 A2 V 10" 1 5 3x10"9 0.40 A3 L 10" 1 5 2x10"9 0.26 A4 L 10" 1 5 5x10" 9 0.66 B. INFILTRATION RATE B1 X 10- 1 5 1x10" 9 0.13 B2 X i o - 1 5 2x10" 9 0.26 B3 X 10" 1 5 5x1 0" 9 0.66 B4 X 10" 1 5 7x10" 9 0.92 B5 V • 10" 1 5 1x10" 9 0.13 B6 V 10-'5 2x1 0" 9 0.26 B7 V 10" 1 5 5x10- 9 0.66 B8 V 10" 1 5 7x10" 9 0.92 TABLE 3.3 SIMULATION RESULTS Q * Q WT max REF. AWT AQ COMMENTS (kg/sec) (m) RUN (m) (%) 1.2x10"2 0.26 1950 • • • 5.6x10"3 0.12 1053 A1 -897 -53 6.4x10"3 0.14 1270 • • • 7.5x10"3 0.17 1977 • • • 5.4X10" 3 0.12 924 A1 -1026 -55 1.0x10"2 0.22 1543 A1 -407 -17 1.2X10"2 0.26 2000 A1 + 50 0 > Threshold 1.2X10"2 0.26 2000 A1 + 50 0 > Threshold 3.3xl0" 3 0.07 490 A2 -563 -41 4.9X10" 3 0.11 805 A2 -248 -13 6.3x10"3 0.14 1486 A2 + 433 + 13 6.4X10" 3 0.14 1896 A2 + 843 + 14 TABLE 3.3 (Cont'd) SIMULATION RESULTS RUN PROFILE k„ I I* u z Q * Q max REF. # TYPE (m2) (m/sec) (kg/sec) (m) RUN C. UPPER ZONE PERMEABILITY Cl X 10' 1 8 2x10" 1 2 0.26 1 .0x10'5 0.22 1 560 B2 C2 X 10" 1 6 2x10" 1 0 0.26 1.0x10"3 0.22 1 550 B2 C3 V 10" 1 8 5x10" 1 2 0.66 6.2x10-6 0.14 1491 B7 C4 V 10" 1 6 5x10" 1 0 0.66 6.2x10-" 0.14 1488 B7 D. PERMEABLE HORIZONS (100 m thick, k = 10ifc > DI X .10" 1 6 2x10- 1 0 0. 22 1 .0x10"3 0.22 1 472 C2 D2 X 1 0 " 1 6 2x10' , 0 0.22 1.0x10"3 0.23 1316 C2 D3 V 1 0 " 1 6 5x10' 1 0 0.66 6.8x10' * 0.15 1 442 C4 D4 V 1 0 - 1 6 5x1 0- 1 0 0. 66 7.4x10-'" 0.16 1395 C4 E. PERMEABLE FRACTURE ZONES E1 X 10" 16 2x10" 1 0 0.22 1 .1 x 1 0" 3 0.24 1292 C2 E2 X 10" 16 2x10" 1 0 0.22 1 .1xl0- 3 0.24 1 263 C2 E3 X 10" 1 6 2x10- 1 0 0.22 1 . 1x10"3 0.24 1 245 C2 E4 X .10" 1 5 2x1 0- 9 0. 22 1 . 1x10-2 0.24 954 B2 L 7.5x10"3 0.17 1067 A3 AWT AQ (m) (%) COMMENTS + 17 -1 00 + 7 -90 + 5 -1 00 + 2 -90 -78 -1 Elev.=-I800m -234 + 4 Elev.=-100m -46 -10 Elev.=-1800m -93 + 19 Elev.=-100m -258 + 10 !,100m,10ku -287 + 10 !,0.1m,10** -305 + 10 /, 0. 1m, 10'* * M -589 + 10 / , 1.0m,10"£ u -203 + 17 Asymm. TABLE 3.3 (Cont'd) SIMULATION RESULTS RUN PROFILE ku .1 z I* Q Q* W?max REF. AWT AQ COMMENTS # TYPE (m2) (m/sec) (kg/sec) (m) RUN (m) (%) F. ANISOTROPIC PERMEABILITY F1 X 10" 1 6 2X1 0-.10 0.22 1 . 1 x 1 0" 3 0.24 1 1 06 C2 -444 +10 kz=^^u, kx=ku G. GLACIER BETWEEN ELEVATION 1500 m and 2000 m G1 X 10" 1 5 2X10" 9 0.26 3.7X10' 3 0.08 517 M5 -1099 -62 G2 V 10" 1 5 5X10" 9 0.66 4.7X1Q- 3 0.10 752 M7 -829 -18 H. BASAL HEAT FLOW H1 X 10 - 1 5 2x1 0-9 0. 26 1 . 1x10- 2 j' 0. 20 1911 B2 + 368 -10 Mb=20 mW/m2 H2 X 10" 1 5 2x1 0" 9 0. 26 9. 1x10' 3 0. 23 1 1 37 B2 -406 + 9 Hb=]20 mW/m2 H3 V i o - 1 5 5x1 0-9 0. 66 5. 1x10" 3 0. 1 1 1697 B7 + 21 1 -19 Hb=30 mW/m2 H4 V 1 o - 1 5 5x1 0' 9 0. 66 8. 3x1 0-3 0. 18 1 1 52 B7 -334 + 32 Hb=]20 mW/m2 TABLE 3.3 (Cont'd) SIMULATION RESULTS RUN PROFILE 4 /, /* Q Q* WTmnr REF. AWT AQ COMMENTS # TYPE (m2) (m/sec) (kg/sec) (m) RUN (m) (%) I. RADIAL SYMMETRY 11 X 1 0 - ' 5 2x10" 9 0. 26 1 . 1 x 1 0"2 0. 24 1 272 B2 -271 + 1 Planar 12 X 1 0 " 1 5 2x10' 9 0. 26 2.0x102 • 786 11 -486 • Radial 13 V 1 0 - 1 5 5x10" 9 0. 66 6 . 4 X 1 0 " 3 0. 14 1 484 B7 -2 + 2 Planar 14 V 10- 1 5 5x10" 9 0. 66 8 . 0 X 1 0 1 • 1 233 13 -251 • Radial J. ASYMMETRIC RIDGE J1 VX 1 0 " 1 5 3x10- 9 0. 40 1.8x10-2 • 1 488 • • • K. ASYMMETRIC VALLEY K1 X 10" 1 5 2x10" 9 0. 26 1.0x10-2 0. 22 1438 B2 -105 0 Relief=2 km L 6.0x10-3 0. 13 1208 A3 + 62 + 5 Relief=2 km K2 X 1 0 " 1 5 2x10" 9 0. 26 1.0x10"2 0. 22 1366 B2 -177 0 Relief=2 km L 3.3x10 - 3 0. 07 869 • • • Relief=1 km K3 V 10" 1 5 5x10" 9 0. 66 6.0x10"3 0. 13 1494 B7 + 8 0 Relief=2 km L 8.5xl0' 2 0 . 19 1758 A4 -219 + 13 Relief=2 km o TABLE 3.3 (Cont'd) SIMULATION RESULTS RUN PROFILE * I7 I Q # TYPE (m2) (m/sec) (kg/sec) L. NON-UNIFORM INFILTRATION PATTERNS L1 X 1 0 - 1 5 2x10" 9 0.26 9.8x10- 3 L2 V 1 0 - 1 5 5x10- 9 0.66 6.1x10" 3 L3 X i o - 1 4 2x10- 8 0.26 1.2x10" 1 L4 X 10" 1 u 2x1 0" 8 • 9.1x10" 2 L5 V 10" 1 4 2x1 0 " 8 0.26 4.9x10" 2 L6 X i o - 1 4 2x1 0- 8 • 4.1x10" 2 M. SURFACE TEMPERATURE CONDITIONS M1 X i o - 1 5 2 x 1 0 " 9 0 . 2 6 .1 , 0 x 10" 2 M2 X 10" 5 2x1 0 - 9 0 . 2 6 9 . 8 x 1 V 3 M3 V 10" 5 5x1 0 " 9 0 . 6 6 6 . 5 x 1 0 " 3 M4 V i o - 5 5x1 0 - 9 0 . 6 6 6 . 1 x 1 0 - 3 M5 X 10" 5 2 x 1 0 " 9 0 . 2 6 9 . 6 x 1 0 " 3 M6 X i o - 1 s 2x1 0 " 9 0 . 2 6 1 . 0 x 1 0 " 2 M7 . V 10" 1 5 5 x 1 0 " 9 0 . 6 6 5 . 7 x 1 0 " 3 M8 • V 10" t 5 5x1 0- 9 0 . 6 6 6 . 9 x 1 0 " 3 Q * ma x REF. AWT AQ COMMENTS (m) RUN (m) (%) 0. 1 5 1 528 B2 -15 -2 Slope Eff 0. 1 3 1 350 B7 -1 36 -3 Slope Eff 0. 22 1 543 » • • Uni form 0. 20 1 477 L3 -66 -24 Orogr. Ef 0. 1 1 805 • • • Uni form 0. 09 705 L5 -100 -16 Orogr. Ef 0. 22 1478 B2 -65 0 c7/=2°C/km 0. 22 1 569 B2 + 26 -2 C /=8°C/km 0. 14 1 41 5 B7 -71 + 3 G /=2°C/km 0. 1 3 1 536 B7 + 50 -3 G /=8°C/km 0. 21 1616 B2 + 73 -3 r r=5°c 0. 23 1 436 B2 -1 07 0 r r=i5 ° c 0. 1 3 1 581 B7 + 95 -10 r r=5°c 0. 15 1 375 B7 -111 + 10 T =15°C TABLE 3.3 (Cont'd) SIMULATION RESULTS RUN PROFILE * u I z I* Q Q* WTmQX REF. AWT AQ COMMENTS # TYPE (m2) (m/sec) (kg/sec) (m) RUN (m) (%) N. THERMAL CONDUCTIVITY OF ROCK-FLUID COMPOSITE N1 X 10" 1 5 2x10- 9 0. 26 1.0x10" 2 0. 22 1478 B2 -65 0 « M=0.01 N2 V 10" 1 5 7x10" 1 0 0. 92 6.2x10" 3 0. 14 1931 B8 + 35 -3 nu=0.0] N3 X 10" 1 5 2X10" 9 0. 26 9.6x10" 3 0. 21 1731 B2 + 188 -4 X*=3.5W/m°C N4 V 10" 1 S 5x1.0- 9 0. 66 5.6x10" 3 0. 12 1601 B7 + 115 -1 1 X*=3.5W/m°C 0. TOPOGRAPHIC RELIEF (1 km r e l i e f ) 01 V 10" 1 5 2x10" 9 0. 26 2.9x10" 3 0. 07 483 B6 -322 -41 02 V 10" 1 5 2x10"9 0. 26 4.8x10" 3 0. 10 836 B6 + 31 -2 Trunc. 03 X 10" 1 5 9x10 - 1 0 0. 1 2 5.0x10" ? 0. 1 1 936 • • • Trunc. P. BASIN DEPTH P1 V 10- 1 5 5X10" 9 0. 26 5.0x10" 3 0. 11 1590 B7 +94 -21 1 km deep TABLE 3.3 (Cont'd) SIMULATION RESULTS Notes: 1) Slope P r o f i l e ; X = convex, V = concave, L = l i n e a r . 2) !,l00m,10fc ; V e r t i c a l fracture zone at extreme right of system, thickness is 100 m, permeability is 10 times ku. 3) /,0.1m,10** ; Steeply dipping fracture zone outcropping at valley f l o o r , thickness is 0.1 m, permeability i s 10* times ku. 4) k kx; v e r t i c a l and horizontal components of anisotropic permeability. 5) H^ = basal heat flow. 6) C7y = thermal lapse rate. 7) T = valley reference temperature. 8) n u - upper zone porosity. 9) Xs = rock thermal conductivity. o Figure 3.1. Groundwater flow patterns, water table configurations and recharge-discharge p r o f i l e s as a function of slope p r o f i l e ; a) convex Run A1, b) concave Run A2. * Figure 3.2. Influence of i n f i l t r a t i o n r a t i o I on water table configuration, flow patterns and recharge-discharge p r o f i l e s ; a. convex (Runs B4, A1, B1), b. concave (Runs B8, A2, B5). Curve # Profile Relief Figure 3.3. Water table elevation ^ as a function of; * a. i n f i l t r a t i o n r a t i o I , and i b. dimensionless t o t a l flow Q * 113 Figure 3.4. Modified surface topographies; a. convex p r o f i l e s , b. concave p r o f i l e s . 1 1 4 Figure 3.5. Influence of r a d i a l symmetry on water table configurations, flow patterns and recharge-discharge p r o f i l e s (Runs 13 and 14). 115 0 2 4 6 8 10 12 DISTANCE (km) Figure 3.6. Influence of ridge asymmetry on water table configuration, flow patterns and recharge-discharge p r o f i l e (Run J1). 1 16 a. Run K1 1-0-1 eg 1.0-> ' S o -1 o s ' 0 -ec - 1.0-2 -t (km) 0 -NOI ELEVATI -2 --4 -b. Run K2 4 6 8 D I S T A N C E (km) Figure 3.7. Influence of valley asymmetry on water table configuration, flow patterns and recharge-discharge p r o f i l e s ; a. r e l i e f of opposing ridges are equal (Run K1), and b. r e l i e f of opposing ridges are unequal (Run K2). 1 1 7 ~ 1.0n C M I T ' o fl CO i - 0 CO x Run < > LU _1 111 System -2 --4 Homogeneous Perm, horizon 100m thick k=10k„ permeable horizon 2 4 DISTANCE (km) Figure 3.8. Influence of a permeable horizon 100 m thick with permeability 10 times ku located 100 m below the valley f l o o r (Run D4). 118 Run System DISTANCE (km) Figure 3.9. Influence of a steeply dipping fracture zone that outcrops at the valley floor (Run E4). The zone i s 0.1 m thick with permeability 10" times V 119 Figure 3.10. Influence of a glacier extending from an elevation of 2000 m to 1500 m above the valley floor on water table configuration, flow pattern and recharge-discharge p r o f i l e (Run G2). 120 Figure 3.11. Influence of basal heat flow on water table configuration, flow pattern and thermal regime; a. Hb = 60 m.W/m2 (Run B2), and b. H, = 120 mW/m2 (Run H2). F A C T O R V A R I A T I O N TOPOGRAPHY Slope Profile Radial Symmetry Valley Asymmetry Slope Effect on l z Topographic Relief Convex •** Concave Planar •*• Radial Convex * Linear Added Reduce by 0.5 GEOLOGY Bulk Permeabi l i ty Permeable Hor i zon: 100 m below valley Ver t i ca l F r a c t u r e Zone at right boundary B a s i n Depth Rock Thermal Conduct iv i ty Poros i ty 5 -Fo ld increase 100 m thick. k - 1 0 * k „ b - 1 0 0 m. k - 1 0 - k u or b-0.1 m, k - 1 0 4 - k u R e d u c e by 0.5 Increase: 2.5 3.5 W / m * C R e d u c e : 0.10 •*• 0.01 CLIMATE Avai lable Infiltration Alpine G l a c i e r s O r o g r a p h i c E f f e c t on l z Thermal L a p s e Rate Val ley Tempera ture 3 - F o l d decrease E l e v ' n : 1.5 2.0 km A d d e d Increase: 2 8 " C / k m Increase: 5 1 5 *C THERMAL REGIME B a s a l Heat Flow Increase: 60 -*• 120 m W / m 2 R E S P O N S E T O V A R I A T I O N * ^ W T m a x < m > Ao (%) * S lope Prof i le: Convex • C o n c a v e Figure 3.12. Summary of the influence of factors c o n t r o l l i n g wate: (fi). r table elevation (WT x) and t o t a l flow 122 CHAPTER 4 THERMAL REGIMES IN MOUNTAINOUS TERRAIN 4.1 INTRODUCTION It i s widely recognized that topographically-driven groundwater flow can perturb conductive thermal regimes (Sass et al. , 1971; Brott et al. , 1981; Mase et al. , 1 982). Numerical results described in Chapter 3 i l l u s t r a t e how h i g h - r e l i e f topography amplifies the influence of factors c o n t r o l l i n g the rates and patterns of groundwater flow in mountainous t e r r a i n . As a consequence, these factors also amplify the influence of groundwater flow on thermal regimes in mountains. In this chapter, the influence of these factors i s examined using the numerical method outlined in Chapter 2. This method d i f f e r s from those used by previous workers (Smith and Chapman, 1983, 1985; Garven and Freeze, 1984) in the way that the upper boundary condition for f l u i d flow i s described. Previous workers assume that the upper boundary of the flow domain coincides with the water table, even in cases where the water table l i e s below the bedrock surface. As a consequence, hydraulic head i s spe c i f i e d everywhere on the upper boundary and the processes of f l u i d flow and heat transfer acting in the unsaturated zone are ignored. The free-surface approach used in th i s study avoids 1 23 making assumptions regarding the water table configuration and includes a s i m p l i f i e d representation of f l u i d flow and heat transfer in the unsaturated zone. In t h i s way, temperatures at the water table are not sp e c i f i e d but rather, they are calculated in the solution procedure. In t h i s chapter, idealized mountain thermal regimes are modeled for a range of conditions representative of the Western C o r d i l l e r a in North America. Factors to be considered are surface topography, geology, climate, and regional heat flow. Elements of surface topography include slope p r o f i l e , r e l i e f and three-dimensional form. Both planar and axisymmetric mountain forms are considered. The geologic environment controls s p a t i a l v a r i a b i l i t y of permeability, porosity, and thermal conductivity. Climatic factors influence available i n f i l t r a t i o n , the presence and extent of alpine g l a c i e r s , and surface temperature conditions. Regional heat flow i s a c h a r a c t e r i s t i c of the tectonic environment within which the mountain i s located. It i s of p a r t i c u l a r interest to examine the influence of these factors on: (1) the nature and magnitude of advective disturbance of conductive thermal regimes, (2) the. interaction of free- and forced-convection within the mountain massif, (3) heat transfer within fracture zones, and (4) the nature of thermal springs. Table 4.1 summarizes the simulations carried out to assess the influence of each factor and provides a guide to the figures that i l l u s t r a t e each e f f e c t . 4.2 SIMULATION PARAMETERS Simulations are performed by assigning a set of f l u i d flow and thermal parameters (Table 4.2) within a geometry similar to that of Figure 2.1. The importance of surface topography i s i l l u s t r a t e d by considering two extremes in slope p r o f i l e ; one convex and one concave (Figure 4.1). Convex p r o f i l e s are t y p i c a l of glaciated c r y s t a l l i n e t e r r a i n while concave p r o f i l e s are often found in folded mountain belts and at volcanic cones. Parameter values given in Table 4.2, and the range of values tested in the numerical experiments, are chosen to r e f l e c t conditions t y p i c a l of mountainous t e r r a i n . Mean annual p r e c i p i t a t i o n is the basic determinant of available i n f i l t r a t i o n rates. In Chapter 3, i t i s shown that 2x10 _ 8 m/sec defines a reasonable upper l i m i t for available i n f i l t r a t i o n rates while a minimum / of 10" 1 2 m/sec i s assumed to occur in a l l but the most a r i d of climates. Although rock permeabilities can range from less than 10"'9 m2 to in excess of 10"9 m2, the bulk permeability of mountainous terrain considered in this study i s assumed to be less than 10" 1 4 m2. Freeze and Cherry (1979) and Neuzil (1986) suggest that this is a reasonable upper l i m i t for moderately fractured c r y s t a l l i n e and argillaceous rocks. In a l l simulations, a basal low-permeability unit (kfr = 10" 2 2 m2) occupies the lower 2 km of the domain to provide a region of conduction-dominated heat transfer (Figure 4.1). This region i s s u f f i c i e n t l y thick that the majority of simulation results show isotherms near the basal boundary to be sub-parallel to the boundary. As a consequence, the v e r t i c a l conductive heat flux applied at the basal boundary is transferred in a consistent manner to the thermal regime. The remainder of the domain contains a higher-permeability unit where advective heat transfer may dominate. In most simulations, both upper and lower zones have homogeneous and isotropic permeability (ku, k^) and uniform porosity (nu, n^). Although thermal conductivity of the s o l i d matrix (X 5) is uniform throughout the system, varying porosity and saturation produce contrasts in thermal conductivity of the s o l i d - f l u i d composite ( X e ) . Because we consider the steady state problem, porosity only has an in d i r e c t influence on the flow system through i t s impact on thermal conductivity of the s o l i d - f l u i d composite. A basal heat flux (#£,) i s applied on the horizontal base of the system while surface temperature conditions are defined in terms of a reference surface temperature (Tr) and a thermal lapse rate (Gj). Longitudinal and transverse thermal d i s p e r s i v i t i e s (a^, at ) are uniform throughout the system and are held constant for a l l simulations. In the following sections, frequent reference i s made to an advectively-disturbed reference case simulated using the parameters l i s t e d in Table 4.2. Results for the reference case, shown in Figure 4.1c, are used as a standard for comparison with subsequent simulation r e s u l t s . The reference Iz of 2x10"9 m/sec (Table 4.2) i s selected to represent the maximum rate of i n f i l t r a t i o n that might be available in a climate t r a n s i t i o n a l between semi-arid and humid. The reference upper unit permeability (*u) of 10" 1 5 m2 represents r e l a t i v e l y permeable conditions that might reasonably be expected in mountainous t e r r a i n . The reference heat flux (H^) of 60 mW/m2 is representative of a normal conductive regional heat fl u x . Thermal conductivity (X s) i s fixed at 2.5 mW/m2, a value that approaches the lower l i m i t for rocks found in mountainous t e r r a i n . Barry (1981) suggests that atmospheric thermal lapse rates may vary from less than 2 K/km to greater than 8 K/km. In thi s study, temperatures at the bedrock surface are defined using a median lapse rate (G[) of 5 K/km (Table 4.2) and a valley reference temperature (Tr) of 10 °C. 4.3 RESULTS OF NUMERICAL SIMULATIONS 4.2.1 Patterns of Advective Thermal Disturbance It i s recognized that the patterns and magnitude of groundwater flow dictate the character of an advective thermal disturbance. Numerical results described in Chapter 3 indicate that h i g h - r e l i e f mountainous ter r a i n amplifies the impact of factors that exert dominant control over groundwater flow systems; geology, climate and surface topography. In t h i s section, the influence of these factors on the character of advective thermal disturbance i s i l l u s t r a t e d by examining the patterns of groundwater flow and heat transfer in idealized mountainous t e r r a i n . Patterns of groundwater flow are depicted in Figure 4.1 by pathlines (dotted lines) representing the track of a f l u i d p a r t i c l e entering the flow system at a specified point on the bedrock surface. Pathline spacing i s inversely proportional to the flux of f l u i d ( s p e c i f i c discharge) through flowtubes bounded by each pair of pathlines* Recall that pathlines should not be confused with the streamlines generated from contour plots of a suitably defined stream function or v e l o c i t y p o t e n t i a l . The magnitude of / provides a convenient reference for comparing the v a r i a t i o n of f l u i d f lux. For example, pathline spacing in Figures 4.1c and 4.1e varies in response to the fact that f l u i d flux in the convex domain increases approximately f i v e f o l d from a value of 2x10"9 m/sec at the free-surface to about 10'8 m/sec at the discharge area. F l u i d flux i s more uniform in the concave domain (Figure 4.1g), decreasing from 2x10 _ 9 m/sec at the free-surface to about 1.6x10"9 m/sec. Patterns of heat transfer are shown in Figure 4.1 by heatlines (thick dashed lines) while subsurface temperatures are represented by isotherms ( s o l i d l i n e s ) . Kimura and Bejan (1983) suggest the use of heatlines as a means of mapping the transfer of thermal energy from the basal boundary by conduction and advection. As a consequence, each heatline shown in Figure 4.1 is l o c a l l y p a r a l l e l to the direction of net energy flow. In addition, the net energy flow across each heatline is zero. Note that t h i s y i e l d s heatlines normal to the isotherms only in regions where advective heat transfer by groundwater flow is ne g l i g i b l e , or where the dir e c t i o n of f l u i d flow i s normal to the isotherms. Heatlines are everywhere normal to isotherms in the conduction-dominated thermal regimes shown in Figures 4.1a, 4.1d and 4.1f. Heatlines and isotherms shown in Figure 4.1 i l l u s t r a t e the increasing influence of advective heat transfer on conductive thermal regimes as rock permeability i s increased three orders of magnitude. In each case, the available i n f i l t r a t i o n rate i s fixed at the reference value (/z = 2x10" 9 m/sec). Three series of simulation results are shown; the convex domain with normal heat flow (Figures 4.1a to 4.1c), the convex domain with doubled heat flow (Figures 4.1d and 4.1e) and the concave domain with normal heat flow (Figures 4.1 £ and 4.1g). Conductive temperature f i e l d s , obtained for ku equal to 10" 1 8 m2, are shown in the lefthand panel of each series (Figures 4.1a, 4.1d and 4.1f). A weak advective disturbance, characterized by heatlines at an oblique angle to the isotherms, is found when ku is 10" 1 6 m2 (Figure 4.1b). A strong disturbance, shown in the righthand panel of each series, i s found when ku exceeds about 10" 1 5 m2 (Figures 4.1c, 4.1e and 4.1g). In each case, the strong 129 disturbance i s characterized by heatlines subparallel to isotherms in the upper permeable unit. Shaded areas shown in Figures 4.1c, 4.1e and 4.1g indicate areas l i t t l e affected by thermal energy or i g i n a t i n g at the basal boundary. Temperatures within these regions r e f l e c t the temperature of groundwater recharge defined at the upper boundary using a thermal lapse rate. Elsewhere, heatline patterns indicate that the entire basal heat flow is absorbed by the c i r c u l a t i n g groundwater and deflected to exit in the discharge area. The basal heat flow i s e f f e c t i v e l y masked by downward flowing groundwater within the shaded regions. The t r a n s i t i o n from conduction-dominated to advection-dominated thermal regimes i s accompanied by l i t t l e change in patterns of groundwater flow within each simulation s e r i e s . Patterns of f l u i d flow d i f f e r somewhat between the high- and low-permeability cases in each series because lower water table elevations increase the thickness of the unsaturated zone where only v e r t i c a l f l u i d flow i s assumed. Patterns of f l u i d flow are also l i t t l e affected by changes in basal heat flow. Doubling the basal heat flow applied to the convex domain, from 60 to 120 mW/m2, causes an approximate doubling of temperatures in conduction-dominated and fluid-upflow regions of each domain (compare Figures 4.1a and 4.1c to 4.1d and 4.1e). Warmer . temperatures cause reduced f l u i d v i s c o s i t y and f l u i d density, a small increase in f l u i d flux (about 10 percent, on average), and reduced water table elevation. These changes have l i t t l e effect on patterns of f l u i d flow and heat transfer. The pattern of the thermal disturbance found for the concave slope p r o f i l e (Figure 4.1g) d i f f e r s from that of the equivalent convex case (Figure 4.1c) because surface topography has a pervasive influence on the pattern of groundwater flow. 4.3.2 Magnitude of Advective Thermal Disturbance The magnitude of an advective thermal disturbance i s characterized by the degree that conductive thermal regimes are cooled, or heated, by groundwater c i r c u l a t i o n . Contour plots of temperature residuals quantify the magnitude and sp a t i a l v a r i a t i o n of cooling and heating. Temperature residuals are calculated by subtracting temperatures obtained for a given simulation from those of the corresponding conductive case, at each node in the f i n i t e element mesh. For example, Figure 4.2a i s derived by subtracting results shown in Figures 4.1b from those shown in Figure 4.1a while Figure 4.2c i s derived using results shown in Figures 4.1d and 4.1e. Values of temperature residual equal zero along the upper boundary because the same boundary conditions are applied in each case. Temperature f i e l d s are affected only s l i g h t l y by a tr a n s i t i o n in thermal conductivity at the water table. This effect i s included in the computation of temperature 131 residual by using the same water table configuration in each pair of simulations. Figures 4.2a through 4.2d i l l u s t r a t e the contour plots of temperature residual obtained for the advectively disturbed thermal regimes shown in Figures 4.1b, 4.1c, 4.1e and 4.lg, respectively. In each case, advective heat transfer cools most of the thermal regime and heats only small regions near the valley f l o o r . Maximum cooling i s found beneath the mountain summit within, or near, the basal low-permeability unit. The degree of advective cooling i s controlled by upper zone permeability, i n f i l t r a t i o n rate, slope p r o f i l e and basal heat flow. In t h i s study, a steady uniform heat flux (rather than a fixed temperature) i s specified at the basal boundary. This allows cooling, caused by advective heat transfer in the upper permeable zone, to propagate to the base of each system by conduction within the basal low-permeability unit. If temperatures were sp e c i f i e d on the basal boundary, temperature residuals on the boundary would equal zero and maximum cooling would be found in central regions of the domain. Modeling r e a l i s t i c systems using this approach would require a knowledge of temperatures at depths of several kilometers. In cases where deep temperature data are unavailable, i t seems more appropriate to adopt the approach used here in defining the basal boundary condition using estimates of regional heat flow. 132 Figures 4.2a and 4.2b i l l u s t r a t e the increased cooling that occurs when ku i s increased tenfold from 10" 1 6 m2 to 10" 1 5 m2. Increased groundwater flux absorbs thermal energy supplied at the basal boundary with less heating of the domain. As a consequence, reduced values of temperature residual are found within a smaller region of heating in the higher-permeability case. Similar results are obtained, but not shown, when groundwater flux i s reduced by decreasing I one order of magnitude to 2x1tV 1 0 m/sec. In t h i s l a t t e r case, the maximum elevation of the water table i s only 40 meters above the valley f l o o r . Temperature residuals shown in Figure 4.2c for the doubled heat flow case (H^ = 120 mW/m2, shown in Figure 4.1e) are calculated with respect to a conductive temperature f i e l d with doubled heat flow, similar to that of Figure 4.1d. Temperatures in both doubled and normal heat flow cases are similar in regions of f l u i d downflow (shaded regions of Figures 4.1c and 4.1e). This s i m i l a r i t y evolves because c i r c u l a t i n g groundwater has the capacity to absorb a l l the thermal energy applied at the base of the system, even in the doubled heat flow case, with l i t t l e change in temperature within upper regions of the domain. Comparing Figures 4.2b and 4.2c i l l u s t r a t e s the increased cooling present in the doubled heat flow case. Given these res u l t s , i t appears impossible to i d e n t i f y the magnitude of the basal heat flow using temperature data c o l l e c t e d in shallow boreholes located at elevations above the valley f l o o r when ku exceeds about 10" 1 5 m2. •Comparing temperature residuals in Figure 4.2d with those of Figure 4.2b highlight the greater cooling that occurs in the domain with convex slope p r o f i l e . Although f l u i d flux i s similar in both convex and concave domains (the same I i s applied at the free-surface), the longer free-surface segment of the convex domain yields a t o t a l flow of groundwater about twice that of the concave domain. This provides a greater volume of f l u i d to absorb the basal heat flux, leads to greater cooling throughout the convex domain, and produces a smaller region of heating. At f i r s t glance, the. large d i s p a r i t y in the size of regions of heating and cooling shown in Figure 4.2 might be attributed to the elevation-dependent decrease in temperature imposed on the upper boundary through the thermal lapse rate. The influence of the thermal lapse rate is examined by simulating a convex domain with the reference conditions modified to maintain a constant surface temperature of 10 °C. Warmer groundwater recharge reduces the degree of advective cooling by about 10 °C but has only a small influence on the patterns of temperature residual and the size of regions of cooling and heating. Regions of cooling predominate in the maps of temperature residual because the basal heat flow i s readily absorbed, with l i t t l e heating, by the r e l a t i v e l y large flux of f l u i d through each system. Temperature residuals calculated for a number of simulations indicate that maximum heating i s found, for a speci f i e d and topographic r e l i e f of 2 km over 6 km, in systems with a bulk permeability of about 10" 1 6 m2 (Figure 4.2a). Reducing groundwater flux by reducing topographic r e l i e f , or permeability, produces larger regions of heating and reduced values of temperature residual. Regions of heating and cooling are approximately equal when a weak thermal disturbance is predicted (ku between 10" 1 8 m2 and 10' 1 7 m 2). Regions of heating extend over at least 50 percent of the lower-relief domains considered by Smith and Chapman (1983, 1985) and Garven and Freeze (1984) when a weak disturbance i s predicted. 4.3.3 Influence of Water Table Configuration The previous discussion indicates that water table configurations have l i t t l e influence on the pattern and magnitude of advective thermal disturbance. Results shown in Figure 4.3 highlight the strong contrast in thermal regimes that can develop under conditions where matching water table configurations are predicted. Matching water table configurations and patterns of groundwater flow are obtained when the r a t i o Iz/ku i s constant. In Figure 4.3 Iz/ku equals 4x10 s s e c ~ 1 i r r 1 . The significance of this r a t i o i s discussed in d e t a i l in Chapter 3. Figure 4.3a is obtained by setting ku equal to 5x10_ 1 5 m2 and Iz to 2x10"9 m/sec. Reducing the value of /_ by a factor of five to 4x10 _ 1 0 m/sec (&,. is also reduced by a factor of fi v e to 1 0 - 1 5 m 2), causes a f i v e f o l d decrease in groundwater flux. This reduced f l u i d flux creates a weaker thermal disturbance and a warmer thermal regime than that found for the h i g h e r - i n f i l t r a t i o n case of Figure 4.3a. Unique values of f l u i d flux cannot be associated with a particular water table configuration, therefore, the position of the water table provides l i t t l e d i r e c t insight into the magnitude of an advective thermal disturbance. Although water table configurations have l i t t l e d i r e c t influence on advective thermal disturbance, they provide insight into the parameters that do exert dominant control over groundwater flux and advective heat transfer; permeability and i n f i l t r a t i o n rate. Simulation results described in Chapter 3 indicate that groundwater flux increases almost l i n e a r l y with increasing k u , in cases with matching water table configurations. When the water table l i e s below the bedrock surface, f l u i d flux at the free-surface equals the available i n f i l t r a t i o n rate and an approximate linear relationship can be assumed between / and k u . The minimum value of Iz/ku required for the water table to coincide with the bedrock surface of the convex domain defines the s o l i d diagonal l i n e shown in Figure 4.4. This l i n e i s termed the water table threshold and i s a ch a r a c t e r i s t i c of a spec i f i e d domain geometry. Points p l o t t i n g below the water table threshold shown in Figure 4.4 indicate i n f i l t r a t i o n and permeability conditions that y i e l d predicted water table elevations below the bedrock surface. Under such conditions, groundwater flux within the domain can be characterized by 1 z . Values of Iz and ku p l o t t i n g above the water table threshold indicate when the water table coincides with the bedrock surface. Under these conditions, groundwater flux for the spec i f i e d domain can be characterized by ku because Iz is unlikely to bear a direct relationship to groundwater flux within the domain. Triangles shown in Figure 4.4 are la b e l l e d to indicate the combinations of 7, and ktl used to obtain the results z u i l l u s t r a t e d in Figures 4.1 and 4.3. The dashed l i n e joining points l a b e l l e d 4.3a and 4.3b in Figure 4.4 indicates that the approximate linear relationship between / and ku also produces matching water table configurations that need not coincide with the bedrock surface. In t h i s case, the maximum water table elevation i s 460 m (Figures 4.3a and 4.3b). Increasing the basal heat flow, or simulating a concave slope p r o f i l e , produces water table elevations lower than those predicted by the reference convex case (compare Figures 4.1e and 4.1g to Figure 4.1c). Therefore, greater i n f i l t r a t i o n rates are required to maintain water table elevations everywhere at the bedrock surface. As a consequence, water table thresholds for the concave and doubled heat flow cases, are shifted to the l e f t of that found for the reference convex case (Figure 4.4). The hydrologic implications of Figure 4.4 provide the basis for deriving insights into the onset of advective thermal disturbance discussed in the following section. 4.3.4 Onset of Advective Thermal Disturbance Previous workers have defined the advective threshold in a variety of ways; Norton and Knight (1977) use a threshold f l u i d mass flux, Smith and Chapman (1983) use a threshold permeability, while Domenico and Palciauskas (1973) and van der Kamp (1982) use a threshold value of a Peclet number that expresses the r a t i o between heat transfer by advection and conduction. In thi s study, advective thresholds are defined on the basis of permeability and f l u i d flux. Simulation results presented by Norton and Knight (1977), by Smith and Chapman (1983), and in thi s study indicate that conductive thermal regimes develop when f l u i d flux i s less than a threshold of about 10" 1 1 m/sec. When the water table l i e s below the bedrock surface, t h i s threshold f l u i d flux can be characterized by / because applying values of available i n f i l t r a t i o n rate in excess of 10" 1 1 m/sec causes a s i g n i f i c a n t disturbance of conduction-dominated thermal regimes (shaded region in Figure 4.4). When the water table coincides with the bedrock surface, the threshold f l u i d flux must be defined on the basis of a threshold k u . This value i s estimated from Figure 4.4 as the value of ku required to produce a water table coincident with the bedrock surface when / equals about 10" 1 1 m/sec. A ku of 4x10" 1 8 m2 provides the threshold f l u i d flux for the convex case with normal heat flow and marks the 138 upper l i m i t of conduction-dominated heat transfer (shaded region in Figure 4.4). Results obtained by Smith and Chapman (1983) indicate a greater advective threshold for their low-relief t e r r a i n of 1 km over 40 km. The dotted l i n e shown in Figure 4.4 represents the water table threshold estimated from their unpublished results and indicates a permeability-controlled threshold of about 3x10" 1 7 m2. Because lower hydraulic gradients develop in regions with reduced topographic r e l i e f , greater permeability i s required to provide the minimum groundwater flux for s i g n i f i c a n t advective heat transfer. This dependence on topographic r e l i e f means that estimated permeability-controlled thresholds are only v a l i d for a s p e c i f i e d domain geometry. The close proximity of the water table thresholds shown on Figure 4.4 suggest that water table thresholds for the mountain-scale r e l i e f considered here are l i t t l e affected by changes in basal heat flow and slope p r o f i l e . As a consequence, the threshold permeability used to define the onset of an advective disturbance varies less than half an order of magnitude in response to variations in slope p r o f i l e or basal heat flow (Figure 4.4). The numerical results described above suggest that a knowledge of available i n f i l t r a t i o n rate, or permeability, cannot be used in i s o l a t i o n to determine the onset and magnitude of advective thermal disturbance. Knowing the approximate position of the water table w i l l a i d in establishing the parameter that provides greatest insight into the magnitude of groundwater flux and the magnitude of the thermal disturbance. Consider the plot of Iz against ku shown in Figure 4.4. If the threshold Iz is applied in high-permeability t e r r a i n , groundwater recharge i s minimal and the water table l i e s deep below the bedrock surface. Because the corresponding groundwater flux i s small, the advective thermal disturbance is n e g l i g i b l e . Thermal disturbance is also n e g l i g i b l e in low-permeability t e r r a i n , despite a higher water table elevation. As I i s increased above i t s advective threshold, the water table' i s found at higher elevations within the system and groundwater flux increases to cause greater thermal disturbance. An upper l i m i t for the magnitude of the disturbance i s obtained, for a sp e c i f i e d value of k u , when I plots at or above the appropriate water table threshold (Figure 4.4). This l i m i t occurs because the water table coincides with the bedrock surface and further increase in / has no effect on the groundwater flow system. In such cases, the magnitude of the advective disturbance is characterized by the bulk permeability of the mountain massif. Figure 4.4 provides insight into the conditions where advective thermal disturbance is un l i k e l y . Maximum groundwater flux can be expected in hig h - r e l i e f t e r r a i n of high permeability with elevated regional heat flow and humid climate. Values chosen in thi s study approach the upper l i m i t for topographic r e l i e f (2 km over 6 km) and basal heat flow (120 mW/m2). These conditions promote maximum rates of groundwater flow and suggest that a minimum permeability, below which advective disturbance of regional-scale thermal regimes i s unlikely, is about 10" 1 8 m2. This lower l i m i t i s l i t t l e affected by variations in slope p r o f i l e and basal heat flow. As topographic r e l i e f is reduced, however, the lower l i m i t for advective disturbance is found at greater values of permeability. 4.3.5 Free-Convection in Mountainous Terrain Elder (1967) describes a numerical experiment designed to examine the influence of external f l u i d flow on patterns of free-convection within a permeable rectangular slab heated from below. In modeling t h i s system, he showed that external f l u i d flows imposed on the upper and v e r t i c a l boundaries of the slab could ob l i t e r a t e the free-convection c e l l s that would otherwise develop. Figure 4.5 i l l u s t r a t e s the results of a similar experiment, performed as part of t h i s study, where the nature of the external f l u i d flow i s e x p l i c i t l y defined by the character of the regional groundwater flow system. Figure 4.5a i l l u s t r a t e s the thermal regime and patterns of f l u i d flow associated with development of free-convection c e l l s in mountainous t e r r a i n , given the geologic conditions shown in Figure 4.5d. A thin low-permeability horizon (100 m thick with permeability ) separates two permeable units; an upper unit with ku equal to 10" 1 5 m2 and a lower unit (equivalent to Elder's permeable slab) with fcj equal to 5x1 tV 1 5 m2. Pathlines shown in the righthand panel of Figure 4.5a depict two free-convection c e l l s that develop in the lower permeable unit when a basal heat flow of 120 mW/m2 is applied and the upper and lower units are convergent heatlines shown in Figure 4.5a highlight the region of upwelling f l u i d in the free-convection c e l l s . Heatlines pass through the low-permeability horizon to be swept in the d i r e c t i o n of topographically-driven flow in the upper permeable zone. Therefore, the basal heat flux has l i t t l e influence on temperatures d i r e c t l y beneath the mountain summit (shaded region of Figure 4.5a). The rectangular geometry of the lower permeable layer allows ca l c u l a t i o n of a Rayleigh number Ra to confirm that free-convection should be expected under the conditions stipulated. Ra for the system shown in Figure 4.5a i s calculated using the d e f i n i t i o n of Cheng (1978) and the parameter values shown below; isolated by setting ki to 10~ 1 9 m2. Upwarped isotherms and p2g$cfk,Hzr (4.1) Ra = M X * where g g r a v i t a t i o n a l constant = 9.8 m/sec2 dynamic v i s c o s i t y of water (70°C) = 4x10'* Pa«s P density of water (70°C) = 978 kg/m3 C s p e c i f i c heat capacity of water = 4186 J/kgK H thickness of lower permeable unit = 1800 m AT = v e r t i c a l temperature difference = 90 K k[ = permeability = 5X10~ 1 5 m2 0 = thermal expansivity of water = 5x10-" K - 1 \ e = thermal conductivity of s o l i d - f l u i d composite = 2.3 W/mK The re s u l t i n g value of Ra i s 173. Because th i s value exceeds the generally accepted threshold of 40 (Cheng, 1978), i t i s reasonable that the numerical model should predict the weak c e l l u l a r free-convection indicated within the lower permeable unit (Figure 4.5a). Increasing the permeability of the intervening unit ki to 10" 1 7 m2 increases groundwater flux between the upper and lower permeable units and increases the t o t a l flow through the system. Increased l a t e r a l flow in the lower permeable unit displaces the region of upwelling to the right, increases cooling beneath the mountain summit and causes warmer temperatures beneath the valley floor (Figure 4.5b). Further increasing kt to 10' 1 6 m2 completely o b l i t e r a t e s the free-convection c e l l s to produce a strongly disturbed thermal regime and patterns of groundwater flow similar to those of Figure 4.1g. Under the conditions simulated here, free-convection c e l l s are e a s i l y disturbed by topographically-driven f l u i d flow because groundwater flux in the permeable upper unit i s about 150 times greater than the buoyancy-driven f l u i d flux of the lower permeable unit. Clearly, free-convection c e l l s are less l i k e l y to be disturbed as the permeability of the upper unit i s reduced or as topographic r e l i e f i s reduced. Increasing k{ by one order of magnitude to 5x10" 1 4 m2 causes more vigorous free-convection within the lower permeable layer (Figure 4.5c) and produces f l u i d flux similar in magnitude to that of the upper permeable layer (approximately 10 - 9 m/sec). Using equation (4.1) and the updated parameters shown below, Ra for thi s system i s 524. M = dynamic v i s c o s i t y of water (45°C) = 6x10-° Pa«s p = density of water (45°C) = 990 kg/m3 AT = v e r t i c a l temperature difference = 40 K kf = permeability = 5 x l 0 - 1 4 m2 Total flow of f l u i d through the system and the pattern of f l u i d flow i s l i t t l e changed from that of Figure 4.5a. Increased buoyancy-driven f l u i d flux in the lower permeable layer decreases the li k e l i h o o d that free-convection c e l l s w i l l be displaced, or obliterated, as the permeability of the intervening horizon is increased. Enhanced rates of buoyancy-driven f l u i d flow focus transfer of the basal heat flow within the region of f l u i d upwelling. This y i e l d s s i g n i f i c a n t cooling throughout the lower permeable layer, except at the top of the region of upwelling (Figure 4.5c). Here, temperatures are similar to those of the weakly convecting case (Figure 4.5a). Temperatures in the upper permeable unit d i f f e r l i t t l e from those found in the weakly convecting case, except in the v i c i n i t y of the valley f l o o r . As a consequence, the character of the thermal regime within the lower permeable unit i s unlikely to be detected using thermal data co l l e c t e d in shallow boreholes. Additional work is underway to examine the interplay of free-convection c e l l s and topographically-driven f l u i d flow within more complex regional flow systems where igneous intrusions and permeable fracture zones may be encountered. This work i s expected to y i e l d useful insights regarding thermal regimes at active volcanic centers and the o r i g i n of epithermal ore deposits. 4.3.6 Influence of Mountain Topography Surface topography has a pervasive influence on the pattern and magnitude of groundwater flow in mountainous t e r r a i n , hence on the pattern of advective heat transfer. The influence of slope p r o f i l e on groundwater flow systems and thermal regimes has been described in previous sections. In the following paragraphs, the implications of topographic symmetry and three-dimensional form are addressed. Asymmetry in surface topography produces an advective thermal disturbance that r e f l e c t s the influence of asymmetry in patterns of groundwater flow. Three mountain valleys are shown in Figure 4.6, each with a convex slope p r o f i l e on the l e f t and a l i n e a r slope p r o f i l e on the r i g h t . Mirror symmetry i s assumed at the v e r t i c a l boundaries in each case. Figure 4.6a shows an asymmetric valley with opposing summits of equal elevation. Figure 4.6b shows an asymmetric valley with the elevation of the righthand summit reduced to one-half that of the lefthand summit. Figure 4.6c represents a symmetric system with a convex ridge adjacent to a wide valley with f l a t topography. The strongly disturbed thermal regimes shown in Figure 4.6 are obtained by assigning the reference conditions given in Table 4.2. Asymmetry in surface topography can cause a strong warping of isotherms near the valley floor (Figures 4.6a and 4.6b) that might be mistaken for the thermal signature of a free-convection c e l l or a permeable fracture zone; both features can cause upflow of heated f l u i d . The pronounced warping of isotherms is centered on groundwater flow divides located within regions of c l o s e l y spaced sub-vertical heatlines. Groundwater flow divides represent imaginary boundaries separating the d i s t i n c t regions of groundwater upflow associated with each slope p r o f i l e . The position of the divide r e f l e c t s a balance between the t o t a l flow of groundwater beneath each slope p r o f i l e . In an ideal symmetric case, the divide i s v e r t i c a l and located at the valley f l o o r . This is the case in Figure 4.6c and the previous Figures 4.1 and 4.4. Total flow through the linear slope p r o f i l e i s about 60 percent of that of the convex p r o f i l e , therefore, both the groundwater divide and the region of upwarped isotherms are displaced to the right in Figure 4.6a. As r e l i e f of the l i n e a r slope i s reduced, f l u i d flow through the linear slope i s reduced. Both the divide and the upwarping of isotherms are displaced further to the right in Figure 4.6b. Figures 4.6a and 4.6b indicate that the groundwater flow divide may intersect the bedrock surface upslope of the valley f l o o r in regions of asymmetric surface topography. As a consequence, chemical and thermal signatures of a subsurface heat source located beneath the lefthand convex slope may be expressed in springs discharging upslope from the valley floor on the opposing li n e a r slope. In addition, samples co l l e c t e d at the valley f l o o r may provide l i t t l e information regarding conditions beneath the linear slope p r o f i l e . As permeability i s reduced, the warping of isotherms becomes less pronounced. The position of the groundwater flow divide, however, i s unchanged. In the preceding discussion, each domain represents a v e r t i c a l section through a linear ridge of i n f i n i t e extent. Although f u l l y three-dimensional systems cannot be simulated with the current model, insights into the influence of other topographic forms can be gained by simulating axisymmetric conical features. The r a d i a l equivalent of the planar domain shown in Figure 4.6c i s an isolated mountain massif with convex slope p r o f i l e and axis of symmetry beneath the ridgetop. Because the horizontal component of f l u i d flux decreases with increasing r a d i a l distance, f l u i d flux in the r a d i a l case decreases from the reference value of /„ with increasing r a d i a l distance. This t r a n s i t i o n is reflected in the plot of f l u i d mass flux normal to the upper boundary of each system shown in Figure 4.7a. / i s the same in each case, therefore, f l u i d mass flux is the same along each free-surface. Elsewhere, f l u i d discharge in the r a d i a l case is less than that of the planar case. This difference i s most readily seen near the break-in-slope, where most discharge occurs. The f l u i d flux p r o f i l e shown in Figure 4.7a indicates that only a small percentage of the t o t a l flow discharges near the righthand boundary. This suggests that, a f l a t valley floor with a width in excess of about two times the horizontal distance covered by the adjacent mountain slope w i l l e f f e c t i v e l y i s o l a t e groundwater flow systems and thermal regimes on opposing sides of the v a l l e y . For example, Figure 4.7 shows a valley half-width of 6 km (for a t o t a l width of 12 km) and a horizontal slope length of 6 km. Despite the contrast in topographic symmetry, temperatures in the r a d i a l case are only s l i g h t l y warmer than those of the planar case. The difference in temperature obtained in each case i s shown in the contour plot of Figure 4.7b. The greatest temperature difference occurs in central regions of the system, beneath the break-in-slope, where greatest contrast in groundwater flux i s encountered (Figure 4.7a). In t h i s region, temperatures in the r a d i a l system are elevated about 8 °C above those of the planar system. L i t t l e difference between the planar and r a d i a l cases i s indicated near the right and l e f t boundaries of the system because the contrast between the patterns and magnitude of groundwater flow in each case is minimal. 4.3.7 Permeable Fracture Zones and Thermal Springs The numerical method used in th i s study allows permeable fracture zones to be embedded, as discrete e n t i t i e s , within the surrounding rock mass. In contrast to previous approaches (Lowell, 1975; Sorey, 1978; K i l t y et al. , 1979; Goyal and Kassoy, 1980; Bodvarsson et al. ,1982), f l u i d flux within the fracture is dictated by conditions c o n t r o l l i n g the regional flow system rather than by a spe c i f i e d f l u i d source at depth or by a sp e c i f i e d uniform f l u i d flux in the fracture zone. An important difference l i e s in the fact that t h i s approach allows f l u i d to enter, or leave, the fracture zone at any point along i t s length. Previous studies often assume that the fracture wall i s impervious. In t h i s modeling study, fracture zones are represented by a series of connected planar segments with uniform thickness (b) and a homogeneous equivalent porous media permeability {kj-). Because an approximate li n e a r r e l a t i o n s h i p exists between permeability and f l u i d flux, various combinations of b and kj- produce the same patterns of groundwater flow in systems with matching water table configurations and matching fcyfc; when b i s less than about 100 m. The product kj'b (expressed in units of m 2«m) i s termed the transmissivity of the fracture zone and i s used in discussing the simulation r e s u l t s . Note that this d e f i n i t i o n of transmissivity d i f f e r s from the product of hydraulic conductivity times thickness usually adopted in isothermal approaches. Figures 4.8a and 4.8b include a steeply dipping fracture zone, extending from the valley floor to the basal boundary. Fracture width b and permeability kj- are assumed constant everywhere along the fracture zone. The transmissivity of the fracture zone kfb i s assumed to equal 10"'&U (m 2»m) and might reasonably correspond to several combinations of b and kj-; b of 1 m and kj- of ^0'>^ku, b of 10 m and kj of 10 3«& U, or b of 100 m and kj- of 10 2«* U. Figure 4.8a i l l u s t r a t e s the influence of the fracture zone described above, with ku of 1 0 - 1 5 m2, on groundwater flow and heat transfer within the asymmetric topography of Figure 4.6a. The t o t a l flow of groundwater increases by 75 percent to y i e l d cooler temperatures everywhere in the domain. Patterns of f l u i d flow and heat transfer are greatly modified because eighty percent of the t o t a l f l u i d flow, and the entire basal heat flow, i s captured by the fracture zone. Rapid f l u i d flow in the fracture zone, and. the surrounding rock mass, yiel d s almost isothermal conditions along the fracture and a spring temperature elevated only 3 °C above the ambient surface temperature of 10 °C. Temperatures are everywhere cooled below those of the corresponding conduction-dominated case and the basal heat flow is e f f e c t i v e l y masked everywhere along the upper boundary (shaded region in Figure 4.8a). Within the low-permeability basal unit, f l u i d flux in the fracture zone i s 6 orders of magnitude less than that of shallower fracture segments because inflow from the basal unit is minimal. As a consequence, the conductive pattern of heat flow within the basal unit i s undisturbed. Above the basal unit, inflow from the surrounding rock mass causes f l u i d flux in the fracture zone to increase almost l i n e a r l y (about one order of magnitude); reaching a maximum at the bedrock surface. Decreasing the permeability of the surrounding rock mass to 10' 1 6 m2, while retaining the transmissivity of the fracture zone at 1 0 9 , * U (m 2«m), reduces f l u i d flux in both the fracture zone and the surrounding rock mass. This e f f e c t , in turn, y i e l d s a a warmer thermal regime (Figure 4.8b) and a spring temperature 13 °C warmer than that of the higher-permeability case. The pattern of upwarped isotherms shown in Figure 4.8b d i f f e r s considerably from those of K i l t y et al. (1979) and Sorey (1978), who predict high temperature gradients at shallow depths in the fracture zone. This difference i s attributed to the fact that previous workers assume an impervious fracture boundary and a uniform f l u i d flux within the fracture zone. In t h i s study, the a b i l i t y for f l u i d to enter, or leave, the fracture zone allows a non-uniform pattern of f l u i d flux to develop. F l u i d flux in the fracture zone increases about one order of magnitude between an depth of 2 km and the valley floor in the system shown in Figure 4.8. As a consequence, a more uniform temperature gradient develops along the fracture zone. A conduction-dominated thermal regime is found when ku equals about 10" 1 8 m2. This result i s independant of the magnitude of the transmissivity of the fracture kfb because the f l u i d flux in the flow system i s dominated by the permeability of the surrounding rock mass. Although rates of f l u i d flow are i n s u f f i c i e n t to cause a perceptible change in the temperature f i e l d , the volume of f l u i d flow in the fracture zone i s s u f f i c i e n t to cause a spring temperature elevated 5 °C above the ambient of 10 °C. Figure 4.9 i l l u s t r a t e s the va r i a t i o n of spring temperature as a function of upper unit permeability (ku), for both the normal and doubled heat flow cases. Peaks in spring temperature (25 °C for H^ equal to 60 mW/m2 and 43 °C for Hb equal to 120 mW/m2) are found when rock mass permeability equals about 10" 1 6 m2. This value of ku also produces maximum advective heating in the unfractured convex slope p r o f i l e (compare temperature residuals shown in Figures 4.2a, 4.2b and 4.2c). This behaviour r e f l e c t s the fact that the ov e r a l l character of the thermal regime i s controlled by regional groundwater flow through the rock mass. Temperature residuals shown in Figure 4.2a indicate that a region of heating develops near the valley floor that contributes to development of a thermal spring. Increasing ku above 10* 1 6 m2 causes reduced spring temperatures (Figure 4.9) because the rapidly c i r c u l a t i n g groundwater absorbs a l l thermal energy applied at the base of the system with l i t t l e increase in temperature. Decreasing ku below 10" 1 6 m2 also causes reduced spring temperatures (Figure 4.9) because f l u i d flow to the fracture is reduced and a greater portion of the basal heat flow i s transferred by conduction through the surrounding rock mass. Changing the e f f e c t i v e length of the fracture zone modifies i t s influence on the thermal regime. The e f f e c t i v e length is increased by considering the influence of permeable cross-fractures that intersect the o r i g i n a l fracture zone. The resultant increase in f l u i d flux throughout the domain causes a greater portion of the basal heat flow to be captured by the fracture zone. As a consequence, spring temperatures are warmer and recharge areas are cooler. For example, the length of the fracture zone shown in Figure 4.8b i s increased by introducing a second fracture zone, 2 km in length, beneath the convex slope (within the upper permeable unit) to intersect the o r i g i n a l fracture at right angles where both fractures meet the top of the basal unit. This increased length y i e l d s a greater spring temperature of 76 °C. Reducing the length of the fracture zone to 1 km causes only a s l i g h t reduction in spring temperature because eliminating inflow from deeper 153 regions of the upper permeable unit has only a small effect on f l u i d flux in the fracture zone. The transmissivity of the fracture zone k f b dictates the degree of disturbance caused by the fracture zone. In the system shown in Figure 4.8b, the fracture zone exerts i t s maximum influence when k f b equals about I 0 a , * t t (m 2«m). Increasing k j - ' b to 1 0 5 « & u (m 2«m) has l i t t l e e f f e c t on the regional flow system and causes only a 1 °C increase in spring temperature (open c i r c l e plotted in Figure 4.9). Reducing k f b to 1 0 3 , ^ U (m 2«m) has a greater influence, causing a 7 °C decrease in spring temperature (Figure 4.9). The fracture zone r e s t r i c t s f l u i d flow in the surrounding rock mass when k f b i s less than 10" - A: w (m 2»m). At values of k j - ' b in excess of this value, f l u i d flux in the system i s r e s t r i c t e d only by the permeability of the surrounding rock mass. The results shown in Figures 4.8 and 4.9 are representative only of the idealized geometries and conditions tested. These results should only be used as a guide to the character of thermal regimes influenced by permeable fracture zones because patterns of thermal disturbance and spring temperatures are strongly controlled by a variety of factors that are d i f f i c u l t to map and quantify. These factors include; fracture position, orientation, transmissivity of the fracture zone, length, and intersection with other fractures. In addition, the three-dimensional nature of intersecting fracture zones and rugged mountainous t e r r a i n makes detailed extrapolation from two-dimensional simulations d i f f i c u l t . Figure 4.9 suggests, however, that an optimal range of bulk permeability exists (ku = 10' 1 7 m2 to 1 0 - 1 5 m2) wherein maximum spring temperatures are most probable. It should be noted that l o c a l i z e d shallow heat sources, not considered in t h i s study, could promote the development of high-temperature springs even when ku l i e s outside t h i s range. Results presented in this section suggest that i t w i l l be d i f f i c u l t to predict accurately the temperature of thermal springs in any but the simplest of cases. 4.4 IMPLICATIONS FOR HEAT FLOW STUDIES Estimates of regional heat flow are often made by measuring v e r t i c a l temperature gradients in shallow boreholes (less than 200 m deep), correcting the results for the influence of surface topography and c a l c u l a t i n g the regional heat flow using estimates of thermal conductivity (Sass et al . , 1971). In mountainous t e r r a i n , an additional correction may be required to account for advective disturbance of the thermal regime. Contour plots of temperature gradient r a t i o provide insight into the magnitude and s p a t i a l v a r i a t i o n of errors in temperature gradient introduced by advective heat transfer. This method of presentation should be of p a r t i c u l a r interest to those attempting to interpret temperature gradients measured in boreholes. Temperature gradient r a t i o s are calculated by normalizing temperature gradients computed at each v e r t i c a l f i n i t e element boundary with respect to temperature gradients calculated in the corresponding conduction-dominated case. Figure 4.10 shows contour plots of gradient r a t i o calculated for the temperature f i e l d s shown in Figures 4.1b, 4.1c, 4.1e and 4.1g. Note that Figure 4.10 i l l u s t r a t e s the potential magnitude of errors introduced only by advective heat transfer. Heat flow estimates must also be corrected for the influence of surface topography, variations in thermal conductivity and the influence of thermal transients. Where exact corrections for these effects cannot be made, errors in heat flow estimates may exceed those indicated in Figure 4.10. Where the gradient ratios shown in Figure 4.10 equal 1.0, temperature gradients are unaffected by the advective disturbance and measured temperature gradients match those of the corresponding conduction-dominated case. Shaded areas shown in Figure 4.10 highlight regions where measured temperature gradients may provide estimates of uncorrected regional heat flow with an error of less than 25 percent. Temperature gradients are only s l i g h t l y affected at depth, even in systems with strong advective cooling, because the geometry of the boundary-value-problem i s designed to ensure that the pattern of conductive heat flow in the basal low-permeability unit i s undisturbed. In more complex situ a t i o n s , temperature gradient r a t i o s at depth may affected to a much greater degree. Where the gradient ra t i o s d i f f e r from 1.0, the basal heat flow is e f f e c t i v e l y masked by groundwater flow. Gradient ratios less than 1.0 indicate regions where temperature gradients measured in shallow boreholes would y i e l d underestimates of uncorrected regional heat flow. In terr a i n with r e l i e f of 2 km over 6 km, such conditions should be expected when the bulk permeabilty of the mountain massif exceeds about 10" 1 6 m2 and humid climate causes high recharge rates. Mase et al. ( 1 982) suggest that such conditions contribute to the regional low in heat flux (less than 30 mW/m2) identifed in permeable volcanic rocks of the Cascade Range in C a l i f o r n i a . Negative gradient ratios (hatched regions in Figures 4.10b, 4.10c and 4.l0d) indicate regions where temperature inversions and negative temperature gradients might be measured in shallow boreholes when bulk permeability equals about 10" 1 5 m2. Terrain with greater ku produces expanded regions of negative gradient. S i g n i f i c a n t negative temperature gradients are reported in the upper several hundred metres of boreholes on the concave flank of Mount Hood, Oregon (Steele and Blackwell, 1982). If the concave domain shown in Figure 4.l0d i s a reasonable approximation of Mount Hood, i t can be inferred that the permeability of the volcanic p i l e l i k e l y exceeds 10" 1 6 m2. In regions of reduced topographic r e l i e f , negative temperature gradients 157 w i l l only be found (on a regional-scale) in ter r a i n with greater bulk permeability. Brott et al. (1981) describe temperature inversions induced by regional groundwater flow within fractured basalts of the Snake River Plain. Topographic r e l i e f i s approximately 1 km over 400 km and bulk permeability for the basalts generally exceeds 10" 1 5 m2. Negative temperature gradients are unlikely in hi g h - r e l i e f t e r r a i n where ar i d climates y i e l d reduced rates of groundwater recharge. The gradient r a t i o map of Figure 4.10a suggests that thermal data c o l l e c t e d in shallow boreholes located on the mountain flank are most l i k e l y to provide good indications of the magnitude of the basal heat flux when ku i s less than about 10~ 1 6 m2 (or Iz i s less than 2X10" 1 0 m/sec). In te r r a i n with ku (or Iz) only one order of magnitude greater, shallow temperature data i s unlikely to y i e l d good estimates of regional heat flow (Figure 4.10b). Because ku and Iz are d i f f i c u l t to estimate to a precision better than one or two orders of magnitude, i t w i l l be d i f f i c u l t to v e r i f y that temperature data c o l l e c t e d from shallow boreholes provides a good estimate of regional heat flow. It should be noted, however, that the influence of thermal transients and uncertainties in understanding the thermal conductivity structure of the subsurface further reduce the l i k e l i h o o d that r e l i a b l e estimates of basal heat flux can be made using data obtained from shallow boreholes. 158 Gradient r a t i o maps indicate that uncorrected regional heat flow may be overestimated by a factor in excess of 2.0 near the valley floor (Figure 4.10). In strongly disturbed regimes i t appears necessary to d r i l l at least 2 km below the valley floor to obtain reasonable estimates of regional heat flow. Heat flow measurements made in rugged mountainous terr a i n are, by necessity, often r e s t r i c t e d to locations near the floor of narrow mountain v a l l e y s . Heat flow data coll e c t e d from shallow boreholes (less than 200 m deep) in mountain valleys of the Coast Mountains, B r i t i s h Columbia y i e l d consistent values of heat flow with plausible magnitude (about 80 to 100 mW/m2) without correcting for possible advective heat transfer (Lewis et al. , 1985). This result suggests that regional groundwater flow in the Coast Mountains causes, at most, a weak advective disturbance. Furthermore, t h i s implies that bulk permeability in the Coast Mountains is unlikely to exceed about 10" 1 6 m2. It should be re c a l l e d that an equivalent porous medium i s assumed to provide an adequate representation of the fractured rock mass, with the exception of major through-going fracture zones. As a consequence, regions of reduced bulk permeability correspond to regions of decreased fracturing rather than to regions of decreased matrix permeability. Souther (1975) notes that high-temperature phenomena such as b o i l i n g springs, fumaroles and mud pots are absent in the Coast Mountains of B r i t i s h Columbia; despite elevated 159 heat flow (Lewis et al. , 1985). Sixteen thermal springs, with temperatures seldom in excess of 65 °C, are scattered throughout a region of 1.3x10" km2 (Souther and Halstead, 1973; Souther, 1975). Because these springs invariably issue from fractures in c r y s t a l l i n e rock, the associated thermal regimes are presumed to be strongly controlled by permeable fracture zones. The sparse d i s t r i b u t i o n of thermal springs suggests that conditions favoring elevated spring temperatures, are found only in a few l o c a l i z e d regions. In the absence of l o c a l heat sources, the location of the observed thermal springs may indicate regions where bulk permeability approaches 10" 1 6 m2. Elsewhere, bulk permeability l i k e l y f a l l s outside the range of ku (10" 1 7 m2 to 10" 1 5 m2) required for maximum spring temperatures. As noted previously, the plausible heat flow values reported for t h i s region suggest that strong advective disturbance i s absent at a regional scale and bulk permeability i s l i k e l y less than 1 0" 1 6 m2. Interpreting heat flow data can be complicated by the presence of alpine g l a c i e r s mantling a mountain massif. Figure 4.11 i l l u s t r a t e s the influence of an alpine glacier on the temperature f i e l d in a convex domain with permeability 10" 1 5 m2 and i n f i l t r a t i o n rate 2X10" 9 m/sec. The valley reference temperature (TR) i s reduced to 5 °C to represent a cooler climate. Using the approach described in Chapter 3, recharge beneath the gl a c i e r i s limited to that derived from sub-glacial melting by thermal energy 160 o r i g i n a t i n g at the basal heat source. Downslope of the g l a c i e r margin, rates of groundwater recharge equal the available i n f i l t r a t i o n rate on the free-surface. This t r a n s i t i o n , from reduced recharge beneath.the glacier to high recharge downslope of the g l a c i e r margin, produces elevated subsurface temperatures and a reduced advective disturbance beneath the g l a c i e r . Borehole temperature data co l l e c t e d near glacier margins in high-permeability t e r r a i n may be d i f f i c u l t to interpret because the distorted groundwater flow system disturbs the patterns of heat transfer (Figure 4 . 1 1 ) . The impact of alpine glaciers is reduced as permeability i s reduced because less i n f i l t r a t i o n i s accepted downslope of the g l a c i e r margin while rates of groundwater recharge generated by sub-glacial melting are l i t t l e a f f ected. It should be noted that similar results would be obtained where extensive units of low-permeability rocks or permafrost are found in upper regions of the domain. Results presented in t h i s section i l l u s t r a t e the use of numerical modelling to assess the disturbance of conductive thermal regimes by groundwater flow. A free-surface approach should be considered when confronted with a deep water table because temperature conditions at the water table are d i f f i c u l t to define without e x p l i c i t l y including the influence of heat transfer in the unsaturated zone. Such situations are most l i k e l y to be found in regions with high-permeability rocks and a r i d climate. A free-surface approach may be unneccessary when temperatures at the water table can be assumed approximately equal to those at the ground surface. Such conditions might be found in low-permeability t e r r a i n where a water table close to the bedrock surface produces a thin unsaturated zone. Al t e r n a t i v e l y , high f l u i d flux found in high-permeability t e r r a i n with humid climate and deep water table may allow temperatures at the bedrock surface to be extrapolated v e r t i c a l l y downward to an assumed water table. If the purpose of a modeling excercise i s to ide n t i f y regions where b o i l i n g might occur, a free-surface approach should be considered to ensure reasonable pressure conditions within the mountain massif. 4.5 CONCLUSIONS 1. Recognizing the pattern and magnitude of groundwater flow i s fundamental to defining the pattern and magnitude of advective disturbance of conductive thermal regimes. A f i n i t e element model of f l u i d flow and heat transfer provides insight into the character of advective heat transfer in mountainous t e r r a i n . A free-surface approach eliminates the dependence of the solution results on d i f f i c u l t - t o - e s t i m a t e water table configurations by computing water table elevations as part of the solution process. Such an approach i s most useful when assessing thermal regimes in high-permeability t e r r a i n with a r i d climate. Here, the water table may l i e deep below the bedrock surface. Conventional approaches that make a priori assumptions regarding the water table configuration require temperatures at the ground surface to be extrapolated to a water table that defines the upper boundary of the groundwater flow system. In more humid climates with lower permeability, water table elevations approach the bedrock surface and a free-surface approach may be unneccessary. 2. The position of the water table y i e l d s insight into the parameters that characterize the magnitude and onset of an advective disturbance. Where the water table l i e s below the bedrock surface, an advective threshold can be defined in terms of the available i n f i l t r a t i o n rate. Conductive thermal regimes are expected when available i n f i l t r a t i o n rates are less than about 10~ 1 1 m/sec. Conditions causing low available i n f i l t r a t i o n rates and deep water table elevations are most l i k e l y to be found in regions with a r i d climate. Where the water table coincides with the bedrock surface, the available i n f i l t r a t i o n rate no longer characterizes f l u i d flux within the domain. In such cases, an advective threshold s p e c i f i c to the topographic r e l i e f of the domain can be defined on the basis of bulk permeability. Numerical results suggest that the permeability-controlled advective threshold i s about 10' 1 8 m2 in t e r r a i n with r e l i e f of 2 km over 6 km. The permeability threshold becomes proportionately greater as r e l i e f i s reduced. Advective disturbance of thermal regimes i s unlikely, on a regional 163 scale, where bulk permeability i s less than 10" 1 8 m2. 3. High-relief topography amplifies the influence of surface topography on patterns of groundwater flow and advective heat transfer. D i f f e r i n g slope p r o f i l e s and summit elevations on opposing valley walls y i e l d an asymmetry in advective heat transfer that complicates the interpretation of thermal data sets. Thermal regimes beneath opposing valley slopes are e f f e c t i v e l y isolated when separated by a f l a t valley floor that is wider than the base of each adjacent mountain massif. 4. Numerical results indicate that special conditions are required in the variation of permeability within the mountain massif i f free-convection c e l l s are to develop. Otherwise, free-convection c e l l s may be o b l i t e r a t e d by topographically-driven groundwater flow. 5. Fracture zones can exert a profound influence on the pattern and magnitude of groundwater flow and advective heat transfer. There i s a non-uniform variation of f l u i d flux within the fracture zone that r e f l e c t s the character of the regional flow system. The strong interaction of f l u i d flow in the fracture zone and the surrounding rock mass produces an optimal range of permeability (between about 10" 1 7 m2 and 10" 1 5 m2) where the temperature of thermal springs may reach a maximum. In lower-permeability t e r r a i n , f l u i d flux i s too slow to cause s i g n i f i c a n t advective heat transfer. As a consequence, spring temperatures approach ambient a i r temperature. In higher-permeability t e r r a i n , rapid f l u i d flux absorbs a l l thermal energy originating at the base of the domain with l i t t l e increase in temperature. These conditions also contribute to spring temperatures approaching ambient a i r temperature. 6. Advective heat transfer causes an ov e r a l l cooling of a mountain massif with small regions of heating near the valley f l o o r . High-permeability and large available i n f i l t r a t i o n rates promote greater cooling and smaller regions of heating. In te r r a i n with permeability in excess of about 10" 1 5 m2, the entire basal heat flow is absorbed by c i r c u l a t i n g groundwater and deflected to exit at the valley f l o o r with l i t t l e increase in subsurface temperatures. Such conditions mask the character of the basal heat flow and preclude using temperature data c o l l e c t e d in shallow boreholes to estimate conductive regional heat flux. 7. V e r t i c a l temperature gradients may be overestimated by at least a factor of 2.0 where measured in mountain v a l l e y s . In high-permeability t e r r a i n (in excess of about 10" 1 5 m 2), negative temperature gradients are predicted at median elevations on the mountain flank and within the mountain massif. Temperature data c o l l e c t e d from shallow boreholes located on the mountain flank may y i e l d reasonable estimates of v e r t i c a l gradients and regional heat flow when bulk permeability i s less than about 10" 1 6 m2. 165 8. The numerical model used in th i s study provides a basis for explaining the general character of thermal conditions observed in boreholes and at thermal springs located in mountainous t e r r a i n . Detailed examination of thermal anomalies at s p e c i f i c s i t e s requires developing a three-dimensional equivalent of the existing two-dimensional model. \ 166 REFERENCES Barry, R.G., Mountain weather and climate, Methuen, 1981. Bodvarsson, G.S., S.M. Benson, and P.A. Witherspoon, Theory of the development of geothermal systems charged by v e r t i c a l f a u l t s , Jour. Geophys. Res., 87, 9317-9328, 1982. Brott, C.A., D.D. Blackwell, and J.P. Ziagos, Thermal and tectonic implications of heat flow in the Eastern Snake River Plain, Idaho, Jour. Geophys. Res., 86, 11709-11734, 1981. Cheng, P., Heat transfer in geothermal systems, Adv. in Heat Transfer, Vol. 14, 1-105, 1978. Domenico,P.A. and V.V. Palciauskas, Theoretical analysis of forced convective heat transfer in regional ground-water flow, Geol. Soc. Amer., B u l l e t i n 84, 3303-3814, 1973. Elder, J.W., Steady free convection in a porous medium heated from below, J. F l u i d Mech., 27, 29-48, 1967. Freeze, R.A., and J.A. Cherry, Groundwater, Prentice-Hall Inc., 1979. Garven, G. and R.A. Freeze Theoretical analysis of the role of groundwater flow in the genesis of stratabound ore deposits: 2. Quantitative r e s u l t s , Am. Jour. S c i . , 284, 1125-1174, 1984. Goyal, K.P. and D.R. Kassoy, Fault zone controlled charging of a liquid-dominated geothermal reservoir Jour. Geophys. Res, 85, 1867-1875, 1980. K i l t y , K., D.S. Chapman, and C.W. Mase, Forced convective heat transfer in the Monroe Hot Springs geothermal area, Jour, of Vole, and Geoth. Res., Vol. 6, 257-272, 1979. Kimura, S., and A. Bejan, The "heatline" v i s u a l i z a t i o n of convective heat transfer, J. of Heat Transfer, Trans, of the ASME, Vol. 105, 917-919, 1983. 167 Lewis, T.J., A.M. Jessop, and A.S. Judge, Heat flux measurements in Southwestern B r i t i s h Columbia: the thermal consequences of plate tectonics, Can. Jour. Earth S c i . , 22, 1262-1273, 1985. Lowell, R.P., C i r c u l a t i o n in fractures, hot springs, and convective heat transport on mid-ocean ridge crests, Geophys. J. Roy. Astr. S o c , 40, 351-365, 1975. Mase, C.W., J.H. Sass, A.H. Lachenbruch, and J.R. Munroe, Preliminary heat-flow investigations of the C a l i f o r n i a Cascades, U.S. Geol. Survey, Open-File Report 82-150. 1982 Neuzil, C.E., Groundwater flow in low-permeability environments, Water Resour. Res., 22(8), 1163-1196, 1986. Norton, D., and J. Knight, Transport phenomena in hydrothermal systems: Cooling plutons, Amer. Jour, of S c i . , 277, 937-981, 1977. Sass, J.H., A.H. Lachenbruch, R.J. Munroe, G.W. Greene, and T.H. Moses, J r . , Heat tflow in the Western United States, J . Geophys. Res., 76(26), 6376-6413, 1971. Smith, L., and D.S. Chapman,The influence of water table configuration on the near-surface thermal regime, Jour, of Geodynamics, 4, 183-198, 1985. Smith, L., and D.S. Chapman, On the thermal e f f e c t s of groundwater flow 1. Regional scale systems, Jour. Geophys. Res., 88(1), 593-608, 1983. Sorey, M.L. Numerical modeling of l i q u i d geothermal systems, U.S. Geol. Surv. Prof. Paper 1044-D, 1978. Souther, J.G. Geothermal potential of Western Canada, Proc. 2nd. U.N. Symp. on Development and Use of Geothermal Resources, San Francisco, C a l i f . , V ol. 1, 259-267, 1975. Souther, J.G., and Halstead, E.C., Mineral and thermal waters of Canada, Canada Dept. of Energy, Mines and Resources, Paper 73-18, 1973. 168 Steele, J . L., and D.D. Blackwell, Heat flow in the v i c i n i t y of the Mount Hood Volcano, Oregon, in Geology and geothermal resources of the Mount Hood area Oregon, Oregon Dept. of Geol. and Min. Ind., Special Paper 14, pp 31-42, 1982. van der Kamp, G., Interactions betwee heat flow and groundwater flow - A review, Proj. 109-17, Waterloo Res. Inst., Univ. of Waterloo, Waterloo, Ont., 1982. 169 TABLE 4.1 SIMULATION SUMMARY AND GUIDE TO ILLUSTRATIONS - CHAPTER 4 FACTOR FIGURE TOPOGRAPHY Slope P r o f i l e 4.1 Valley Asymmetry 4.6 Radial Symmetry 4.7 GEOLOGY Bulk Permeability 4.1, 4.3 Permeable Fracture Zones 4.8 CLIMATE Available I n f i l t r a t i o n 4.1, 4.3 Alpine Glaciers 4.11 Surface Temperature NS THERMAL REGIME Basal Heat Flow 4.1 Free-Convection 4.5 Thermal Springs 4.8, 4.9 ADVECTIVE THERMAL DISTURBANCE Onset 4.1,4.4 Pattern 4.1,4.3 Magnitude 4.2, 4.10 Not es: NS Not Shown 170 TABLE 4.2 TYPICAL SIMULATION PARAMETERS - CHAPTER 4 Fluid Flow Parameters permeability of basal unit 1.0x10"22 m2 * permeability of upper unit 1.0x10"15 m2 * v e r t i c a l i n f i l t r a t i o n rate 2.0x10"9 m/s Thermal Parameters Hb * basal heat flow 60.0 mW/m2 Gi * thermal lapse rate 5 K/km Tr * reference surface temperature 10 °C nb porosity of basal unit 0.01 nu porosity of upper unit 0.10 Xs s o l i d thermal conductivity 2.50 W/mK x' f l u i d thermal conductivity 0.58 W/mK Xv vapor thermal conductivity 0.024 W/mK Cf s p e c i f i c heat capacity of water 4186.0 J/kgK s saturation above water table 0.0 al longitudinal thermal d i s p e r s i v i t y 100.0 m at transverse thermal d i s p e r s i v i t y 10.0 m Note: * Denotes parameters changed in simulation se r i e s . Groundwater pathline _,—• ^* Heatline ~ — 5 Isotherm ( ° C ) Water table KM Region not penetrated by basal heat flux o Basal low-permeability unit DISTANCE (km) 2 4 DISTANCE (km) F i g u r e 4.1 P a t t e r n s of groundwater flow and heat t r a n s f e r i n convex and concave topography s i m u l a t e d wi th the r e f e r e n c e c o n d i t i o n s of T a b l e 4.2 (except where n o t e d ) ; a . *« = 10" 1 8 m 2 b . ku = 10" 1 s m2 c . k = u 10" 1 5 m 2 ( r e f e r e n c e c a s e ) , d . 10" 1 8 m 2 and Hb = 120 mW/m 2, e. k = u 10" 1 5 m 2 and Hb = 120 mW/m 2, f . ku = 10" 1 8 m2 g. ku - 10-1 5 m 2 • DISTANCE (km) DISTANCE (km) Figure 4,2. Contour plots of temperature residual (°C); a. ku -• 10" 1 6 m2 and »b m 60 mW/m2, b. ku -= 1 0 " 1 5 m2 and Hb = 60 mW/m2, c. ku -= 10" 1 s m2 and Hb - 120 mW/m2, d. ku -= 10" 1 5 m2 and Hb " 60 mW/m2. 0 2 4 6 D I S T A N C E (km) Figure 4.3. Matching water table configurations in convex topography with Hb = 60 mW/m2; a. ku = 5X 1 0 - 1 5 m2 and Iz = 2X10" 9 m/sec, b. ku = 10" 1 5 m2 and Iz = 4x10- 1 0 m/sec. 174 10 - 9 10 H 01 -10 £10 H 10 - i i - 1 2 10 (4.1c,4.1e,4.1g) ;(4.1a,4.1d,4.1f) 10 - 1 9 10 • 1 4 ku(m2) A (2a) Figure number showing simulation results for plotted Iz, k u i I Conduction-dominated: Convex slope profile, H b = 60 mW/m 2 Water Table Thresholds Concave slope profile Convex slope profile Convex slope profile Relief mW/rrr km/km 60 2/6 120 2/6 60 2/6 .f' Smith & Chapman (1983) 60 1/40 Max. water table elevation = 460 m ^ • ' ^ Convex slope profile 60 2/6 Figure 4.4. Influence of i n f i l t r a t i o n Iz and permeability ku on water table elevations and advective thresholds. 175 F i g u r e 4 . 5 . Thermal regimes and p a t t e r n s of groundwater flow i n mixed f r e e - and f o r c e d - c o n v e c t i o n s c e n a r i o s w i t h Hb = 120 mW/m 2, Iz = 5x10" 9 m/sec and ku = 1 0 " 1 5 m 2 ; a . ki • -- 1 0 -1 9 m 2 and h • = 5 x 1 0 "1 s m 2 , b . = 1 0 " 1 7 m 2 and k l • = 5 x 1 0 -1 5 c . ki • = i o -1 9 m 2 and k l • = 5 x 1 0 -1 4 m 2 d . g e o l o g i c c o n d i t i o n s for s i m u l a t i o n s . 176 (km) Figure 4.6. Thermal regimes beneath mountain valleys simulated with reference conditions of Table 4 . 2 ; a. opposing summits have equal r e l i e f and d i f f e r i n g slope p r o f i l e , b. linear summit has r e l i e f one half that of the convex summit, c. convex summit adjacent to a f l a t p l a i n . 177 2 4 6 8 10 12 DISTANCE (km) Figure 4.7. Comparison of thermal regimes in r a d i a l and planar symmetry; a. f l u i d flux normal to the upper boundary, b. temperature difference contours (°C). 178 Figure 4.8. Influence of a steeply dipping fracture zone with uniform kfb = 10 a'* u (m 2»m) on thermal regimes within the asymmetric topography of Figure 4.6a; a. ' Jfc = 10" 1 5 m2, b. ku = 10- 1 6 m2, 46 O UJ rx r-< CE LU Q. HI tr CO 42 38 -34 -30 -26 -fj 22 Z 18 -O H b = 120 mW/m k, «b = 10 5 »k u (m2.m) j. • kf -b - 10 4 -k u (m2.m) o k f -b = 1 0 3 « k u (m 2 Tn) j H b = 60 mW/m 2 * k f 'b = 1 0 4 « k u ( m ^ m ) / V 14 -10 10 -19 -18 10 10 10 10 -15 10 k u (mz) Figure 4.9. Spring temperature as a function of upper zone permeability ku, transmissivity of the fracture zone kj--br and basal heat flow . 180 DISTANCE (km) DISTANCE (km) [ j 0.75 - temperature gradient ratio ^ 1.25 [ ' -J Temperature gradient ratio < 0.0 Figure 4.10. Contour plots of temperature gradient r a t i o ; a. ku = 10- 1 6 m2 and Hb = 60 mW/m2, b. ku = 10" 1 5 m2 and Hb = 60 mW/m2, c. ku = 10" 1 5 m2 and ,Hb = 120 mW/m2 d. ku = 10" 1 5 m2 and Hb = 60 mW/m2 181 Figure 4.11. Influence of a glacier mantling a convex mountain from 1500 m to 2000 m above the valley f l o o r . 182 CHAPTER 5 SUMMARY OF CONCLUSIONS A f i n i t e element model i s developed to simulate steady groundwater flow and heat transfer through v e r t i c a l sections of t e r r a i n with mountainous topography. A free-surface approach eliminates the dependence of the solution results on d i f f i c u l t - t o - e s t i m a t e water table configurations by computing water table elevations as part of the solution process. Conventional approaches that make a priori assumptions regarding the water table position require temperatures at the ground surface to be extrapolated to a water table that defines the upper boundary of the groundwater flow system. This d i f f i c u l t y i s avoided by assuming one-dimensional advective heat transfer by v e r t i c a l f l u i d flow in the unsaturated zone. When the water table l i e s below the bedrock surface, the rate of groundwater flow in the unsaturated zone i s controlled by the available i n f i l t r a t i o n rate. This parameter, viewed as a percentage of the mean annual p r e c i p i t a t i o n rate, defines the maximum rate of groundwater recharge possible for a given set of climatic conditions. This method i s most useful when simulating mountainous ter r a i n where the water table l i e s deep below the bedrock surface. Numerical results suggest that such conditions are most l i k e l y in high-permeability t e r r a i n (greater than 10" 1 5 m2) with a r i d climate, or where 183 groundwater recharge is r e s t r i c t e d by extensive alpine g l a c i e r s . Simulation results indicate that slope p r o f i l e , rock permeability, i n f i l t r a t i o n rate, g l a c i e r s and basal heat flow exercise a strong influence on water table configurations and the rates of groundwater flow. For example, a three-fold increase in permeability (or a three-fold decrease in available i n f i l t r a t i o n rate) can cause at least 500 m decline in water table elevation and a 45 percent change in t o t a l flow through a system with r e l i e f of 2 km over 6 km. Variations of similar magnitude are found to result when; a) two extremes of slope p r o f i l e are compared (concave and convex), b) alpine g l a c i e r s act to r e s t r i c t groundwater recharge, d) thin permeable fracture zones and horizons are incorporated in the flow system and d) two extremes of basal heat flow are compared (30 to 120 mW/m2). Recognizing the influence of these factors on rates and patterns of advective heat transfer i s fundamental to assessing the disturbance of conductive thermal regimes by groundwater flow. The high topographic r e l i e f of mountainous terrain amplifies the impact of slope p r o f i l e and thin permeable zones on patterns of groundwater flow and advective heat transfer * For example, asymmetry in ridge topography alone can cause s i g n i f i c a n t displacement of upland groundwater divides that, in turn, may enhance inter-basin groundwater flow. Furthermore, lowland groundwater divides displaced by asymmetry in valley topography can lead to uncertainties in defining the source of chemical or thermal signatures found in groundwater samples obtained from springs and shallow boreholes. Thin permeable fault zones and horizons (thickness of 0.1 m and permeability 10' times the surrounding rock or thickness of 100 m and permeability 10 times the surrounding rock) exert a strong influence on groundwater flow patterns. This influence further enhances the p o s s i b i l i t y of inter-basin groundwater flow and complicates accurate interpretation of chemical sampling results and borehole temperature measurements. The strong interaction of f l u i d flow in the fracture zone and the surrounding rock mass produces an optimal range of permeability (between about 10* 1 7 m2 and 10" 1 5 m2) where the temperature of thermal springs may reach a maximum. In t e r r a i n with bulk permeability values outside t h i s range, spring temperatures approach ambient a i r temperature. Simulation results indicate that simply i d e n t i f y i n g the position of the water table in mountainous t e r r a i n w i l l provide l i t t l e insight into the nature of the groundwater flow system and the character of the advective thermal disturbance. This result r e f l e c t s the fact that the water table should be viewed as an internal c h a r a c t e r i s t i c of the flow system, rather than an important c o n t r o l l i n g factor. The position of the water table, however, does y i e l d insight into the parameters that characterize the magnitude and 185 onset of an advective disturbance. Where the water table l i e s below the bedrock surface, an advective threshold can be defined in terms of the available i n f i l t r a t i o n rate. Conductive thermal regimes are expected when available i n f i l t r a t i o n rates are less than about 10~ 1 1 m/sec. Conditions causing low available i n f i l t r a t i o n rates and a deep water table are most l i k e l y found in regions with a r i d climate. In more humid climates where the water table coincides with the bedrock surface, the available i n f i l t r a t i o n rate no longer characterizes f l u i d flux within the domain. In such cases, an advective threshold can be defined on the basis of bulk permeability that i s s p e c i f i c to the topographic r e l i e f of the domain. Numerical re s u l t s suggest that the permeability-controlled advective threshold i s about 1 0 - 1 8 m2 in t e r r a i n with r e l i e f of 2 km over 6 km. The permeability threshold becomes proportionately greater as r e l i e f i s reduced. The numerical model used in t h i s study provides a basis for explaining the character of thermal conditions observed in boreholes and at thermal springs located in mountainous t e r r a i n . In addition, the numerical results provide insight into the possible errors associated with making heat flow measurements in mountainous t e r r a i n . Active c i r c u l a t i o n of groundwater in t e r r a i n with permeability in excess of about 1 0 - 1 5 m2 e f f e c t i v e l y masks the character of the basal heat flux. Such conditions preclude using temperature data c o l l e c t e d in shallow boreholes to estimate conductive regional heat flux and to i d e n t i f y underlying geothermal systems. Errors associated with estimating regional heat flux in advectively disturbed thermal regimes are quantified using a temperature gradient r a t i o . This r a t i o expresses the deviation of v e r t i c a l temperature gradients from those of the corresponding conductive case. V e r t i c a l temperature gradients may be overestimated by at least a factor of 2.0 where measured in mountain v a l l e y s . In high-permeability t e r r a i n (in excess of about 1 0 " 1 5 m 2), negative temperature gradients are predicted at median elevations on the mountain flank and within the mountain massif. Temperature data coll e c t e d from shallow boreholes located on the mountain flank may y i e l d reasonable estimates of v e r t i c a l gradients (and regional heat flow) when bulk permeability i s less than about 1 0" 1 6 m2. CUMULATIVE REFERENCE LIST Baca, R.G., R.C. Arnett, and D.W. Langford, Modelling f l u i d flow in fractured-porous rock masses by f i n i t e element techniques, International Journal for Numerical Methods in Fluids, 4, 337-348, 1984. Barry, R.G., Mountain weather and climate, Methuen, 1981. Bear, J., Dynamics of f l u i d s in porous media, American Elsevier, 1972. Bear, J., Hydraulics of Groundwater, McGraw-Hill, 1979. Birch, F., Flow of heat in the Front Range, Colorado, B u l l . Geol. Soc. Amer., 61, 567-630, 1950. Black, G.L., D.D. Blackwell, and J.L. Steele, Heat flow in the Oregon Cascades, in Geology and Geothermal Resources of the Central Oregon Cascades Range, ed. by G.R. Priest and B.F. Vogt, Oregon Dept. of Geol. and Min. Ind., Special Paper 15, 69-76, 1983. Blackwell, D.D., A transient model of the geothermal system of Long Valley Calder, C a l i f o r n i a , Jour. Geophys. Res., 90, 11229-11242, 1985. Blackwell. D.D., and J.L. Steele, A summary of heat flow studies in the Cascade Range, Geoth. Resources Council Trans., 7, 233-236, 1983. Bodvarsson, G.S., S.M. Benson, and P.A. Witherspoon, Theory of the development of geothermal systems charged by v e r t i c a l f a u l t s , Jour. Geophys. Res., 87, 9317-9328, 1982. Bodvarsson, G.S., and K. Pruess, Modeling studies of geothermal systems with a free water surface, Stanford Ninth Annual Workshop on Geothermal Resources Engineering, 1983. Bortolami, G.S., B. R i c c i , G.G. Susella, and G.M. Zuppi, Hydrogeochemistry of the Corsaglia Valley, Maritime Alps Piedmont I t a l y , Jour, of Hydrol., 44, 57-79, 1979. 188 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 , Idaho, Jour. Geophys. Res., 86, 11709-11734, 1981. Bryson, R.A., and F . K . Hare, Climates of North America, E l s e v i e r , 1974. Buntebarth, B., Geothermics, Springer-Verlag, 1984. Chapman, D.S., and L. Rybach, Heatflow anomalies and their interpretation, Jour, of Geodynamics, 4, 3-37, 1985. Cheng, P., Heat transfer in geothermal systems, Adv. in Heat Transfer, Vol. 14, 1-105, 1978. Domenico,P.A. and V.V. Palciauskas, Theoretical analysis of forced convective heat transfer in regional ground-water flow, Geol. Soc. Amer., B u l l e t i n 84, 3303-3814, 1973. Elder, J.W., Steady free convection in a porous medium heated from below, J. F l u i d Mech., 27, 29-48, 1967. Faust, C.R., J.W. Mercer, S.D. Thomas, and W.P. Balleau, Quantitative analysis of ex i s t i n g conditions and production strategies for the Baca Geothermal System, New Mexico, Water Resour. Res., 20(5), 601-618, 1984. Fox, F.M., The Simplon tunnel, Minutes of the Proceedings of the Inst, of C i v i l Engineers, 168, 61-86, 1907. Freeze, R.A., and J.A. Cherry, Groundwater, Prentice-Hall Inc., 1979. Freeze, R.A., and P.A. Witherspoon, Theoretical analysis of regional groundwater flow: 2. Ef f e c t of water-table configuration and subsurface permeability v a r i a t i o n , Water Resour. Res., 3, 623-634, 1967. Frind, E.O., Simulations of long-term transient density-dependant transport in groundwater, Adv. Water Resourc, 5, 73-88, 1 982. Garven, G. and R.A. Freeze, Theoretical analysis of the role of groundwater flow in the genesis of stratabound ore deposits: 2. Quantitative r e s u l t s , Am. Jour. S c i . , 284, 1125-1174, 1984. Gosnold, W.D., Heat flow and ground water flow in the great plains of the United States, Jour, of Geodynamics, 4, 247-264, 1985. Goyal, K.P. and D.R. Kassoy, Fault zone controlled charging of a liquid-dominated geothermal reservoir Jour. Geophys. Res, 85, 1867-1875, 1980. Halstead, E.C., Groundwater investigation, Mount Kobau, B r i t i s h Columbia, Canada Inland Waters Branch, Dept. of Energy Mines and Resources, Tech. B u l l . 17, 1969. Hennings, F., On the question of long railway tunnels Construction, v e n t i l a t i o n and operation, B u l l , of the Internl. Railway Congress Assoc., 24(10), 943-983, 1910 Home, R.N., and M.J. O'Sullivan, Numerical modeling of a desaturating geothermal reservoir, Numerical Heat Transfer, 1, 203-216, 1978. Huyakorn, P.S., and G.F. Pinder, Computational methods in subsurface flow, Academic Press, 1983. Ingebritsen, S.E., and M.L. Sorey, A quantitative analysis of the Lassen hydrothermal system, North Central C a l i f o r n i a , Water Resour. Res., 21(6), 853-868, 1985. Jamier, D., Etude de la f i s s u r a t i o n , de 1'hydrogeologie et de l a geochemie des eaux profondes des massifs de l ' A r p i l l e et du Mont Blanc, PhD. d i s s e r t a t i o n , Faculte des Sciences, Universite de Neuchatel, Switzerland (in French), 1975. 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. Keays, R.N., Construction methods in the Moffat Tunnel, Amer. Soc. C i v i l Eng., 62, 63-112, 1928. Keenan, J.H., F.G. Keyes, P.G. H i l l , and J.G. Moore, Steam Tables, 162 pp., John Wiley, 1978. K i l t y , K., D.S. Chapman, and C.W. Mase, Forced convective heat transfer in the Monroe Hot Springs geothermal area Jour, of Vole, and Geoth. Res., Vol. 6, 257-272, 1979. 190 Kimura, S., and A. Bejan, The "heatline" v i s u a l i z a t i o n of convective heat transfer, J. of Heat Transfer, Trans, of the ASME, Vol. 105, 917-919, 1983. Lahsen, A., and P. T r u j i l l o , The geothermal f i e l d of E l Tatio, Chile, Proc. 2nd. U.N. Symp. on Development and Use of Geothermal Resources, San Francisco, C a l i f . , Vol. 1, 157-175, 1975. Lauscher, F., Wettweite typen der hohenabgangigkeit des Niederschlags, Wetter U. Leben., 28, 80-90, 1976. Lewis, T.J., A.M. Jessop, and.A.S. Judge, Heat flux measurements in Southwestern B r i t i s h Columbia: the thermal consequences of plate tectonics, Can. Jour. Earth S c i . , 22, 1262-1273, 1985. Lowell, R.P., C i r c u l a t i o n in fractures, hot springs, and convective heat transport on mid-ocean ridge crests, Geophys. J. Roy. Astr. S o c , 40, 351-365, 1975. Majorowicz, J.A., F.W. Jones, H.L. Lam, and A.M. Jessop, T e r r e s t r i a l heat flow and geothermal gradients in r e l a t i o n to hydrodynamics in the Alberta Basin, Canada, Jour. Geodynamics, 4, 265-283, 1985. Martinec, J., H. Oeschager, U. Schotterer, and U. Siegenthaler, Snowmelt and groundwater storage in an alpine basin, in, Hydrological Aspects of Alpine and High Mountainous Terrain, ed. J . W. Glen, IAHS Publ. 138, 169-175, 1982. Mase, C.W., J.H. Sass, A.H. Lachenbruch, and J.R. Munroe, Preliminary heat-flow investigations of the C a l i f o r n i a Cascades, U.S. Geol. Survey, Open-File Report 82-150. 1982 Mears, F., The eight mile Cascade Tunnel, Great Northern Railway Part I I . Surveys, construction methods and a comparison of routes, Amer. Soc. C i v i l Eng., Trans., 96, 915-1004, 1932. Neuman, S.P., and P.A. Witherspoon, F i n i t e element method of analyzing steady seepage with a free-surface, Water Resour. Res., 6(3), 889-897, 1970. Neuzil, C.E., Groundwater flow in low-permeability environments, Water Resour. Res., 22(8), 1163-1196, 1986. Norton, D., and J. Knight, Transport phenomena in hydrothermal systems: Cooling plutons, Amer. Jour, of S c i . , 277, 937-981, 1977. Patankar, S.V., Numerical heat transfer and f l u i d flow, McGraw H i l l , 1980. Paterson, W.S.B., The physics of g l a c i e r s , Pergamon, 1981. Prats, M., The effect of horizontal f l u i d flow on thermally induced convection currents in porous mediums, Jour. Geophys. Res., 71(2), 4835-4837, 1966. Province of B r i t i s h Columbia, Groundwater observation wells of B r i t i s h Columbia, B r i t i s h Columbia Water Resources Service, Water Investigations Branch, 1974. Raithby, G.D., and K.E. Torrance, Upstream weighting differencing schemes and their application to e l l i p t i c problems involving f l u i d flow, Comput. Fluids, 2, 191-206, 1974. Reader, J.F., and Fairbank, B.D., Heat flow in the v i c i n i t y of the Meager Volcanic Complex, Southwestern B r i t i s h Columbia, Geoth. Resour. Council Trans., Vol. 7, 535-539, 1983. Ross, B., A conceptual model of deep unsaturated zones with negligible recharge, Water Resour. Res., 20(11), 1627-1629, 1984. Ruffner, J.A., and F.E. B l a i r , Climates of the States with current tables of normals 1941-1970 and means and extremes to 1975, National Ocean and Atmospheric Adimin., Vols. I and II, 1975. Sass, J.H., A.H. Lachenbruch, R.J. Munroe, G.W. Greene, and T.H. Moses, J r . , Heat tflow in the Western United States, J. Geophys. Res., 76(26), 6376-6413, 1971. Sauty, J.P., A.C. Gringarten, H. Fabris, D. Thiery, A. Menjoz, and P.A. Landel, Sensible energy storage in aquifers, 2, F i e l d experiments and comparisons with theoreti c a l r e s u l t s , Water Resour. Res., 18(2), 245-252, 1982. Schadt, M.H., Les resultats s c i e n t i f i q u e s du percement du tunnel du Simplon, B u l l e t i n Technique de la Suisse Romande, 125-178, 1905. S c h a r l i , U., and L. Rybach, On the thermal conductivity of low-porosity c r y s t a l l i n e rocks, Tectonophysics, 103, 307-313, 1984. Schermerhorn, V.P., Relations bewteen topography and annual p r e c i p i t a t i o n in Western Oregon and Washington, Water Resour. Res., 3(3), 707-711, 1967. Sklash, M.G., and R.N. Farvolden, The role of groundwater in storm runoff, Jour, of Hydrol., 43, 45-65, 1979. Slaymaker, H.O., and L.J. Zeman, Influence of a l t i t u d e and continentality on watershed hydrology in the Coast Mountains of B r i t i s h Columbia, Proc. Canadian Hydrology Symposium, National Research Council, 1975. Smart, C.C., The hydrology of the Castleguard Karst, Columbia I c e f i e l d s , Alberta, Canada, A r c t i c and Alpine Research, 15(4), 471-486, 1985. Smith, L., and D.S. Chapman,The influence of water table configuration on the near-surface thermal regime, Jour, of Geodynamics, 4, 183-198, 1985. Smith, L., and D.S. Chapman, On the thermal e f f e c t s of groundwater flow 1. Regional scale systems, Jour. Geophys. Res., 88(B1), 593-608, 1983. Smith, M.W., Microclimatic influences on ground temperatures and permafrost d i s t r i b u t i o n , Mackenzie Delta, Northwest T e r r i t o r i e s , Canadian Journal of Earth Science, 12, 1421-1438, 1975. 193 Sorey, M.L., Evolution and present state of the hydrothermal system in Long Valley Caldera, Jour. Geophys. Res., 90, 11219-11228, 1985a. Sorey, M.L., Types of hydrothermal convection systems in the Cascade Range of C a l i f o r n i a and Oregon, Proc. Workshop on Geothermal Resources of the Cascade Range, ed. M. Guffanti and L.J.P. Muffler, 63-67, 1985b. Sorey, M.L. Numerical modeling of l i q u i d geothermal systems, U.S. Geol. Surv. Prof. Paper 1044-D, 1978. Souther, J.G. Geothermal potential of Western Canada, Proc. 2nd. U.N. Symp. on Development and Use of Geothermal Resources, San Francisco, C a l i f . , Vol. 1, 259-267, 1975. Souther, J.G., and Halstead, E.C., Mineral and thermal waters of Canada, Canada Dept. of Energy, Mines and Resources, Paper 73-18, 1973. Spalding, D.B., A novel f i n i t e difference formulation for d i f f e r e n t i a l equations involving both f i r s t and second derivatives Int. J. Numer. Methods Eng., 4, 551-559, 1972. Steele, J. L., and D.D. Blackwell, Heat flow in the v i c i n i t y of the Mount Hood Volcano, Oregon, in Geology and geothermal resources of the Mount Hood area Oregon, Oregon Dept. of Geol. and Min. Ind., Special Paper 14, pp 31-42, 1982. Storr, D., and H.L. Ferguson, The d i s t r i b u t i o n of p r e c i p i t a t i o n in some mountainous Canadian watersheds, in, D i s t r i b u t i o n of P r e c i p i t a t i o n in Mountainous Areas, Geilo Symposium, Norway, Proc. World Meteor. Assoc., WMO/OMM No. 326, Vol. II, 243-263, 1972. Thompson, W.T., How and why to di s t i n g u i s h between mountains and h i l l s , Prof. Geographer, 16, 6-8, 1964. van der Kamp, G., Interactions betwee heat flow and groundwater flow - A review, Proj. 109-17, Waterloo Res. Inst., Univ. of Waterloo, Waterloo, Ont., 1982. V e r r u i j t , A., Theory of Groundwater Flow, Macmillan Press, 1970. Walsh, J.B., and E.R., Decker, Effect of pressure and saturating f l u i d on the thermal conductivity of compact rock, Jour. Geoph. Res., 71, 3053-3061, 1966. Waring, G.A., Thermal springs of the United States and other countries of the world - a summary, U.S. Geol. Survey Prof. Paper 492, 1965. Watson, J.T.R, R.J. Basu, and J.V. Sengers, An improved representative equation for the dynamic v i s c o s i t y of water substance, Jour. Phys. Chem. Ref. Data, 9(3), 1255-1279, 1981. Woodbury, A.D., and L. Smith, On the thermal e f f e c t s of three-dimensional groundwater flow, Jour. Geophys. Res., 90, 759-767, 1985. Zablocki. C.J., R.I. T i l l i n g , D.W. Peterson, R.L. Christiansen, G.V. K e l l e r , and J.C. Murray, A deep research hole at the summit of an active volcano, Kilauea, Hawaii, Geoph. Res. Lett., 1(7), 323-326, 1974. APPENDIX I NOMENCLATURE width of fracture zone L sp e c i f i c heat capacity of f l u i d L 2 / t 2 T conduction dispersion tensor for L 2/T f l u i d in rock matrix conduction dispersion tensor for L 2/T f l u i d in fracture zone thermal lapse rate T/L gravitational constant L / t 2 basal heat flux M/t 3 equivalent freshwater head L available i n f i l t r a t i o n rate L/T i n f i l t r a t i o n r a t i o dimensionless hydraulic conductivity L/t reference hydraulic conductivity L/t permeability of basal unit L 2 permeability of fracture f i l l i n g L 2 permeability tensor for porous L 2 medium permeability of upper unit L 2 horizontal distance between L v e r t i c a l boundaries c h a r a c t e r i s t i c length for L individual f i n i t e element pathline length L mass flux normal to upper boundary M/L 2t porosity dimensionless porosity of basal unit dimensionless porosity of fracture material dimensionless 196 u Pe P Q AC * Q ai Ra s S T T. W WT, max WT max x x, outward normal to free-surface porosity of upper unit Peclet Number f l u i d pressure t o t a l mass flow percentage change in t o t a l mass flow dimensionless t o t a l flow magnitude of f l u i d flux f l u i d flux vector f l u i d flux normal to the free-surface f l u i d flux along fracture zone horizontal component of f l u i d flux v e r t i c a l component of f l u i d flux Rayleigh Number coordinate d i r e c t i o n p a r a l l e l to fracture element degree of saturation temperature reference surface temperature c h a r a c t e r i s t i c t r a n s i t time along pathline width of boundary value problem normal to page elevation of highest point on water table change in elevation of highest point on water table horizontal coordinate x-coordinate at l e f t boundary x-coordinate at right boundary v e r t i c a l coordinate (elevation) dimensionless dimensionless dimensionless ML/t 2 M/t dimensionless dimensionless L/T L/T L/T L/T L/T L/T dimensionless L dimensionless T T L L L L v e r t i c a l coordinate at the upper L boundary of domain for s p e c i f i e d x position v e r t i c a l coordinate at the water L table for specified x position longitudinal thermal d i s p e r s i v i t y L transverse thermal d i s p e r s i v i t y L thermal expansivity for water 1/T thermal conductivity tensor for ML/t 3T sol i d - v a p o r - f l u i d composite thermal conductivity tensor for ML/t 3T s o l i d f l u i d thermal conductivity ML/t 3f vapor thermal conductivity ML/t 3T f l u i d density M/L3 f l u i d density at s p e c i f i e d M/L3 reference temperature r e l a t i v e f l u i d density M/L3 orientation of tangent to water radians table measured from horizontal dynamic f l u i d v i s c o s i t y M/Lt 198 APPENDIX II DISTINCTION BETWEEN POINT OF DETACHMENT AND HINGE POINT ON A SEEPAGE FACE Conventional approaches to solving free-surface groundwater flow problems assume that a single point, the exit point, marks the boundary between the free-surface and the seepage face. Two conditions occur at t h i s point; the free-surface deviates from the seepage face and the upper l i m i t of discharge on the seepage face i s defined. These conditions occur at a single point for free-surface problems similar to the triangular earth dam shown in Figure I I . 1. In t h i s isothermal example, a homogeneous isotropic hydraulic conductivity K i s assumed for a l l panels shown in Figure II.1. The horizontal base CD i s impermeable while segments BC and DE are constant head boundaries with head d i f f e r e n t i a l A/i = hj - h2* The form of free-surface AB, the length of the seepage face AE and the position of the exit point A are controlled by th i s head d i f f e r e n t i a l , the geometry of the dam, the hydraulic conductivity K and the pattern of i n f i l t r a t i o n applied on the upper surface of the dam. In a s i m p l i f i e d mountain groundwater flow problem, the v e r t i c a l lefthand boundary of Figure II.1a becomes a symmetry boundary (Figure II.1b). At a high uniform i n f i l t r a t i o n rate (IzQ) the water table coincides with the ground surface across the flow system. In th i s case, a free-surface i s absent and the entire upper surface can be considered a seepage face. The hinge point (HP) at point E on the seepage face marks the boundary between recharge and discharge and f u l f i l l s one condition of an exit point by defining the upper l i m i t of discharge. The point of detachment A i s undefined in Figure II.1b because a free-surface i s absent under conditions where the water table i s everywhere at the bedrock surface. Reducing the i n f i l t r a t i o n rate to 199 (7 Z) causes a free-surface to develop (Figure I I . 1c) with a point of detachment (POD) at A and a hinge point (HP) at E. In t h i s case, a single point cannot be defined that f u l f i l l s the d e f i n i t i o n of an exit point and recharge may occur on the seepage face between the POD and HP. Further reducing the i n f i l t r a t i o n rate decreases the separation between POD and HP (A and E). In low-relief topography, the separation becomes ne g l i g i b l e and the usual exit point d e f i n i t i o n i s v a l i d . In h i g h - r e l i e f mountainous topography, the separation between POD and HP can be substantial and must be considered in the numerical formulation. Bear (1972, pg. 272) notes that mapping free-surface groundwater flow problems into the hodograph plane i s a useful method for examining flow conditions at the boundaries. While the boundary between the free-surface and the seepage face i s i n i t i a l l y unknown in the physical plane, in the hodograph plane i t i s completely defined. This mapping procedure provides an a n a l y t i c a l argument that supports the concept of separated POD and HP on seepage faces in mountainous t e r r a i n . D etails on methods of mapping from the physical to hodograph planes can be found in Bear (1972) and V e r r u i j t (1970). Figure II.1d is the hodograph representation of the physical system of Figure I I . 1c. In the hodograph plane, v e r t i c a l and horizontal components of f l u i d flux (qx, qz) at each point on the boundary of the physical system form the hodograph. The outline of the hodograph i s defined as the locus of points marking the d i s t a l ends of s p e c i f i c discharge vectors o r i g i n a t i n g at point C (see vectors F 1, F2 and F3 shown in Figures I I . 1c and II.1d). In the mapping process, however, the s p a t i a l relationships between adjacent points in the physical plane are no longer defined in the hodograph plane. Along the horizontal impermeable boundary CD, qz i s zero and q% increases with increasing distance from the o r i g i n to a maximum value equal to the hydraulic 200 conductivity K at point D. Si m i l a r l y , along the v e r t i c a l impermeable boundary CB, qx i s zero and the magnitude of qz increases to a maximum with absolute value equal to the v e r t i c a l i n f i l t r a t i o n rate / at B. In the hodograph plane, the free-surface is described by a c i r c u l a r arc with radius (K - I )/2 (Bear, 1972). The highest point on the free-surface occurs at the intersection with the v e r t i c a l flux axis (point B) while the lowest point occurs at the POD (point A) where the free-surface intersects the upper boundary of the domain. Three f l u i d flux vectors F 1, F 2 and Fg are shown in Figure I I . 1c with their corresponding hodograph representation in Figure II.1d. In the physical plane, vector F^ i s directed outward indicating discharge across the seepage face, while vector F^ i s directed inward and indicates recharge. At the hinge point, vector F 2 is p a r a l l e l to the upper boundary in the physical plane and perpendicular to the corresponding l i n e AD in the hodograph plane (Figure II.Id). Because F 2 i s p a r a l l e l to the upper physical boundary, f l u i d flux normal to the boundary at point E i s zero. Therefore, point E i s the hinge point that marks the boundary between recharge and discharge on the seepage face. Increasing the v e r t i c a l i n f i l t r a t i o n rate from I 2 to Izj causes the HP to move upslope from E2 to Ej and the POD to move upslope from A2 to Aj in the physical plane of Figure II.1e. The HP remains fixed in the hodograph plane (points Ej and E2), because t h i s point always marks the point of zero flux normal to the upper boundary in the hodograph representation (Figure 11.1 f ) . Separation between POD and HP increases with increasing i n f i l t r a t i o n in the physical plane (Figure II.1e). Although the physical distance between POD and HP cannot be defined in the hodograph plane, i t seems l o g i c a l to assume that increasing separation in the hodgraph plane corresponds to increasing separation in the physical plane (Figure II.1e). In the simple systems shown in Figures I I . 1c and II.1e, the horizontal and v e r t i c a l components of f l u i d flux are e x p l i c i t l y defined at points A, B, C, D and E in the hodograph plane (Figures II.1d and II. 1 f ) . Hydraulic head solutions obtained with the f i n i t e element method are used to estimate the free-surface configuration and calculate boundary fluxes. Boundary fluxes calculated with the f i n i t e element model correspond well to those defined in the hodograph plane. 202 z . Figure I I . 1. Mapping water table configurations and flow patterns from physical to hodograph planes; a. triangular dam, b. idealized f u l l y saturated mountain p r o f i l e , c. mountain free-surface problem with i n f i l t r a t i o n rate 7 , d. hodograph of (c), e. mountain free-surface problem with i n f i l t r a t i o n rates I z ] and I z 2 , f. hodograph of (e). 

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-0052356/manifest

Comment

Related Items