UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Limnology of the Meromictic Island Copper Mine pit lake Fisher, Timothy Simon Richmond 2002

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

Item Metadata

Download

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

Full Text

L I M N O L O G Y OF THE MEROMICTIC ISLAND COPPER M I N E PIT L A K E by TIMOTHY SIMON RICHMOND FISHER B.E., The University of Canterbury, 1993 M.E. , The University of Canterbury, 1995 A THESIS SUBMITTED IN PARTIAL F U L F I L M E N T OF THE REQUIREMENTS FOR THE D E G R E E OF DOCTOR OF PHILOSOPHY IN THE F A C U L T Y OF G R A D U A T E STUDIES (DEPARTMENT OF CIVIL ENGINEERING) We accept this thesis as conforming to the requited standard THE U N I V E R S I T Y ^ BRITISH C O L U M B I A April 2002 ©Tim S.R. Fisher, 2002 In p resen t ing this thesis in partial fu l f i lment of t h e requ i remen ts f o r an advanced degree at the Univers i ty o f Brit ish C o l u m b i a , I agree that t h e Library shall make it f reely available f o r re ference and study. I fu r ther agree that pe rmiss ion f o r extens ive c o p y i n g o f th is thesis f o r scholar ly pu rposes may be g ran ted by the head o f my d e p a r t m e n t or by his o r her representat ives. It is u n d e r s t o o d that c o p y i n g or pub l i ca t i on o f th is thesis fo r f inancial ga in shall n o t b e a l l o w e d w i t h o u t m y w r i t t e n permiss ion . D e p a r t m e n t of Civil £n#j tue&n'/iCPJ The Univers i ty o f Brit ish C o l u m b i a Vancouver , Canada Date 25 fifril 2O0Z DE-6 (2/88) ii ABSTRACT The BHP Billiton owned Island Copper Mine pit near Port Hardy, Vancouver Island, BC, Canada, was flooded in 1996 with seawater and capped with fresh water to form a meromictic pit lake (maximum depth 350 m and surface area 1.71 km ). The pit lake is being utilized as a passive treatment system for acid rock drainage (ARD). Conductivity-temperature-depth profiles (July 1997 to March 2001) document the evolving structure of the pit lake. Intensive measurements with two thermistor "chains" and a meteorological station (December 1999 to November 2000) enabled physical processes to be investigated. The pit lake has developed into three distinct layers: a brackish and well mixed upper layer; an intermediate layer circulated by ARD plumes; and a thermally convecting lower layer. The upper halocline has risen due to the injection of buoyant ARD into the base of the intermediate layer and entrainment of the halocline by the impinging ARD plume. The lower halocline fluctuated seasonally due to the balance between ARD plume erosion and lower layer thermal convection. The true salinity was determined from the Practical Salinity Scale 1978 plus a time and layer dependent correction for non-seawater ions sourced from the ARD and submerged pit walls. The plume stirred intermediate layer was analogous to the filling box model (Baines and Turner, 1969), which explains the observed salinity and counter-stable temperature gradients. Lake length and width internal waves were identified at the upper halocline and their periods agreed with predictions. These internal waves appeared to degenerate by damping, not by non-linear steepening as predicted by theory. Conceptual models used to model the long-term evolution of the pit lake identified upwelling and/or erosion of the intermediate layer as potential failure modes. Clear from the conceptual models was the importance of maintaining the density stratification for the long-term stability of the meromictic structure. Detailed modeling with a one-dimensional process model DYRESM replicated some details of the pit lake, but did not adequately model wind mixing, nor did it allow for plume entrainment of the upper halocline. iii TABLE OF CONTENTS Abstract ii Table of Contents iii List of Tables vi List of Figures vii List of Abbreviations xi Preface xii Acknowledgements xiii 1 Development of the Island Copper Mine Meromictic Pit Lake for the Treatment of Acid Rock Drainage 1 1.1 Introduction 1 1.2 Methods 5 1.3 Results 6 1.4 Discussion 9 1.5 Conclusions 14 2 Determination of Salinity and Density 30 2.1 Introduction 30 2.2 Application of Theory 32 2.3 Methods 36 2.4 Results 37 2.5 Discussion 41 2.6 Conclusions 44 3 Plume Discharge and Filling Box Behaviour 56 3.1 Introduction 56 3.2 Literature Review 58 3.3 Methods 62 iv 3.4 Results 64 3.5 Discussion 67 3.6 Conclusions 71 4 Internal Waves at the Upper Halocline 83 4.1 Introduction 83 4.2 Theory 86 4.3 Site Description and Method 92 4.4 Results 93 4.5 Discussion 98 4.6 Conclusions 101 5 Modeling the Long-term Evolution 117 5.1 Introduction 118 5.2 Lake Models 120 5.3 Upwelling Model 124 5.4 Results 127 5.5 Discussion 130 5.6 Conclusions 134 6 Conclusions and Recommendations 145 6.1 Conclusions 145 6.2 Recommendations 148 References 150 Appendices A Water Balance 161 A.l Introduction 161 A.2 Method 162 A.3 Results 163 A. 4 Plume Entrainment 164 B Island Copper Meteorology 173 B. l Introduction 173 B.2 Regional Weather 173 B.3 Island Copper Weather 175 B.4 Converting Island Copper Wind Speed 175 B.5 Port Hardy Maximum Annual Wind Frequency Analysis 176 B.6 Comparison of Island Copper and Port Hardy Wind 177 B. 7 Conclusions 180 C Conductivity Temperature Depth (CTD) Measurements 194 C. l Introduction 194 C.2 OS200 Performance 194 C.3 OS200 Operation 195 C. 4 OS200 Calibration 196 D Thermistor Chain Mooring Design and Thermistor Specification and Calibration 201 D. l Introduction 201 D.2 Thermistor Chain Design 202 D.3 Thermistors Chain Deployment and Recovery 203 D.4 Thermistor Instruments 203 D. 5 Thermistor Calibration 204 E Meteorological Measurements 218 E. l Introduction 218 E.2 Meteorology Station Setup 218 E.3 Meteorology Station Instrumentation 219 E.4 Meteorology Station Calibration 220 E.5 Port Hardy Airport Meteorological Station 220 F Photographs 222 VI LIST OF TABLES 1.1 Island Copper pit lake volume and area 15 2.1 Comparison of density from PSS78 and EOS80 to density measured directly 46 2.2 Alternative salinity measurements of November 29 2000 pit lake samples 47 2.3 Coefficient of haline contraction 48 2.4 Coefficient of thermal expansion 49 2.5 Major chemical constituents and excess concentration 50 3.1 Classification of entrainment regime by Cotel and Breidenthal (1997) 73 4.1 Observed internal waves compared to Defant and Merian predictions 103 4.2 Parameterization of Island Copper pit lake to wind disturbance 104 5.1 Annual upper halocline depths from CTD profiles 135 A. l Halocline depths and changes in position 167 A.2 Interface volume fluxes 168 C. l Summary of CTD profiles using the OS200 unless noted otherwise 198 D. l Specification for Thermistors 206 D.2 Total Instrument Summary (working) 206 D.3 Central Mooring Thermistor Chain - Instrument Deployment Summary 207 D.4 South Injection System Thermistor Chain - Instrument Deployment Summary 210 D. 5 Thermistor Calibration Summary 212 E. l Meteorological station instrumentation and deployment periods 221 vii LIST OF FIGURES 1.1 Oblique aerial photograph of the Island Copper Mine site and pit lake 16 1.2 Island Copper pit lake bathymetry and cross-section 17 1.3 Vertical profiles showing salinity, temperature, density and dissolved oxygen 18 1.4 Weekly average north and south injection system flow 19 1.5 True salinity and seawater (Cl%o) salinity trends 20 1.6 Temperature trends for upper, intermediate and lower layers 21 1.7 Dissolved oxygen trends for upper, intermediate and lower layers 22 1.8 pH trends for upper, intermediate and lower layers 23 1.9 Copper concentration trends for upper, intermediate and lower layers 24 1.10 Zinc concentration trends for upper, intermediate and lower layers 25 1.11 Time series of salinity profiles through the upper halocline and lower 26 halocline 1.12 Change in intermediate layer volume versus total injection system flow 27 1.13 Lowest measured dilutions in north injection system plume on 9 Oct 1998 ... 28 1.14 Dilution of the south injection system plume on 2 - 3 March 1999 29 2.1 Island Copper pit lake cross-section and profiles of temperature and salinity 51 2.2 Density versus salinity for Island Copper samples 52 2.3 Density of int. layer diluted with injection system and double distilled water 53 2.4 Non-seawater salinity versus days since pit lake was flooded 54 2.5 Salinity trends for upper, intermediate and lower layers 55 3.1 Island Copper pit lake bathymetry and cross-section showing thermistor chain and meterological station locations 74 3.2 Sketch of the development of stratified environment due to buoyant fresh/cool source discharging into a saline/warm receiving environment 75 3.3 Daily average temperature of upper layer, upper interface and top of vm intermediate layer 76 3.4 Intermediate layer profiles of salinity, potential temperature and sigma-0 77 3.5 Intermediate layer density gradient by least squares regression of CTD data . 78 3.6 Daily average temperature of the intermediate layer and total injection 79 system flow and temperature 3.7 Lower layer profiles of salinity, potential temperature and sigma 80 3.8 Plume characteristics (g0', B, ht and tn,) for north and south injection systems 81 3.9 The non-dimensional height versus non-dimensional density gradient in the intermediate layer for winter compared to filling box model 82 4.1 Daily average temperature contours at the upper halocline 105 4.2 Bi-monthly profiles of daily average temperature through the upper halocline 106 4.3 Daily average thermistor temperatures at the central mooring and over the south injection system 107 4.4 Upper thermocline hourly average temperatures at the central mooring and over the south injection system compared to hourly average wind energy and direction and 30 minute north and south injection system flows 108 4.5 Spectrum from thermistors in the halocline at the central mooring and over the south injection system for period 1-18 January 2000 109 4.6 Windrose for 1-18 January 2000 110 4.7 Spectrogram of thermistors in the halocline at the central mooring and over the south injection compared to wind energy and direction I l l 4.8 Detailed temperature at the central mooring and south injection system compared to wind and south injection system flow 112 4.9 Upper thermocline hourly average temperatures compared to hourly average wind speed and direction and 30 minute SIS and NIS flows 113 4.10 Windrose for 10 June - 12 July 2000 114 4.11 Spectrum from thermistors in the halocline at the central mooring and over the south injection system versus for period 10 June - 12 July 2000 115 4.12 Detailed upper thermocline temperatures compared to 15 min average wind energy and direction 116 ix 5.1 Frequency distribution of Island Copper annual maximum hourly wind 136 speeds 5.2 Model 1 results for Run 1 with both NIS and SIS discharging to intermediate layer 137 5.3 Model 1 results for Run 2 with SIS discharging to intermediate layer and NIS to the upper layer 138 5.4 Model 1 predictions of return period for the wind required to cause upwelling and risk of upwelling 139 5.5 Model 2 results for Run 2 with SIS discharging to intermediate layer and NIS to the upper layer 140 5.6 Model 2 results for Run 1 with both NIS and SIS discharging to intermediate layer 141 5.7 Upper layer DYRESM temperature compared to observed temperature 142 5.8 Upper layer DYRESM salinity compared to observed halocline position 143 5.9 Intermediate layer DYRESM temperature compared to observed temperature and injection system flow and temperature 144 A. l Schematic diagram of mixing at the upper interface 169 A.2 Depth of lower halocline, injection system inflow and lower halocline 170 volumetric flux A.3 Depth of upper halocline, injection system inflow and gross and net upper 171 halocline volumetric flux A. 4 Relationships between injection system flow and upper layer volumetric 172 flux, upper layer net volumetric flux and lower layer volumetric flux B. l Southwest British Columbia and region map showing the Island Copper 181 Mine B.2 Island Copper hourly temperature, relative humidity, short-wave radiation, wind speed and direction and rainfall 182 B.3 Island Copper hourly wind in u and v components 183 B.4 Port Hardy hourly wind in u and v components 184 B.5 Frequency distribution of Port Hardy annual maximum hourly wind speeds .. 185 B.6 Hourly wind speed and direction at Island Copper compared to Port Hardy... 186 X B.7 u and v wind component at Island Copper compared to Port Hardy 187 B.8 Correlation between Port Hardy and Island Copper hourly wind speed 188 B.9 Correlation between Port Hardy and Island Copper hourly u wind component 189 B. 10 Correlation between Port Hardy and Island Copper hourly v wind component 190 B.l 1 Probable wind direction at Island Copper compared to Port Hardy 191 B.12 Wind spectrum, coherence and phase for Island Copper and Port Hardy 192 B. 13 Comparison of Island Copper and Port Hardy extreme wind speeds 193 C. l Island Copper pit lake bathymetry showing locations for CTD casts 200 D. l Design drawing of central mooring thermistor chain in elevation view 214 D.2 Design drawing of south injection system thermistor chain mooring in plan view 215 D.3 Design drawing of south injection system thermistor chain mooring showing cross-section AA 216 D.4 Design drawing of south injection system thermistor chain mooring showing cross-section BB 217 F.l Oblique aerial of Island Copper Mine and Rupert Inlet (August 1993) 222 F.2 Flooding the Island Copper Mine pit with seawater (June 1996) 223 F.3 Island Copper Mine pit post closure (9 April 1996) 224 F.4 Island Copper Mine pit during flooding (17 Julyl996) 224 F.5 Island Copper Mine pit at end of flooding (29 July1996) 225 F.6 North injection system outfall 225 F.7 South injection system outfall 226 F.8 South injection system sedimentation pond and outfall intake 226 F.9 UBC "Terra" meteorological station located near the south injection system 227 F.10 OS200 conductivity temperature depth profiler 227 F . l l South injection system and central mooring thermistors 228 F.l2 Deploying south injection system thermistor chain 228 xi LIST OF ABBREVIATIONS ARD acid rock drainage CMTC central mooring thermistor chain EOS80 1980 Equation of State (UNESCO 1981) NIS north injection system PSS78 1978 Practical Salinity Scale (Lewis and Perkin 1978, Perkin and Lewis 1980) SIS south injection system SISTC south injection system thermistor chain UWL datum used for pit lake depths representing the design ultimate water level of 2.44 m above mean sea level (1008 ft level in pit) xii PREFACE The main body of the thesis is comprised of five chapters, each of which has been developed as a paper intended for journal publication. Each chapter has a different focus regarding the limnology of the lake and therefore stands alone. Repetition in the introduction and methods of each chapter has been minimized by references to the appropriate chapter. The numbering of pages, figures, equations and tables within each paper have been changed for thesis continuity. As well the acknowledgements and references of individual papers have been collated. The material presented in each of the papers are my original ideas and credit to co-authors at the time of publication will reveal the advice and critical review of my supervisor and collaborators. The thesis abstract, conclusions and recommendations summarize the collective work. Chapter 1 is an introductory paper, which gives background on the Island Copper pit lake and details the development of the lake structure and water quality. Chapter 2 establishes the unique temperature-salinity-density relationship for each layer in the pit lake. Chapter 3 details the filling box behaviour of the ARD plumes discharged to the intermediate layer, and Chapter 4 focuses on the upper halocline with regard to the plume discharge and wind generated internal waves. Chapter 5 develops models to better understand the long-term stability of the pit lake. Appendices cover details not appropriate for the papers, but essential to a thorough thesis: water balance, weather and instrumentation details. Xlll ACKNOWLEDGEMENTS I would like to thank Chris Hanks, Brian Welchman and Eric Schwamberger of BHP Billiton for financial and logistic support, as well as access to the site and data. BHP's Annual Environmental Assessment Data Reports 1998-2000 were an invaluable resource from which water quality data was complied. I would also like to thank Steve and Trudy Lacasse and Andy Hanke of BCL Biotechnologies Ltd. for their field assistance and hospitality, and acknowledge their sampling and laboratory analysis of water quality data. Clem Pelletier, George Poling, Deborah Muggli, Marc Wen and Shane Uren of Rescan Environmental Services Ltd. supplied invaluable technical assistance by way of their own plume tracer studies, bathymetric mapping and submersible inspection of the outfall systems. The University of British Columbia provided me with a University Graduate Fellowship. Thanks are extended to Roger Pieters for his general support and assistance with instrumentation and to Mike Wilton for his initial research. Field help and moral support was provided by Dan Walker, Colin Rennie and Violeta Martin. Stephen Pond provided help with the moorings and field trips and provided careful review of manuscripts. Thanks to committee members Ken Hall, Peter Ward and Noboru Yonemitsu and also to Craig Stevens for their reviews. My gratitude to Greg Lawrence for introducing me to such an interesting topic and as my supervisor giving me the independence to do what I wanted while still providing thoughtful advice and constructive criticism. Special thanks to my family for their support, and to my wife Anna, for following me to Canada and for her love and never ending patience. 1 CHAPTER 1 Development of the Island Copper Mine Meromictic Pit Lake for the Treatment of Acid Rock Drainage ABSTRACT The physical structure and water quality of the meromictic pit lake created at the Island Copper Mine near Port Hardy, Vancouver Island BC in 1996, has developed into three distinct layers: a brackish and well mixed upper layer; a plume stirred intermediate layer; and a thermally convecting lower layer. To date the upper halocline has risen due to the injection of buoyant acid rock drainage (ARD) into the base of the intermediate layer. An upper layer "equilibrium" depth is expected in the near future reflecting a balance between meteorological forcing from above and rising ARD plumes from below. The lower halocline depth fluctuates seasonally: deepening in the winter due to high ARD flows, which erode the lower halocline; and rising in the summer when thermal convection in the lower layer dominates over reduced intermediate layer stirring. The proposed treatment of acid rock drainage by metal-sulphide precipitation using sulphate-reducing bacteria is dependent on anoxia developing in the intermediate and lower layers. By March 2001, the lower layer was close to anoxia with dissolved oxygen (DO) = 0.03 mg L"1, whereas, the intermediate layer with DO = 2.5 mg L"1 still contained some oxygen. 1.1 INTRODUCTION The Island Copper Mine pit lake is both a practical solution to an environmental problem and a unique limnological experiment. A meromictic (permanent density stratification) pit lake has been purposefully created for the treatment of acid rock drainage 2 (ARD). ARD results from the chemical and biological oxidation of reactive sulphide minerals when exposed to air and water. The pit lake contains three distinct layers and two buoyant plumes, which makes it a truly remarkable waterbody. Pit lakes - Pit lakes form when abandoned mine pits fill with rain water, groundwater, or are deliberately flooded. Depending on the mining environment, pollution in the form of ARD can be a concern, and pit lakes are often polluted with heavy metals and high acidity (Geller et al. 1998; Watkins, 2000). Vertical chemical gradients are often observed in pit lakes and may cause meromixis. The Brenda Mines pit lake, Canada, is meromictic from chemical stratification (Stevens and Lawrence 1998), and pit lakes in Merseburg-Ost, Germany are meromictic from salt laden groundwater inflow (Bohrer et al. 1998). Meromixis in the Island Copper mine pit lake was caused by flooding the pit with seawater and capping it with fresh water. The Island Copper mine pit lake has two sharp haloclines with salinity changes of about 20 at the upper interface and 3 at the lower interface. Thus meromictic conditions are similar to evaporative basins such as Mahoney Lake with salinities of 4-40%o (Northcote and Hall 1983; Ward et al. 1990), but not as extreme as Mono Lake with salinities of 76-87%o (Jellison and Melack 1993; Romero and Melack 1996; Jellison et al. 1998). Anoxic basins such as Powell Lake (Sanderson et al. 1986), Lake Nitinat (Ozretich 1975) and Sakinaw Lake (Perry and Pederson 1993), all in British Columbia, and marine anoxic basins such as the Black Sea and Cariaco Trench (Millero and Sohn, 1992) are related because of their seawater origins. Common to the anoxic monimolimnion (lower unmixed layer) of these systems are sulphate-reducing bacteria, which are anticipated to become part of the passive treatment system in the Island Copper pit lake. The bathymetry of pit lakes resulting from hard rock mining is characterized by high aspect ratios (Stevens and Lawrence 1998). The Island Copper mine pit was certainly deep, the bottom was 350 m below sea level making it among the lowest surface depressions on earth, and the resulting aspect ratio (depth divided by equivalent diameter) for the pit lake is 0.24. Such high aspect ratios are found naturally only in crater lakes (Martin 1985; Chikita et al. 1993, 1995; Crawford and Collier 1997). 3 Unique to the Island Copper pit lake are two buoyant plumes used for the disposal of ARD. In lakes, bubble plumes used for artificial destratification have been researched (Wiiest et al. 1992; Schladow and Fisher 1995) and thermal plumes observed (Colomer et al. 2001). However, the pit lake situation is more like an ocean outfall discharge (Fischer et al. 1979; Wood et al. 1993), but within a confined space like the "filling box" experiments of Baines and Turner (1969) and Kumagai (1984). In filling box experiments a plume rises within a fixed volume until the top is reached, there it spreads horizontally and downward to replace the flow entrained into the plume as it rises. Island Copper Mine History and Closure Plan - The Island Copper Mine, located near Port Hardy on Vancouver Island, is owned by BHP Billiton. The mine site (Figure 1.1; Figure F.l) is adjacent to Rupert Inlet, which connects via Quatsino Sound to the Pacific Ocean. Mining operations, which began in October 1971 and ended in December 1995, produced mostly copper with lesser amounts of molybdenum, silver and gold (Aspinall 1995; BHP 1996). A submarine tailings disposal system discharged mine tailings into Rupert Inlet for 23 years without significant deleterious effect to marine life (Ellis et al. 1995). Waste rock was dumped into Rupert Inlet forming the 260 ha beach dump and on land around the open pit. During the 1980s, ARD was detected from the land waste rock dumps (BHP 1996). Traces of ARD seepage from the beach dump are undetectable away from the face of the dump (Ellis et al. 1995). Tests suggest that the on land dump material will produce moderately acidic leachate and release elevated concentrations of sulphate and metals over the long-term (Dagenais 1996). The Island Copper Mine pit was flooded from 15 June until 23 July 1996 with seawater from Rupert Inlet, and capped with fresh water to form a meromictic lake (Figures F.2 - F.5). The resulting lake has a surface area of 171 ha, maximum depth of ~ 350 m, average depth of 141 m and long axis of 2.35 km (Figure 1.2). The lake reached its design water level - the ultimate water level (UWL) of 2.44 m above mean sea level - during the 1998/99 winter. Depths in the pit lake are measured as positive below the UWL. Lake volume and surface area are detailed in Table 1.1. BHP's Closure Plan for the Island Copper Mine is to use the open pit as a passive treatment system for ARD (BHP 1996). ARD is collected from land waste rock dumps and 4 discharged via two injection systems into the intermediate layer of the pit lake, where it rises as buoyant plumes through the more dense salt water (Figure 1.2). The upper 'brackish' water layer forms a cap over the intermediate layer, preventing the plumes from surfacing and impeding oxygen transfer to the intermediate layer. A distinct lower layer has formed below the injection system discharges. The ARD is currently treated in the intermediate layer by dilution within the large volume and by neutralization by the natural alkalinity in the seawater. It was proposed that removal of heavy metals would be provided by metal-sulphide precipitation via anaerobic sulphate-reducing bacteria (BFfP 1996), once anoxic conditions were established in the lake, although anoxia has not developed at the time of writing. The two injection systems were installed in October 1996. The south injection system (SIS) and north injection system (NIS), see Figure 1.2, drain different areas of the waste rock dump. ARD is collected by ditches and routed to sedimentation ponds (Figure F.5). It is then discharged into the pit lake by static head via two high-density-polyethylene outfall pipelines at a depth of -220 m (Figures F.6 and F.7). The two diffusers consist of 3 sets of 200 mm discharge ports at quarter points around the circumference, spaced along a 14 m pipe section, and the 428 mm open end of the pipe. Inspection of the injection systems by submersible remotely operated vehicle revealed that both diffusers were lying on the pit walls at an angle of 37° from horizontal (Rescan, 2001). The SIS diffuser has only three ports open (out of 12) plus the end of pipe. Furthermore, both diffusers were only discharging from the first set of diffuser ports, which indicates the intrusion of saline pit lake water was causing the diffusers to be partially blocked (Wilkinson, 1985). An old conveyer tunnel discharges drainage into the intermediate layer at a depth of 154 m (UWL). The conveyer tunnel flow is estimated to be 4 % of the total injection system flow, based on a comparison of catchment areas. There are no groundwater inflows to the pit. The upper layer receives fresh water from direct precipitation and drainage from the pit walls. Outflow is only from the surface of the lake (upper layer), through the beach dump as a diffuse subsurface discharge, and eventually into Rupert Inlet 700 m away. Outline - Decommissioning of mine pits and the treatment of ARD are topical subjects, and in the broader context, disposal of waste in stratified fluids is an important 5 environmental issue. This chapter documents background information and details field methods. General results regarding the development of the lake's physical structure and water quality are given. In the discussion, the results are interpreted and issues facing the passive treatment system are explored. 1.2 METHODS Field Measurements and Analysis - Routine water quality data have been collected since the pit lake was flooded. Temperature was obtained by reversing thermometers and water samples obtained by Niskin samplers analyzed for chlorinity, dissolved oxygen (modified Winkler titration), dissolved heavy metals (flame atomic absorption spectrophotometry), pH, alkalinity, and turbidity. Samples were taken at the west station (Figure 1.2) and were the layer averages of one or more samples taken at different depths. Water level was monitored manually with a staff gauge until January 1999, and by an ultrasonic water level meter since then. Injection system inflow measurement started on 27 November 1997 using acoustic flow meters at each of the injection system outfall pipes. Conductivity-temperature-depth (CTD) profiling has been carried out on 18 occasions: July, August, September and November 1997; April, August and December 1998; March, June, September and December 1999; February, March, June, July, August and November 2000; and March 2001. The CTD instrument, an Ocean Sensors' Model OS200 averaged 99 raw measurements every -0.5 s. Laboratory calibration of the CTD was performed before and after each trip. Field verification of salinity was provided by bottle samples analyzed in a Guildline Model 8400 Autosal and by occasional use of supplementary CTD's. The temperature accuracy and resolution were 0.02 °C and 0.001 °C, and for salinity they were 0.03 and 0.001, respectively (see Appendix C). A depth uncertainty of 0.23 m was accomplished by comparing against in situ thermistors and a Sea Bird Electronics SBE19 CTD. A depth resolution of 0.15 m was achieved in interface areas by slowing the CTD descent rate. Salinity - In the pit lake, salinity is the dominant contributor to density differences, so accurate determination of salinity is essential to understanding the physical mixing processes and the long-term stability of the stratification. The accurate determination of true salinity 6 (defined as mass of dissolved salts in one kilogram of solution) is not a trivial matter. Wilton (1998) used the Practical Salinity Scale 1978 (PSS78) (Perkin and Lewis 1980) to describe the pit lake salinity from CTD measurements. However, true salinity in the pit lake is significantly different from PSS78 salinity, as the chemical composition of the originally seawater pit lake has changed due to the addition of ARD, dissolution of ions from the pit walls, and other chemical changes. To calculate the true salinity from CTD measurements a time and layer dependent correction is applied to PSS78 salinity. A third measure of salinity, the seawater salinity determined from chlorinity measurements, is also important as it indicates the salinity contribution from the original seawater. Chapter 2 is dedicated to the issue of salinity determination and details the salinity correction relationships. 1.3 RESULTS The pit lake has evolved into a stably stratified three-layer system. The CTD profile taken 6th December 1999 was typical of other dates, with true salinity, temperature, density and dissolved oxygen being approximately constant in each layer (Figure 1.3). The layers were separated by strong haloclines. ARD flow to the intermediate layer via the NIS and SIS was highly seasonal (Figure 1.4). Flows increase through the fall, reach a peak in December and or January, then drop through the spring to minimal summer values. The annual daily average flows [and peak daily flow] for the NIS and SIS were 7,190 m3 day"1 [75,180 m3 day" '] and 4,650 m3 day"1 [71,980 m3 day"1], respectively. The ARD salinity (PSS78) measured in the SIS pond from 1997 to 2000 had a mean of 0.86, but was variable with a standard deviation of 0.23. Salinity - The pit was initially flooded with seawater from Rupert Inlet. A brackish upper layer formed on the surface of the lake due to precipitation, runoff and additional fresh water (pumped from the nearby Marble River until 15 March 1998). The salinity of the upper layer decreased rapidly at first, due to dilution from the Marble River supply, and more slowly since 15 March 1998 reaching 4.08 in March 2001 (Figure 1.5a). There was a slight seasonal variation in upper layer salinity; winter precipitation diluted the upper layer reducing the salinity, whereas, during the summer (when there was relatively little precipitation) mixing across the upper halocline of intermediate layer water causes the upper layer salinity to increase. 7 The intermediate layer salinity decreased to 26.62 in March 2001 from the original seawater salinity estimated to be 28.87 (Chapter 2), due to the relatively fresh injection system discharge (Figure 1.5b). The most rapid decreases occur in the winter when injection system flows were the greatest. Diversion of the NIS from the intermediate layer to upper layer on 29 June 2000 reduced the dilution of the intermediate layer salinity during the 2000/01 winter. The lower layer salinity increased to 29.78 in March 2001, which is greater than the original pit lake seawater salinity (Figure 1.5b). The seawater salinity, which represents the concentration of original seawater ions (particularly Cl), was diluted over this period confirming the exchange of fluid between the intermediate and lower layers. The effective diffusion rate to cause this decrease in salinity is lxlO"6 m s"1 (based on a lower halocline average salinity gradient of 0.3 nf'), which is 2000 times the molecular diffusivity of salt (NaCl). In all layers there was a growing difference between the true and seawater salinity due to the addition of non-seawater ions, which were primarily ions dissociated from the pit walls and from ARD (Chapter 2). In the lower layer the concentration of non-seawater ions from these sources is exacerbated by lower pH, smaller volume to surface area ratio and the greater availability of loose rock (Chapter 2). In addition, the lower layer may be biochemically active due to impending anoxia and accumulation of biomass from the upper layer. Temperature - The temperature of the upper layer followed the seasonal cycling of temperature with summer highs up to 24 °C and winter lows down to 4 °C (Figure 1.6a). The intermediate and lower layers remained between 12 and 13 °C (Figure 1.6b). The temperature of the intermediate layer was generally decreasing with seasonal variation due to winter influx of cold ARD and summer solar heating. The temperature in the lower layer rose asymptotically to -12.9 °C, probably as a result of geothermal heating (Figure 1.6b). Dissolved Oxygen - The dissolved oxygen (DO) concentration of the upper layer remained near saturation and fluctuated seasonally with water temperature (Figure 1.7). The intermediate and lower layer DO concentrations decreased to 2.48 mg L"1 and 0.03 mg L"1, respectively, by 29 November 2000. The lower layer was approaching anoxia. The 8 intermediate layer had 80 times the oxygen concentration of the lower layer and does not appear to be headed for anoxia in the near future. pH - The pH of the upper layer fluctuates seasonally (Figure 1.8). The variations in the pH of the intermediate and lower layers has not been so regular, they have exhibited mean [and standard deviation] values of 6.63 [0.15] and 6.31 [0.16], respectively, since 1997. Copper and Zinc - Copper and zinc concentration increase with depth (Figures 1.9 and 1.10). Zinc concentrations in the upper layer decrease each summer coinciding with fertilization of the pit lake. The decrease is attributed to biological removal by microorganisms (BFLP, 2000) by adsorption or active uptake. Haloclines - The upper halocline rose 10.2 m from 18.05 m (UWL) in July 1997 to 7.85 m (UWL) in March 2001. The halocline rate of rise was slowing with time and the halocline has become thinner (see Figure 1.11). The upper halocline was remarkably thin with an average thickness of 0.83 m, since January 1999. The interface has been rising due to the discharge of ARD to the intermediate layer, which increases that layer's volume and displaces the upper layer upwards. The halocline rose rapidly during the winter when ARD flows were high. During the summer the halocline rose less and even dropped slightly in the summer of 2000. The lower halocline had an average thickness of 8.1 m and its position oscillates seasonally, by as much as 9 m (Figure 1.11). It deepens in the fall/winter and rebounds in the spring/summer. The halocline deepened considerably in the winters of 1997/98 and 1998/99, but less so for the winters of 1999/2000 and 2000/2001, corresponding to the observed variation in ARD flow (Figure 1.11). The unusual shape of the lower halocline, particularly the observed double step, was probably due to the variable flow at the NIS and SIS diffusers, which was further complicated by the steep diffuser angle and the influence of salt water blocking on the diffuser hydraulics. The change in intermediate layer volume can be determined from the variation in upper and lower halocline depths and the lake bathymetry. Effective halocline depths were determined by replicating the haloclines (Figure 1.11) as step functions. The step was described by the salinity of the layers above and below the halocline, with the position of the step calculated by conservation of salinity (Appendix A). Figure 1.12 shows that the 9 intermediate layer volume changes (_>/) were proportional to the total injection system flow (QA), as was to be expected. However, the coefficient of proportionality is 1.59, so the expansion of the intermediate layer is caused not only by the input of injection system water, but also to processes linked to the injection system flow, i.e. entrainment. Note negative Qi occurs at times when the injection system flow was minimal and was due to the upper halocline deepening and/or the lower halocline rising, which must be caused by entrainment processes (see Section 1.4) in the upper and lower layers, respectively. 1.4 DISCUSSION DO levels - The intermediate layer DO decreased less rapidly than the lower layer DO, primarily because the intermediate layer receives oxygen via the ARD injection. In addition some DO is mixed downward across the upper halocline, and the lower layer experiences a greater sediment oxygen demand (Wilton, 1998). To hasten the depletion of oxygen in the intermediate and lower layers, the pit lake has been fertilized with liquid nitrate and phosphate each summer since 1997. There is little natural occurring nutrient input to the lake and the fertilization program was initiated to boost primary productivity in the surface euphotic zone. Fallout of organic debris to the intermediate and lower layers of the lake accelerated the depletion of oxygen in these layers (BHP 2000). A secondary benefit of fertilization is the reduction of zinc concentrations in the upper layer during the summers due to adsorption by algae (BHP 2000). In order to reduce the oxygen input to the intermediate layer, the NIS (containing -60 % of the annual total injection system flow) was diverted from the intermediate layer to the upper layer on 29 June 2000. The diversion of the NIS should not significantly affect upper layer water quality if biological removal of dissolved zinc continues to occur, and will actually improve the upper layer alkalinity (BHP, 2000). Metal concentrations - Treatment by metal-sulphide precipitation using sulphate reducing bacteria is desirable to control long term metal concentrations. Sulphate reducing bacteria are not yet active in the intermediate or lower layers, as sulphide has yet to be detected (BHP 2001) and metal concentrations have not dropped. Without metal-sulphide precipitation, dilution will only remain effective in the short-term, as future metal 10 concentrations will asymptotically approach the inflow metal concentrations (Wilton, 1998). The temporal change in metal concentration without biological treatment options can be explored using the assumptions of a continuously stirred reactor problem and a conservative substance, yielding: -£-=1. f c N I (1.1) where C is the concentration in the reactor, C,n is the inflow concentration, and Co is the initial reactor concentration. The time, t, is non-dimensionalized by the retention time, tr, which is 127 years for the "reactor" or intermediate layer volume divided by the SIS flow rate. Current intermediate layer concentrations of copper and zinc are 0.1 mg L"1 and 0.9 mg L"1 (Figures 1.9 and 1.10), respectively. The inflow concentrations of copper and zinc to the intermediate layer from the SIS are 0.79 mg L"1 and 5.93 mg L"1 (average of 1998 to 2000 values), respectively (BHP 1999, 2000, 2001). After 100 years, assuming no biological treatment and constant inflow concentrations, Equation (1.1) predicts that the intermediate layer copper and zinc concentrations will rise to 0.49 mg L"1 and 3.62 mg L"1, respectively. The upper layer has a retention time of only 1.6 years based on a 7 m depth and inflows from the NIS, intermediate layer (at a rate equal to the SIS inflow) and direct precipitation. Thus the concentration of metals in the upper layer will quickly converge to the volumetric average of inflows from the NIS and intermediate layer, diluted by direct precipitation. The NIS has copper and zinc concentrations of 0.02 mg L"1 and 1.88 mg L"1 (average of 1998-2000 values), respectively (BHP 1999, 2000, 2001). Given current intermediate layer copper and zinc concentrations, the expected upper layer copper and zinc concentrations are 0.03 mg L"1 and 0.70 mg L"1, respectively. Indeed upper layer zinc concentrations are observed to rise each year, but are controlled by fertilization and algal absorption each summer. After 100 years assuming no biological treatment, the upper layer concentrations would be 0.11 mg L"1 for copper and 1.43 mg L"1 for zinc. The permitted maximum concentrations of 0.05 mg L"1 for copper and 1 mg L"1 for zinc, are specified at offshore buoys in Rupert Inlet. Dilution within the groundwater discharge and at the Rupert Inlet foreshore should always ensure that these requirements are met. However, it is desirable that the surface water of the pit lake should also meet these 11 requirements. Therefore, long-term treatment requires metal-sulphide precipitation or failing that continued fertilization of the lake. Upper layer thickness - The rise of the upper halocline has been a concern because a thinner upper layer makes the lake more susceptible to destratification from wind induced upwelling. Maintaining the upper layer is important to the passive treatment system as it prevents the injection system plumes from surfacing and impedes oxygen transfer to the intermediate layer. The rate of rise of the halocline has decreased with time, which partially reflects that the 1999/2000 and 2000/2001 winters in the region were drier than previous years and the reduced inflow to the intermediate layer in 2000/2001 due to removal of the NIS discharge. However, it also indicates that wind and penetrative convection are increasingly eroding the halocline from above. An equilibrium upper layer depth is expected to be reached when upper layer mixing processes (wind and penetrative convection) can erode the halocline at a rate equal to its rate of rise from injection system inflow and entrainment of the upper halocline by the impinging plume (Figure A.l). The equilibrium upper layer depth is a conceptual idea representing a stable inter-annual depth. Short-term changes in halocline position are expected due to the relative strengths of competing processes. In support of an impending upper layer equilibrium depth is the upper layer salinity trend, which had been dropping due to direct precipitation, but has leveled off at 4 - 4.5 since 1999. This leveling off can only be attributed to the transport of saline intermediate layer water into the upper layer by downward erosion of the upper halocline. In fact the summer increases in upper layer salinity of 0.05, 0.1 and 0.25 (see Figure 1.5) for 1998, 1999 and 2000, respectively, are direct evidence of increased downward erosion with rise of the upper halocline. The results of a water balance indicate the increasing effectiveness of upper layer mixing processes at eroding the upper halocline as the upper layer thickness decreases (Appendix A). Diversion of the NIS to the upper layer provides the additional benefit of reducing the rate of rise of the upper halocline, thereby deepening the upper layer and preserving the intermediate layer salinity for longer, both of which reduce the long-term risk of upwelling. 12 Plume entrainment is an additional mixing process that may act to raise the halocline. The injection system discharges essentially fresh water creating a buoyant plume, which entrains the more saline intermediate layer water as it rises. If the injection system plume rises to the upper halocline, mixing will occur as the plume impinges against the upper halocline and is redirected horizontally. Plume entrainment is important for not only determining the upper halocline depth, but also as a source of oxygen to the intermediate layer affecting anoxia and the potential for metal-sulphide precipitation. The results of a water balance indicate that plume entrainment of the halocline is occurring at a rate of 0.26(24, which is supported by theoretical energy considerations (Appendix A). Indeed there is experimental and theoretical support for plume entrainment in the studies of Baines (1975), Kumagai (1984) and Cotel et al. (1997). Intermediate layer dynamics - The discharge of plumes into the intermediate layer results in slight temperature, salinity and density gradients with depth in that layer. These observations are consistent with the "filling box" mechanism, whereby a buoyant plume discharged into a confined region causes the receiving environment to become stratified (Baines and Turner 1969). This stratification reduces plume entrainment of the upper halocline by limiting the plume height of rise especially for low injection system flows. Indeed, during low injection system flows the plume was observed by dye tracer experiments to be trapped at depth by the intermediate layer density gradient, and for higher flows to rise to the upper halocline before spreading horizontally (Rescan 1999a, 1999b; Muggli et al. 2000). These dye tracer experiments were performed on the plumes in October 1998 at the NIS for a flow 1,800 m3 d"1 (25 % of NIS average), and in March 1999 at the SIS for a flow of 15,500 m3 d"1 (334 % of SIS average). For the low flow at the NIS the plume rose to a depth of 90 m (UWL) and slowly spread in a 10 m deep layer horizontally across the lake (Figure 1.13). Whereas, at the higher flow at the SIS the plume rose rapidly to the upper halocline, and spread quickly in a 60 m deep layer horizontally across the lake reaching the pit walls within 24 hours (Figure 1.14). The front of the horizontally spreading plume was estimated to be traveling at 0.036 m s"1 after two hours and 0.023 m s"1 after six hours. Lower halocline - At the lower halocline, the large seasonal oscillation is believed to result from the balance of mixing energies across the halocline. The large depth changes of the lower halocline, compared to the upper halocline, result from smaller density difference 13 rather than stronger mixing processes. Turner (1973) summarized the effect of entrainment as being the transport of fluid from the quiescent fluid into the stirred fluid. Thus for deepening of the lower halocline (or upward entrainment of fluid) there must be greater stirring in the intermediate layer than the lower layer. Indeed, the intermediate layer is strongly stirred by the injection system plumes, particularly because the injection system diffusers have some ports facing downwards (Muggli et al. 2000, Rescan 2001). Rapid deepening of the halocline is correlated with periods of high ARD flow (Figures 1.11). The rate of entrainment of the lower halocline is estimated at 0.33<2A by subtracting the injection system inflow of QA and estimated plume entrainment at the upper halocline of 0.26QA from the intermediate layer volume change of \.59QA-The observed summer rise of the lower halocline (and downward entrainment of fluid) is somewhat more difficult to explain, because it implies that the lower layer is stirred more vigorously than the intermediate layer during these periods. This behaviour is contrary to the previously made assumption that the lower layer is "quiescent" (Wilton, 1998). The lower layer may actually be well mixed as the uniformity of potential temperature, salinity and density imply, by high Rayleigh number (Ra) thermal convection, see below. Double diffusion may also be an important mechanism at the lower halocline (Fisher and Lawrence 2001). At the lower halocline, warm saline fluid underlies fresher colder fluid, creating the situation where salinity and temperature are making opposing contributions to the vertical density gradient, the right condition for double diffusive convection. Double diffusive convection is known to transport heat and salt (Turner 1965, Huppert 1971), and cause the upward migration of interfaces (Turner 1968, Huppert and Linden 1979). Lower layer - The behaviour of the lower layer may be analogous to the classic thermal convection problem of fluid between two plates, uniformly heated below and cooled above. The lower layer is estimated (from the observed temperature change 1998-2000) to be heated geothermally at approximately 0.023 W m" , which is consistent with geothermal heating rates of 0.02-0.04 W m"2 for British Columbian fiords (Hyndman, 1976). Note the initial high rate of heating is attributed to the gradual release of heat stored in the pit walls from solar radiation during the summer of 1996. Cooling from above is provided by the 14 decrease in temperature ~ 0.8 °C across the lower halocline. The stability of the fluid is described by Ra, Dv (1.2) where the coefficient of thermal expansion, a = -2 x 10"4T ]; the depth, d = 38 m, is the average lower layer height (volume/surface area); thermal molecular diffusivity, D = 1.4 x 10"7 m2 s"1, and the kinematic viscosity v = 1.3 x 10"6 m2 s"1. If AT = 0.8 °C, then Ra = 4.7 x 1014. This far exceeds the critical Rayleigh number for transition to convective motion (Ra > 1700) and for transition to turbulence (Ra > 14,000 Pr° 6 = 56,000 for Prandtl number = v/D = 10) (Turner, 1973). Thus the lower layer is turbulent and has the ability to erode the lower halocline, particularly during the summer when the ARD flow is negligible, resulting in a rise in the halocline. 1.5 CONCLUSIONS The Island Copper pit lake has evolved into three distinct layers: a brackish and well mixed upper layer; a plume stirred intermediate layer; and a thermally convecting lower layer. To date the upper halocline has risen over 10 m, due to the injection of buoyant acid rock drainage (ARD) into the base of the intermediate layer. An upper layer "equilibrium" depth is expected in the near future reflecting a balance between meteorological forcing from above and rising ARD plumes from below. The lower halocline depth fluctuates seasonally: deepening in the winter due to high ARD flows, which erode the lower halocline; and rising in the summer when thermal convection in the lower layer dominates over reduced intermediate layer stirring. The proposed treatment of acid rock drainage by metal-sulphide precipitation using sulphate-reducing bacteria is dependent on anoxia developing in the intermediate and lower layers. The lower layer is close to anoxia (DO = 0.03 mg L"1); whereas, with a DO = 2.5 mg L"1, the intermediate layer has a long way to go. Without metal-sulphide precipitation metal concentrations in the pit lake are likely to increase, although they may be controlled by fertilization via biological removal. Table 1.1. Island Copper pit lake volume and area (adapted from Rescan 1999a). Depth interval *' Layer Cumulative Area Volume Volume (m) (m3) (m3) (m2) 0 12.192 2.08E+07 2.08E+07 1.71E+06 12.192 24.384 1.95E+07 4.03E+07 1.60E+06 24.384 36.576 1.85E+07 5.88E+07 1.52E+06 36.576 48.768 1.76E+07 7.64E+07 1.44E+06 48.768 60.96 1.67E+07 9.31E+07 1.37E+06 60.96 73.152 1.59E+07 1.09E+08 1.30E+06 73.152 85.344 1.50E+07 1.24E+08 1.23E+06 85.344 97.536 1.42E+07 1.38E+08 1.16E+06 97.536 109.728 1.32E+07 1.51E+08 1.08E+06 109.728 121.92 1.22E+07 1.64E+08 1.00E+06 121.92 134.112 1.12E+07 1.75E+08 9.19E+05 134.112 146.304 1.01E+07 1.85E+08 8.28E+05 146.304 158.496 8.94E+06 1.94E+08 7.33E+05 158.496 170.688 8.02E+06 2.02E+08 6.58E+05 170.688 182.88 7.24E+06 2.09E+08 5.94E+05 182.88 195.072 6.52E+06 2.16E+08 5.35E+05 195.072 207.264 5.74E+06 2.21E+08 4.71E+05 207.264 219.456 4.78E+06 2.26E+08 3.92E+05 219.456 231.648 3.80E+06 2.30E+08 3.12E+05 231.648 243.84 3.03E+06 2.33E+08 2.49E+05 243.84 256.032 2.44E+06 2.35E+08 2.00E+05 256.032 268.224 1.91E+06 2.37E+08 1.57E+05 268.224 280.416 1.39E+06 2.39E+08 1.14E+05 280.416 292.608 9.47E+05 2.40E+08 7.77E+04 292.608 304.8 6.93E+05 2.40E+08 5.68E+04 304.8 316.992 4.56E+05 2.41E+08 3.74E+04 316.992 329.184 1.81E+05 2.41E+08 1.48E+04 329.184 350 3.59E+04 2.41E+08 2.94E+03 * l Bathymetry was measured when pit lake was 2.44 m below ultimate water level. 16 Figure 1.1. Oblique aerial photograph of the Island Copper Mine site and pit lake looking east. Points of interest are: (a) pit lake; (b) on land waste rock dumps; (c) the beach waste rock dump; (d) Rupert Inlet; and (e) remaining mine buildings, sludge thickening tanks and docking facilities. 17 System 1 1 1 1 : r 607000 607500 608000 608500 609000 upper layer (mixed) ( D 0 250 500 750 1000 1250 1500 1750 2000 Distance (m) Figure 1.2. Island Copper pit lake (a) bathymetry and (b) cross-section AA. Horizontal datum is WGS-84 and depths are referenced to the ultimate water level (UWL), which is 2.44 m above mean sea level (adapted from Muggli et al. 2000). 18 (a) DO I DO (mg/L) 4 T 10 15 20 S, T(°C), atG<g/mJ) 25 30 DO (mg/L) 4 10 15-215 220 (b) DO: I 225 230 235-(c) DO + \—<—I H — ' — h 0 5 10 15 20 25 30 S, T(°Q, at0<g/m3) Figure 1.3. Vertical profiles of the (a) full depth, (b) upper interface and (c) lower interface showing salinity (S), temperature (T), density as a, = p(S,T,P-0)-1000 and dissolved oxygen (DO) for 6 December 1999. 19 Figure 1.4. Weekly average ARD flow to the intermediate layer via south injection system [—] and north injection system ["']. The NIS was diverted to the upper layer 29 June 2000. 25.0 20.0 15.0 '2 03 00 10.0 -i 1 1 1 1 1 i 1 r Jul-96 Jan-97 Jul-97 Jan-98 Jul-98 Jan-99 Jul-99 Jan-00 Jul-00 Jan-01 Jul-01 Date 25.0 -\ 1 1 1 1 r Jul-96 Jan-97 Jul-97 Jan-98 Jul-98 Jan-99 Jul-99 Jan-00 Jul-00 Jan-01 Jul-01 Date Figure 1.5. True salinity (bold symbols) and seawater salinity trends for (a) upper layer [ and (b) intermediate layer [A,A] and lower layer [•,•]. 21 Figure 1.6. Temperature trends for (a) upper layer [•] and (b) intermediate layer [_] and lower layer [•]. 22 Figure 1.7. layer[•]. Dissolved oxygen trends for upper layer [•], intermediate layer [A] and lower 23 9.0 Jul-96 Jan-97 Jul-97 Jan-98 Jul-98 Jan-99 Jul-99 Jan-00 Jul-00 Jan-01 Date Figure 1.8. pH trends for upper layer [•], intermediate layer [A] and lower layer [•]. 24 0.7 Jul-96 Jan-97 Jul-97 Jan-98 Jul-98 Jan-99 Jul-99 Jan-00 Jul-00 Jan-01 Date Figure 1.9. Copper concentration trends for upper layer [•], intermediate layer [A] and lower layer[•]. 25 1.0 Jul-96 Jan-97 Jul-97 Jan-98 Jul-98 Jan-99 Jul-99 Jan-00 Jul-00 Jan-01 Date Figure 1.10. Zinc concentration trends for upper layer [•], intermediate layer [A] and lower layer[•]. 26 00 2 C/3 GO Jul97 Jcn98 Jul98 Jcn99 Jul99 JcnOO JulOO JcnOl JulOl Figure 1.11: Time series of salinity profiles through the (a) upper halocline and (b) lower halocline with surface water level [—] and effective halocline depth ['"]. Thicker lines represent the extent of the halocline. The bars adjacent to (b) indicate the depths of the north injection system (NIS) and south injection system (SIS), (c) Total injection system flow (QA) averaged over the period between CTD casts. 27 Figure 1.12. Change in intermediate layer volume (Qi), determined from upper and lower effective halocline depths, versus total injection system flow (QA)- The regression line Qi=1.59 Q A - 7950 has a correlation coefficient (r) of 0.89. 2 8 Figure 1.13. Lowest measured dilutions in north injection system plume on 9 October 1998 for mean flow rate of 1,800 m3 d"1. North-south cross-section during 6 hours of continuous rhodamine dye injection (from Rescan 1999a, Muggli et. al. 2000). 29 30 CHAPTER 2 Determination of Salinity and Density Abstract The mine pit at Island Copper Mine, Vancouver Island BC, was flooded with seawater in 1996, to create a meromictic lake for the treatment of acid rock drainage (ARD). The original seawater has changed in chemical composition, such that the Practical Salinity Scale 1978 (PSS78) is no longer accurate. Non-seawater ions are increasing with time and are sourced primarily from ARD and dissociation of ions from the pit walls. A time variant salinity correction to PSS78 was developed and is different for each of the three layers within the lake. The salinity correction is based on the density-derived salinity, which was back-calculated from laboratory density measurements using the seawater Equation of State 1980 (EOS80), whose accuracy was verified for pit lake waters. Other estimates of salinity such as PSS78, chlorinity, drying and chemical composition proved to be insufficient or unsatisfactory. If the salinity correction is included, the lower layer is found to be more saline than the original seawater. Future changes to the chemical constituents of the pit lake will necessitate ongoing adjustments to the salinity correction. The work illustrates the difficulty of determining the salinity and density in a biochemically active system such as a pit lake. It is recommended that verification of salinity and density equations be a regular first step in the investigation of biochemically active lakes. 2.1 Introduction Island Copper pit lake - The BFLP Billiton owned Island Copper Mine pit near Port Hardy, Vancouver Island, BC, Canada, was flooded in 1996 with seawater and capped with fresh water to form a meromictic pit lake. The pit lake is being utilized as a passive treatment system for acid rock drainage (ARD) (BHP 1996). ARD is caused by the chemical and biological oxidation of reactive sulphide minerals when exposed to air and water. Tests 31 suggested that the on-land waste rock dump material would produce moderately acidic leachate and release elevated concentrations of sulphate and metals over the long-term (Dagenais 1996). The passive treatment system collects ARD from the on-land dumps and discharges it via two injection systems at a depth of -220 m into the pit lake. These discharges rise as buoyant plumes through the more dense salt water (Figure 2.1). The ARD is presently treated in the intermediate layer by dilution and by neutralization utilizing the natural alkalinity of the seawater. The upper 'brackish' water layer forms a cap over the intermediate layer, to prevent the plumes from surfacing and to restrict oxygen transfer into the intermediate layer. It was also anticipated that treatment would be provided by metal-sulphide precipitation via anaerobic sulphate-reducing bacteria (BHP 1996), once anoxic conditions were established in the lake (Chapter 1). Background information on the pit lake and passive treatment system was provided in Chapter 1. Of significance to the current chapter is the observed development of three distinct layers within the pit lake: a brackish and well mixed upper layer; a plume stirred intermediate layer; and a thermally convecting lower layer (Figure 2.1). In November 2000 the dissolved oxygen concentrations of the intermediate and lower layers were 2.58 mg L"1 and 0.03 mg L"1, respectively. The pH of the intermediate and lower layers dropped soon after flooding, but from 1997 to 2000, the pH of the intermediate and lower layers showed similar small variations with mean [and standard deviation] values of 6.63 [0.15] and 6.31 [0.16], respectively. The pH of the upper layer (@ 3 m) had a mean of 7.41 and varied seasonally between - 6.9 and ~ 8.7. The alkalinity of the upper layer has decreased to near zero. The intermediate layer alkalinity was steady with a mean in year 2000 of 1.90 mg CaC03 L"1, while the lower layer alkalinity has increased to a mean in year 2000 of 3.27 mg CaC03 L"1. Of the heavy metals of environmental concern: cadmium, copper, manganese and molybdenum, all increased with depth, whereas, zinc concentrations were highest in the intermediate layer which was attributed to the ARD inflow. Salinity and Density - Salinity is the measure of dissolved salts in the water, and can be used with temperature and pressure to determine the density of the water. In the pit lake salinity is the dominant contributor to density differences, so accurate determination of 32 salinity is essential to understanding the physical mixing processes and the long-term stability of the stratification. Wilton (1998) used PSS78 and Equation of State 1980 (EOS80) to describe the pit lake salinity and density, which was a reasonable first assumption as the pit lake water originated from seawater. However, the addition of ARD, dissociation of ions from the pit walls and perhaps other biochemical changes have changed the lakes' chemical composition. If these chemical changes are significant enough to invalidate PSS78 and EOS80, then new and unique methods for salinity and density will be required for the pit lake. The objective of this chapter is to determine whether PSS78 and EOS80 adequately describe the salinity of the pit lake, and to determine methods to correct for the differences. The chapter reviews literature on salinity and density in seawater and lakes. The emphasis is on lakes with unusual water chemistry, and the methods used to develop unique salinity and density equations. Our methods of analysis are detailed. In the results, the alternate methods for determining salinity are compared and the density method is selected as the preferred method. EOS80 is tested and found to accurately describe density for pit lake waters. Equations for salinity correction are then developed. The discussion focuses on the sources of the salinity correction and the limitations of the density method. We conclude by summarizing our findings for the Island Copper pit lake and explaining the source of non-seawater ions. 2.2 Application of Theory Salinity - True salinity is defined as the mass of dissolved salts in one kilogram of solution (%0). Direct determination of salinity is very difficult due to the complexity of seawater. Drying and weighing the solution is inaccurate due to the loss of components, unless the International Commission (Forch et al. 1902) method is used, but it is difficult and slow (Sverdrup et al. 1942). Finding salinity via analysis of the chemical composition is difficult too, because of the analysis precision required to provide the necessary accuracy (American Public Health Association 1995). Thus, salinity is determined from indirect (or practical) methods. In oceanography, salinity was traditionally determined from chlorinity, a measure the halide content (CI, mostly chloride) by silver nitrate titration, using an established relationship between CI and salinity (Sverdrup et al. 1942). Other indirect 33 methods are based on physical properties such as conductivity, density, sound speed or refractive index (Millero and Sohn 1992, American Public Health Association 1995). At present the determination of salinity from conductivity is usually used, particularly for seawater, as its relative chemical composition is nearly constant (Millero and Sohn, 1992). PSS78 was introduced to calculate the salinity of seawater from measurements of conductivity, temperature and pressure (Lewis and Perkin 1978, Perkin and Lewis 1980). PSS78 low salinities are based on the dilution of standard seawater (with nearly constant chemical composition) with pure water. The Practical Salinity Scale is a unit-less quantity and does not require a unit symbol. PSS78 assumes that all waters with the same conductivity ratio have the same salinity. The conductivity ratio is the measured conductivity divided by the conductivity of standard seawater (S=35, T=15 °C). If the chemical composition is sufficiently different from standard seawater, the salinity and subsequent calculations for density will be in error. Non-electrolytes such as SiC>2 are of particular concern because they do not affect the conductivity. Deviations from PS.S78 are known to occur in the natural environment, such as estuaries, anoxic basins, hydrothermal vents and evaporated basins, where processes such as precipitation, dissolution, evaporation, freezing and oxidation can change the chemical composition of the seawater (Millero and Sohn 1992). Density derived salinity can be back-calculated from an equation of state, such as EOS80 (American Public Health Association 1995). Millero and Sohn (1992) examined estuarine waters by mixing river water with seawater and found that density derived salinities were in better agreement with their estimates of the true salinities than the PSS78 values. An advantage of the density method is that all chemical components are included (including organic compounds and non-electrolytes). The method is dependent on the equation of state accurately describing the true salinity-density relationship; i.e. EOS80 must be unaffected by the change in chemical composition. While the salinity relations are well understood for seawater, they are more complicated for lakes because the relative proportions of major ions vary from lake to lake and sometimes even vertically within a lake (WUest et al. 1996). Relationships between in situ conductivity and salinity must be established for each unique water volume in order to 34 utilize field measurements from conductivity-temperature-depth (CTD) probes. A variety of methods have been used. McManus et al. (1992) experimentally determined temperature-conductivity-salinity for Crater Lake, where salinity was determined from chemical constituents. Wuest et al. (1996) developed a conductivity-salinity relationship for Lake Malawi via chemical composition for salinity and summing of individual species conductivities. Determining conductivity from chemical composition is limited to natural fresh water with total dissolved solids less than about 2.5 g L"1 (American Public Health Association 1995). Hall and Northcote (1986) found for meromictic Mahoney Lake that chemical composition changed with depth and time, and that samples were unstable, thus it became necessary to make T-C plots in the field. Jellison and Melack (1993) and Romero and Melack (1996) developed a direct temperature-salinity-density relationship for meromictic Mono Lake, from salinity measurements by drying and weighing and direct density measurements. Most recently Millero (2000a) applied a salinity correction to PSS78 salinity, where the salinity correction was based on estimates of true salinity back-calculated using EOS80. With regard to mine pit lakes, Schimmele (1998) found that high concentrations of dissolved substances influenced the measurement of conductivity and density measurements were higher than those calculated from CTD measurements using various lake equations of state. Stevens and Lawrence (1998) assumed salinity was proportional to temperature-corrected conductivity following Wuest et al. (1992), for a pit lake with salinity ranging from 0.6 to 1. To overcome the difficulties of working through salinity in pit lakes with temporally variable chemical stratification, Boehrer (2001) tested an ultrasonic sensor to directly measure density with promising results. Density - Density for oceanographic waters is defined by the EOS80 as a function of temperature (IPTS68), salinity (PSS78) and pressure (UNESCO 1981). Hill et al. (1986) extended EOS80 for salinities less than 2. The function results from an empirical fit to data for standard seawater and is non-linear. While seawater is remarkably consistent in composition, there are slight variations in chemical composition between oceans, which are not important if relative differences at a given locality are of interest. Millero (2000b) found changes to the composition of deep water due to the breakdown of plant material, adding 35 alkalinity, carbon dioxide, silica and nitrate can affect the calculated densities at a given PSS78 salinity. In large scale and inter ocean studies these differences can be important. Millero (2000b) suggests two alternative approaches to correct the density for varying chemical composition by: (1) using EOS80 with a true salinity corrected from PSS78 salinity by accounting for the mass of added dissolved solids; or (2) EOS80 with a correction for the change in density from added solutes. Determining the density of lake water from salinity has been routine, provided a good estimate of true salinity can be established, because the equation of state formulations are relatively insensitive to changes in chemical composition. McManus et al. (1992) and Wiiest et al. (1996) used the limnological equation of state (Chen and Millero 1986) to calculate density after estimating the true salinity. McManus et al. (1992) found computed densities to be in good agreement with direct measurements of density. Wiiest et al. (1996) found the coefficient of haline expansion J3 (ionic species) determined from chemical composition and theoretical considerations for Lake Malawi to be only 1.028 times greater than the value for /? derived from seawater EOS80. Millero (2000a) compared density from molal volume data to density from seawater EOS80 using chemical composition salinity, for Lake Baikal, Lake Malawi and Lake Tanganyika and confirmed that the densities were in good agreement (within ~ 0.01 kg m" ). Furthermore, for Lake Tanganyika these calculated densities were in excellent agreement with density measurements made using a vibrating densimeter. Spigel and Priscu (1996) determined the density of a chemically stratified Antarctic lake using PSS78 salinity and EOS80, but found it necessary to modify EOS80 for hypersaline waters (S > 42). Island Copper Approach - Based on the literature review, the density method was chosen as the most promising technique for estimating the true salinity. Other methods for determining salinity such as PSS78, chlorinity, drying, and chemical composition were not expected to be accurate, due to the high salinity and non-standard chemical composition of the pit lake water. Nevertheless, all these methods were employed to confirm this hypothesis. Following Millero (2000a), the true salinity (Srrae) as estimated by the density-derived salinity can be related to the conductivity salinity (Spss) by, 36 where the salinity correction, AS, is an empirical function of the difference between Sprue and Spss- The density method is only viable if the pit lake water behaves according to EOS80. The approach was applied to Island Copper pit lake samples to develop equations for AS using laboratory analysis. Then, AS was applied to SPSS field data to establish the true salinity of the pit lake. Alternately, a salinity correction could be developed for salinity determined by another reliable field method, such as chlorinity. 2.3 Methods To investigate the salinity and density relationships, samples were taken on 29 November 2000 from the upper layer (5 m), intermediate layer (100 m), lower layer (300 m and 350 m), the south injection system and the surface of Rupert Inlet. Alternative measures of salinity were investigated for the 29 November 2000 samples and are detailed below. Of these conductivity and chlorinity salinity have been measured routinely since the pit lake was flooded. (i) Spss the conductivity salinity, was determined from PSS78 using measurements made with a Guildline Model 8400 Autosal (accuracy 0.005 plus any sampling errors). (ii) Sci the chlorinity salinity was determined using silver nitrate titration (accuracy 0.02 plus any sampling errors) and the 1966 Joint Panel for Oceanography Table and Standards definition of SCL = 1.80655 (CI %c). (iii) SFR the filtrable residue was the fraction of the sample remaining after passing through a V W R Brand 696 glass filter followed by evaporation and drying at 103 ° C for one hour. (iv) SFFR the fixed filtrable residue was the fraction of the sample remaining after passing through a V W R Brand 696 glass filter followed by evaporation, drying and igniting at 550 ° C for one hour. The SFR and SFFR measurements were done at a commercial laboratory, and the stated salinity accuracy was 0.01 for injection system water, 0.05 for the upper layer and 0.7 for the pit lake and Rupert Inlet samples. 37 (v) Sdiem the salinity from chemical composition was based on routine commercial laboratory analysis of eight anions (by ion chromatography analysis), 28 cations (by inductively coupled plasma atomic emission) and silica (Si02). The accuracy of the chemical composition measurements was ~5 %. The Cl measure was improved using the chlorinity measurement. Bicarbonate was estimated from the alkalinity. (vi) SEOS the density-derived salinity was back-calculated from density measurements using EOS80. Density and temperature measurements were made with an Anton Paar DMA 5000 Oscillating U-Tube Density Meter, which had density and temperature accuracies of 0.01 kg m" and 0.01 °C, respectively. Measurements were made at atmospheric pressure. Based on the literature review findings and results from the above comparison, SEOS was pursued as the most promising method for finding the true salinity of the pit lake. To develop an equation for AS the conductivity-temperature-salinity relationship was tested. Samples were stored in a refrigerator and filtered through an 11 pm Whatcom Grade 1 Qualitative filters prior to analysis. Sub-samples for each layer sample (29th November 2000) were created to achieve a range of concentrations by diluting with double distilled water and evaporating. Each sub-sample was measured for conductivity and temperature while the temperature was varied over the expected temperature range; 4-24 °C upper layer and 10-15 °C for the intermediate and lower layers. The initially cold sub-sample, in a 250 ml sealed jar, was stirred constantly and allowed to warm slowly by heating from the ambient room temperate and hand warmth. Measurements were made with a Guildline 9540 Digital Platinum Resistance Thermometer (accuracy 0.01 °C) and Radiometer Copenhagen CDM210 Conductivity Meter with CDC566T probe (accuracy was ± 0.2 % of reading ±3 on the fourth significant figure). Density measurements were made of each sub-sample and dilution to determine SEOS-2.4 Results Salinity - PSS78 and EOS80 do not accurately predict the density of the Island Copper pit lake water (Table 2.1). Density calculated with EOS80 using Spss was less than density measured directly with the Anton Paar by 0.23, 0.28 and 1.07 for the upper, intermediate and lower layers, respectively (for 29th November 2000). Therefore, the 38 addition of ARD and other chemical changes to the pit lake have significantly altered the chemical composition of the lake from the original seawater. The true salinity needs to be established and to this end the six alternative salinity methods are examined. Salinity was determined by PSS78, chlorinity, filtrable residue, fixed filtrable residue, chemical constituents and from density (methods are detailed in the methods section). Samples were taken on 29th November 2000 from the injection system and pit lake at 5 m, 100 m, 300 m and 350 m. A sample from Rupert Inlet was included as a seawater control. The wide range of salinity values obtained using different salinity methods illustrates the difficulty of determining salinity (Table 2.2). The differences between Spss and Sci confirm that pit lake water was no longer seawater, because for seawater these methods would be in agreement (Lewis and Perkin 1981). Neither PSS78 or chlorinity salinity methods can be trusted as both methods are based on standard seawater. Note, Spss values for the pit lake were greater than Sci, which indicates that non-seawater ions cause an increase in conductivity that was not compensated for in chlorinity. Determination of salinity by drying and weighing at low temperatures, SFR, and high temperatures, SFFR, produced quite different results that tend to bracket SPSS. The inaccuracies of these methods was expected because low temperature drying is insufficient to evaporate all water and at temperatures necessary to drive off the last tracers of water other components are vaporized (Millero and Sohn, 1992). Salinity from chemical constituents was greater than SPSS and Sci for pit lake waters (except for the intermediate layer were an error in sodium concentration is suspected), whereas it was slightly lower for Rupert Inlet seawater, which again indicates the presence of additional non-seawater sourced ions in the pit lake. Inaccuracies preclude the use of Schem as a method for determining salinity, as the error on constituent measurements (except CI) was 5 %. Also, the equivalence ratio (the sum of cation charges over the sum of anion charges) was an average of 0.86, which indicates the analytical inaccuracies and/or that major ions were not included. The remaining and most promising salinity method was to estimate the density-derived salinity, SEOS, by back-calculation from direct measurements of density using EOS80. SEOS and SPSs should be the same for Rupert Inlet seawater, but were found to differ by 0.23. 39 The salinity of coastal and estuarine waters using PSS78 should be good to ± 0.04, Millero and Sohn (1992). The larger difference may be due to the nature of the sample, which was collected from the surface of Rupert Inlet at the Island Copper wharf, where the seawater has been diluted and altered chemically by local runoff and subject to biological growth. PSS78 measures only the ionic salinity whereas density-derived salinity measures total salinity (ionic, non-ionic, organic). The method was tested on northeast Pacific Ocean water (at Station P) and SEOS was within 0.01 of SPSS, which confirms the suitability of the density method. The differences between SEOS and Spss were 0.30, 0.37 and 1.40 for the upper, intermediate and lower (300 m) layers, respectively. True salinity, if it was indicated by SEOS-, was greater than SPss and the difference was different for each layer and increased with layer depth. The implication was that PSS78 was not sufficient for calculating salinity in the pit lake. Density- Millero (2000a) and Wiiest et al. (1996) found that changing chemical composition did not adversely affect the equation of state. To confirm this result for our own samples we experimentally determined the coefficients of thermal expansion (a) and haline expansion (J3) for pit lake waters, and compared them to values calculated from EOS80. A simplified density equation relates changes in temperature (7) and salinity (5) and Grand fito changes in density (p), where pm is a mean density (Turner, 1973). The coefficient of haline contraction was investigated using sub-samples prepared by diluting and evaporating 29th November 2000 samples from the upper, lower and intermediate layer. Additional sub-samples were made from the intermediate layer sample by diluting with injection system water. Salinity was determined in two ways: SEOS by measuring density (@ 15 °C) and back calculation using EOS80 and Sdu from the known volumetric dilution assuming SEOS of the undiluted sample to be the correct starting value. SEOS and Sdu were in good agreement and the slope dp/dS was the same for all samples (Figure 2.2). PEOS and pM were then calculated for each sub-sample using Equation (2.2) and p = pji+c<r+/3s) where (2.2) 40 SEOS and Sdu- PEOS varied by a maximum of 0.15 % across all sample types indicating the relative insensitivity of B to changes in salinity and chemical composition (Table 2.3). Good agreement was found between j3Eos and /_•/, although there was a tendency for fiEos to be slightly greater than /3dU. In fact the maximum difference between (3Eos and fya of 2.66 % occurred for the Rupert Inlet control sample. The greatest source of error was in the volumetric dilution and evaporation of samples used to determine /?„/. Next, changes to the coefficient of thermal expansion were investigated and the results are summarized in Table 2.4. Density was measured for changes in temperature from 14-15 °C and an experimental c^ .vp was calculated according to Equation (2.2). This was compared to aEos calculated directly from EOS80. The coefficient of thermal expansion decreased with increased salinity for both c x^p and CCEOS (Table 2.3). cCeXp was found to be in excellent agreement with aEos- The maximum difference between CCEOS and cXeXp was 0.56 % for the intermediate layer sample. These experiments on a and (3 show EOS80 to be independent of the changes in chemical composition as experienced in the pit lake. Thus, EOS80 can be reliably used to back-calculated salinity from measured density, and was the best of the methods investigated for determining salinity. Salinity Correction - The remaining task was to develop an empirical equation for the salinity correction, AS, based on conductivity-temperature-salinity experiments. This equation provides the necessary correction to Spss to enable density to be computed from CTD field data. AS was calculated from Equation (2.1) assuming SEOS to be the true salinity. = SEOS ~ SPSS (2-3) Correcting Spss rather than direct calculation of SEOS = f(C,T), was chosen as Spss = f(C,T,P) is a multi order polynomial that would be difficult to replicate. PSS78 also accounts for in situ pressure effects on salinity that are of order 0.1 over the 350 m depth of the pit lake. A unique equation for AS in each layer was determined using a linear two independent variable equation of the form AS = a + b.T + C.SPSS- The coefficients a, b and c were chosen by method of least squares. 41 Upper layer: AS = 4.5482 - 1.00 SPSS for 4 < T< 20, 3 < S < 6 (2.4) Int. layer: AS = -0.182 - 0.00530 T+ 0.0204 S W 5 for 11 < T< 13, 20 < S < 27 (2.5) Lower layer: AS = 0.0420+ 0.00206 T + 0.0447 S m for 12 < T < 14, 27 < S < 30 (2.6) The standard deviation of AS was 0.04, 0.05 and 0.06 for the upper, intermediate and lower layers, respectively. AS increases with salinity for all layers and to a lesser extent with temperature in the intermediate and lower layers. Higher order polynomials and second order 'TS" terms achieved only slight reduction in the standard deviation of AS so the simpler linear equation was chosen. Applying Equations (2.4), (2.5) and (2.6) for 29 November 2000 pit lake conditions the salinity correction was 0.29, 0.29, 1.34 for upper, intermediate and lower layers, respectively. 2.5 Discussion Salinity Correction - A limitation of the salinity correction equations is that they are based on samples on a given date diluted with double distilled water or evaporated. In reality the chemical constituents of the pit lake are changing due to the injection of ARD and changes to water chemistry. As an example of the effect of ARD, intermediate layer water was diluted with an equal amount of injection system water and the resulting density was 1.28 kg m3 greater than if it had been diluted with double distilled water (Figure 2.3), which is the basis of our AS equations. Fifty percent dilution of the intermediate layer by ARD is expected to take -60 years. An alternative approach could have been to base AS equations on dilution with injection system water, but it would still not account for mixing between layers, and chemical changes to the layers and injection system water. To apply the salinity correction to data prior and after November 2000 we need to know how AS has been changing with time. The original pit lake water was seawater, so AS - 0 at that time. The original salinity was estimated as 28.87 ± 0.10 by extrapolating Spss and Sci data for the intermediate and lower layers back to the time of flooding. Salinity by PSS78 and chlorinity should be in close agreement at this time as the water was seawater (Lewis and Perkin, 1981). 42 Examining the temporal difference between SPSS and Sci provides an insight to how the pit lake salinity has evolved and AS has changed. Sci can be thought of as a tracer for the original seawater, since Sci is based on the halide content (mostly CI) and the injection system water contains very little CI. Spss on the other hand responded to the conductivity of all ions in the pit lake. Thus the difference Spss less Sci is postulated to represent the increase in non-seawater sourced ions, which contribute to conductivity and Spss but not to Sci- The difference has increased with time in a roughly linear fashion (Figure 2.4). Therefore the salinity correction, which results from non-seawater sourced ions, must also be increasing with time. So the salinity correction determined for November 2000 is applicable only at that time. Source of non-seawater ions - The chemical constituents of pit lake water were compared to Rupert Inlet seawater to identify non-seawater sourced ions (Table 2.5). The chemical concentrations of Rupert Inlet seawater were multiplied by a ratio of pit lake chloride to Rupert Inlet chloride to allow for the dilution of pit lake seawater. These discounted concentrations were subtracted from the chemical concentrations of pit lake samples to determine the chemical differences between the pit lake water and seawater. The lower layer, for example, had excess: sulphate, calcium, [100 - 1000 mg L"1]; sodium [100-10 mg L"1]; magnesium, silica, manganese, and strontium [10-1 mg L"1]; and zinc and copper [1-0.1 mg L"1]. The excess chemical constituents of the lower layer were also found in high concentrations in the south injection system ARD. The match was weaker for the upper and intermediate layer, but still apparent. Calcium sulphate was the dominant compound. The comparison was limited by the accuracy of the chemical analysis, the variability of lake and injection system chemistry, and by the influence of other chemical processes. However, there is sufficient similarity between the ARD chemical composition and the pit lake excess concentration, to say that the ARD or the rock type producing ARD seems to be an important source of non-seawater ions. The salinity correction for the lower layer was considerably greater than for the intermediate layer even though it was not directly receiving ARD. Ions from the intermediate layer can not be transported or diffused to the lower layer without decreasing the lower layer's salinity. Thus most of the non-seawater ions in the lower layer must have 43 their origin in the lower layer. Dissociation of ions from the pit walls may have caused the increase of non-seawater sourced ions observed in the lower layer, and would be greater in the lower layer for several reasons. The pH of the lower layer at -6.3 was 0.3 lower than the intermediate layer (Chapter 1). The lower layer has a lesser volume to pit wall area ratio than the intermediate layer, 35 m compared to 140 m. Compounding these factors was the abundance of fine rock in the lower layer from: waste rock left in the bottom of the pit during the final year of mining operations; fine tailing material scoured out of the water quality ponds during flooding; rock scoured out of the pit walls during flooding; and waste rock from the northwest dump that was pushed into the pit lake to control its ARD (Brian Welchman pers. com.). The fine rock decreases the "effective" volume to surface area ratio of the lower layer. The sub-oxic conditions and decomposition of biomass may also elevate the concentrations of non-seawater sourced ions in the lower layer. Rescan (1999a) undertook biochemical studies in October 1998, at a time when the lower layer was approaching anoxic conditions (dissolved oxygen under 0.9 mg L"1), and found the lower layer chemical composition to be quite different to the overlying intermediate layer. Under these reducing conditions lower layer concentrations of alkalinity, ammonia and nitrite were greater than the intermediate layer. Silicate concentrations increased in the lower layer due to remineralization of organic matter and/or release via weathering of soluble minerals from pit walls. Also, observed was an increase in manganese concentrations due to redox cycling, and concurrent increases in cadmium, cobalt and copper. Increases in metal concentrations were observed above bottom sediments where there was a concurrent decrease in dissolved oxygen. Time Dependent Salinity Correction - The salinity correction is a constant correction to Spss for each layer. Clearly these corrections are not appropriate for dates closer to flooding when the water was more similar to seawater. We propose as a first approximation a salinity correction, ASt, for each layer, which increases linearly with time from the end of flooding until November 2000. 44 AS, = t AS (2.7) 1590 where t is the time in days since flooding ceased (23 July 1996). Thus our best estimate of salinity from in situ CTD measurements is the sum of SPSS and AST. Note, there is a unique salinity correction for each of the three layers in the pit lake (Equations 2.4, 2.5 and 2.6). Within the interfaces AS, is proportionally applied based on the AS, for the adjacent layers. EOS80 can then used to calculate density. Care should be applied in the future when applying this correction, as the chemical input from ARD and the pit walls should decrease and major chemical changes are anticipated if the intermediate and lower layers become anoxic. Salinity calculated with the time variant salinity correction, Spss + ASt, gives an estimate of lower layer salinity that is slowly increasing with time (Figure 2.5). The lower layer salinity as of November 2000 was 1.34 greater than the Spss estimate and 0.83 greater than the original seawater of the pit lake. Upper and intermediate layer salinities are still decreasing, but at lesser rate than that predicted from SPSS or Sci-2.6 Conclusions The addition of ARD and chemical changes to the pit lake seawater have sufficiently changed the chemical constituents so that PSS78 no longer provides an accurate estimate of salinity. Salinity determination by chlorinity, drying and weighing, and chemical constituents were also unsatisfactory. The density-derived salinity proved to be the most accurate estimate of salinity. Based on the density-derived salinity a correction to PSS78 salinity was developed for application to field CTD measurements. The salinity correction was found to be a function of Spss and temperature, and a unique equation for each layer was required. Observations suggest that non-seawater sourced ions are increasing with time, and are sourced primarily from the ARD inflow and dissociation of ions from the pit wall. So the ST - S PSS + AS{ (2.8) 45 salinity correction is only applicable to the date of samples used for its development. To compensate, a time variant salinity correction was developed. If the salinity correction is included the lower layer is found to be more saline than the original seawater. Future changes to the ARD and chemical changes to the pit lake may invalidate the salinity correction, particularly for the lower layer where AS is the largest and the water chemistry the most complex. This work illustrates the difficulty of determining the salinity in a biochemically active system such as a pit lake. It is recommended that salinity and density equations be verified for lakes where the salinity has the potential to dominate or even moderately influence density relationships. The density method is likely to provide a good estimate of salinity, because equation of state formulae are less sensitive to changes in chemical composition than conductivity methods such as PSS78. 46 Table 2.1. Comparison of density from PSS78 (SPSs) and EOS80 (a,iEos) to direct density measurements ((Jt.measured) for samples at 20 °C, where cr,= p(S,T,0) -1000. Sample Depth (m) Sample Date Spss <Jt,EOS (kg m"3) Gtuneasured (kgrn"3) ACT, = Ap (kgrn"3) 5 5-Jun-2000 4.080 1.310 1.524 0.214 29-Aug-2000 4.321 1.493 1.703 0.210 29-Nov-2000 4.233 1.426 1.653 0.227 15-Mar-2001 3.888 1.165 1.403 0.238 100 5-Jun-2000 26.409 18.224 18.510 0.286 29-Aug-2000 26.389 18.209 18.469 0.260 29-Nov-2001 26.364 18.189 18.465 0.276 15-Mar-2001 26.302 18.143 18.415 0.273 300 5-Jun-2000 28.389 19.728 20.817 1.089 29-Nov-2000 28.349 19.698 20.765 1.067 15-Mar-2001 28.347 19.697 20.767 1.070 47 Table 2.2. Alternative salinity measurements of November 29 2000 pit lake samples. Sample Spss Sci SFR SFFR Schein SEOS (%o) (%o) (%c) (%c) (%c) Injection System Water 0.97 - 1.817 1.68 2.31 2.23 Upper Layer (5 m) 4.23 3.95 4.87 3.22 4.36 4.55 Intermediate Layer (100 m) 26.36 26.21 28.94 22.33 25.74* 26.70 Lower Layer (300 m) 28.35 27.68 32.04 25.61 28.94 29.70 Lower Layer (350 m) 28.35 27.71 31.85 26.14 28.96 29.80 Rupert Inlet 28.45 - 30.63 25.68 28.26 28.68 * see note Table 2.5. 48 Table 2.3. Coefficient of haline contraction determined by dilution experiments, fya, compared to EOS80 values, fiEos, where DDW is double distilled water and SIS is south injection system water. Sample PEOS Difference (IO"4) (io-4) (•%) Upper Layer (5 m) diluted with DDW 7.515 7.424 1.21 Intermediate Layer (100 m) diluted with DDW 7.507 7.521 -0.18 Intermediate Layer (100 m) diluted with SIS 7.504 7.367 1.83 Lower Layer (300 m) diluted with DDW 7.510 7.365 1.93 Rupert Inlet diluted with DDW 7.511 7.311 2.66 49 Table 2.4. Coefficient of thermal expansion determined by experiments (14-15°C) compared to EOS80 values. Sample CCEOS 1 Q - 4 T - i 1 0 - 4 T - l Difference % Upper Layer (5 m) -1.546 -1.545 0.06 Intermediate Layer (100 m) -1.968 -1.957 0.56 Lower Layer (300 m) -2.010 -2.013 0.15 Rupert Inlet -1.991 -1.996 0.25 50 Table 2.5. Major chemical constituents and excess concentration (CE = Ciayer - (C l i a y e r / C1RI).CRO where increases are bold. Chemical Constituent Concentration, C (mg L"1) Excess Concentration, CE (mg L"1) Upper Layer Int. Layer Lower Layer (300m) Lower Layer (350m) Rupert Inlet South Injection System Upper Layer Int. Layer Lower Layer (300m) Lower Layer (350m) Chloride (CI) 2189.00 14824.00 15724.00 15728.00 16141.00 5.00 - - - -Sulphate 605.00 2340.00 3056.00 3225.00 2340.00 1794.00 287.61 190.92 776.44 944.88 Bromide 8.00 59.00 50.00 64.00 71.00 0.50 -1.63 -6.21 -19.17 -5.18 Phosphorus 3.00 50.00 25.00 50.00 50.00 0.50 -3.78 4.08 -23.71 1.28 Nitrite 1.33 22.13 13.28 22.13 22.13 0.22 -1.67 1.81 -8.28 0.57 Nitrate 0.33 9.85 9.85 9.85 13.14 1.12 -1.45 -2.21 -2.95 -2.95 Fluoride 0.50 10.00 5.00 10.00 10.00 2.90 -0.86 0.82 -4.74 0.26 Silica 1.80 4.70 10.00 9.90 2.60 23.80 1.45 2.31 7.47 7.37 Sodium 1169.00 7355.00* 8481.00 8402.00 8641.00 38.30 -3.03 -580.99* 63.20 -17.87 Calcium 169.00 452.00 893.00 780.00 298.00 345.00 128.58 178.31 602.70 489.63 Magnesium 134.00 909.00 1038.00 1031.00 1051.00 59.50 -8.55 -56.25 14.15 6.90 Potassium 40.70 246.00 260.00 255.00 313.00 4.20 -1.75 -41.46 -44.92 -49.99 Strontium 1.84 3.84 8.41 6.74 4.82 1.51 1.19 -0.59 3.71 2.04 Boron 0.70 3.05 2.85 2.29 2.96 0.38 0.30 0.33 -0.03 -0.59 Manganese 0.11 1.38 4.27 3.59 0.00 6.92 0.11 1.38 4.27 3.59 Zinc 0.12 0.70 0.50 0.42 0.04 9.13 0.12 0.67 0.46 0.38 Copper 0.01 0.03 0.15 0.15 0.01 1.24 0.00 0.03 0.15 0.15 Aluminum 0.05 0.14 0.27 0.05 0.13 23.30 0.03 0.02 0.14 -0.08 * The laboratory measurement of intermediate layer sodium is suspected to be erroneous, as the intermediate layer sodium concentration and Schem are ~ 600 mg L" 1 lower than anticipated without any obvious explanation. 52 1026 1008 -I 1 1 1 1 1 10 15 20 25 30 35 Salinity Figure 2.2. Density versus salinity for Island Copper samples diluted and evaporated with double distilled water (DDW) and south injection system (SIS) water. The density-derived salinity, SEos, and salinity computed by volumetric dilution/ evaporation, Sda, are plotted. Figure 2.3. Density of intermediate layer water diluted with injection system water [A] compared to intermediate layer water diluted with double distilled water [ A ] . 1 -0.2 J Days Figure 2.4. Non-seawater salinity indicated by SPSS less Sa versus days since pit lake flooded, for upper layer [•], intermediate layer [A] and lower layer [•]. Figure 2.5. Salinity trends for (a) upper layer, (b) intermediate layer and (c) lower layi showing SA SPSS [A] and SPSS+ AS, [*]. 56 CHAPTER 3 Filling Box Behaviour of a Plume Discharged into a Meromictic Pit Lake Abstract The mine pit at Island Copper Mine, Vancouver Island, BC, was flooded with seawater in 1996, to create a meromictic lake for the treatment of acid rock drainage (ARD). The pit lake has developed into three distinct layers: a brackish and well mixed upper layer; a plume stirred intermediate layer; and a thermally convecting lower layer. The intermediate layer is analogous to the filling box model (Baines and Turner, 1969), as ARD is injected at depth and rises as a buoyant plume within the confined receiving environment. The gradients of salinity, temperature and density observed in the intermediate layer are consistent with the filling box model. However, application of the filling box model is limited because of the complexities of the pit lake: the injection system diffuser geometry and variable performance, bottom heating from the warmer lower layer which causes convective mixing in the bottom of the intermediate layer, and the intermittency of the injection system flow which tends to stratify the top of the intermediate layer. Furthermore, during the summer, heating of the top of the intermediate layer and minimal injection system flow cause the plume to be trapped within the intermediate layer. 3.1 Introduction The filling box model has been experimentally investigated (Baines and Turner 1969, Baines 1975, Kumagai 1984) and used to model the behaviour of the atmospheric boundary layer (Baines and Turner, 1969) and deep sea basins (Mannis, 1973) and lakes (Killworth and 57 Carmack, 1979). The filling box model explains how a source of buoyancy in a confined receiving environment circulates fluid and controls the gradients of the buoyant source (temperature or salinity) and density in the receiving environment. The Island Copper pit lake is presented as an example of filling box behaviour and raises some limitations of the model. The pit lake at Island Copper Mine was developed as part of the mine closure plan for the treatment of acid rock drainage (ARD) (BHP 1996; Chapter 1). The mine pit was flooded from 15 June until 23 July 1996 with seawater from Rupert Inlet and capped with fresh water to form a meromictic lake. ARD is collected from on-land waste rock dumps and discharged via two injection systems into the intermediate layer of the pit lake, where it rises as buoyant plumes through the more dense salt water (Figure 3.1). The upper 'brackish' water layer forms a cap over the intermediate layer, preventing the plumes from surfacing and impeding oxygen transfer to the intermediate layer. A distinct lower layer has formed below the injection system discharges. Background information on the development of the pit lake and passive treatment system is detailed in Chapter 1. Of significance to the current chapter are details specific to the discharge of ARD to the intermediate layer. The south injection system (SIS) and north injection system (NIS), see Figure 3.1, drain different areas of the waste rock dump. ARD is collected by ditches and routed to sedimentation ponds. It is then discharged into the pit lake by static head via two high-density-polyethylene outfall pipelines (with diameters of 2 feet (0.6 m) for the NIS and 3 feet (0.9 m) for the SIS) at a depth of -220 m (UWL). The two diffusers consist of 3 sets of 200 mm discharge ports at quarter points around the circumference, spaced along a 14 m section pipe, and the 428 mm open end of pipe. Inspection of the injection systems by a submersible remotely operated vehicle revealed that the SIS diffuser has only three ports open (out of 12) plus the end of pipe, while the NIS has all ports open (Rescan, 2001). Both diffusers were lying at an angle of -35° below horizontal. Furthermore, both diffusers were only discharging from the first set of diffuser ports, which indicates that the intrusion of saline pit lake water was causing the diffusers to be partially blocked. The degree of blockage may be flow dependent since at higher flows the diffuser may purge (Wood et al. 1993) and more ports may become operational. 58 The NIS discharge was diverted to the upper layer of the pit lake on the 29 June 1999 in order to reduce the oxygen input to the intermediate layer, reduce the rate of upper halocline rise, improve the upper layer alkalinity and to reduce the long-term risk of upwelling (Chapter 1). An old conveyer tunnel discharges drainage into the intermediate layer at a depth of 154 m (UWL). The conveyer tunnel flow is estimated to be 4% of the total injection system flow, based on a comparison of catchment areas. In this chapter the physical limnology is investigated focusing on the intermediate layer and plume dynamics with particular focus on similarities with the filling box model. The chapter reviews plume and filling box theory appropriate to the pit lake. Field methods are introduced and the results of conductivity-temperature-depth (CTD) profiling and thermistor chain measurements describing the pit lake physical structure are given. The discussion focuses on interpretation of the results and comparison to the filling box model predictions. 3.2 Literature Review Axisymmetric Plumes - The mean properties of axisymmetric plumes in a homogeneous environment where developed by Morton, Taylor and Turner (1956) by dimensional analysis and experiments. The centreline time averaged velocity (w), the radius (b) where the velocity has fallen to lie of w for an assumed Gaussian shape and buoyancy (g'=g(pP-Pa)lp a), are functions of the buoyancy flux at the source (B), height above the nozzle (z) and the entrainment coefficient (a= 0.1). B = Q g ; M , 6 (3.2) b=-az 5 l (3.3) 5 f 18 aB^' w = — 6a y5x z , 59 g=Yn 5 (5n\3 18, 1 l (3.4) a4z5 The initial plume buoyancy is given by g0 / = g[p0-pa]/pa- Q is the initial plume discharge and g is gravitational acceleration. The different densities in use are the initial plume (Po), plume (pp) and the ambient reference densities (pa) taken as the average density of the receiving environment. These plume equations are valid for a buoyant jet (discharge with momentum and buoyancy) for z » l,„, where the characteristic length scale lm = M ViB~Vl and momentum flux M = QV, with V as the mean discharge velocity. Equations 3.1-3.4 assume a virtual origin (which does not coincide with the port location) and that z is in the zone of established flow (where the flow is self-similar and plume-like). For jets the virtual origin and start of the zone of established flow are ~ 4 and ~ 10 port diameters downstream (Fischer et. al. 1979), so for deep discharges the origin can be assumed to be at the port. In the case of a stratified environment a terminal height of rise (ht) is reached when the density of the plume fluid equals that of the ambient environment. ht =3.8 JL (3.5) Z?4 3 where the variation of density within the receiving environment is characterized by e' (Fischer et al. 1979) £ , = z]Ldp^ (3-6) Pa dz where P2 is the local density (subscript 2 to indicate the intermediate layer of the pit lake). Potential density must be used for pi in Equation (3.6) to remove pressure effects so that s represents the stability. If in situ density is used the neutral case appears very stable since pi decreases with height as pressure drops with height. N = (e'g)'0'5 is the natural (or buoyancy or Brunt Vaisala) frequency. 60 Filling Box Model - Baines and Turner (1969) proposed the "filling box" model to describe the effects of continuous convection from a buoyant discharge on the properties of the confined receiving environment. In the filling box experiment a plume rises, entraining fluid from its receiving environment. The plume arrives at the top of the fluid container with buoyancy and spreads out laterally to form a layer of light fluid. The continuous plume now entrains light fluid from this layer and hence arrives at the top of the container even lighter. The more recently arrived plume thus spreads out above the plume layer, displacing the latter downwards. The downward flow of fluid replaces the fluid entrained into the plume from the ambient environment and itself is entrained back into the plume. In this way, the fluid in the confined receiving environment is circulated and a stratified profile is produced (see Figure 3.2). The filling box model assumes the plume spreads out instantly across the top of the region (of height H, radius R) and becomes part of the non-turbulent environment at that level. The buoyancy of the receiving environment (#2') changes at a fixed height due to the vertical motion given by dt dz 2 a n d -R2U = b2w <3-8) where gi'- g(pi-pa)lpa and U is the downward vertical velocity of the receiving environment (provided b « R). The filling box time scale, tp,, is the time for the first front of plume fluid to descend the depth of the container and is similar to a circulation time. 1 (3.9) 1 ( 7TR6 ^ ya4BH2 The dependence of the buoyancy of the receiving environment on time can be eliminated if g2'is assumed to be linear in t. This is a reasonable assumption since the result of steady heating (or addition of freshwater) for a long time should be a uniform decrease in density through the whole of the environment. Baines and Turner's (1969) solution was 61 based on dimensional analysis of B and H and scaling of the physical variables z and t, expressed in their simplest non-dimensional form as C,- z/H and r= t/fa, (3.10a) 82 (AnfaAH5 where (3.10b) f=g 3 (3.3-0.84c;-0.062c;2)-2.4 where/is a shape function determined from algebra and integration. The solution is valid for large non-dimensional time t » fa. The buoyancy of the receiving environment changes linearly in time at all heights, while the shape of the density profile remains fixed. The density distribution is independent R, but R enters into the scaling of fa. The filling box model is valid for aspect ratios, H/R less than 1, because for greater aspect ratios the plume inertia overcomes the buoyancy and rapid overturning occurs instead of the filling box process (Baines and Turner, 1969). Worster and Huppert (1983) developed an approximate analytical expression for the time dependence of the density profile, which tends to the Baines and Turner solution at the limit t —» 00. Killworth and Turner (1982) found that the asymptotic state for time-varying buoyancy to be qualitatively very similar to that for a steady plume with buoyancy equal to the maximum value of the time-varying buoyancy. Plume Entrainment of the Halocline - Baines (1975) and Kumagai (1984) extended the filling box model for the case of a plume discharged to a lower layer capped by an upper layer of lesser density. Upper layer fluid is entrained into the lower layer through the end of the plume as it impinges against the density interface. The addition of less dense upper layer fluid to the horizontally spreading plume, modifies the density stratification of the receiving environment (3.9) such that/2 becomes f (3.11) ( i - / 0 1/3 3.3-0.84 -0.062-i - / 0 ( i - / 0 ' •2.4 62 where /? = B*/(B*+B) and B* = Q.Egn is the buoyancy entrained from the upper layer. QEis the entrainment flux and §12' = g{pi-p\)l(h is the buoyancy difference between the plume (effectively the lower layer) and the upper layer (p\). Baines (1975) modification of Baines and Turner's (1969) model is not large and need not be included if entrainment is small; in fact if B* « : B and Equation (3.11) becomes (3.10b). Baines (1975) observed that the entrainment was confined to a region about the size of the plume cross-section, and established that QE/WD2 O C Ri3'2, a result that was confirmed by Kumagai (1984). The Richardson number, Ri = gn d/v2, where v and d are the characteristic velocity and width of the incident flow. Fernando (1991) in an extensive literature review, found the exponent, °= in Ria, to be dependent on conditions and that it varied from -1/2 to -2. Cotel and Breidenthal (1997) proposed a model for entrainment that includes Ri dependence plus Reynolds and Schmidt numbers and the vortex persistence. Plumes are non-persistent as they are intermittent in nature and their entrainment follows that of discrete vortex rings. In Regime III the baroclinic vorticity, responsible for engulfing fluid from the upper layer, is limited as even small eddies are stratified and entrainment follows RF (as per Baines (1975) and Kumagai (1984)). As Ri increases the interface becomes flat. Regime IV is define by Ri > Sc (Schmidt number, Sc = vID, where is the D molecular diffusivity, so that Sctemp -10 and Scsait -1000) and molecular diffusion limits entrainment, which follows Sc'xl2RiA. Cotel et al. (1997) observed a dramatic change in flow character at Ri = 10 for a thermally stratified system, from a highly convoluted interface at lesser values to gentle undulations with no overturning vortices at Ri > 10. Further increase in Ri results in Regime V, which is defined by Ri > Rem (Reynolds number, Re = wb/v, where v is the kinematic viscosity). For Regime V the stratification is so dominant that mass can only be transported across the interface by diffusion with gradients enhanced by surface renewal and the transport rate follows Sc2l3Rim. 3.3 Methods Fieldwork - Conductivity-temperature-depth (CTD) profiling has been carried out on 18 occasions: July, August, September and November 1997; April, August and December 1998; March, June, September and December 1999; February, March, June, July, August and 63 November 2000; and March 2001. The CTD instrument, an Ocean Sensors' Model OS200 averaged 99 raw measurements every -0.5 s. Laboratory calibration of the CTD was performed before and after each trip. Field verification of salinity was provided by bottle samples tested in a Guildline Model 8400 Autosal and by occasional use of supplementary CTDs. The temperature accuracy and resolution were 0.02 °C and 0.001 °C, and for salinity they were 0.03 and 0.001, respectively (Appendix C). A depth uncertainty 0.23 m was accomplished by comparing against in situ thermistors and a Sea Bird Electronics SBE19 CTD (Appendix C), and a depth resolution of 0.15 m was achieved in interface areas by slowing the CTD descent rate. Thermistors were deployed from December 6, 1999 until November 29, 2000 in five consecutive deployments. One thermistor chain was located over the south injection system diffuser (SISTC) and a second thermistor chain was located near the West Station and was named the central mooring (CMTC) (Figure 3.1). The results of only the CMTC are reported here as it had a better depth coverage. Instrumentation was concentrated on the upper layer, upper interface and lower interface where temperature changes were expected. Instruments were moved around to optimize their usage (Appendix D). As an example for 3 February until 13 March 2000 instruments were located at depths of 1, 2, 3, 4, 5 x , 6, 7T, 7.5 T , 8 T , 8.5 T , 9T, 9.5T,10T, l l x , 12, 15, 20, 30 T , 100 x , 200w, 220 w , 222 T , 224 T , 227.5 T , 231 T , 235 w and 250 x m. For 8 June - 29 November 2000 five fewer instruments were used. Brancker TR1000F [T] had an accuracy of 0.02 °C and sampling rate of 30 s. TSKA Wadars [w] had an accuracy of 0.02 °C and sampling interval of 10 min. Onset Stowaways had an accuracy of 0.2 °C and sampling interval of 10 min. Brancker XL200 [x] had an accuracy of 0.1 °C for temperature and sampling interval of 10 min. Calibration details are given in Appendix D. Injection system inflow measurement started on 27 November 1997 using acoustic flow meters at each of the injection system outfall pipes. Water level was monitored manually with a staff gauge until January 1999, and by an ultrasonic water level meter since then. Water levels were used to correct CTD and thermistor depths to the UWL datum. Salinity - In the pit lake, salinity is the dominant contributor to density differences, so accurate determination of salinity is essential to understanding the physical mixing processes 64 and the long-term stability of the stratification. The accurate determination of true salinity (defined as mass of dissolved salts in one kilogram of solution) is not a trivial matter. Wilton (1998) used the Practical Salinity Scale 1978 (PSS78) (Perkin and Lewis 1980) to describe the pit lake salinity from CTD measurements. However, true salinity in the pit lake is significantly different from PSS78 salinity, as the chemical composition of the originally seawater pit lake has changed due to the addition of ARD, dissolution of ions from the pit walls, and other chemical changes. To calculate the true salinity from CTD measurements a time and layer dependent correction is applied to PSS78 salinity. Chapter 2 discusses methods of salinity determination and the salinity correction relationships. 3.4 Results Development of three layers - The pit lake has evolved into a stably stratified three-layer system. The CTD profile taken 6th December 1999 was typical of other dates, with salinity, temperature, potential density and dissolved oxygen being approximately constant in each layer (Figure 1.3). Layers were sharply defined by the upper and lower haloclines, which separated the upper from intermediate, and intermediate from lower layers, respectively. The lower halocline forms just below the level of the injection system discharges. Chapter 1 detailed general trends for salinity and temperature, these are briefly summarized with values given for 15 March 2001. The upper layer salinity has decreased to 4.08 while the temperature exhibits seasonal variation. The intermediate layer salinity [26.62] and temperature [12.02 °C] have been slowly decreasing due to the injection of ARD to that layer. Meanwhile, the lower layer salinity [29.78] has increased due to the dissociation of ions from the pit walls and the lower layer temperature [12.90] has increased due to geothermal heating. Upper Layer - The upper layer salinity was uniform for all CTD casts. Thermistor data provides a detailed view of the temperature changes. The upper layer becomes temporarily thermally stratified in the spring and summer until wind mixes it, and remains well mixed in the fall and winter due to cooling by penetrative convection (Figure 3.3). Of interest was that the top of the intermediate layer continued to warm in the fall from solar heating, as it was protected from surface cooling and penetrative convection by the upper halocline. 65 Intermediate Layer - Subtler than the general trends, but crucial to understanding the hydrodynamics of the intermediate layer, were the slight gradients of temperature, salinity and density with depth observed in the intermediate layer. The potential temperature (6) was used to eliminate the effect of adiabatic temperature change. The in situ temperature (7) increases by -0.01 °C for a 100 m depth increase. The potential density (= pg [S, 6, P = constant]) referenced to atmospheric pressure was expressed for convenience as <Jg = pe -1000 [kg m" ]. Equation (3.13) demonstrates that potential density decreases with height (negative or stable potential density gradient) due to increasing temperature (positive temperature gradient for 6> 0) and decreasing salinity (S) (negative salinity gradient). do Q dpfl dz dz d6 adS dz dz 1 dp I dp (3.13) where a - p =——— p d6 pdS d) («) The intermediate layer a and /3 were -1.96xl0"4 °C _ 1 and 7.5xl0"4, respectively, from measurements described in Chapter 2. Height (z) was positive upwards. Salinity decreased with height throughout the intermediate layer and gradients were more negative during winter (Figure 3.4a). The top of the intermediate layer (down to 90 m) was heated during the summer by solar radiation (Figure 3.4b), consistent with the loss of visible light at - 90 m (Rescan, 2001). The potential temperature gradient from 100 m to 200 m was weakly negative in winter and weakly positive in summer. The gradients were at the limits of our CTD accuracy so there was some instrument noise and nonsensical patterns. The resulting potential density profiles (Figure 3.4c) were stable as changes in salinity dominate density; a one unit change in salinity changes the density by -5 times that caused by a one unit change in temperature. The stability (e) (Equation 3.6) can be divided into the contribution from changing temperature (Eg) and salinity (Es), which are calculated from Equation 3.13 terms (i) and (ii), respectively. By definition the stability is positive if the density profile is stable. The stability for the intermediate layer (30-200 m) was strongly seasonal (Figure 3.5), which was consistent with observations made from the raw profiles. The stability was dominated by the contribution to density from salinity and was at all times stable. The stability due to salinity was strongest in winter and weakest in summer, whereas the stability due to temperature was 66 positive in the summer and negative in the winter. The stability was smaller if only the bottom (100-200 m) of the intermediate layer was considered. Still the salinity effected remained positive and seasonal, while temperature effect was near zero. Thermistor data provides a detailed picture of the intermediate layer temperature structure and its dependence on injection system flow (Figure 3.6). The intermediate layer had a negative (unstable) temperature gradient from December 1999 to February 2000 when flows were high and ARD temperatures were less than 12 °C. During this period the negative temperature gradient weakened concurrent with a decrease in injection system flow. The layer became essentially uniform with a temperature of = 12 °C in March and April 2000. The top of the intermediate layer warmed from May to October to a depth of -90 m (UWL). The positive (stable) temperature gradient during this time was a result of solar heating, rather than warm ARD inflow, as flows were minimal from July to mid-October. With the onset of higher injection system flows from mid-October to November the layer begins to mix starting at the bottom of the intermediate layer. Mixing is indicated by more uniform temperatures; the top of the layer cools, while the bottom of the layer warms. The thermal stratification above 25 m (UWL), but below the upper interface, remained until the last day of record 29 November 2000, indicating that the plume had yet not penetrated this high. However, by 15 March 2001 the date of the next CTD profile the layer was back to a winter structure with negative temperature and negative salinity gradients (Figures 3.4 and 3.5). Lower Layer - Salinity, potential temperature and crgweTe more uniformly distributed in the lower layer than in the intermediate layer (Figure 3.7). Potential temperature has no gradient and neither does potential density. The absence of a density gradient suggests a well mixed layer. Chapter 1 proposed thermal convection as an explanation for the observed erosion of the lower halocline from below. The lower layer is geothermally heated from below and cooled from above. The Rayleigh number of 4.7 x 1014 was far larger than the critical Rayleigh number required for turbulent convection. The apparent increase in salinity at the top of the lower layer is not believed to truly represent an increase in salinity, but rather it is a manifestation of the changing chemical composition observed between the intermediate and lower layers (Chapter 2). There was a 67 unique salinity correction for each of the three layers in the pit lake and within the interfaces the salinity correction was proportionally applied based on the change in salinity between defined interface shoulders (Appendix A). This simple correction scheme does not adequately correct changing water chemistry at the top of the lower layer and a slight over-calculation of the salinity results. 3.5 Discussion Plume Dynamics - Driving the intermediate layer dynamics is the plume buoyancy flux (Equation 3.1). The plume buoyancy (g0') was calculated for November 1997 to March 2000 from densities of the intermediate layer (via CTD records) and ARD (calculated with EOS80 from estimates of SIS salinity and temperature). The ARD salinity measured in the SIS pond was variable and without any discernable trend, and a mean salinity of 0.86 was used based on 1997-2001 water quality measurements, which had a standard deviation of 0.23. The Port Hardy Airport air temperature plus 2.57 °C was used for ARD temperature, based on its good correlation [correlation coefficient = 0.94] with a thermistor in the SIS pond during 2000 (Appendix D). gd fluctuates seasonally due to ARD temperature and decreases each year as a result of the gradual decrease in intermediate layer salinity (Figure 3.8a). The buoyancy flux (B) shows large seasonal variance mostly due to the variation in injection system flow (Figure 3.8b). The filling box model (Equation 10) predicts the density gradients to increase for large buoyant flux (g2 °= Bm), and during the winter when the buoyancy flux was large the observed density gradients increased (Figures 3.4c and 3.5). The variable buoyancy flux affects the plume properties (Equations 3.2-3.4), the plume height of rise (Equation 3.5) and the filling box characteristics such as the density of the receiving environment (Equation 3.10) and filling box time scale (Equation 3.9). Assuming for simplicity that each diffuser behaves like an axisymmetric plume in a linearly stratified environment (values extrapolated from CTD casts, see Figure 3.5), then the predicted plume terminal height of rise varies from zero to heights exceeding the intermediate layer depth (Figure 3.8c). Heights exceeding the intermediate layer depth were not possible because of the strong density step at the upper halocline forces the plumes to spread horizontally under the upper halocline. These predictions are both higher and lower 68 than the height of rise observed during dye tracer studies (Muggli et al. 2000), which emphasizes the difficulty of modeling the diffusers. Dye tracer experiments on 9 October 1999 at NIS (B = 0.004 m 4 s"3) observed ht= 130 m compared with predicted h,= 160, and 3 March 2000 at SIS (B = 0.034 m4 s"3) observed h,> 210 m compared with predicted h,= 203 m. The plume and the density gradient are intrinsically linked, as the plume creates the density gradient by the filling box process, and in turn is influenced by it. However, solar heating also exerts an influence by increasing the density gradient at the top of the intermediate layer, thus trapping plumes at lower depths than they would otherwise have reached. Plume Entrainment - Plume entrainment of the upper halocline can be analyzed within the Cotel and Breidenthal (1997) framework based on Richardson, Reynolds and Schmidt numbers (Table 3.1). The plume initial buoyant flux, B, is calculated using Equation (3.1) from daily flow measurements and density estimated from CTD data. The maximum and average B for the high flow months of November 1999 to January 2000 have been used. A homogeneous receiving environment is assumed for simplicity and w is calculated via Equation (3.2) using g0'= 0.19 m s"2, and via Equation (3.3) b = 25.2 m using z = 210 m. Ri is then calculated with the measured 1999/2000 winter gx{ of 0.17 m s"2. By Cotel and Breidenthal's (1997) classification, plume entrainment at the upper halocline will be in Regime V (Ri > Rem and Ri < Scsan -1000) with limited entrainment following Sc'2/3Ri~ l / 4 . However, for high plume flows Ri decreases and entrainment increases to Regime III (Ri < Rem) and follows Ri~ 3 / 2 (Table 3.1). The plume initial buoyancy flux to transition from Regime V —» Regime III can be calculated by substituting Equations (3.2) and (3.3) into Ri - Re1'4, making the assumption of a homogeneous receiving environment. D „' 4 / 3 , /V 3 B = g n 1 v' '6aN 4 ( 57t] I5 J [l8aj 2 Z (3.12) The plume initial buoyant flux for the transition from Regime V to Regime HI is 0.084 m 4 s"3 for gn' = 0.17 m s"2 and z = 210 m. Buoyancy fluxes greater than 0.084 m4 s"3 occur for only - 1.5% of days (Figure 3.8b). A stratified receiving environment would 69 reduce w, increase Ri and make Regime V more likely. Thus for the current stratification entrainment is minimal and the Baines (1975) modification of Baines and Turners (1969) model need not be included. Transition to the higher entrainment Regime III may occur in the long-term due to weakening of the stratification and decrease in g\{. This transition would immediately increase the entrainment by an order of magnitude and precipitate destratification by erosion of the upper layer. Note, while the entrainment flux is small compared to the incident flow, it may still be significant compared to the initial plume flow. Baines (1975) experimental relationship QE°cRiy2 (equivalent to Cotel and Breidenthal's (1997) Regime III) extended up to Ri = 25 and at this value QEITTWD2 = 0.001, which corresponds to QEIQ = 0.7. Entrainment by the plume at the upper interface was inferred from observations as QEIQA = 0.3 (Appendix A), where QA is the total injection system inflow. The Baines' (1975) prediction and observed entrainment do not match because the injection system plumes and halocline have higher Ri and perhaps different entrainment behaviour (i.e. Cotel and Breidenthal's (1997) Regime V). Plume entrainment was evident from observations that the thermocline over the SIS was higher and thinner than that found in the centre of lake during periods of high flow (Chapter 4). Filling Box Model - The hypothesis that the intermediate layer of the pit lake behaves like the filling box is not unreasonable given the physical similarities between it and filling box experiments. Indeed the observed stratification is consistent with the filling box process, where a fresh cold plume discharged into saline warm water from below will produce a positive (stable) salinity gradient and positive (unstable) temperature gradient, provided there is sufficient flow and that salinity dominates the density. Variables such as container shape, variable buoyancy flux and multiple sources can be accounted for in the filling box model. In reality development of a specific filling box model for the intermediate layer is too difficult because of the complexities of the diffuser geometry and functionality; some ports are not open and others may be blocked by the intrusion of intermediate layer water, which may be dependent on flow conditions. In addition, solar heating and heat transfer from adjacent layers influence the density gradient of the intermediate layer. Therefore, the filling box model is only used as a conceptual model for qualitative comparison. 70 An indicator of filling box behaviour is the filling box time scale (Figure 3.8d). The filling box time scale drops to ~ 4 days during winter months, which is much less than the seasonal time scale, but close to the time scale for frontal weather systems on the BC coast (Chapter 4; Appendix B). Therefore, during the winter when fa is small and plumes rise the full height of the intermediate layer, it might be reasonable to compare the density gradients observed in the intermediate layer to the filling box asymptotic solution of Equation (3.10). The filling box time scale increased dramatically and became more irregular after 29 June 1999, as the NIS was diverted on this date so the injection system flow was only from the SIS, which had minimal flow in the summer and was generally more intermittent than the NIS. The filling box model compares reasonably to the observed winter density gradients for 0.5 < C < 0-9, but the shape of the filling box density distribution is unlike the observed profiles and large differences are observed at the top and bottom of the layer (Figure 3.9). Convection at the bottom of the intermediate layer caused by heating from the lower layer, which was -0.8 °C warmer than the intermediate layer, is proposed to explain the difference at the bottom of the intermediate layer. Baines and Turner (1969) experimentally observed the increasingly sharp gradient in density distribution near the source as predicted by the filling box model, but noted for real applications, such as the atmospheric boundary layer, there would be a well mixed layer close to the boundary from convective or mechanical mixing. They proposed a well mixed layer of depth SH, on top of which would be a layer of depth (l-S)H with a stable density distribution formed by plumes, where £was the fraction of the heat flux that goes into heating the lower layer. Baines and Turner (1969) calculated that for S= 0.1 the density gradient would change by only 1% at the most compared to the previous asymptotic solution. The deviation of the observed profiles from the filling box model at the top of the intermediate layer is somewhat harder to explain. Convection is expected in the winter when the overlying upper layer is cooler than the intermediate layer, but it would act to increase vertical mixing and reduce g{ not to enhance the density stratification as observed. Entrainment of upper layer fluid by the plume impinging against the upper interface would draw fresher less dense water into the top of the intermediate layer. Baines (1975) and 71 Kumagai (1984) accounted for this plume entrainment in the filling box models and predicted a curve for the receiving environment density similar to the observations shown in Figure 3.9. In reality the intermittency of the ARD flow, which is largely filtered out of the weekly average buoyancy flux (Figure 3.8b), may be the cause for the strong stratification at the top of the intermediate layer. ARD flows are peaky due to the rapid runoff of surface water from the waste rock dumps following storm events. The peak of the storm hydrograph has a duration of hours, rather than the days required for the filling box time scale. Plumes with higher buoyant fluxes are more buoyant at a given height (3.4) and thus rise higher in a stratified ambient (3.5). These plumes will collect at the top of the intermediate layer because subsequent plumes with less buoyancy flux are less buoyant and will not rise as high. Entrainment by the plume of upper layer fluid increases with buoyant flux (Chapter 1), so the density stratification at the top of the intermediate layer may be further enhanced by the addition of entrained brackish water from the upper layer. Baines and Turner (1969) also noted that the shape of the container can have a profound effect on the density distribution of the receiving environment. For a wedge shaped region (not unlike the intermediate layer) with a line plume located at its vertex the density distribution changes to be nearly linear with height. Furthermore the time to achieve the asymptotic state is reduced. The complexities of the pit lake limit the applicability of the filling box model. These complexities include: the injection system diffuser geometry and variable performance; bottom heating of the intermediate layer from the warmer lower layer; solar heating of the top of the intermediate layer; variability of the injection system flow and unique bathymetry. A process based lake model would be better adapted at dealing with these issues. In Chapter 5 a one-dimensional layer model DYRESM, which includes a plume routine, is used to simulate the lake plume behaviour. 3.6 Conclusions The intermediate layer is similar in appearance to the filling box model (Baines and Turner, 1969) as ARD is injected at depth and rises as a buoyant plume within the confined 72 receiving environment. The filling box model analogy explains the gradients of salinity, temperature and density observed in the intermediate layer. Plume entrainment of upper layer fluid is not important to the filling box model because entrainment is relatively small. However, application of the filling box model is limited because of the complexities of the diffuser, bottom heating from the warmer lower layer which causes the bottom of the intermediate layer to mix, and the intermittency of the injection system flow which tends to stratify the top of the intermediate layer. Furthermore, during the summer heating of the top of the intermediate layer and minimal injection system flow cause the plume to be trapped within the intermediate layer. Detailed modeling is required to resolve this complicated behaviour and in Chapter 5 a one-dimensional layer model DYRESM is used to simulate the plume behaviour. 73 Table 3.1: Classification of entrainment regime by Cotel and Breidenthal (1997). Parameter NIS SIS Average*1 Max.*1 100 year*2 Average*1 Max.*1 100 year*2 Q ( m V ) 0.131 0.798 1.970 0.128 0.476 1.130 B ( m V 3 ) 0.025 0.151 0.374 0.024 0.091 0.215 w (m s"1) 0.199 0.367 0.492 0.198 0.307 0.409 Re 5.0xl06 9.2xl06 1.2X107 5.0xl06 7.7xl06 l.OxlO7 Re'/4 47 55 58 47 53 56 Ri 104 31 17 105 44 24 Regime V III III V III III *' From November 1999 to Jan 2000 daily flow measurements. * 2 Hydrological estimate of 100 year return period flow BHP (1996). 74 (a) 350 4 . 1 1 1 1 — i 1 r— 0 250 500 750 1000 1250 1500 1750 2000 Distance (m) Figure 3.1. Island Copper pit lake (a) bathymetry and (b) cross-section AA showing thermistor and meteorological station locations. CMTC and SISTC show the location of the central mooring and south injection system thermistor "chains". Horizontal datum is WGS-84 and depths are referenced to the ultimate water level (UWL) (adapted from Muggli et al. 2000). 75 Figure 3.2. Sketch of the development of stratified environment due to buoyant fresh/cool source discharging into a saline/warm receiving environment, (a) Profiles of the environment and plume (p) density (p), salinity (S) and temperature (T). (b) Motions of the plume and environment fluid. The position of the first front of descending environmental fluid and the resulting profiles are given at two times -1] [—] and t2 [—] (adapted from Baines and Turner, 1969). 76 Date Figure 3.3. Daily average temperature of upper layer, upper interface and top of intermediate layer for 9 December 1999 until 28 November 2000 with vertical breaks indicating thermistor deployments. Thermistors were spaced at 1 m intervals in the upper layer and 0.5 m intervals in the upper interface and depths are corrected for fluctuating lake surface level. The upper halocline shoulders [—] are based on 7 CTD casts undertaken during the period. 77 Figure 3.4. Intermediate layer profiles of CTD data for (a) salinity, (b) potential temperature and (c) sigma-9. 78 -l.E-07 Jul-97 Jan-98 Jul-98 Jan-99 Jul-99 Jan-00 Jul-00 Jan-01 Jul-01 Dates Figure 3.5. Intermediate layer density gradient (e) [—] and contribution from temperature (ee) [ ] and salinity (es) [—], by least squares regression of CTD data from 30-200 m. 79 (a) 200 Dec Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec (b) 4 ~ 3 m o ft ! x 10 1 — 1 1 1 1 I I I I A A 1 f V> [v V ^ --\7 \ 1V \ A / l ^ t t 1 \ ' " I n I v\f v^/v—' v... 1 1 1 1 1 i . I . I i . l 22 12 I O J H Dec Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec Figure 3.6. (a) Daily average temperature of the intermediate layer for 9 December 1999 until 28 November 2000 with breaks between thermistor deployments. Thermistors were spaced at 12 m, 15 m, 20 m, 30 m, 100 m and 200 m. (b) Total injection system flow [—] and temperature [••••]. 80 Figure 3.7. Lower layer profiles of CTD data for (a) salinity, (b) potential temperature and (c) sigma-0. 81 0.22 0.18 0.05 300 150 I 1 — 1 1 1 1 1 (c) \ l AAAA fl \A-Jul97 Jan98 Jul98 Jan99 Jul99 JanOO JulOO Figure 3.8. Plume characteristics for north injection system ["] and south injection system [—] as weekly averages: (a) initial plume buoyancy (g0'); (b) buoyancy flux (B); (c) height of rise (ht) predicted for axisymmetric plumes in a linearly stratified intermediate layer compared to the observed heights of rise for tracer experiments at NIS [A] and SIS [•] (Muggli et al. 2000) and height of intermediate layer [ ]; and (d) filling box time scale (tfb). 82 Figure 3.9. The non-dimensional height vs non-dimensional density gradient in the intermediate layer for winter profiles compared to filling box model (Equation 3.10). 83 CHAPTER 4 Internal Waves at the Upper Halocline Abstract The mine pit at the Island Copper Mine near Port Hardy, Vancouver Island BC, was flooded in 1996 with seawater and capped with an upper layer of fresh water to form a meromictic lake. The lake is being used as a passive treatment system for acid rock drainage (ARD), which is discharged to the intermediate layer at a depth of 220 m as a buoyant plume. The upper halocline, which separates the upper layer from the intermediate layer, was exceptionally sharp with a salinity change of 22.5 over less than 1 m. The plume actively impinges against the halocline, lifting and eroding it during the winter and fall when there is sufficient flow. Spectral analysis was used to identify horizontal modes of the vertical mode-1 lake-length seiche and their periods matched those predicted from the Defant procedure. Additional waves were attributed to vertical mode-1 lake-width seiches. Winds were sufficiently variable in direction to generate seiches in both directions. The mixed layer behaviour was consistent with Spigel and Imberger's (1980) Regime 4. Seiches persisted for at 0(100) wave periods (many days) and degenerated by damping, not by non-linear steepening and the development of solitons as predicted from the results of Horn et al. (2001). The damping period was overestimated suggesting that the parameterization for wave dampening based on boundary mixing over the lake bottom, was not sufficient for the deep pit lake with small amplitude waves. 4.1 Introduction The Island Copper pit lake is a natural laboratory for studying environmental fluid mechanics and limnology because of the strong salinity stratification and interaction with 84 buoyant plumes. The focus of this chapter is on observations of internal gravity waves on the unusually strong halocline of the pit lake. This extreme case is useful as a validation of the existing theory of internal gravity waves. The pit lake is classified using the mixed layer regimes of Spigel and Imberger (1980) and the wave degeneration regimes of Horn et al. (2001) and compared to observations. The influence of the injection system plumes on the halocline is also investigated. Island Copper Mine - The Island Copper Mine was flooded with seawater in June-July 1996 and capped with fresh water, thus forming a meromictic (permanently stratified) lake for the treatment of acid rock drainage (ARD) (BHP 1996; Chapter 1). ARD is collected from waste rock dumps and injected into the lake via two outfalls - the south injection system (SIS) and North Injection System (NIS) - at a depth of ~ 220 m (Figure 1.2). Three distinct layers have developed (Figure 1.3). The upper layer [S = 4.08 in March 2001] is important to the development of the passive treatment system because it prevents the two injection system plumes from surfacing and impedes oxygen transfer to the intermediate layer. The intermediate layer [S = 26.62 in March 2001] was observed to circulate as a "filling box" (Baines and Turner, 1969) during the winter when injection system flows were high and the plumes rose the full height of the layer (Chapter 3). In the summer, the plumes are less buoyant due to reduced injection system flow, and become trapped at a terminal height within the intermediate layer (Muggli et al. 2000). The lower layer was observed to be homogeneous in temperature and salinity as a result of turbulent thermal convection from geothermal heating (Chapter 1; Chapter 3). The upper halocline (separating the upper and intermediate layers) has been closely monitored due to its importance to the long-term stability of the stratification. The halocline has been rising due to the injection system discharge into the intermediate layer causing the upper layer to become thinner as the only outflow is from the surface of the lake. An equilibrium upper layer thickness may be formed if upper layer mixing processes (wind and penetrative convection) are able to erode the halocline at a rate equal to its rate of rise due to injection system inflows (Chapter 1). As a result of these competing mixing processes the upper halocline was exceptionally sharp with a salinity change of 22.5 over less than 1 m. The average peak Brunt-Vasiala frequency was 0.87 rad s"1 for seven measurements Dec 1999 - Nov 2000, which corresponds to a period of 7.2 s. 85 The pit was constructed out of benches that are 12.2 m high and wide, with some double height benches that were 24.4 m high and 12.2 m wide. The pit walls deviate from this structure for approximately a third of the perimeter due to waste rock from the north-west dump pushed into the pit lake, waste rock placed on the south wall of the pit lake prior to closure and from slope failures and weathering. A bench coincides with the pit lake's ultimate water level (UWL). Internal waves at the halocline, which was 7.5 - 8.5 m below UWL (during December 1999 - November 2000), are therefore exposed to essentially a vertical rock wall or steeply sloping rubble surface. So strong reflection is expected and energy loss to shoaling should be minimal. If the halocline coincided with a pit wall bench then much greater shoaling and energy dissipation would be expected. Winds - The pit lake is situated adjacent to Rupert Inlet, which is a coastal fjord reaching west via Quatsino Sound to west coast of Vancouver Island and the Pacific Ocean (Figure B.l). The prevailing wind patterns along the west coast of Canada are south-east to south-west from October through April from anti-clockwise air flow around the Aleutian Low, and north-west from May through September due to clockwise flow around the North Pacific High (Thomson, 1981). The greatest departures from the prevailing winds originate from eastward moving high and low synoptic systems, which result in a roughly three day cyclic variation in wind speed. The cyclonic lows in particular intensify into storms that invade the coast between fall and late spring. These storms characteristically extend hundreds of kilometres and have durations of days, and cause strong and rapidly changing winds and heavy precipitation. The anti-clockwise airflow typically results in southeast to southwest winds that reverse in direction as the tail of the cyclonic depression passes. Coastal winds are accelerated and redirected orographically by the coast mountains and by Queen Charlotte and Johnstone Straits to be more southeasterly. In the vicinity of the Island Copper mine winds are further rotated to be more easterly (and westerly) by the Vancouver Island Range and Rupert Inlet and Holberg Inlets which, bisect north Vancouver Island in a west-east axis. Wind patterns are detailed in Appendix B. The pit lake is enclosed on the west, north and east sides by these steep pit walls, which consist of up to seven benches (~ 85 m high), while the south side is open to the flat beach dump area and Rupert Inlet. 86 Outline - The objective of this chapter is to investigate the internal wave climate in the Island Copper pit lake with particular regard to the strong salinity stratification, interaction with the buoyant plumes and the wave degeneration. A theoretical background on internal waves is introduced and the field methodology detailed. Then results from temperature, wind and injection system flow measurements are presented for winter and summer conditions. In the discussion the theoretical predictions are compared to observations. The conclusions summarize the internal wave climate at the Island Copper pit lake and the effectiveness of theoretical predictions. 4.2 Theory Internal Wave Climate - Waves are generated when a system is disturbed and a restoring force exists to bring the system back to its undisturbed state. Inertia causes the system to repeatedly overshoot the undisturbed state. The system only returns to the undisturbed state through damping of the wave energy. In stratified water bodies internal waves are generated by disturbances such as wind, and the restoring force provided by gravity acting on the density difference in the fluid. The thermal stratification of lakes and the ocean create such density differences and internal waves are commonly observed as temperature fluctuations at the thermocline. Internal waves are crucial for mixing stratified water columns as they induce shear in the interior of the water body and at lake boundaries. In large lakes the earth's rotation is important and Kelvin and Poincare waves are observed, while for smaller lakes internal gravity waves are common (Home and Goldman, 1994). Internal waves range in frequency from the inertial frequency to the buoyancy frequency (Imboden and Wiiest, 1995). The relative importance of the earth's rotation is given by the Burger Number (S) defined as: (4.1) S = fr r where g&P hh2 8 = and h'= P2 87 where g is the reduced gravity of the internal density difference Ap between the upper layer density p\ and lower layer density pi. K is the reduced depth calculated from the upper layer depth hi and the lower layer depth hj. The horizontal length scale r is taken as the square root of the halocline area. The inertia! frequency/is given by 2Qs,m<t> where Q, = 7.27 x IO"5 s"1 and 0is the latitude. The Burger number can be thought of as the ratio of the time it takes for the lake to rotate about its axis to the travel time of an internal wave across the lake, or alternately as the ratio of internal Rossby radius of deformation (RR) to the lake radius. If S « 1, then the earth's rotation dominates the wave dynamics and one has planetary or Rossby waves. If S » 1, then internal oscillations take on the characteristics of simple gravitation seiches (Imberger, 2001). When S is 0(1) both gravity and rotation are important. For example, in Lake Kuttara, a Japanese crater lake with r = 1.2 km and S = 2.9 - 3.7, Poincare type internal waves (gravity waves modified by rotation) were observed (Chikita et al. 1993). For Island Copper (winter 1999/2000 Ap = 17 kg m.3, hx = 8.45 m, Have =h{ + h2 = 141 m, r = 737 m and latitude = 50°41') the Burger number was -14, so gravitational internal waves in the form of internal seiches are expected to dominate over rotational waves. Internal gravity waves - Wind blowing across a lake piles up water in the surface layer at the downstream end of the lake. When the wind relaxes the tilted interface starts to oscillate. The solid barriers of lake boundaries reflect the internal waves giving rise to standing internal waves called seiches. Internal seiches are organized into vertical and horizontal harmonic waves with periods Tnjn associated with the horizontal (n) and vertical (m) modes. The mode describes the number locations (or nodes) of zero wave amplitude. The simplest vertical mode (m = 1) is appropriate for the pit lake as it corresponds to a two-layer system. Longitudinal oscillation periods Tn and wavelengths An, for a basin length L can be approximated by the generalized Merian formula: T, 2L (4.2) n 4 7 ^ n 2L (4.3) 'n n 88 Interna] waves are not harmonic in the strictest sense (Imboden and Wiiest, 1995) and a better approximation of first vertical mode periods is given by the Defant procedure (Mortimer, 1979). The Defant procedure was originally developed for the prediction of periods of surface seiches in basins of irregular shapes (Defant, 1960), but was adapted for internal seiches in a system of two homogeneous layers (Mortimer, 1979). The procedure describes the "basin" bathymetry using cross-sections perpendicular to the thalweg ("valley" line along the lake bottom) and specific inputs include upper and lower layer vertical cross-section areas, horizontal interface areas between cross-sections and the distance between cross-sections. The stratification is described by the upper layer depth and density difference. Seiching periods at various modes are those for which the "basin" internal waves satisfy continuity. The procedure also predicts positions of nodes, relative interface displacements and average horizontal layer velocities. Lemmin and Mortimer (1986) applied Mortimer's modified Defant procedure to eight lakes, and the Defant-predicted basin-mode periods were in good agreement (within ± 10%) of periods identified in energy spectra. Mixed Layer Behaviour - The mixed layer Richardson number (Ri) quantifies the stability of the thermocline from buoyancy to the destabilizing effect of wind and is related to the Wedderburn number (W) which determines the likelihood of upwelling (W < 1) for a two layer system (Thompson and Imberger 1980; Monismith 1985). u*L L The friction velocity u* is a velocity scale related to the wind stress T= pu* and is derived using the wind-speed measured 10 m above the water (U\o) and the ratio of air density (pa) to water density (p\). The drag coefficient (Co) varies as a function of wind speed and wave age (Donelan 1998), but is often taken as a constant value of 1.3 x 10"3 (Fischer et al. 1979). 89 w* = 1 (4.6) LD UIQ If the wind blows evenly over the water surface and manifests itself as an evenly distributed force in the surface layer, the slope of the thermocline is given by RiA (Imberger 1985). The wind must blow for a time 7V4 to tilt the interface from initially horizontal to a slope Ri'1 (Spigel and Imberger 1980; Stevens and Lawrence 1997b). For small amplitude seiches the endwall deflection of the fundamental vertical wave can be determined as follows: = L_ (4.7) ? e Ri Spigel and Imberger (1980) classified the response of a lake to wind in terms of four regimes with respect to thermocline deepening behaviour. The dependence is described in terms of a mixed layer Richardson number (Ri) and the aspect ratio of the mixed layer thickness to length. Regime 1 [Ri < 1] is for all practical purposes homogeneous. For Regime 2 [1 < Ri < (LIlhfi.iHIiH-hi))* ] deepening occurs rapidly via large interfacial displacements which cause interfacial shear and Kelvin Helmhotz billows and even upwelling. Regime 3 [(L/2hi).(H/(H-hi)f < Ri < (L2IAh2).(HI(H-hx))] develops with increased wind or with reduced g' h\. Buoyancy is still strong enough that entrainment and billowing have minor effects on the density. Internal wave amplitudes increase, but the interface still remains sharp. Most the energy for deepening comes from stirring effects introduced at the surface. In Regime 4 [Ri > (L2IAh2).(HI(H-h\))], buoyancy completely dominates all processes resulting in internal waves with short periods and small amplitudes that are hardly affected by entrainment. Billowing is not noticeable at all and the interface remains sharp. Spigel and Imberger (1980) classified 11 "lakes" by this scheme and 10 fell into Regimes 2 and 3; only the one laboratory case (Wu 1973,1977) was classified as Regime 4. Thus it is of interest to report on "natural" lakes that fall into Regime 4, as does the Island Copper Pit Lake, see below. Wave Degeneration - Spigel and Imberger (1980) proposed a time scale for wave damping (Tj) for Regime 3 and 4 lakes, based on energy dissipation in the boundary layer, 90 (4.8) where 5Ab Yd = and V u max11 471V 2 where Yd the decay modulus, V is lake volume, At, is lake bottom surface area, 8 is boundary-layer thickness, umax is the maximum velocity due to the wave [u,nax = u* Tt/(4h])], v is kinematic viscosity and e - 0.06 m is the value of the roughness suggested by Spigel and Imberger (1980) based on a smooth river channel. Spigel and Imberger (1980) found good agreement between Td and the observed rates of damping in Lake Windermere. However, Lemmin (1987) found it overestimated the damping period for Balderggersee, and suggested that circulation induced shear is more important in wave damping than boundary mixing, at least for small, deep lakes with low wind energy input. Recently, Horn et al. (2001) identified additional mechanisms for the degeneration of large-scale interfacial gravity waves, and defined regimes by comparing timescales over which each of these degeneration mechanisms acts. Td is used to define Regime I (viscous damping) found in lakes with small amplitude waves and large Wedderburn numbers. Most lakes actually fell into Regime II (solitons), where non-linear steepening transferred energy to higher frequency solitons, which lose energy by shoaling. Regime II occurs when the steepening timescale, Ts, is less than the Td and is typical of Wedderburn numbers 1-10. The steepening timescale is given by Equation (4.9), where t]0 is the maximum displacement of the interface and a is the non-linear coefficient based on the wave speed C = (g'h'f' and layer thicknesses. L (4.9) ar/0 where iJFh7^ -h{) 2h\hi 91 Higher regimes occur due to supercritical flow and Kelvin-Helmhotz billowing at Wedderburn number < 1. Each lake has its own regime diagram defined by timescale solutions, plotted on a Wedderburn number versus ratio of thermocline depth to total depth. Plume Disturbance - The discharge of ARD via two injection systems creates two plumes in the intermediate layer. During the winter, when the plumes have sufficient buoyancy flux to overcome the stratification in the intermediate layer, they will rise the full height of the intermediate layer and impinge against the upper halocline. A boil will be created where the elevated surface provides the pressure to transform the vertical flow to horizontal flow (Wood et al. 1993). The disturbance of the plume impinging against the interface may cause gravity waves to radiate away from the impinging zone. Indeed experiments of jet impinging against free surfaces cause surface waves (MacLatchy, 1998). Papantoniou and List (1989) observed intermittencies (likelihood of occupation of plume derived fluid at a position of the flow) of 90% at the centreline of an axisymmetric plume and the occurrence of "large-scale structure in the far field of buoyant jets". Fisher (1995) tracked fronts or pulses of plume fluid from the axisymmetric plume through the surface impinging zone to the horizontal surface spreading field. The frequency (fp) of such events can be estimated by dividing the plume velocity by a characteristic width, and following Fischer et al. (1979) we use the maximum time-averaged velocity (wm) and the velocity half-width (bw) to get l 1 (4.10) _w m _4.7ff 3 1 _47ff3 t p ~ bw ~ I O.lz ~ 4 z3 z3 where buoyant flux B = Qg'0 is defined by the initial buoyancy difference between the plume and ambient fluid (g'0 - g(pi-p0)lp2 where p0 is the initial plume density) and the plume discharge flow Q. The distance from the plume discharge, z, can be taken as the height of the intermediate layer. 92 4.3 Site Description and Method Meteorological Station - Pit walls on the west, north and east sides are presumed to cause wind sheltering and differential wind on the surface of the lake. Wind sheltering effects have been observed at Brenda pit lake, where pit walls extended 60-120 m high around the ~ 700 m diameter lake (Steven and Lawrence 1998, Hamblin et al. 1999). In addition, Valdivia (2001) and Laval et al. (2002) have reported differential wind on the surface of lakes caused by local topography. Indeed variable winds were observed on the surface of the lake, with greater sheltering near the pit walls and at the north-west end of the lake. As there was such variability in wind conditions across the lake the advantages of deploying the single meteorological station on a raft on the lake surface were nullified. An onshore meteorological station ensured less risk to equipment and greater data security. The meteorological station was position on the south side of the pit ~5 m above the elevation of the water surface (Figure 3.1; Appendix E). The meteorological station was equipped with an R.M Young propeller vane anemometer for measuring the wind speed and direction, Vaisala temperature and relative humidity sensors mounted in a R.M. Young radiation shield and Li-cor pyranometers for short-wave solar radiation. A Terra Sciences data logger recorded measurements every 15 min and represented point measurements except for the speed, which was a 15 min, average. The wind measured at 4.4 m was converted to a height of 10 m with a scaling factor of 1.29, assuming a neutral surface layer and applying the "law of the wall" equation following Stull (2000). A stand-alone tipping bucket rain gauge measured precipitation. A "high" wind was defined as a wind exceeding an hourly average wind speed of 7.7 m s"1, which was 1% of the record. Details of the meteorological station are given in Appendix B and Appendix E. Thermistor Chains - Two thermistor "chains" were deployed from December 6, 1999 until November 29, 2000 in five consecutive deployments. The south injection system thermistor chain (SISTC) was located over the south injection system diffuser and the central mooring thermistor chain (CMTC) was located near the West Station close to the centre of the lake (Figure 3.1). The upper halocline was intensively monitored with 7 Brancker TRIOOOs placed at half-meter intervals on each thermistor chain. The Brancker TRIOOO's had an accuracy of 0.02 °C. The CMTC TRIOOO's had internal thermistors (response time 30s) 93 and sampled every 30s, whereas the SISTC TRIOOO's had external thermistors (response time 2-3 s) and a sampling rate of 20 or 30 s was used. The depth accuracy of the thermistors at the upper halocline was -0.1 m. Flow in each of the injection systems was measured with acoustic flow meters at 30 min intervals (Appendix A). Lake water level was measured daily with an ultrasonic water level meter. These water depths were used to correct thermistor depths to the ultimate water level (UWL) datum, where the UWL is 2.44 m above mean sea level. 4.4 Results Upper thermocline mean properties - The thermocline was relatively sharp with temperature changes occurring over a depth of one metre (Figure 4.1). The upper layer water was coldest on 22 January 2000 at 4.3 °C, it warmed to 20.5 °C on 10 August 2000, before cooling again to 8.2 °C on 29 November 2000 (end of measurements). The thermocline was counter-stable during the winter when the upper layer temperature was less than the intermediate layer temperature. The density stratification was maintained by the halocline with a salinity difference of ~ 22 across < 1 m depth. The intermediate layer was stirred over its full depth by injection system flows during the winter, but at other times when the injection system stirring was weaker, the top of the intermediate layer (down to 90 m UWL) heated from solar energy. After the upper layer starts to cool on 10th August, the interface continued to warm until 3rd September 2000 (reaching 20.8 °C). Likewise, the top of the intermediate layer continued to warm until mid November (11m was 16.8 °C, 12m was 15.7 °C and 15 m was 14.0 above an intermediate layer background temperature of 12.2 °C). Warming of water in and below haloclines has been observed for other meromictic lakes (Northcote and Hall 1983; 2000) and is the principle behind solar ponds. The halocline protects the top of the intermediate layer from cooling, which occurs at the surface of the lake and is propagated downward by penetrative convection and wind mixing. A comparison of temperature-depth profiles at the CMTC and SISTC shows differences in mean thermal structure at the halocline. The upper thermocline was higher and thinner over the south injection system than at the central mooring for December 1999 to March 2000 (Figure 4.2). The thermocline over the south injection system thins to < 0.5 m 94 in February-March, whereas the thermocline at the central mooring was 0.5-1 m thick. During the remainder of the year the thermocline was the same at both locations. The differences seen in the thermocline structure can be explained by the concurrent presence of the SIS plume. In Chapter 3 it was shown that the plume rose the full height of the intermediate layer from December 1999 until April 2000, after which the injection system flow decreased and the thermal stratification at the top of the intermediate layer was able to develop. The residual summer stratification was sufficient to trap the injection system plume at depth during the fall, when injection system flows again increased. Resolving upper thermocline internal waves - Internal waves travel along the pycnocline (change in density), which is primarily determined by the halocline (change in salinity). The upper layer salinity of 4 and intermediate layer salinity of 26.6 results in a density difference (@ 12 °C) of 17.5 kg m" . The seasonal temperature change of the upper layer has a relatively minor affect on density, as upper layer temperatures of ~ 4 °C in the winter and ~ 20 °C in the summer cause the density difference to change by - 0.6 kg m"3 (-3%) and + 1.4 kg m" (+8%), respectively. Thermistors measuring temperatures were used as a surrogate for the more important salinity because measurement of temperature is simpler and more reliable than salinity (via conductivity). Temperature was a good indicator of the halocline from December to March as both the upper and intermediate layers were well mixed during this time and the thermocline and halocline coincide. However, using temperature fluctuations as an indicator of internal wave activity was ineffective in May and November when the upper layer warms and cools past the intermediate layer temperature and the temperature gradient disappears. In the late spring, summer and fall, the thermocline broadened downward as the top of the intermediate layer warmed and does not necessarily coincide with the halocline, thus care was required interpreting the thermistor data. To resolve internal wave activity in the halocline, thermistors were located at 0.5 m intervals over 2-3 m. Despite these efforts there were times when no thermistors were located in the thermocline, especially at the thinner south injection system thermocline. During the December to March period, the only thermistors found completely in the thermocline were at SISTC 8.39 m(25 December 1999 - 18 January 2000), CMTC 8.43 m (1 January - 15 April) and CMTC 8.93 m (6 December 1999 - 15 February 2000) (Figure 4.3). 95 The thermocline rose past these instruments due to the injection system discharge to the intermediate layer, which raises the halocline. A concurrent but smaller rise in lake water level (and lifting of the thermistor chain) occurs, but was not sufficient to keep the thermistors in the thermocline, as water drains from the surface of the lake. Thermistors, outside the thermocline could not resolve disturbances caused by internal waves because the temperature gradients within the layers were very slight and the resulting temperature fluctuations were small. Thus, detailed internal wave analysis focuses on the period of 25 December 1999 to 18 January 2000 by default, as this was the only period where thermistors were found concurrently in the thermocline at both the SISTC and CMTC. Fortunately, this period coincides with the most and least active windy and rainy periods during the entire experiment. Standard practice is to convert temperature fluctuations at a thermistor to internal wave displacements using the temperature gradient established from adjacent thermistors (Emery, 1998). Three instruments in an area of relatively consistent temperature gradient are required. In our case only one or at the most two instruments were in the sharply defined thermocline, so the temperature gradients and displacements could not be adequately resolved. To use thermistors from outside the thermocline would underestimate the temperature gradient and overestimate the displacement. Spectral analysis using displacement better reflects internal wave energy than using temperature fluctuations, which depend not only on the wave amplitude but also the temperature gradient. However, because of the inaccuracies in estimating displacements we were forced to use temperature fluctuations for spectral analysis. As a result the comparison of spectral density magnitudes is compromised. Nevertheless, the spectra of temperature time series do correctly diagnose frequencies (Lemmin, 1979). Upper thermocline internal waves during winter - The winter period of 15 December 1999 - 18 January 2000 was examined in detail as thermistors were located in the thermocline and periods of both stormy and calm weather were encountered. Figure 4.4 shows temperatures, wind and injection system activity during this period, although only hourly temperatures and wind are shown for reasons of presentation. Starting the period was typical December weather with several cyclonic low-pressure systems that produced high 96 winds events on the 12, 15 and 17 December and two SIS flow peaks of 338 L s"1 on 12 December and 258 L s"1 on 16 December. An extended period of calm sunny weather followed, from 20-29 December, caused by an atypical high-pressure system. During this period of calm weather the mean hourly wind speed was 0.94 m s"1, peak hourly wind speed was 2.59 m s"1 and there was a total of 0.5 mm of rain. The weather returned to the pattern of successive cyclonic low-pressure systems in the new year, producing six separate high wind events from 3-8 January and two SIS moderate flow peaks of 160 L s"1 on 4 January and 210 Ls"1 on the 9 January. The cyclonic systems in both December and January were distinguished by the oscillation of wind direction between westerly (~ 90°) and easterly (~ 270°) as the fronts pass. The weather improved during 10-18 January, with low-moderate winds and decreasing flow at the SIS, while flow at the NIS, which was less peaky, held steady. Temperatures in the thermocline exhibit considerable high frequency variation and this activity increases during periods of high wind and injection system flow (Figure 4.4). Greater temperature fluctuations were observed at the SISTC due to the sharper thermocline and perhaps its location away from the centre of the lake. Horizontal mode-1 seiches would produce smaller amplitude fluctuations in the centre of the lake because of the central location of the wave node. At times the observed temperature fluctuations were limited by contact with upper or intermediate layer water, indicating that internal waves displaced upper layer fluid down to the level of the thermistor or the intermediate layer fluid up to the thermister. The spectra of thermistors SISTC 8.39 m and CMTC 8.43 m were compared for the period of 1- 18 January 2000 where they coexisted in the thermocline (Figure 4.5). The greatest spectral density occurs at long periods due to seasonal and frontal weather patterns and changes to the interface position (and the position of thermistors in the interface due to water level changes). At the halocline over the south injection system, the spectral density was 10 - 100 times greater than that measured at the central mooring. The greater spectral density at the SIS may result from larger displacements due to the location of measurement relative to the mode-1 seiche. The sharper temperature gradient at the SIS thermocline would also result in a larger temperature fluctuation for a given displacement. 97 Numerous large peaks were observed in the spectra (Figure 4.5) and these correspond to vertical mode-1 seiches of various horizontal modes. The lowest frequency wave occurred at 53.4 min, which matches the Defant procedure prediction of 52.6 min and was similar to the Merian formula prediction of 64.2 min for horizontal mode one lake-length internal waves. Spectral peaks at higher frequencies correspond well to Defant procedure predictions of higher horizontal modes (Table 4.1). However, not all spectral peaks were accounted for by the various horizontal modes of the lake-length internal waves. The outstanding spectral peaks were found to correspond well with predicted lake-width internal waves. The lowest frequency spectral peak not attributed to a lake-length horizontal mode has a frequency of 24.8 min, which matched the lake-width horizontal mode-1 internal waves predicted by the Defant procedure of 24.4 min and was similar to the Merian formula prediction of 29.4 min. Wind measured at the site was predominantly from the east and west quadrants and was largely absent from north and south quadrants, consistent with the west-east alignment of Rupert Inlet and the Vancouver Island range. The large average winds were from the northeast, southeast, southwest and west (Figure 4.6). Therefore, winds were available with cross lake components and sufficient magnitude to cause lake-width seiches. The spectrogram of thermistor SISTC 8.39 m and CMTC 8.43 m reveals the frequency versus time dependency of the wave climate (Figure 4.7). In response to the initial input of energy from wind events 12-17 December at the CM, internal waves L l (lake-length mode 1), L2 (lake-length mode 2) and WI (lake-width mode 1) were generated (although the response at the SIS was obscured because the thermistor was not yet fully in the thermocline). These waves become damped after the wind stops, but were still visible in the spectrogram. L l actually appears stronger during 28-30 December because the thermistor was closer to the middle of the thermocline (Fig 4.6), as there was no new input of wind energy. The wind events during 3-10 January 2000 input new energy and cause an energetic response across the spectrum but particularly at internal wave frequencies. After the 10 January 2000 energy was lost, except at internal wave frequencies, which persist and gain energy from smaller wind events. Two half days of halocline water temperatures, wind and south injection system inflows are detailed in Figure 4.8. The 29th December 1999 was a calm day at the end of the 98 10-day period of good weather. Correspondingly, the halocline was relatively inactive, except for at SISTC 8.35 m where the remnant L l internal wave dominates. In contrast, the 8th January 2000 was during stormy weather with interna] wave generation in progress. The peak average hourly wind was 7.0 m s"1 and the south injection system flow was 13,734 m3d"1. During this period the halocline was very energetic, as evident from the temperature fluctuation at a range of frequencies, particularly at the SIS. Internal waves L l at SISTC 8.35 m and L2 at CMTC 8.35 m were obvious. Upper thermocline internal waves during summer - The summer period of 10 June - 12 July 2000 is presented as a comparison to winter conditions (Figure 4.9). During this time there was minimal injection system flow into the pit lake. Winds were consistently diurnal in nature and five high wind events were recorded. The winds were sea-land breezes as they peak each day in the late afternoon and drop to near zero at night. The direction of these winds was variable (Figure 4.10) perhaps due to the concurrent existence of synoptic weather patterns. The diurnal wind forcing shows up as daily temperature fluctuations at the halocline thermistors (Figure 4.9). The apparent magnitude of these fluctuations seems small, but the temperature gradient through the halocline was small (Figure 4.1), thus the internal wave signal was obscured. The diurnal forcing was observed in the temperature spectra at 1.2 x 10"5 s"1 (Figure 4.11). Lake-length and lake-width internal waves were evident (Figure 4.11), but their periods were slightly longer than in winter due to the reduced summer g' and h' (see Equation 4.1). The magnitude of the spectral densities at both the SISTC and CMTC were similar, which indicates that the differences in magnitude of the spectral densities observed during the winter were probably a result of the greater temperature gradient (from winter plume erosion) at the SIS, and not its location away from the centre of the lake. 4.5 Discussion Seiche Periods - Lake-length and lake-width internal wave periods were well predicted by Defant's procedure (Table 4.1). For example the observed winter horizontal mode-one lake-length internal wave (Ti) 53.4 min was in good agreement with Defant prediction of 52.6 min, but less well matched with the Merian formula estimate of 64.2 min. 99 Higher modes of both lake-width and lake-width internal waves were also reasonably predicted. As the wind period to cause tilting of the interface was only -15 min (7V4), short temporal changes to wind speed and/or direction will disturb the interface. Ri, W and ge for winter 1999/2000 conditions are given in Table 4.2. For the maximum observed wind of 13.7 m s"1 and high wind of 7.7 m s"1, the internal wave amplitudes are predicted to be 0.238m and 0.075 m, respectively. So disturbances while being frequent are always small. Mixed Layer Behaviour - In terms of Spigel and Imberger's (1980) mixed layer classification, Regime 3 occurs at Ri > (LI2h\).(HI(H-h{))Vl = 134 and Regime 4 at Ri > (l}lAhi2).{HI(H-h\)) = 18,000. For the year of observation the pit lake was classified as Regime 4 for 98.3% of the time, and as Regime 3 only when Ri decreases in response to strong winds. The pit lake was consistent with Spigel and Imberger's (1980) Regime 4 descriptions as buoyancy dominates resulting in small amplitude, short period internal waves and the interface remains sharp. Wave Persistence - With regard to the wave degeneration the time-scales Td and Ts can be calculated for the maximum and high wind events using winter 1999/2000 conditions. Td is expected to be long as the lake seiche amplitude is small and the lake is deep; Yd is small from a thin boundary layer 5 due to small umax and because AtJV « 1. Td was calculated to be 24 days for the maximum wind generated wave and increased to 76 days for the "high" wind generated wave. Whereas, Ts was only -13 hrs for the maximum wind generated wave and increased to 43 hrs for the high wind generated wave. Thus, Horn et al.'s (2001) classification suggests that solitons are produced and dominate wave degeneration, which occurs within several days. This is an unexpected prediction given the high Wedderburn number and low amplitude waves, which are expected to be well described by linear approximations and thus unlikely to steepen due to non-linearities. Ts is small because pit lake conditions result in an unusually large a. The wave-speed c0 = (gh{)A is an order of magnitude higher than thermally stratified lakes due to the strong salinity contribution to density difference, and (hi - h\) is large compared the case h\ - hi where the non-linear coefficient vanishes (see Equation 4.9). The observed degeneration of internal waves was not well predicted by Horn et al.'s (2001) theory as the soli ton time-scale was too short and the damping time-scale too long. 102 Waves generated by "high" winds were predicted to degenerate by the development of solitons in less than two days (following Horn et al. (2001)), but observations showed the persistence of waves for greater than 10 days (>270 wave periods). There was no evidence of the development of solitons, rather wave amplitudes were slowly reduced with time suggesting wave degeneration by wave damping. Spigel and Imberger's (1980) parameterization for wave dampening is based on boundary mixing over the lake bottom, which seems unrealistic for the pit lake as the predicted decay periods are too long. Additional work needs to be done to determine the role of mixing at the side walls and from damping due to intermediate layer turbulence. High frequency internal waves on the upper halocline were predicted, but were not able to be observed and this possibility remains an outstanding issue to be investigated. 103 Table 4.1. Observed internal waves compared to Defant and Merian predictions*1. Observed Lake-length Internal Waves Lake-width Internal Waves (min) mode Defant Merian mode Defant Merian (min) (min) (min) (min) 53.4 1 52.6 64.2 30.7 2 29.0 32.1 24.8 1 24.4 19.4 19.4 3 20.0 21.4 15.4 4 15.3 16 14.5 2 13.4 14.6 13.4 5 12.7 12.8 11.8 3 9.5 9.8 10.6 6 10.6 10.7 9.5 4 7.6 7.3 *' Defant and Merian calculations are based on winter 1999/2000 conditions at the upper halocline which were Ap = 17 kg m3 and hi = 8.45 m. The pit lake has Have =h\ + h2 = 141 m. Table 4.2. Parameterization of Island Copper pit lake to wind disturbance. Parameter Average Wind Typical Peak Wind Maximum Wind U 1 0 (m s"1) 2.53 7.70 13.7 u* (m s"1) 0.0032 0.0097 0.0173 Ri 135,000 14,700 4,630 W 520 56 18 Se.(m) 0.008 0.075 0.238 5(m) - 0.058 0.187 Yd - 4.76e-4 1.52e-3 T d - 76 days 24 days a - 0.19 0.19 T s - 43.0 hrs 13.5 hrs 105 Dec Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec Date Figure 4.1. Daily average temperature contours of the upper halocline at the central mooring and from thermistors at 0.5 m vertical separation. The upper halocline [ ] and its shoulders [—] are based on 7 C T D casts undertaken during the period. Figure 4.2. (a) Profiles of daily average temperature through the upper halocline at bi-monthly intervals at the central mooring [—] and south injection system [—] from thermistors at 0.5 m vertical separation, (b) South injection system [—] and north injection system [—] flows. 107 (a) Dec Jan Feb Mar Apr May Jun Jul A u g Sep Oct Nov Dec Dec Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec Figure 4.3. Daily average thermistor temperatures at the (a) central mooring in the upper layer at 2.93 m [—], in the thermocline at 6.93 m [ ], 7.43 m [—], 7.93 m [—], 8.43 m [—], 8.93 m [—] and in the intermediate layer at 19.93 m [ ]. (b) South injection system at the corresponding depths which are 0.04 m lower. Figure 4.5. Spectrum of CMTC 8.43 m [—] and SISTC 8.39 m [—] thermistors versus (a) frequency and (b) period 0-60 min for period 1-18 January 2000. 110 607000 608000 609000 Figure 4.6. Wind for 1-18 January 2000 and Island Copper pit lake. The windrose is centred at the location of the meteorological station and is constructed from 15 minute average winds binned in 36 directions plotted as percentage wind direction and average wind speed in each direction. The pit lake is sheltered by steep pit walls consisting of 12 m high benches as indicated by lines around the west-north-east sides of the pit. 112 29-Dec-1999 8-Jan-2000 12 F l O OJ 1 8 CD C L I 6 h-4 12 1 8 CD E L I 6 i — 4 60 CO CN 40 CN -o 20 CZ 5 (a) CM i i i i i i i i i t i i T i i i i i i i i i i i i i i i Jbl_SIS a i i i i i ( i i i i I i i i i i i i i i i i i (c) Win j <so1id> and SISJ Flow <d|ash> 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 9.85m 9.35m 8.35m 8.85m 8.35m 7.85m (d) CIV 1 1 1 i i > >— 1 1 . f t ) sis ! i 1 ! y i l • i (f) Wir cj <so!id> and SIS Flow <dc sh> ; - H — ' \ 8.35m 8.85m 1.35m '.85m 00:00 02:00 04:00 06:00 08:00 10:00 12:00 00:00 02:00 04:00 06:00- 08:00 10:00 12:00 Time (hours) Time (hours) Figure 4.8. (a) CMTC and (b) SISTC temperature and (c) wind and SIS flow for 12 hours on 29 December 1999, and (d) CMTC and (e) SISTC temperature and (f) wind and SIS flow for 12 hours on 8 January 2000. 114 8OT0D SHED earn Figure 4.10. Wind for 10 June - 12 July 2000 and Island Copper pit lake. The windrose is the centred at the location of the meteorological station and is constructed from 15 minute average winds binned in 36 directions plotted as percentage wind direction and average wind speed in each direction. The pit lake is sheltered by steep pit walls consisting of 12 m high benches as indicated by lines around the west-north-east sides of the pit. 115 (a) 10 10 CO c SO o 0 Q. CO 10 10 10 r m ! — n Frequency (s ) (b) 30 Period (min) Figure 4.11. Spectrum of [—] and SISTC 8.39 m [ ] thermistors versus (a) frequency and (b) period 0-60 min for period 10 June - 12 July 2000. 116 (a) 01/09 01/10 01/11 01/12 Date (mm/dd) Figure 4.12. Upper thermocline temperatures at the (a) central mooring (sampled at 30 s) and (b) south injections system (sampled at 20 s), compared to (c) 15 min wind energy [—] and direction [- ] for the 9-12 January 2000. 117 CHAPTER 5 Modeling the Long-term Evolution Abstract The meromictic pit lake at the Island Copper Mine is being used for the treatment of acid rock drainage (ARD). Maintaining a permanent stratification is crucial to the performance of the passive treatment system because the ARD is expected to be generated for over 100 years. Two conceptual models were used to evaluate the long-term evolution of the pit lake and potential failure modes. Model 1, a two layer continuously stirred reactor box model, predicts that upwelling is the likely failure mode because the density difference between the upper and intermediate layer decreases with time, which decreases the wind speed required for upwelling. Extreme wind frequency analysis was used to estimate the likelihood of failure. In Model 2 the upper layer deepens in response to the decreasing density difference between the two layers. This model predicts that upwelling is not likely and that erosion of the intermediate layer is the expected failure mode. Both models demonstrate that maintaining the density difference between the upper and intermediate layers extends the life of the stratification. This extension can be achieved by diverting the north injection system discharge from the intermediate layer to the upper layer. The one-dimensional lake model DYRESM was also used. It was able to resolve the short-term and annual trends such as mixed layer behaviour and plume dynamics. However, the model does not yet simulate the mixed layer well enough to balance the annual heat budget of the upper layer. For this reason multi year modeling was not possible. DYRESM does not include entrainment of the upper halocline by the impinging plume. Plume entrainment of the halocline provides an additional source of dilution to the intermediate layer that was estimated at 30% of the injection system flow (Chapter 1; Appendix A). Plume entrainment 118 of the halocline is anticipated to increase with time as the density difference between upper and intermediate layers decreases. Therefore, plume entrainment needs to be added to the DYRESM before long-term simulations are attempted. 5.1 Introduction Island Copper Pit Lake - The Island Copper pit lake was flooded with seawater and capped with fresh water to form a meromictic (permanently stratified) lake for the treatment of acid rock drainage (ARD) (Chapter 1). The proposed passive treatment system utilizes metal-sulphide precipitation via anaerobic sulphate reducing bacteria (BHP 1996). Maintaining the meromictic structure of the pit lake is crucial as the permanent stratification prevents the transport of oxygen into the lower layers and the ARD injection system plume from reaching the surface of the lake. Estimates of how long ARD will be generated vary widely. Li (1991) indicated that the Island Copper north dump had the potential to generate ARD for 650 years. Lister (1994) predicted that on land waste dumps at Island Copper will generate ARD for 99 to 551 years. Based on Morin and Hutt's (1996) sulphate measurements the beach dump waste rock may produce ARD for up to 120 years (pers. com. G. Poling). Dagenais (1996) also found the north dump material had a high potential to generate acid. The treatment for ARD must be long-term. Thus it is important to predict the long-term structure and stability of the stratification as it is crucial to the performance of the passive treatment system. The stratification of the pit lake is expected to weaken with time due to the discharge of ARD (essentially fresh) to the intermediate layer, which will decrease the salinity and density difference across the upper haloline. Destratification of the lake may eventually occur due to turnover (mixing of the lake) initiated by upwelling. Upwelling occurs when the wind is sufficiently strong to tilt the haloline, such that the halocline and intermediate layer water come to the surface at the upwind end of the lake. Alternately, as the stratification weakens, the upper layer may deepen until it reaches the lower haloline, and/or the upper and intermediate layers might turnover (if density increases to the upper layer from fall/winter cooling overcome the stabilizing density difference due to salinity). An additional complication is the north injection system (NIS) and south injection system (SIS) plumes, 119 which impinge against the upper halocline and entrain fluid into the intermediate layer, because as the density difference decreases the plume entrainment will increase (Chapter 3; Appendix A). Lake Modeling - Modeling of meromictic lakes or pit lakes or lakes with plumes is briefly reviewed. The simplest type of model is a continuously stirred reactor (CSR) box model where the entire lake is considered a well-mixed volume, or as two boxes for the epilimnion and hypolimnion (Imboden and Schwarzenbach 1985; Chapra 1997). The two box approach can be enhanced to model mixed layer deepening by environmental fluxes and is called mixed-layer modeling (Ward et al. 1990; Stevens et al. 1994). The set of box models includes filling-box models, which can be applied to turbulent buoyant discharges (i.e. plume discharges) into a confined region (Baines and Turner 1969; Kumagai 1984; Chapter 3). More complicated is a one-dimensional diffusion model whereby the transport processes are modeled analogously to diffusion. Determining appropriate turbulent diffusion rates is the difficulty with this type of model. Diffusivity is determined from semi-empirical equations (Imboden and Wiiest 1995), in situ measurements and calibration to observations. Diffusion type models have been applied to pit lakes (Stevens and Lawrence 1995 and 1997a; Bohrer et al. 1998) and other meromictic lakes (Sanderson et al. 1968). More sophisticated are one-dimensional physical process models, which quantify all the important physical processes and their interactions. DYRESM is a widely used model of this type and is detailed in Fischer et al. (1979) and Imberger and Patterson (1981, 1990). DYRESM has been applied to meromictic lakes (Romero and Melack 1996; Jellison et al. 1998), pit lakes (Hamblin et al. 1999) and to model bubble plumes (Schladow and Fisher 1995; Lindenschmidt and Hamblin). Two and three-dimensional models are available but are less well suited for long-term modeling due to computing restrictions. The use of higher dimension models is premature, since the one-dimension assumption in good in small lakes, particularly if the lake is strongly stratified. Island Copper Modeling - Considerable effort has been expended modeling the Island Copper pit lake. The majority of modeling was done prior to the flooding of the lake, as part of the design for the passive treatment system. Initially two and three layer box 120 models with mixed layers were used (Allen 1991, 1992 and 1993), but a diffusion type model was chosen for the final design (Dunbar and Hodgins 1990; Rescan 1994 and Dunbar 1995). Wilton (1998) found that even after two years there were differences between the predictions of the final design model and observed conditions. The final design model was deficient in the following areas: heat flux into the lower layer; unstable salinity and potential density profiles for the intermediate layer; and no mixing across the upper interface resulting from the impinging buoyant plume. Also not captured in the final design model was the thinning of the upper layer due to the injection system discharge into the intermediate layer (Fisher et al. 1999, Fisher and Lawrence 2000). Wilton (1998) used a two layer box model and predicted elevated metal concentrations in the long-term for the worst case scenario of no biochemical treatment. Rescan (1999a) modeled the lake out of concern that oxygen concentrations were not decreasing as predicted by the final design model. The new model used in situ measurements of biological oxygen demand and predicted that intermediate water oxygen concentrations would cycle seasonally between summer concentrations of zero and low-levels of < 0.2 mg L"1 in winter. Outline - Modeling should start at the simplest level and build with complexity once the simple model has been exploited (Imboden and Schwarzenbach 1985). The purpose of modeling the pit lake is to understand the its long-term evolution and stability, rather than the short-term details, so simple conceptual models are appropriate. Three models are used in the present study. Model 1 is a simple CSR model that represented the upper and intermediate layers of the pit lake as two fixed volume boxes. Model 2 is an extension of Model 1, which allows for deepening of the upper layer with time in response to reduced density difference across the upper interface. The third model is DYRESM, which replicates most aspects of the lake physics. These models are detailed and the results of the simulations evaluated. Model predictions of the wind speed to cause upwelling are coupled with frequency analysis to predict the risk of destratification by upwelling. Finally the limitations of each model are explored and recommendations for future modeling efforts given. 5.2 Lake Models Model 1 - A CSR model is the simplest model for a natural water body. It is appropriate for lakes that are fully mixed (Chapra 1997) and can be extended to lakes with 121 distinct layers. These criteria are met for the upper and intermediate layers, as each layer is discrete and well mixed (Chapter 1). Wilton (1998) linked two CSRs to model the upper and intermediate layers as a combined system. The lower layer is less critical to the long-term stability of the pit lake and was excluded from the model. Long-term changes to salinity and heavy metals were modeled. Wilton's (1998) model assumes that the mixing time is negligible compared to retention time and salinity is conservative. Furthermore, the upper layer depth is assumed to be constant, which means upper and lower layer volumes are constant and the inflow into each layer equals outflow. A simple mass balance gives the change in salinity with time (dSldt) from the mass of salt in {qSin) less the mass of salt out (qS), VdS = {qSin-qS)dt (5.1) where q is the inflow ( = outflow) and V is the layer volume. Integrating the above expression gives S i„ 1-A ^ " M l J T, e R where the original salinity is So and the retention time is TR = Vlq. Note, that at large time (t » TR) that S = Sin. Model 2 - The CSR two layer box model was modified to allow for changing upper layer depth. Model 1 was limited by the assumption of a constant upper layer depth. In reality the halocline may deepen because upper layer mixing processes will become more effective as the stratification weakens. The upper layer depth is critical to the calculation for upper layer salinity and upwelling wind speed and duration. Therefore a more dynamic model that accounts for the change in upper layer depth will be a better analogue of the pit lake and allow for better prediction of the failure mode. Note, both Model 1 and 2 assume the lake water level remains fixed, whereas the observed (1999-2001) seasonal fluctuation in water level was ~ 1 m, but the assumption is reasonable for inter-annual modeling. 122 An equilibrium upper layer depth is expected when upper layer mixing processes erode the halocline at a rate equal to its rate of rise from injection system inflow and entrainment of the halocline by the impinging plume (Chapter 1; Appendix A). The depth can be found from the energy balance, where the change in the potential energy of stratification is equal to the amount of mixing energy that has been utilized. The work done (mechanical energy ME per unit time) changes the potential energy by entraining fluid and mixing it through the upper layer (per unit area); lifting entrained fluid of density difference (Pi-P\) by half the height of the upper layer (h\) (P2 -Pl)gh.(^L+ QA +QE^_d(ME) (5-3) 2 v dt A Adt where p\ is the density of the upper layer, pi the density of the intermediate layer, g is gravity and A is the area of the upper halocline. When the upper layer is in equilibrium (dh\ldt = 0) the flux of fluid entrained into the upper layer equals the ARD inflow (QA) and plume entrainment (QE) of the halocline. Equation (5.3) predicts for a given QA+QE and d(ME)/dt, that the upper layer will deepen (dh\ldt > 0) as the density difference (pz-pi) decreases. Note that QE equals 2Ck(Pi-p§)l(Pi-p\) QA, where Q is the fraction of the plume power available and po is the density of the injection system discharge (Appendix A). The validity of Equation (5.3) is examined by a brief review of the literature on mixed layer deepening. Important to the energy balance is the type of mixing and how the erosion rate changes as the lake evolves. Spigel and Imberger's (1980) classification places the pit lake's mixed layer in regime 4 where buoyancy dominates all processes (Chapter 4), so internal waves have short periods and small amplitudes, and interfacial shear is not important for mixed layer deepening. Rather, mixing is by surface stirring effects such as wind and penetrative convection and these are expected to deepen the mixed layer very slowly. The deepening of a mixed layer by stirring has been widely investigated by grid experiments and the rate of non-dimensional deepening ranges from Ri'1 to RF115 (Nokes 1988; Fernando 1991). In the context of mixed-layer deepening in lakes and if interfacial shear can be neglected, the model of Kraus and Turner (1965), experiments of Wu (1973) 123 and Kranenburg (1984) and review of Imberger and Patterson (1990) suggest a relationship of the form, dhi _ CfU (5.4) dt Ri where Ri = g% /' u2 with g'= g{jh-P\)lpi, CF is the fraction of the mean kinetic energy flux, expressed as u the root mean (over the depth of the mixed layer) square velocity that arrives at the interface. CF may be itself a function of both upper layer depth and of the shear velocity available at the interface, but has been modeled as a constant with good results (Spigel et al. 1986). Rearranging Equation (5.4) gives {p2-p\)ghi^L = CFp2u3 y J ' J ) dt which describes how reduced h\ or pi-p\ causes ah increase in rate of deepening for a constant amount of mixing. Note that the mixing energy CFP ii3 = d(ME)/Adt, so that Equation (5.5) has the same form as Equation (5.3). Similarly for the mixed layer of a lake, the change in potential energy of stratification can be expressed at the sum of penetrative convection and wind induced mixing (Fischer et al. 1979, Ward et al. 1990), / N , , an, _ 3 „ 3 K?-V) \P2 ~ Pi )Sh = P\.CkU f + PyC*U* where Uf is the velocity of a falling penetrative convection plume and u* is the wind induced shear velocity. Ck and C* are the penetrative convection and wind efficiency factors, respectively. Again the form of Equation (5.6) is the same as Equation (5.3), except for the inclusion of injection system inflow, as the right-hand-side of Equation (5.6) is d(ME)/Adt. Ward et al. (1990) solved Equation (5.6) for Q and C* from wind and temperature data and observed changes in the stratification. DYRESM's mixed layer algorithm essentially solves the same energy balance (Antenucci and Imerito 2001). Model 2 solves Equation (5.3) by using an estimate of mechanical energy, rather than by detailed modeling of the mixed layer processes such as wind and penetrative convection. 124 The available mixing energy was estimated from field observations of the change in halocline position for known injection system flow. In this simplistic way the calculation of mixing energy is already calibrated. A limitation of the model is that the mixing energy of the one year is assumed to be representative of all future years. Model 2 solves Equation (5.3) at annual time-steps while simultaneously accounting for the changes to salinity via Equation (5.1) in a similar manner to Model 1. DYRESM - The conceptual models are useful and improve our understanding of the lake physics, but their predictive ability is limited by their simplicity. A more sophisticated model like DYRESM accounts for a multitude of environmental fluxes and lake physical processes Fischer et al. (1979), Imberger and Patterson (1981, 1990). DYRESM is a one-dimensional hydrodynamics model for predicting the vertical distribution of temperature, salinity and density in lakes and reservoirs. The model is constructed of horizontal layers of uniform properties. Importantly, DYRESM has been applied to many lakes and Hamilton and Schladow (1997) claim that the hydrodynamic sub-model performs so well that calibration is not required. Algorithms for submerged buoyant discharges have recently been added to the model (Antenucci and Imerito 2001). Each submerged discharge is treated as an axisymmetric buoyant jet (Fischer et al. 1979), which entrains fluid from each layer until it reaches a level of neutral buoyancy when it detrains into that layer. Before attempting to make long-term predictions with DYRESM, it is first necessary to validate the model. Validation is performed for the period of December 1999 to November 2000 against water column thermistor and CTD measurements (Appendix C-E). The full water column (all three layers) was modeled. DYRESM does not model double diffusion at the upper and lower interfaces; entrainment of the upper halocline by the impinging plume, enhanced mixing in the intermediate layer by plume advective currents and geothermal heating in the bottom layer. 5.3 Upwelling Model Upwelling - Destratification of the lake may eventually occur due to turnover initiated by upwelling. Upwelling occurs when the wind is sufficiently strong to tilt the interface, such that the halocline and intermediate layer water come to the surface at the upwind end of 125 the lake. The Wedderburn number indicates the potential for upwelling (Thompson and Imberger 1980; Monismith 1985), W=$Ar where u2 =Cd^-U2 Lu2 * Pi * where L is the fetch, u* is the shear velocity, Cd is a drag coefficient (Fischer et al. 1979 recommended Q = 1.3 x 10~3), pa is the density of air and U the wind speed at a height of 10 m. Theoretically, upwelling of the halocline will occur when W < 1 and partial upwelling may occur when W ~ 3 (Monismith 1986). The wind speed required to cause upwelling (Uu) can be determined by rearranging Equation (5.7) and setting the Wedderburn Number equal to one. u = Uh A U \ L paCd The duration (Tu) this wind speed needs to be sustained to cause upwelling has been found to be T\/4 (Spigel and Imberger 1980, Stevens and Lawrence 1997b), where T\ is the fundamental wave period defined by the Merian formulae as (g'\h2/{\+h2))U2 where I12 is the lower layer thickness. If the density difference decreases with time then Uu will decrease and Tv will increase. From a probabilistic perspective, a lesser wind is more likely, but as the required duration increases the wind becomes less likely. A fetch length of 2350 m corresponding to the lake-length (ESE) axis (Figure 1.2) was used to calculate the wind speed for upwelling. To utilize the frequency analysis (described in the next section) it is necessary to convert the annual upwelling wind speed to an hourly annual upwelling wind speed. Resio et al. (2002) gave the following empirical relationship to scale extreme wind speeds from an hourly average period to another averaging interval (in seconds), taken in our case as Tu. 126 U U _ u = 1.277 + 0.296 tank U,hr 0.9log1Q for t < 3600 (5.10) Risk of Upwelling - The likelihood of upwelling in the pit lake in any year is the probability of obtaining the annual upwelling wind speed in that year. A frequency distribution for maximum annual hourly wind speed was determined from the Port Hardy Airport meteorological station, which is 12 km away and has 48 years of record (Appendix B). Wind observations at Island Copper (December 1999 to November 2000) were reasonably well correlated to Port Hardy, particularly the east-south-east (ESE) wind that dominates at both sites (Appendix B). The wind scaling factor of 0.825 was determined by the ratio of wind speed at Island Copper to that at Port Hardy, for the five maximum Island Copper wind events occurring in the year of simultaneous observations (Appendix B). An Island Copper frequency distribution was calculated by multiplying the Port Hardy distribution (Appendix B, Equation B.3) by the wind scaling factor. This approach was reasonable because all five maximum Island Copper wind events were all ESE winds, as were the 48 winds used to construct the Port Hardy frequency distribution. The ESE wind blows straight down the long axis of the pit lake (Figure 1.2). The resulting Island Copper maximum annual hourly wind (Uic) frequency distribution is shown in Figure 5.1 and described by the equation ( ( TR =21.52 23-U IC 0.825 1.346+-U IC 0.825x300 (5.H) J where TR is the return period 1/P and P is the probability that the wind speed required for upwelling is exceeded in any given year. The Island Copper maximum annual wind frequency distribution predicts a 100-year return period hourly wind of 18.7 m s"1. Note, the maximum hourly wind observed at Island Copper for December 1999 - November 2000 was 10.5 m s"1. Wind speeds greater than 18.7 m s"1 are assigned a zero probability as they can not be reliably estimated and because the wind frequency distribution trends asymptotically upwards suggesting that higher winds are not possible. The return period of the hourly 127 annual upwelling wind speed occurring was calculated for each year. The likelihood of upwelling is best expressed by the risk, R = l-(l-P])x(l-P2)x---x(l-Pn), which is the probability that the wind speed required for upwelling is exceeded at least once in n successive years (Viessman and Lewis 1996). 5.4 Results Model 1 - The model used monthly time steps. Initial salinities were 4 and 27 for the upper and intermediate layers, respectively. The upper layer depth was 7 m. Upper and intemediate layer volumes were based on detailed bathymetry (Rescan, 1999a). Injection system inflows had a salinity of 1.1 and flow rates were based on the monthly averages from 3 years of monitoring. Rainfall inflows to the upper layer were based on hydrological estimates (BHP, 1996). Inflows and outflows were determined by conservation of volume: the upper layer entrains intermediate layer fluid equal to the injection system inflow (plume entrainment of the upper halocline was ignored in this model); and the upper layer outflow equals the combined inflow from rain and entrainment from the intermediate layer. Layer densities were calculated using Equation of State 1980 (UNESCO 1981) from the model salinity and intermediate layer temperatures of 12 °C and seasonally varying upper layer temperature (4-20 °C). Run 1 of the model was for both north and south injection systems discharging to the intermediate layer (Figure 5.2). The upper layer salinity was predicted to increase in the short-term due to mixing with the more saline intermediate layer water, then decrease to the proportional average of the rainwater and intermediate inflow salinity. The model predicted that the intermediate layer salinity asymptotically decreased with time, due to the salinity of the injection system. The difference in salinity (and density) between the upper and intermediate layers decreased with time, and thus the stability of the stratification also decreased. Run 2 of the model was for NIS diverted to the upper layer (Figure 5.3). Under these conditions the difference in salinity (and density) between the upper and intermediate layers decreased more slowly. The models predicted, for both Run 1 and 2, that the wind speed required for upwelling decreased with time as the salinity difference (and hence g') decreased. Seasonal 128 temperature change in the upper layer caused annual variation in the upwelling wind speed (Figure 5.2 and 5.3). An annual upwelling wind speed was defined as the average of November to March values each year, as all maximum wind events occurred in these months (Appendix B). The benefits of diverting the NIS to the surface layer for the stability of the lake are clear, as the annual wind speeds to cause upwelling at 100 years change from ~ 16 m s"1 to ~ 23 m s"1, for Run 1 to 2, respectively (Figures 5.2 and 5.3). The risk of upwelling exceeds 0.5 in 43 years for Run 1 and 96 years for Run 2 (Figure 5.4). Model 2 - Model 2 solves Equation (5.3) at annual time-steps while simultaneously accounting for the changes to salinity via Equation (5.1). The mechanical energy was calculated for 2000/01, when the upper halocline remained at approximately a constant depth (dh\/dt=0) (Table 5.1). The near constant depth for 2000/01 indicates an equilibrium upper layer depth existed during that year, for that density difference with the configuration of the SIS discharging to the intermediate layer and the NIS discharging to the upper layer. The mechanical energy needed to maintain this equilibrium depth was calculated to be 864 J m" year"1 based on hx = 7 m, pi-px = 17.71 kg m"3 [S2 = 27, T2 = 12.35 °C, Si = 4 and Tx = 12 °C] and (QA +QE)IA = 1.42 m year"1 [from QA =1.759 x 106 m3 year"1 (Table A.2: March 2000 -March 2001), QEIQA = 0.26 (Appendix A) and halocline area = 1.6 x 106 m2]. These values were used to start the model. Run 1 and Run 2 were configured as per Model 1. Run 2 is detailed first, as data from this configuration was used to determine the mechanical energy. Run 2 (Figure 5.5) found the salinity of the upper and lower layers to decay exponentially (in a similar manner to Model 1 Run 2). The upper layer deepened and at 100 years had a depth of 24 m. The increased upper layer depth caused the wind speed for upwelling to increase with time and at 100 years the wind speed required for upwelling was 79 m s"1 blowing for Ti/4 = 38 min. Thus upwelling was improbable during the first 100 years. Run 1 behaved quite differently to Run 2 due to the re-configuration of the NIS discharge to the intermediate layer (Figure 5.6). The upper layer depth decreased to a minimum depth of 4.3 m after 2 years before deepening rapidly. For Run 1, there was greater injection system discharge to the intermediate layer, that can only be overcome by an upward displacement of the halocline, which results in a shallower upper layer equilibrium 129 depth. The upper layer salinity rose quickly as high salinity intermediate water was transported into the shallower upper layer and the dilution provided by the NIS discharge to the upper layer, as per Run 2, was not available. As a result the density difference decreased rapidly, the upper layer deepened rapidly and the intermediate layer was eroded away after 45 years. Model 2 again demonstrated the benefits of diverting the NIS discharge to the upper layer. DYRESM - The verification of DYRESM used initial conditions of water column temperature, salinity and depth from 6 December 1999; at this time the upper halocline was at a depth of 8.7 m UWL. Lake bathymetry was described using bathymetric data from Rescan (1999a). The model was started on 9 December 1999 using daily inputs of inflow (discharge depth, flow, temperature and salinity), outflow (withdrawal depth and flow) and meteorological conditions. Inflow for the north and south injection systems was measured by flow meters (Chapter 1; Appendix A). The inflow from conveyer tunnel and surface streams were 12% and 23% of the SIS inflow, respectively, based on a comparison of catchment areas. Ouflow from the upper layer was taken as the annual daily average of inflow. The temperature (Appendix D) and salinity (Chapter 1) of the SIS were measured and these values were used for all inflows. The dominant surface stream (draining the west dump and entering the pit lake upper layer from the south-west) and flow into the conveyer tunnel both appear to contain ARD. Meteorological inputs included short-wave radiation, cloud cover, air temperature, vapour pressure, wind speed and precipitation. The Island Copper meteorological station measured short-wave radiation, air temperature, relative humidity and wind speed from 9 December 1999 to 29 November 2000 (Appendix B). The model used cloud cover and simulated surface water temperature to calculate long-wave radiation. Cloud cover was estimated from daily noon observations of weather conditions at Port Hardy Airport (Appendix E). Input vapour pressure was calculated from relative humidity and air temperature (Antenucci and Imerito 2001). The light extinction coefficient (r}A) was set as 0.37 m"1 based on an average secchi depth of 6.2 m; assuming that the secchi depth was equivalent to 10% transparency (Wetzel and Likens 2000) and transparency = exrX-TfaXsecchi) 130 The model predicted the upper layer temperatures to be cooler in the winter and warmer in the summer than observed temperatures (Figure 5.7). The upper layer minimum temperature predicted by the model was 1.10 °C compared to the 4.37 °C observed, whereas, the upper layer maximum temperature predicted by the model was 25.45 °C compared to the 20.54 °C observed. The upper layer was uniformly mixed in the winter and spring for both the model and observations. However, the model predicted a strong seasonal thermocline in the summer that was not observed, indicating that the model underestimates wind mixing. The simulated depth of the upper thermocline (Figure 5.7) and halocline (Figure 5.8) rose more rapidly than was observed. The model reasonably predicted the complex thermal structure of the intermediate layer (Figure 5.9). Mixing from the submerged discharges (NIS, SIS and conveyor tunnel) was evident in the simulation until the end of February, but continued in the observations until the end of April. The cause may be the lower light extinction coefficient of 0.15 m"1 chosen for Figure 5.9, which was done to better replicated the depth of solar heating in the summer. The consequences of a higher light extinction coefficient were that solar heating started earlier and temperatures were higher than observed. In the simulation the thermal stratification (at the top of the intermediate layer) started in February compared to May in reality. In the model the plume mixed water was warmer than observations by ~ 0.1 °C. The model captured the mixing caused by the large plume flow in mid-October. 5.5 Discussion Conceptual Models - Model 1, the simplest conceptual model with two fixed volume boxes, predicted upwelling as a cause of failure. The limitations of Model 1 were made clear in Model 2, which allowed for a variable upper layer depth determined from the balance of upper layer mixing versus injection system inflow and plume entrainment of the halocline. Model 2 predicted a stable lake structure for the NIS discharged to the upper layer, and failure by upper layer deepening for the NIS discharged to the intermediate layer. A limitation of the model was that the mixing energy of the one year was assumed to be representative of all future years, also that the mixing energy available for mixing remains constant, whereas it may decrease as the upper layer depth increases (Imberger and Patterson, 131 1990). The importance of upper layer depth determination suggests a model with a more rigorous mixed layer representation is required and for this reason DYRESM was evaluated. DYRESM - The DYRESM simulation produced an upper layer that was too cool in winter, too warm in summer, did not mix enough in summer and did not deepen as much as was observed. These are all symptoms of insufficient mixing in the upper layer. Increased upper layer mixing would entrain intermediate layer water that was warmer (relative to the upper layer) in the winter and cooler (relative to the upper layer) in the summer. Thus both the interface position and temperature deficiencies would be remedied by greater upper layer mixing. The DYRESM simulation of upper layer mixing was actually worse than was apparent in Figure 5.7, because plume entrainment of the upper halocline, which would lift the halocline by an additional 30% (Chapter 1; Appendix A), was not modeled (see further discussion below). The simulation of upper layer mixing was not improved by varying the model configuration [settings for Figures 5.7, 5.8 and 5.9] such as time-step [daily meteorological forcing], non-neutral stability [off], boundary layer thickness [= 0], layer thickness [minimum = 0.2 m, maximum = 2 m], light extinction coefficient and diffusion volume fraction [= 0]. The model was found to be relatively insensitive to changes to short-wave radiation and cloud cover, which were tested for the range of measurement uncertainty expected for these inputs. Increasing wind speeds and/or increasing the mixing coefficients improved the simulation of summer upper layer temperature, but did not improve the interface position or winter low temperatures. DYRESM did a reasonable job at simulating the plume dynamics and solar heating of the intermediate layer. The light extinction coefficient is a constant in DYRESM (unless CAEDYM the coupled biochemical model is used). However, in the pit lake the transmissivity (amount of radiation passed through water) varies with depth. In particular, the transmissivity was drastically reduced at the upper halocline. So exact replication of the solar heating was not possible. The simulated intermediate layer was ~ 0.1 °C warmer than observed, which is attributed to the model not accounting for entrainment of the upper interface by the impinging plume. The rate of plume entrainment has been established as ~ 0.26 times the 132 injection system flow (Chapter 1; Appendix 2). Fluid entrained from the upper layer at 5 °C with a volume equal to 30% of the injection system flow (which from December 1999 -April 2000 was 2.06 x 106 m3) will decrease the intermediate layer temperature by 0.07 °C. Interestingly the simulated temperature increases with depth (counter-stable) during the high injection system flow period and decreases with depth (stable) when the flow was minimal and solar heating occurs. The simulated salinity gradient always increased with depth, which maintains a stable water column. This simulated behaviour was consistent with observations and the filling box model, discussed in Chapter 3. Due to the problems with upper layer mixing the annual heat budget does not close. At the end of the simulation period the upper layer temperature was 6.64 °C compared to the observed temperature of 8.05 °C. When run into a second year the upper layer water temperatures drop to below 0 °C, at which point the model aborts because it does not simulate ice-cover. Thus long-term simulation has not been possible with DYRESM. Improving the DYRESM mixed layer model with particular regard to increasing the wind mixing is a current project of the model developers at Centre for Water Research (pers. com. Jason Antenucci). Plume Entrainment of the Upper Halocline- Model 1 and DYRESM do not account for entrainment of the upper halocline by the impinging plume. Plume entrainment of the upper halocline was estimated at 30% of the injection system inflow for 1997 to 2001 conditions (Chapter 1; Appendix A). Model 2 accounts for plume entrainment of the halocline using <2E = 2Ck(P2-po)l{pi-P\)QA (Appendix A). The relationship predicts that decreasing Pi-p\ results in increasing QE, and at the limit of Pi-pi = 0, QE —> °°. In Model 2 Run 2, as Pi-p\ got smaller and QEIQA increased to 2-4 in the years prior to failure. The effect of increasing QE was to dilute the intermediate layer more rapidly which further reduces Pi-p\ and increases QE. Thus the plume hastened the reduction in density difference and failure of the stratification. Alternate plume entrainment models express QE «= Ri'3/2, where Ri = gr>lw2 with b and w being some characteristic length and velocity of the plume (Baines 1975; Kumagai 1984; Cotel and Breidenthal 1997). These models concur that QE increases for a reduction in density difference, which is incorporated in Ri. 133 Double diffusion - Double diffusion may be occurring at the upper and lower halocline. At the lower halocline, warm saline fluid was under fresher colder fluid, creating the situation where salinity and temperature were making opposing contributions to the vertical density gradient, the right condition for double diffusive convection. Similar conditions existed at the upper halocline during the winter, when more saline intermediate layer was warmer than the fresher upper layer. Double diffusive convection is known to transport heat and salt (Turner 1965, Huppert 1971), and cause the upward migration of interfaces (Turner 1968, Huppert and Linden 1979). Double diffusion has been observed in meromictic lakes that are geothermally heated (Newman 1976; Sanderson et al. 1986; Wuest et al. 1992). An argument for double diffusive convection across the lower halocline of the Island Copper pit lake was made by Fisher and Lawrence (2001), but no direct observations have been made. Wuest et al. (1992) found that fluxes due to double diffusion were minor compared to shear induced turbulence. In the pit lake turbulent entrainment processes are thought to occur at both haloclines (Chapter 1), and these are expected to be more important for the transport of heat and salt than double diffusion. Management Scenarios - The modeled scenarios have only allowed for variation in NIS discharge to the upper or intermediate layers. The NIS may in the future be discharged to the upper layer for most of the year, but be diverted back to the intermediate layer during the fall "first flush" of ARD, which contains high concentrations of oxidation products (ARD more acidic with higher concentrations of heavy metals). The conceptual models illustrated the importance of maintaining the density difference between the upper and intermediate layers, as diverting the NIS to the upper layer kept the layer fresher and increased the long-term salinity of the intermediate layer. The lake could be manipulated in other ways to maintain the density stratification. Fresh water could be added to the upper layer to lower its salinity, as was done originally with fresh water pumped from the Marble River to help establish the upper layer (Chapter 1). The salinity of the intermediate layer could be replenished by adding seawater, while simultaneously removing intermediate layer water to Rupert Inlet (not to the upper layer as this would increase its salinity). These management options need thorough investigation in terms of the long-term physical stratification and stability and the environmental consequences of any discharge. 134 5.6 Conclusions Model 1, the two layer CSR box model, provided some insight into the dilution of the upper and intermediate layers and the potential for combining model results with wind frequency analysis to predict the risk of lake overturning by upwelling. However, model 2 exposed the limitations of the fixed upper layer depth assumption of model 1. Model 2 incorporated a mixing layer and predicted that the upper layer would deepen with time due to the reduced density difference between upper and intermediate layer. In the case of Run 2 (the NIS discharged to the upper layer), the deeper upper layer made upwelling impossible. For Run 1 (with both NIS and SIS discharge to intermediate layer) the upper layer initially thinned to an equilibrium depth of 4.3 m, which was followed by a rapid decrease in density difference and consequential deepening of the upper layer to the bottom of the intermediate layer in 45 years. While the conceptual models may not be accurate in their ability to predict the time of destratification, they both clearly illustrate the importance of maintaining the density difference between upper and intermediate layers. Maintaining the density difference can be improved by discharging the NIS to the upper layer (provided water quality considerations are satisfied). DYRESM was used to simulate the pit lake hydrodynamics for December 1999 to November 2000. The model was able to resolve the short-term and annual trends such as mixed layer behaviour and plume dynamics. However, the model did not simulate the mixed layer well enough to balance the annual heat budget of the upper layer. For this reason multi year modeling was not possible. DYRESM does not include entrainment of the upper halocline by the impinging plume. Plume entrainment of the halocline provides an additional source of dilution to the intermediate layer that was estimated at 30% of the injection system flow for current conditions (Chapter 1; Appendix A). Plume entrainment of the halocline is anticipated to increase with time as the density difference between upper and intermediate layers decrease. Therefore, plume entrainment needs to be added to the DYRESM before long-term simulations are attempted. 135 Table 5.1. Annual upper halocline depths from CTD profiles. Date Upper Halocline Depth*1 (mUWL) Annual Change to Upper Halocline Depth (m) 24 July 1997 16.92 -08 June 1998 * 2 12.81 4.11 09 July 1999 9.31 3.50 06 July 2000 * 3 7.56 1.74 15 March 2001 7.64 -0.08 * 4 *' Halocline depths estimated from method described in Appendix A.. * 2 Estimated by averaging upper layer depths from 03 April 1998 and 14 August 1998. * 3 Prior to 29 June 2000 both the SIS and NIS were discharging to the intermediate layer, but after this date the NIS was diverted to the upper layer. * 4 Note year not complete. 136 Maximum Annual Wind Speed at 10m (m/s) Figure 5.1. Frequency distribution of Island Copper annual maximum hourly wind speeds adapted from Port Hardy Airport wind frequency distribution for 1953 to 2000. Curve is Equation (5.11). 137 Figure 5.2. Model 1 Run 1 (both NIS and SIS discharging to intermediate layer) showing upper and intermediate layer salinity (monthly), wind speed for upwelling (monthly) and duration required for upwelling (annual). 138 Time (years) Figure 5.3. Model 1 Run 2 (SIS discharging to intermediate layer and NIS to the upper layer) showing upper and intermediate layer salinity (monthly), wind speed for upwelling (monthly) and duration required for upwelling (annual). 139 Time (years) Figure 5.4. Model 1 predictions of return period for the wind required to cause upwelling [—] and risk of upwelling [—]. 140 Figure 5.5. Model 2 Run 2 (SIS discharging to intermediate layer and NIS to the upper layer) showing upper and intermediate layer salinity (annual), wind speed for upwelling (annual) and duration required for upwelling (annual). 141 40 35 4-30 4-S Intermediate 5 + Upper Layer Depth Wind Speed ' I 1 1 1 1 I 20 40 60 80 Time (years) 100 120 100 4- 90 4- 80 4- 70 4- 60 50 £ 40 + 30 4- 20 + 10 n. T3 c 140 Figure 5.6. Model 21 Run 1 (both NIS and SIS discharging to intermediate layer) showing upper and intermediate layer salinity (annual), wind speed for upwelling (annual) and 'duration required for upwelling (annual). 142 (a) Dec Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec Date Figure 5.7. Upper layer temperature for (a) DYRESM (run icpit7a) with r\A = 0.37 m"1 and (b) observed from thermistors at 0.5-1 m spacing. Dotted lines are 98% halocline shoulders from 6 CTD casts. IH 3 Figure 5.8. Upper halocline salinity simulated by DYRESM (run icpit7a) with T | A = 0.37 m' . Dotted lines are 98% halocline shoulders from 6 CTD casts. 144 Figure 5.9. (a) DYRESM (run icpit7j) with T | A = 0.15 m"1. (b) Observed daily temperature from thermistors spaced at 12 m, 15 m, 20 m, 30 m, 100 m and 200 m. (c) Injection system flow [—] and temperature [—]. 145 6 CONCLUSIONS AND RECOMMENDATIONS 6.1 Conclusions The Island Copper Mine pit lake is both a practical solution to an environmental problem and a unique limnological experiment. As such, my motivation has been twofold. Firstly, to improve the understanding of the meromictic pit lake to optimize its performance as a passive treatment system so that the technology can be applied to other sites around the world. Secondly, in the broader context, to investigate fundamental aspects of environmental fluid mechanics applicable to lakes and any stratified environment receiving a pollutant discharge. In this regard the pit lake is rich in environmental fluid mechanics including limnology, meromixis, buoyant plumes, filling box behaviour, thermal convection and double diffusion. The Island Copper pit lake has evolved into three distinct layers: a brackish and well mixed upper layer; a plume stirred intermediate layer; and a thermally convecting lower layer (Figure 1.2). To date the upper halocline has risen over 10 m, due to the injection of buoyant acid rock drainage (ARD) into the base of the intermediate layer. An upper layer "equilibrium" depth is expected in the near future reflecting a balance between meteorological forcing from above and rising ARD plumes from below. The lower halocline depth fluctuates seasonally: deepening in the winter due to high ARD flows, which erode the lower halocline; and rising in the summer when thermal convection in the lower layer dominates over reduced intermediate layer stirring. The proposed treatment of acid rock drainage by metal-sulphide precipitation using sulphate-reducing bacteria is dependent on anoxia developing in the intermediate and lower layers. The lower layer is close to anoxia (DO = 0.03 mg/L); whereas, with a DO = 2.5 mg/L, the intermediate layer has a long way to go. Without metal-sulphide precipitation metal concentrations in the pit lake are likely to increase, although they may be controlled by fertilization via biological uptake. 146 The addition of ARD and chemical changes to the pit lake seawater have sufficiently changed the chemical constituents so that 1978 Practical Salinity Scale (PSS78) no longer provides an accurate estimate of salinity. Salinity determination by chlorinity, drying and weighing, and chemical constituents were also unsatisfactory. The density method proved to be the most accurate method for the estimation of salinity because equation of state formulae are less sensitive to changes in chemical composition than conductivity methods such as PSS78. Based on the density method a salinity correction to PSS78 salinity was developed for application to field CTD measurements. The salinity correction was found to be a function of Spss and temperature, and a unique equation for each layer was required. Observations suggest that non-seawater sourced ions are increasing with time, and are sourced primarily from the ARD inflow and dissociation of ions from the pit wall. To compensate, a time variant salinity correction was developed. The true salinity (including correction) of the lower layer was greater than the original seawater. Future changes to the ARD and chemical changes to the pit lake may invalidate the salinity correction, particularly for the lower layer where AS is the largest and the water chemistry the most complex. This work illustrates the importance of verifying the salinity and density in biochemically active systems such as mine pit lakes. The intermediate layer is similar in appearance to the filling box model (Baines and Turner, 1969) as ARD is injected at depth and rises as buoyant plumes within the confined receiving environment. The filling box model analogy explains the gradients of salinity, temperature and density observed in the intermediate layer. Plume entrainment of upper layer fluid is not important to the filling box model because entrainment is small. However, application of the filling box model is limited because of: the complexities of the diffuser; bottom heating from the warmer lower layer, which causes the bottom of the intermediate layer to mix; and the intermittency of the injection system flow, which tends to stratify the top of the intermediate layer. Furthermore, during the summer heating of the top of the intermediate layer and minimal injection system flow causes the plume to be trapped within the intermediate layer. The halocline over the south injection system (SIS) was higher and thinner than at the central mooring at times when the plume rose the full height of the intermediate layer 147 (similar behaviour is expected at the north injection system). The plume lifts the interface at the SIS as the vertical plume flow impinges against the interface and is redirected horizontally. Plume entrainment of the interface causes it to be thinner at the SIS than at the central mooring. Greater temperature fluctuations were observed at the SIS during the winter due to the sharper thermocline at that time. Spectra analysis was used to identify the vertical mode-1 lake-length seiche and higher horizontal modes, which match the Defant prediction well. Additional spectral peaks were attributed to vertical mode-1 lake-width seiches and its higher horizontal modes. Winds are sufficiently variable in direction to generate both types of waves. Internal waves were predicted to degenerate by the development of solitons within a few days, but observations showed the persistence of waves for many days and 0(100) wave periods. There was no evidence of the development of solitons, rather wave amplitudes were slowly reduced with time, suggesting wave degeneration by wave damping. The parameterization for wave dampening is based on boundary mixing over the lake bottom, which seems unrealistic for such a deep lake and the predicted decay periods are too long. Additional work is required to determine the role of sidewall mixing and intermediate layer turbulence on wave damping. Maintaining a permanent stratification is crucial to the performance of the passive treatment system because the ARD is expected to be generated for over 100 years. Two conceptual models were used to evaluate the long-term evolution of the pit lake and potential failure modes. Model 1, a two layer continuously stirred reactor box model, predicts that upwelling is the likely failure mode because the density difference between the upper and intermediate layer decreases with time, which decreases the wind speed required for upwelling. Extreme wind frequency analysis was used to estimate the likelihood of failure. In Model 2 the upper layer deepens in response to the decreasing density difference between the two layers. This model predicts that upwelling is not likely and that erosion of the intermediate layer is the expected failure mode. Both models demonstrate that maintaining the density difference between the upper and intermediate layers extends the life of the stratification. This extension can be achieved by diverting the north injection system discharge from the intermediate layer to the upper layer. 148 The one-dimensional lake model DYRESM was applied to the pit lake. It was able to resolve the short-term and annual trends such as mixed layer behaviour and plume dynamics. However, the model does not yet simulate the mixed layer well enough to balance the annual heat budget of the upper layer. For this reason multi year modeling was not possible. DYRESM does not include entrainment of the upper halocline by the impinging plume. Plume entrainment of the halocline provides an additional source of dilution to the intermediate layer that was estimated at 26% of the injection system flow. Plume entrainment of the halocline is anticipated to increase with time as the density difference between upper and intermediate layers decreases. Therefore, plume entrainment needs to be included in any model before long-term simulations are attempted. 6.2 Recommendations A good understanding of the physics of the Island Copper mine pit lake has been developed, but more detailed work in some areas is required. In collaboration with Craig Stevens (NIWA, New Zealand) additional work has been undertaken to quantify the mixing processes via direct measurements of turbulent dissipation using temperature microstructure. Specifically measurements of turbulent dissipation at the side walls versus the middle of the lake made in July 2000 will improve the understanding of side wall mixing and its role in wave degeneration, plus the influence of wind on mixing at the upper halocline. A second set of measurements February 2002 focused on the plume with the goal of quantifying plume turbulent dissipation and plume entrainment of the upper halocline. A secondary goal was to establish whether double diffusion was active at the upper halocline, as the winter conditions of warm saline water underlying cool fresh water are conducive to double diffusive convection. Analysis of these data will be carried out in the future. Thermal convection in the lower layer driven by geothermal heating was hypothesized to explain the well mixed temperature and salinity profiles through that layer. Lower interface transport processes are dependent on this mechanism, so direct measurements of the thermal convection by temperature microstructure or tracer experiment would be insightful. Again conditions are right for double diffusive convection at the lower halocline and further work is required to confirm its presence and quantify its importance. Important to the long-term stratification of the pit lake is the role of plume entrainment at the 149 upper halocline. It is apparent from measurements, energy arguments and scaling arguments that that plume entrainment exists and will become more important as the density difference between the upper and intermediate layer decreases. Plume entrainment under strongly stratified conditions requires more research and is ideally suited to laboratory experiments. Additional work on the long-term predictive modeling is required to better estimate the stability of the pit lake and its viability as a treatment system for acid rock drainage. Such a model should include plume entrainment as well as lake physical mixing processes. Also requiring modeling is the biochemical system to establish the likelihood of anoxia and sulphate reducing bacteria for the treatment of acid rock drainage by metal-sulphide precipitation. This model once developed could be used as a management tool to evaluate the impact of any changes and to optimize the system. 150 REFERENCES Allen, S.E. 1991. Chemocline Erosion in Proposed Island Copper Mine Lake. Prepared for BHP-Utah Mines. Allen, S.E. 1992. Chemocline Erosion in Proposed Island Copper Mine Lake - Part II. Prepared for BHP-Utah Mines. Allen, S.E. 1993. Chemocline Erosion in Proposed Island Copper Mine Lake - Part III. Prepared for BHP-Utah Mines. American Public Health Association, American Water Works Association, and Water Pollution Control Federation. 1995. Standard methods for the examination of water and wastewater. 19th Edition. 2.46-2.49. Antenucci, J., and A. Imerito. 2000. The CWR Dynamic Reservoir Simulation Model: DYRESM. Centre for Water Research, University of Western Australia. Aspinall, C. 1995. The story of Island Copper. BHP Minerals Canada, Vancouver. Baines, W. D., and J. S. Turner. 1969. Turbulent buoyant convection from a source in a confined region. J. Fluid. Mech. 37: 51-80. Baines, W.D. 1975. Entrainment by a plume or jet at a density interface. J. Fluid Mech. 68: 309-320. BHP Minerals Canada LTD. 1996. Island Copper Mine closure plan - addendum. Prepared for Vancouver Island Mine Development Review Committee, Ministry of Energy, Mines and Petroleum Resources, Victoria, B.C. BHP Minerals Canada LTD. 1999. 1998 Annual environmental assessment report, BHP Minerals Canada Ltd., Island Copper Mine, Rupert Inlet, British Columbia. Prepared for submission to The Provincial of British Columbia, Victoria, B.C. 151 BHP Minerals Canada LTD. 2000. 1999 Annual environmental assessment report, BHP Minerals Canada Ltd., Island Copper Mine, Rupert Inlet, British Columbia. Prepared for submission to The Provincial of British Columbia, Victoria, B.C. BHP Billiton. 2001. 2000 Annual environmental assessment report, BHP Minerals Canada Ltd., Island Copper Mine, Rupert Inlet, British Columbia. Prepared for submission to The Provincial of British Columbia, Victoria, B.C. Bohrer, B., H. Heidenreich, M. Schimmele, and M. Schultze. 1998. Numerical prognosis for salinity profiles of future lakes in the opencast mine Merseburg-Ost. International Journal of Salt Lake Research. 7: 235-260. Boehrer, B. 2001. In-situ measurement of density and stability. Proceedings of the Sixth Workshop on Physical Processes in Natural Waters, July, Girona. 113-117. Chapra, S.C. 1997. Surface water-quality modeling. The McGraw-Hill Companies, Inc. Chen, C.T., and F.J. Millero. 1986. Precise thermodynamic properties for natural waters covering only the limnological range. Limnol. Oceanogr. 31: 657-662. Chow, V.T. 1959. Open Channel Hydraulics. McGraw-Hill Book Company Inc. New York. Chikita, K., Y. Hosogaya, and S. Natsume. 1993. The characteristics of internal waves in a caldera lake induced from field measurements: Lake Kuttara, Hokkaido. Jpn. J. Limnology. 54(3): 213-224. Chikita, K., K. Sakata, and S. Hino. 1995. Transportation of suspended sediment slowly settling in a caldera lake. Jpn. J. Limnology. 56(4): 245-257. Colomer, J., T. Serra, J. Piera, E. Roget, and X. Casamitjana. 2001. Observations of a hydrothermal plume in a karstic lake. Limnol. Oceanogr. 46: 197-203. Cotel, A.J. and R.E. Breidenthal. 1997. A model of stratified entrainment using vortex persistence. Applied Scientific Research. 57: 349-366. Cotel, A.J., J.A. Gjestvang, N.N. Ramkhelawan, R.E. Breidenthal. 1997. Laboratory experiments of a jet impinging on a stratified interface. Experiments in Fluids. 23: 155-160. 152 Crawford, G. B., and R. W. Collier. 1997. Observations of a deep mixing event in Crater Lake, Oregon. Limnol. Oceanogr. 42: 299-306. Dagenais, PJ. 1996. An investigation into the construction, excavation and geochemical history of a waste rock dump and implications for long-term water quality at the Island Copper Mine. Port Hardy, British Columbia. M.A.Sc. Thesis, University of British Columbia. Defant, A. 1961. Physical Oceanography. Pergamon. Donelan, M.A. 1998. Air-Water Exchange Processes. In Physical Processes in Lakes and Oceans, Coastal and Estuarine Studies, American Geophysical Union. 54: 19-36. Dunbar D.S. and Hodgins S.L.M. 1990 An investigation of the evolution of water properties in the Island Copper Mine following sea water flooding. By Seaconsult Marine Research Ltd. Prepared for BHP-Utah Mines Ltd. Dunbar D.S. 1995. Island Copper Mine Closure plan - Numerical Simulation of Water Properties in the Flooded Pit. By Ocean Applied Research Ltd. Prepared for BHP Minerals Canada Ltd. Ellis, D.V., T.F. Pedersen, G.W. Poling, C. Pelletier, and I Home. 1995. Review of 23 years of STD: Island Copper Mine, Canada. Marine Georesources and Geotechnology. 13: 59-99. Environment Canada. 1993. Canadian Climate Normals 1961 - 1990, British Columbia. Canadian Climate Program, Atmospheric Environment Services, Environment Canada. Emery, W.J., and R.E. Thomson. 1998. Data analysis methods in physical oceanography. 1st ed., New York, Pergamon. Fernando, H.J.S. 1991. Turbulent mixing in stratified fluids. Ann. Rev. Fluid. Mech. 23: 455-493. Fischer, H.B., E J . List, R.C.Y. Koh, J. Imberger, and N.H. Brooks. 1979. Mixing in inland and coastal water. Academic Press. 153 Fisher, T.S.R. 1995. Dilution of axisymmetric buoyant jets and surface spreading fields. Masters of Engineering Thesis. University of Canterbury, Christchurch, New Zealand. Fisher, T. S. R., G. A. Lawrence, and M. J. Wilton, 1999. The evolution of the Island Copper Mine pit lake. Proceedings of ACSE International Water Resources Conference (CDRom), Seattle, U.S.A.. Fisher, T. S. R., and G. A. Lawrence. 2000. Observations at the upper halocline of the Island Copper pit lake. Proceedings of 5th International Symposium on Stratified Flows, Edited by G.A. Lawrence, R. Pieters and N. Yonemitsu, Vancouver, British Columbia. 413-418. Fisher, T. S. R., and G. A. Lawrence. 2001. Double diffusion in the Island Copper Mine pit lake. Proceedings of 6th Workshop on Physical Processes in Natural Waters, Edited by Xavier Casamitjana, 27-29 June 2001, Girona, Catalonia, Spain. 107-111. Geller, W., H. Klapper and W. Salomons. 1998. Acidic mining lakes - acid mine drainage, limnology and reclamation. Springer. Jellison, R., and J. M. Melack. 1993. Meromixis in hypersaline Mono Lake, California. 1. Stratification and vertical mixing during the onset, persistence, and breakdown of meromixis. Limnol. Oceanogr..38: 1008-1019. Jellison, R., J. Romero, and J. M. Melack. 1998. The onset of meromixis during restoration of Mono lake, California: Unintended consequences of reducing water diversions. Limnol. Oceanogr. 43: 706-711. Hall, K.J., and T.G. Northcote. (1986) Conductivity - temperature standardization and dissolved solids estimation in a meromictic saline lake. Can. J. Fish. Aquat. Sci. 43: 2450-2454. Hamblin, P. F., C. L. Stevens, and G. A. Lawrence. 1999. Simulation of vertical transport in mining pit lake. J. Hyd. Eng. 125: 1029-1038. Hamilton, D.P. and Schladow, S.C. 1997. Prediction of water quality in lakes and reservoirs: Part 1 - Model Description. Ecological Modelling, 96: 91-110. 154 Henderson, F.M. 1966. Open channel flow. Macmillan Publishing Co., Inc. New York. Hill, K.D., T.M. Dauphinee, and DJ. Woods. 1986. The extension of the PSS-78 to low salinities. IEEE J. Oceanic Eng. 11:109. Home, A J. , and C.R. Goldman. 1994. Limnology. 2nd Edition, McGraw-Hill. New York. Horn, D.A., J. Imberger, and G.N. Ivey. 2001. The degeneration of large-scale interfacial gravity waves in lakes. J. Fluid Mech.. 434: 181-207. Huppert, H.E. & Linden, P.F. 1979. On heating of a stable salinity gradient heated from below. J. Fluid. Mech. 95: 431-464. Huppert, H.E. 1971. On the stability of a series of double diffusive layers. Deep-Sea Res. 18: 1005-1021. Hyndman, R.D. 1976. Heat flow measurements in the Inlets of Southwestern British Columbia. J. Geophysical Research. 81: 337-349. Imberger, J. and Patterson, J. C. 1981. A dynamic reservoir simulation model - DYRESM:5. In "Transport Models for Inland and Coastal Waters". H.B. Fischer (ed). Academic. pp310-361. Imberger, J. and Patterson, J.C. 1990. Physical Limnology. Advances in Applied Mechanics, 27: 303-475. Imberger, J. 2001. Characterizing the dynamic regimes of a lake. Proceedings of Physical Processes in Natural Waters. Edited by X. Casamitjana. Girona, Catalonia, Spain. 77-92. Imboden, D.M and Schwarzenbach, R.P. 1985 Spatial and temporal distribution of chemical substances in lakes: modeling concepts, in Chemical Processes in Lakes, Stumm, W. ed., 1-30. Imboden, D.M., and A. Wiiest. 1995. Mixing mechanisms in lakes. In Physics and Chemistry of Lakes. Edited by A. Lerman, D.M. Imboden and J.R. Gat. 83-138. Killworth, P.D. and J.S. Turner. 1982. Plumes with time-varying buoyancy in a confined region. Geophys. Astrophys. Fluid Dynamics. 20: 265-291. 155 Killworth, P.D. and E.C. Carmack. 1979. A filling-box model of river-dominated lakes. Limnol. Oceanogr. 24: 201-217. Kranenburg, C. 1985. Mixed-layer deepening in lakes after wind setup. Journal of Hydraulic Engineering. I l l : 1279-1297. Kraus, E.B., and J.S. Turner. 1967. A one-dimensional model of the seasonal thermocline U. The genera] theory and its consequence. Tellus, 19: 1104-1121. Kumagai, M. 1984. Turbulent buoyant convection from a source in a confined two-layer region. J. Fluid Mech. 147: 105-131. Laval, B, J. Imberger, and B.R. Hodges. 2002. Modeling surface layer energetics in lakes: spatial and temporal variations. Submitted to Limnol. Oceanogr. LeBlond, P.H., and Mysak, L.A. 1978. Waves in the Ocean. Elsevier Oceanography Series 20. Elsevier Scientific Publishing Company. Amsterdam. Lemmin, U. 1987. The structure and dynamics of internal waves in Baldeggersee. Limnol. Oceanogr. 32(1): 43-61. Lemmin, U., and C H . Mortimer. 1986. Tests of an extension to internal seiches of Defant's procedure for determination of surface seiche characteristics in real lakes. Limnol. Oceangr. 31(6): 1207-1231. Lewis, E.L., and R.G. Perkin. 1978. Salinity: its definition and calculation. J. Geophys. Res. 83(C1): 466-478. Lewis, E.L., and R.G. Perkin. 1981. The practical salinity scale 1978: conversion of existing data. Deep-Sea Research. 28(A): 307-328. Li, M.G. 1991. Chemistry of the drainage from a waste dump at BHP-Utah Mines Ltd, Island Copper mine. Masters Thesis. University of British Columbia. Lindenschmidt, K.-E. and Hamblin, P.F. 1997. Hypolimnetic Aeration in Lake Tegel, Berlin. Wat. Res. Vol 31: 1619-1628. Lister, D. 1994. An assessment of acid rock drainage potential of waste rock and implications for long term weathering of the North Dump at Island Copper Mine, Port Hardy, B.C. Masters Thesis. University of British Columbia. 156 MacLatchy, M.R. 1999. Radially Spreading Surface Flows. Ph.D. Thesis. University of British Columbia. Mannis, P.C. 1973. A filling-box model of the Red Sea. Memoires, Societe Royale des Sciences de Liege. 6: 153-166. Martin, J.M. 1985. The Pavin Crator Lake. Chemical Processes in Lakes, W. Stumm, ed., Wiley, New York. 169-188 McManus, J., R. W. Collier, C. A. Chen, and J. Dymond. 1992. Physical properties of Crater Lake, Oregon: A method for the determination of a conductivity- and temperature-dependent expression for salinity. Limnol. Oceanogr. 37(1): 41-53. Millero, F.J. and M.L. Sohn. (1992) Chemical Oceanography. CRC Press, Inc. Millero, F.J. 2000a. The equation of state of lakes. Aquatic Geochemistry. 6: 1-17. Millero, F.J. 2000b. Effect of changes in the composition of seawater on the density- salinity relationship. Deep-Sea Research I. 476: 1583-1590. Monismith, S.G. 1985. Wind-forced motions in stratified lakes and their effect on mixed layer shear. Limnol. Oceanogr. 30(4): 771-783. Monismith, S.G. 1986. An experimental study of the upwelling response of stratified reservoirs to surface shear stress. J. Fluid Mech. 171: 407-439. Morin, K.A. and N.M. Hutt. 1996. Summary of studies on acidic drainage and metal leaching at Island Copper Mine, Port Hardy, British Columbia. Minesite Drainage Assessment Group. Prepared for BHP Minerals Canada. Mortimer, C H . 1979. Strategies for coupling data collection and analysis with dynamic modeling of lake motion. Hydrodynamics of Lakes; Proceedings of a Symposium, Lausanne, Switzerland. Elsevier. 287-322. Morton, B. R., G. I. Taylor, and J. S. Turner. 1956. Turbulent gravitational convection from maintained and instantaneous sources. Proc. R. Soc. Lond. Ser. A 234: 1-23. Muggli, D.L., C A . Pelletier, G.W. Poling, and E . C Schwamberger. 2000. Injected ARD plume behaviour in a pit lake utilizing in situ dye studies. Proceedings from the Fifth 157 Internationa] Conference on Acid Rock Drainage, Denver, Colorado. Society for Mining, Metallurgy, and Exploration, Inc. 305-318. Nokes, R.I. 1988. On the entrainment rate across a density interface. J. Fluid Mech. 188: 185-204. Newman, F.C. 1976. Temperature steps in Lake Kivu: a bottom heated saline lake. J. Phys. Oceanogr. 6: 157-163. Northcote, T. G., and K. J. Hall. 1983. Limnological contrasts and anomalies in two adjacent saline lakes. Hydrobiologia. 105: 179-194. Northcote, T. G., and K. J. Hall. 2000. Short-term (decadal,annual,seasonal) changes to the limnology of a saline uniVbimeromictic lake: causes and consequences. Verh. Internat. Verein. Limnol. 27. Ozretich R.J. 1975. Mechanism for deep water renewal in Lake Nitinat, a permanently anoxic fjord. Estuarine and Coastal Marine Science, 3: 189-200. Papantoniou, D., and E.J. List. 1989. Large-scale structure in the far field of buoyant jets. J. Fluid Mech. 209: 151-190. Perkin, R.G., and E.L. Lewis. 1980. The practical salinity scale 1978, fitting the data. IEEE J. Oceanic Eng. 5: 9-16. Perry, K.A., and T.F. Pedersen., 1993. Sulphur speciation and pyrite formulation in meromictic ex-fjords. Geochimica et Cosmochimica Acta. 57: 4405-4418. Pieters, R., L.C. Thompson, L Vidmanic, M. Roushorne, J. Stockner, K. Hall, M. Young, S. Pond, M. Derham, K. Ashley, B. Lindsay, G. Lawrence, D. Sebastian, G. Scholter, F. McLaughlin, A. Wiiest, A. Matzinger and E. Carmack. 1999. Arrow Reservoir limnological and trophic status report, year 2 (1998/99). RD 72, Fisheries Branch, Ministry of Environment, Lands and Parks, Province of British Columbia. Rescan. 1994. Preliminary design of a meromictic lake in the flooded pit: a passive sulphide precipitation system for ARD. By Rescan Environmental Services Ltd. Prepared for BHP-Utah Mines Ltd. 158 Rescan Environmental Services Ltd. 1999a. Tracer and biochemical studies to determine the status of meromictic conditions in the flooded Island Copper pit. Report prepared for BHP Minerals Canada Ltd. Rescan Environmental Services Ltd. 1999b. Dye dispersion study of injected plume in pit lake, march 1999. Report prepared for BHP Minerals Canada Ltd. Rescan Environmental Services Ltd. 2001. Inspection of pit lake injection systems. Report prepared for BHP Billiton Ltd. Resio, D., S. Bratos, and E.F. Thompson. 2001. Meteorology and wave climate. In Coastal Engineering Manual, Chapter 2. Coastal and Hydraulics Laboratory, U.S. Army Corps of Engineers. II—2: 1-67. Romero, J. R., and J. M. Melack. 1996. Sensitivity of vertical mixing in a large saline lake to variations in runoff. Limnol. Oceanogr. 41: 955-965. Sanderson, B., K. Perry, and T. Pedersen. 1986. Vertical diffusion in meromictic Powell Lake, British Columbia. Journal of Geophysical Research. 91: 7647-7655 Schimmele, M. 1998. Density measurements in mining lakes. Proceedings of the Third Workshop on Physical Processes in Natural Waters, Magdeburg. Schladow, S. G., and I. H. Fisher. 1995. The physical response of temperate lakes to artificial destratification. Limnol. Oceanogr. 40: 359-373. Spigel, R.H., and J. Imberger. 1980. The classification of mixed-layer dynamics in lakes of small to medium size. J. Phys. Oceanogr. 10: 1104-1121. Spigel, R.H., J. Imberger, and K.N. Rayner. 1986. Modeling the diurnal mixed layer. Limnol. Oceanogr. 31: 533-556. Spigel, R.H., and J.C. Priscu. 1996. Evolution of temperature and salt structure of Lake Bonney, a chemically stratified Antarctic lake. Hydrobiologia. 321: 177-190. Stevens, C.L., Lawrence, G.A. Rogers C K . and Hamblin P.F. 1994. Modeling the thermal stratification of water-filled mine pits. Environmental Fluid Mechanics Group, Department of Civil Engineering, University of British Columbia. Supplied to Environment Canada. June 8, 1994. 159 Stevens, C. L., and G. A. Lawrence. 1997a. The effect of sub-aqueous disposal of mine tailings in standing waters. Journal of Hydraulic Research. 35: 147-159. Stevens, C.L. and Lawrence G.A. 1997b. Estimation of wind-forced internal seiche amplitudes in lakes and reservoirs, with data from British Columbia, Canada. Aquat. Sci. 59, 115-134. Stevens, C. L., and G. A. Lawrence. 1998. Stability and meromixis in a water filled mine pit. Limnol. Oceanogr. 43: 946-954. Stull, R.A. 2000. Meteorology for scientists and engineers. Second Edition. Brooks/Cole Thomson Learning. Sverdrup, H.U., M.W. Johnson, and R.H. Fleming. 1942. The oceans their physics, chemistry and general biology. Prentice-Hall, Inc. Thompson, R.O.R.Y. and J. Imberger. 1980. Response of a numerical model of a stratified lake to wind stress. Second International Symposium on Stratified Flows, IAHR, Edited by T. Carstens and T. McClimans, Trondheim, Norway. 562-570. Thomson, R.E. 1981. Oceanography of the British Columbia Ccoast. Canadian Special Publication of Fisheries and Aquatic Sciences 56, Department of Fisheries and Oceans, Ottawa. Turner, J.S. 1965. The coupled turbulent transports of salt and heat across a sharp density interface. Int. J. Heat Mass Transfer. 8: 759-767. Turner, J.S. 1968. The behaviour of a stable salinity gradient heated from below. J. Fluid. Mech. 33: 183-200. Turner, J. S. 1973. Buoyancy Effects in Fluids. Cambridge University Press. Turner, J.S. 1984. Turbulent entrainment: the development of the entrainment assumption, and its application to geophysical flows. J. Fluid Mech. 173: 431-471. Unesco. 1981. Tenth report of the joint panel on oceanographic tables and standards. Unesco technical papers in marine science 36. Sidney, B.C., Canada, 1-5 September, 1980. Valdivia, F.J.R. 2001. A three-dimensional hydrodynamic and transport model for lake environments. Ph.D. Thesis, University of California, Davis. 160 Viessman, W. Jnr., and G.L. Lewis. 1996. Introduction to hydrology. Fourth Edition. Harper Collins College Publishers. Wetzel, R.G. and G.E.Likens. 2000. Limnological Analysis. 3rd Edition. Springer. New York Ward, P. R. B., K. J. Hall, T. G. Northcote, W. Cheung and T. Murphy. 1990. Autumnal mixing in Mahoney Lake, British Columbia. Hydrobiologia. 197: 129-138. Watkins, T.H. 2000. Hard Rock Legacy. National Geographic Magazine, March 2000. Wilkinson, D.L. 1984. Purging of saline wedges from ocean outfalls. Journal of Hydraulic Engineering. 110: 1815-1829. Wilton, M.J. 1998. The Evolution of the Island Copper Mine Pit Lake. MASc. Thesis in Department of Civil Engineering, University of British Columbia. Vancouver, B.C. Wood, I.R., R.G. Bell and D.L. Wilkinson. 1993. Ocean disposal of wastewater. World Scientific. Singapore. Wu, J. 1973. Wind-induced turbulent entrainment across a stable density interface. J. Fluid Mech. 61: 275-287. Wu, J. 1977. A note on the slope of a density interface between two stable stratified fluids under wind. J. Fluid Mech. 81: 335-339. Wiiest, A., N.H. Brooks, D.M. Brooks. 1992. Bubble plume modeling for lake restoration. Water Resources Research. 28: 3235-3250. Wuest, A., W. Aeschbach-Hertig, H. Baur, M. Hofer, R. Kipfer, and M. Schurter. 1992. Density structure and tritium-helium age of deep hypolimnetic water in the northern basin of Lake Lugano. Aquatic Sciences. 54: 204-218. Wuest, A., G. Pieke, and J. D. Halfman. 1996. Combined effect of dissolved solids and temperature on the density stratification of Lake Malawi. Paleoclimatology of the East African Lakes. Gordon and Breach, Toronto. 161 APPENDIX A Water Balance A.l Introduction The movement of water between layers was investigated by undertaking a water balance on the Island Copper pit lake. CTD data was used to determine the interface position, from which changes in layer volumes and volume fluxes (volume change per day) were calculated. Quantitative results were affected by the accuracy of CTD depth determination. Nevertheless, the water balance illuminates transport mechanisms at the upper and lower haloclines. This analysis extends Chapter 1, which inferred transport mechanisms from the temporal variation in interface position. The analysis was included as an appendix because the uncertainties of depth and flow measurements limit the quantitative interpretation, although trends are still clear, and the work illuminates the physical transport processes. In summary, the upper halocline has been rising due to the injection of ARD to the intermediate layer and entrainment of the upper layer by the impinging plume. However, upper layer mixing processes such as wind and penetrative convection were becoming increasing effective at eroding the upper halocline, due to the reduced upper layer thickness. An upper layer equilibrium thickness is expected when these competing processes balance (Figure A.l). The entrainment from the plume impinging against the upper halocline was -0.29 times the injection system flow, which coincides reasonable well with the entrainment rate predicted from energy considerations. The lower halocline is observed to fluctuate seasonally as a result of the balance between intermediate layer stirring and lower layer thermal convection. 162 A.2 Method The halocline depth was defined as the depth at which an average salinity of the two adjacent layers was found. The upper and lower shoulders of each halocline were defined as the position at which the salinity change exceed threshold values for the first and last time. Threshold values for the salinity change at the upper and lower interfaces were chosen as 0.15 and 0.015, respectively, as these most consistently captured the interface shoulders. The volume of each layer was calculated from halocline depths and lake water levels using an accurate bathymetric survey (Rescan 1999a). Depths are specified relative to the surface water level (WL) or the ultimate water level (UWL) datum, which was 2.44 m above mean sea level. Changes in layer volumes between field trips result from volumetric fluxes between layers i.e. across the upper and lower haloclines and discharge of ARD to the intermediate layer. A layer volumetric flux of 10,000 m3 d"1 over any three month period reflects changes to the upper and lower haloclines of approximately 53 cm and 290 cm, respectively. The large difference results from the smaller surface area of the lower halocline relative to the upper halocline. The uncertainty of CTD depth measurements was 23 cm (Appendix C). The depth uncertainty limits the quantitative water balance results because changes in interface position between CTD casts, particularly at the upper interface, are of similar magnitude. Nevertheless, trends over longer periods are still clear and physical transport processes are illuminated. SIS and NIS flows are used in the ensuing calculations and the accuracy of these measurements requires comment. Flows were measured by Milltronics OCMIII acoustic flow meters, which in fact measure the depth of water in the injection system pipes, and internally calculate flow using Mannings equation (Henderson, 1966). The specified depth accuracy is 1 mm, which is equivalent to 0(1%) of the flow. However, the accuracy of the calculated flow was more dependent on the estimates of Manning roughness (n) and pipe slope (S) programmed into the instrument. The flow metres calculated the flow using the default pipe slopes of 0.025 rather than the actual slopes, which are estimated at 0.0186 for the SIS and 0.0344 for the NIS (pers. com. Shane Uren, Rescan Environmental Services, 163 Ltd.). These slope differences would result in errors at the mean flow of -16 % at the SIS and +15 % at the NIS, which coincidentally sum to an error for the total flow of just -4 %. Verification of flow measurements was not undertaken, as alternate measurements of flow in the injection systems were not possible. A.3 Results Upper and lower halocline depths are given in Table A. l and calculated volume changes are given in Table A.2. The water balance results are plotted in Figures A.2 and A.3 and are now discussed. Lower halocline - In Chapter 1 it was observed that the lower halocline oscillated seasonally, by as much as 7.5 m (Aug 1998 to March 1999), see Figure A.2. The halocline deepened considerably in the winters of 1997/98 and 1998/99, but less so for the winters of 1999/2000 and 2000/2001, which corresponds to the observed decrease in winter injection system flow. The lower halocline rebounded in the spring/summer. Changes in halocline position can be represented as volumetric fluxes across the lower halocline (Figure A.2). Negative volumetric fluxes (decrease in halocline position or upward fluxes of fluid across the halocline) correspond to the winter and high injection system flows. Negative volumetric fluxes were on average 0.66 of the corresponding injection system flow, with a range of 0.04 to 1.68. Positive fluxes (increase in halocline position or downward fluxes of fluid across the halocline) correspond to the summer and low injection system flows. Upper halocline - In Chapter 1 it was shown that the upper halocline (position of maximum dS/dz) has been rising with time, from 18.05 m (UWL) in July 1997 to 7.85 m (UWL) in March 2001, but that the rate of rise had slowed (Figure A.3). A seasonal signal was observed as the halocline rose more rapidly during the winter, and flattened during the summer, with even a slight drop in the summer of 2000. The gross volumetric fluxes at the upper halocline were mostly positive (Figure A.3), which is consistent with the observed rise of the halocline. Periods of large positive gross volumetric fluxes corresponded to high winter injection system flows, which caused the intermediate layer to fill and lift the halocline. The 164 "net" volumetric flux at the upper halocline is defined as the gross volumetric fluxes less the injection system discharge, and represents fluxes across the upper halocline from sources other than direct injection system discharge. The net volumetric flux was positive during the high injection system flow winters of 1997/98 and 1998/99 and was calculated to be 53 % and 21 %, respectively, of the injection system discharge. Entrainment of upper layer fluid into the intermediate layer by the injection system plume is proposed as an explanation for this downward flux, as the timing was concurrent with high injection system flows. During the fall of 1998 and 1999, the net volumetric flux was negative, which implies the downward erosion of the upper halocline by upper layer mixing processes. Through the time series, periods of negative net volumetric fluxes became more common, which indicates the increased importance of downward entrainment by upper layer mixing. This observation supports the Chapter 1 hypothesis of a "equilibrium upper layer depth", achieved when downward erosion of the halocline by upper layer mixing matches the rate of rise of the halocline due to ARD inflow and plume entrainment of the halocline. Relationship to Injection System Flow - In Chapter 1 the intermediate layer volume was found to expand at 1.59QA (r = 0.89) times the total injection system flow (QA). Similar relationships for the volumetric fluxes across the upper and lower halocline were calculated and shown in Figure A.4. The volumetric flux across the upper halocline was 1.26<2A (r = 0.87). The net volumetric flux across the upper halocline and volumetric flux across the lower halocline was estimated to be 0.272A (r = 0.44) and -0.42<2A (r = - 0.51) respectively, although the correlations were weaker. A better estimate of volumetric flux across the upper halocline is 1.26QA - QA = 0.26QA and lower halocline is l.59QA - 1.26QA = 0.33QA, as the relationships with the best correlations are used. These results are consistent with the observations that as the injection system flow increases the upper halocline rises (at a rate greater than that attributed to injection system flow alone) and the lower halocline descends. A.4 Plume Entrainment The expansion of the intermediate layer at a rate greater than that attributed to the injection system flow has been explained in terms of plume entrainment at the upper halocline and stirring of the intermediate layer and direct entrainment at the diffusers for the 165 lower halocline (Chapter 1). The theoretical plume entrainment of the upper halocline is now explored to see if the observed rate of 0.26<2A is reasonable. The phenomenon of plume entrainment was introduced in Chapter 1 and explored theoretically in Chapter 3. Plume entrainment is now quantified from energy arguments by equating the power available in the rising plume to the change in potential energy caused by entrainment. The power (P) available for mixing is obtained by integrating for power across the plume using expressions for the plume velocity (w) distribution and the maximum velocity (wm) and velocity half width (bw) defined in Chapter 3, P = \Ppw3dA ( A " 1 } where w=wm exp dA = 27lxdx which gives P = PPBH (A.2) where pp is the plume density, H is the thickness of the intermediate layer and plume maximum height of rise, and B is the initial buoyancy flux defined in Chapter 3. The rate of change of potential energy of the system is set equal to a fraction (Ck) of the power available in the rising plume. dPE= (A.3) dt The change in potential energy results from mixing of upper layer fluid uniformly through the intermediate layer, dPE ^{p2-Pl)gHQE (A.4) dt 2 where pi is the upper layer density, pi is the intermediate layer density and QE is the entrainment rate. An expression for the plume entrainment is derived by substituting Equations (A.2) and (A.4) into (A.3) and rearranging. Note that the plume density in 166 Equation (A.2) is taken as the density of the intermediate layer, as much dilution has occurred by the time the plume reaches the upper interface. Ck is analogous to the efficiency of penetrative convection plumes in the mixed layer of a lake eroding the thermocline. For penetrative convection Fischer et al. (1979) suggest Q = 0.13, whereas Ward et al. (1990) found Q - 0.20 for meromictic Mahoney Lake, British Columbia. Representative densities are 1003 kg m 3 , 1020 kg m 3 and 1000.7 for the upper (/?/), intermediate (pi) layers and injection system water (po), respectively. The upper layer and injection system waters where assumed to be at 8°C, which is an average fall/winter temperature. Substituting these values into (5) yields a calculated entrainment rate Qe=0.3Qd based on Ck - 0.13 (or Qe=0A5Qd if Ck = 0.20 is used). Thus the entrainment rate (for the plume impinging against the upper halocline) from energy considerations is in reasonable agreement with the observed volumetric change at the upper halocline. In the analysis of plume entrainment we assumed a homogeneous intermediate layer and this is an over simplification. In reality a density gradient exists through the intermediate layer, resulting from the "filling box" mechanism (Chapter 3). This reduces plume entrainment by slowing plume velocities and at lower ARD flows limiting the plume height of rise. Indeed, during low injection system flows the plume was observed by dye tracer experiments to be trapped at depth by the intermediate layer density gradient, and for higher flows to rise to the upper halocline before spreading horizontally (Muggli et al. 2000), see discussion Chapter 1. QE=2Ck (Pz-Po) (P2-P1) QA (A.5) 167 Table A. l : Halocline depths and changes in position. CTD Water Lower Halocline Upper Halocline Date level Depth Depth dz Depth Depth dz (m UWL) (m WL) (mUWL) (m) (m WL) (m UWL) (m) 24-Jul-97 10.06 207.69 217.75 7.45 17.51 28-Aug-97 9.75 207.57 217.32 0.42 7.91 17.66 -0.15 24-Sep-97 8.84 207.13 215.97 1.35 8.21 17.05 0.62 24-Nov-97 7.01 214.79 221.80 -5.83 9.18 16.19 0.86 3-Apr-98 2.59 219.27 221.86 -0.05 10.69 13.28 2.91 14-Aug-98 2.44 215.46 217.90 3.96 10.26 12.70 0.59 3-Dec-98 1.37 219.30 220.67 -2.78 11.33 12.70 0.00 3-Mar-99 -0.46 225.71 225.25 -4.58 10.69 10.23 2.47 9-Jun-99 0.36 221.70 222.06 3.19 8.95 9.31 0.92 22-Sep-99 0.77 218.95 219.72 2.33 8.43 9.20 0.11 6-Dec-99 -0.03 221.55 221.52 -1.80 8.80 8.77 0.43 31-Jan-00 -0.25 220.98 220.73 0.80 8.36 8.11 0.65 13-Mar-00 -0.28 220.51 220.23 0.50 8.14 7.86 0.26 5-Jun-00 -0.04 219.00 218.96 1.27 7.62 7.58 0.28 6-Jul-OO 0.15 218.39 218.54 0.42 7.41 7.56 0.02 29-Aug-OO 0.36 218.20 218.56 -0.02 7.58 7.94 -0.39 29-Nov-OO 0.03 219.86 219.89 -1.32 7.59 7.62 0.32 15-Mar-01 -0.30 220.79 220.49 -0.60 7.94 7.64 -0.02 Table A.2: Interface volume changes. CTD Date Lower Interface Vol. Change Upper Interface Vol. Change ARD Flow Vol.*1 Upper Interface Net Vol. Change (106m3) (106m3) (106m3) (106m3) 24-Jul-97 28-Aug-97 0.151 -0.247 na na 24-Sep-97 0.493 0.977 na na 24-Nov-97 -2.18 1.35 na na 3-Apr-98 -0.087 4.67 3.23 1.44 14-Aug-98 1.43 0.957 0.420 0.537 3-Dec-98 -1.03 -0.0111 1.00 -1.01 3-Mar-99 -1.60 4.04 3.56 0.481 9-Jun-99 1.10 1.52 1.11 0.415 22-Sep-99 0.838 0.194 0.231 -0.0369 6-Dec-99 -0.655 0.703 1.24 -0.533 31-Jan-00 0.281 1.08 1.08 0.00753 13-Mar-00 0.180 0.425 0.544 -0.19 5-Jun-00 0.464 0.471 0.835 -0.364 6-Jul-00 0.155 0.0321 0.0592 -0.0271 29-Aug-00 -0.00572 -0.639 0.0025 -0.641 29-Nov-OO -0.487 0.533 0.305 0.227 15-Mar-01 -0.221 -0.0358 0.557 -0.593 na - injection system flow not measured at this time. *' includes estimate of conveyer tunnel flow as 4% of total injection system flow. 169 Wind Q(ARD) (1) Injection system flow and plume entrainment (2) Wind mixing and penetrative convection deepen the raise the interface by dz=[Q(Plume)+Q(ARD)]/area interface by dz=[Q(Wind Mixing+Q(PC)]/area Figure A . l . Schematic diagram of mixing at the upper interface (Fisher, 1999). 170 (a) 2101 1 1 1 1 1 Jcn98 Jcn99 JcnOO JcnOl Date Figure A .2 . (a) Depth of lower halocline average salinity [—] and shoulders ["], (b) injection system total flow ['"] and average flow for period between field trips [—], (c) lower halocline volumetric flux [—]. Positive volumetric flux values represent a rising halocline and downward net volumetric flux of fluid across the interface. 171 Figure A .3. (a) Depth of upper halocline average salinity [—] and shoulders ["'], (b) injection system total flow [ '] and average flow for period between field trips [—], (c) gross upper halocline volumetric flux [—] and net volumetric flux [""]. Positive volumetric flux values represent a rising halocline and downward net volumetric flux of fluid across the interface. 172 x 10 (a) 5 4 *T 3 T3 m x I 1 i o -2 x 10 Q, UI-ARD' = 0 . 2 6 Q A R D + - 3 8 7 1 , r=0.44 I 0 x 10 0.5 1.5 2 Q A ( mJ d 3 ^ - 1 \ 2.5 3.5 (c) x 10 x 10 Figure A .4 . Relationships between injection system flow (QARD)^ (a) upper interface volumetric flux (Qui), (b) upper interface "net" volumetric flux (QUI-ARD) and (c) lower interface volumetric flux (Qu)- Relationships based on regression by least-squares-fit on y variable with correlation coefficients given. 173 APPENDIX B Island Copper Meteorology B.l Introduction The meteorological station at the Island Copper Mine was deployed from December 1999 to November 2000 to coincide with thermistor chain measurements. The Island Copper meteorological station was located onshore near the south injection system at 50°36' 127°29' and at an elevation of ~8 m above mean sea level (Figure 3.1). Appendix E details the meteorological station deployment, instrumentation and calibration. In this appendix the background to weather on the south coast of British Columbia is given and the Island Copper meteorological data summarized. In particular, the wind data are carefully analyzed to predict the likelihood of extreme winds. The probability of extreme wind events is estimated by wind frequency analysis using the long-term Port Hardy Airport (Appendix E) wind record. The coherence and relationship between Port Hardy Airport and Island Copper winds are examined by comparing concurrent wind measurements during 2000. The resulting scaling relationship can then be used to estimate extreme winds at Island Copper from those observed at Port Hardy Airport. These results are used to predict destratification of the pit lake by upwelling as discussed in Chapter 5. B.2 Regional Weather The pit lake is situated adjacent to Rupert Inlet, which is a coastal fiord reaching west via Quatsino Sound to west coast of Vancouver Island and the Pacific Ocean (Figure B.l). Environment Canada operates a permanent weather station at Port Hardy Airport (Station ID 1026270, location 50°41' 127°22' and elevation of 22 m), which is 12 km north-east of Island Copper and on the east coast of Vancouver Island (Figure B.l). The Port Hardy Airport 174 weather station is referred to as Port Hardy in this appendix. Climate normals at Port Hardy are for daily mean temperatures ranging from 3.0 °C in the winter to 13.9 °C in the summer, mean annual precipitation of 1871 mm and average winds of 3.3 m s"1 from the southeast (Environment Canada, 1993). Cloud cover is high both in summer and winter averaging between 75% and 90% along the B.C. coast (Thomson, 1981). The prevailing wind patterns along the west coast of Canada are south-east to south-west from October through April from anti-clockwise air flow around the Aleutian Low, and north-west from May through September due to clockwise flow around the North Pacific High. The greatest departures from the prevailing winds originate from eastward moving high and low synoptic systems, which result in a roughly three day cyclic variation in wind speed. In particular, cyclonic lows intensify into storms that invade the coast between fall and late spring. These storms characteristically extend hundreds of kilometres and last for days. These storms cause strong and rapidly changing winds and heavy precipitation. The anti-clockwise air flow typically results in south-east to south-west winds that reverse in direction as the tail of the cyclonic depression passes (Thomson, 1981). These storms interact with the coast mountains to cause strong south-easterly winds up the east coast of Vancouver Island (Figure B.l). Winds at the Island Copper mine are modified locally by the Vancouver Island Range, Rupert Inlet and Holberg Inlets (Figure B.l). The Vancouver Island Range bisects north Vancouver Island in a west-east direction. Winds are channeled through Quatsino Sound, which extends 40 km in a south-west direction to the Pacific Ocean, and via Holberg Sound, which extends west-north-west for 40 km. Easterly winds from Queen Charlotte Strait on the east coast of Vancouver Island funnel into these waterways via an 11 km narrowing of the island, which has relatively low relief (maximum elevation -100 m). Island Copper is sheltered from the south by the Vancouver Island Range and from the north by hills exceeding 300 m. Note wind directions follow meteorology convention and are the direction in which the wind comes from. 175 B.3 Island Copper Weather The Island Copper meteorological 15 min data was converted to hourly values by excluding the most extreme value (of the five measurements over an hour - 0:00, 0:15, 0:30, 0:45, 1:00) and averaging the remaining four measurements. Hourly temperatures ranged from 27.9 °C to -3.0 °C (Figure B.2), whereas average daily temperatures ranged from 19.3 °C to -0.8 °C. The hourly relative humidity ranged from 100% to 62% with an average of 87 %. A seasonal trend of increased relative humidity in the winter was observed. Short-wave solar radiation was seasonal, peaking in summer with a maximum hourly value of 1077 W m" 2 . Rainfall was measured from 9 December 1999 to 29 August 2000 only. During the period of observation it rained on 51% of the days. The peak hourly rainfall was 6.5 mm hr"1 and maximum daily accumulation was 35.6 mm. The most dominant wind at Island Copper was an east-south-east (ESE) wind (100-110°), followed by a westerly (-260°) (Figure B.2 and Figure B.3). Both winds are consistent with the alignment of geographical features discussed in Section B.3. These winds will have a significant impact on the lake as they blow down the axis of the lake, which has a bearing of -115° (-295°). For comparison, the most dominant wind at Port Hardy was also from the ESE (110-120°) (see Figure B.4). B.4 Converting Island Copper Wind Speed Measurements of wind speed at Island Copper were made by an anemometer 4.4 m off the ground. These were converted to wind speed at a height of 10 m above the ground for consistency with the standard atmospheric observations and the limnological literature. To convert the wind speed a "law of the wall" log profile and a statically-neutral surface layer were assumed. The wind speed increases logarithmically with height and the shape of the profile is dependent on the surface roughness (z0) (Stull, 2000): u. I M (B-l) w(z) = —In k If a wind speed (JJ\) is known at height (zi) then a new wind speed (U2) at any other height (Z2) can be calculated: 176 u 2 = U i l A h i i A ( B - 2 ) The surface roughness was selected as 0.25 m, which is classified as "rough" for a "landscape of crops of variable height" and "scattered obstacles" (Stull, 2000). The correction from a wind speed measured at 4.4 m (C/4.4) to wind at 10 m (U\Q) from Equation (B.2)is Uio= 1.29£/4.4. B.5 Port Hardy Maximum Annual Wind Frequency Analysis. The wind speed to cause upwelling has been calculated for the future from modeling of the pit lake (Chapter 5). The likelihood of upwelling occurring is the probability of obtaining the upwelling wind speed in that year. The Port Hardy weather station was used to assess the probability of obtaining the upwelling wind speed, as a long record was available. The Port Hardy wind speed and direction are measured at the standard height of 10 m. The wind measurements are two minute averages taken every hour. Hourly average wind values were estimated by averaging two consecutive measurements (to give an estimate of the average over the hour as was done with pit lake data). The annual maximum hourly wind speed was found for each year from 1953 to 2000. The annual maximum hourly wind often persisted for longer than an hour. Four years (1965, 1970, 1982 and 1997) had two storms where the annual maximum hourly wind was achieved and only one record was considered in the frequency analysis. The annual maximum hourly wind speeds were all ESE winds confined to a narrow band between 90-120° with a median of 110°. The annual maximum hourly winds mostly occurred during the winter, as 79% occurred in November to January, whereas 96% occurred in the six months from October to March. These winds are consistent with the general winter weather patterns discussed in Section B.2. The annual maximum hourly wind speeds were ranked in descending order with the rank (m) equal to one given to the largest value. Probabilities (P) were determined using the Weibull plotting scheme, P = m/(n+l), where n is the number of years, and the return period (Tr) is the reciprocal of the probability (Viessman and Lewis, 1996). Return periods up to 177 twice the record length can be estimated, so a 100 year event can be calculated from the Port Hardy record of 48 years. The frequency distribution is shown in Figure B.5 and has been fitted with the curve described by: The estimated 100 year return period annual maximum hourly wind speed at Port Hardy is 22.7 m s"1. This is only marginally greater than the observed 49 year return period maximum hourly wind speed of 22.4 m s"1 due to the shape of the frequency distribution (Figure B.5). B.6 Comparison of Island Copper and Port Hardy Wind To apply the Port Hardy maximum annual wind frequency distribution to the Island Copper site the wind patterns must be reasonably correlated. In this section the wind correlation between the two sites is examined, then a relationship is found to scale the Port Hardy maximum annual wind frequency distribution to winds at Island Copper. Comparisons are made between concurrent measurements at Island Copper and Port Hardy for 9 December 1999 to 30 November 2000. Comparison Techniques - In the comparison analysis linear regression using least squares fits to the data were used. Linear regression assumes independent and dependent variables. Two least square fits can be made by taking either variable to be the independent. The average of these two fits can be used as the predictive equation (Pieters et al, 1999). The goodness of fit was measured by the standard error. The similarity of data sets is measured in terms of the correlation coefficient, r, (B.3) X(*-*Xy-y) (B.4) r = where x and y are the means of samples x and y. Probability of direction was determined by binning the data into 10° increments. The number of occurrences in each bin divided by the total number of events gives the probability (xl00=%) of that wind direction. 178 The periodicity of the wind was investigated using spectral analysis. The spectral density is the measure of the variability at a given frequency, and the area under a spectra curve is the variance. The units of spectral density are (m s"1)2 (c d"1) ~l. The Fourier Series used to construct the spectra is the sum of sine waves each with an amplitude, phase and frequency. The coherence is the fraction of the amplitude that is related. When the coherence = 0.7, half of the variance at that frequency is related at the two stations. The phase is the difference in phase between the two sinusoids representing the variations in the two signals at that frequency. When the phase = 0°, the two signals are in phase (Pieters et al, 1999). Correlation between Island Copper and Port Hardy Wind - The Island Copper and Port Hardy wind record is compared for 1-14 January 2000 in Figure B.6. Wind components in the v and u are plotted in Figure B.7. The wind storm on 3rd January 2000 was the largest at Port Hardy and fifth largest at Island Copper in terms of wind speed, during the observation period. For this period, the correlation is reasonable (r = 0.61) for wind speed, better (r = 0.78) for u, but poor (r = 0.22) for v. Similarly, for the full observation period the wind speeds are reasonably correlated (r = 0.48), but better correlated in u (r - 0.65) and less well correlated in v (r = 0.37), see Figures B.8, B.9 and B.10. Winds less than 1 m s"1 were discarded because these are flukier and the Island Copper anemometer has a theshold speed of 1 m s"1. Using all winds the correlation coefficients were not too different at 0.58, 0.60, 0.34 for wind speed, u and v, respectively. Winds are better correlated in the u component because easterly winds dominate at both Island Copper and Port Hardy (Figures B.ll). The most probable wind direction at Island Copper is 100-110° and 50% of winds occurred in the quadrant from 55-145°. The most probable wind direction at Port Hardy was 110-120° and 43% of winds occurred in the quadrant from 75-165°. These winds occur at the same time as seen in the time series and result in good correlation in the u component. The -10° difference in direction is likely due to local topographical differences. The second most common wind at Island Copper is a westerly wind centered at 260°, and 28% of all winds occurred in the 225-315° quadrant. The westerly winds are not significant at Port Hardy probably because the site is sheltered 179 from the west by the Island's interior mountains. The second most common wind at Port Hardy is 330-340°, and 17% of all winds occur in the 295-25° quadrant. The spectral density of Port Hardy wind speeds is greater than for Island Copper at all 2 2 frequencies (Figure B. 12). In fact, the variance of Port Hardy winds is 6.34 m s" compared to 1.95 m2 s"2 at Island Copper. At both sites the diurnal (1 cycle day"1) wind pattern contributes strongly to the variance. The second and third highest peaks in spectral density for both sites are at periods of 15-60 days and 6-15 days, respectively. The Port Hardy spectral density has a further peak at 2-3 days. The diurnal wind exhibits very strong coherence and the phase is near zero. Elsewhere, the coherence varies from 0.7 at low frequencies to 0.4 at 2 cycles day"1. Phase varies but is centered about zero. Wind Scaling Relationship - The comparison so far has been based on all hourly winds and a reasonable correlation between Island Copper and Port Hardy has been established. A scaling relationship could be based on all hourly winds and use the u component wind as this captures the dominant wind and has the best correlation. Averaging the two regression curves gives Uic = 0.35 + 0.66 Upn, where U/c and UPH are the Island Copper and Port Hardy winds, respectively. More appropriate, would be to derive the wind scaling relationship from maximum winds, since it is the maximum annual wind speed frequency distribution that we want to reproduce at Island Copper, and these may scale differently than all hourly winds. Hourly wind speeds exceeding three times the standard deviation from the mean, which sum to less than 0.1% of the record, were found for both Island Copper and Port Hardy. The maximum wind in any wind event was identified and other high winds in the event were discarded. These maximum winds were matched with the maximum wind occurring for the concurrent wind event at the other station. An estimate of the wind scaling factor is the ratio of the maximum event wind at Island Copper to that at Port Hardy. Note, maximum event wind speeds at the two stations rarely occurred at the same time; e.g. ESE winds tended to peak at Island Copper an hour earlier than for Port Hardy. Thus, it was important to use the maximum event wind speeds, rather than simultaneous wind speeds as these erroneously reduce the wind scaling factor. 180 Pairs of maximum event winds are plotted in Figure B.13. The wind scaling factor ranges from 0.407 to 1.696 with an average of 0.887. The ratio is sensitive to wind direction (at Island Copper). ESE winds that comprise 60% of the maximum event winds tend towards a lower average ratio of 0.68. Winds from the north, west and south quadrants tend towards a higher average ratio of 1.18. The five largest Island Copper maximum event winds are distinct from the rest of the events. These five largest Island Copper maximum winds correspond to Port Hardy events with rank=[14,11,9,2,1] and are all ESE winds. The average ratio of these events is 0.825, which is similar to the average of all maximum wind events. The wind scaling relationship of Uic = 0.825 UPH was chosen to convert the Port Hardy maximum annual wind frequency distribution to one for Island Copper. The scaling factor is based on the ratio of the five largest Island Copper winds compared to concurrent measurements at Port Hardy and is similar to the scaling factor of all maximum events. These five largest wind speed events at Island Copper were produced by ESE winds and this same wind direction was responsible for all 48 annual maximum wind speeds observed at Port Hardy. The use of this wind relationship is further supported by the good correlation between u wind components at Island Copper and Port Hardy. B.7 Conclusions Wind at the Island Copper mine was consistent with synoptic weather patterns modified by local geography. Island Copper and Port Hardy Airport winds were found to be reasonably correlated particularly in the easterly quadrant. A scaling relationship of Uic = 0.825 UPH was determined from extreme winds for the year of simultaneous record. The scaling relationship can be used to convert the Port Hardy maximum annual wind frequency distribution (Equation B.3) to one for Island Copper. 181 122UW 124°W M n92l. M 08cU T3 C ca o JC c o C H ca 6 e o "5b <u CH 1 T3 SH > O o c ca > c SH <D C H C H O U T3 C ca ta) c '$ o JC OO JO e 3 U JC 0 0 '•4—* •c m -*—> 0 0 <U J = "3 o 00 pq 60 -a J = t3 o c u JC o oo JC o J C 60 c ca sc CO > 3 O o c ca > JC T3 C 3 O 00 o c '33 ca 3 a "c I-H 60 S-. <u j=> "o K t-H tJ CO C H 3 P< O <u _> ca a c SH <D C H Cu O U O0 ca <u i oo 183 Figure B.3. Island Copper hourly wind presented as v (north-south component) versus u (west-east component), u and v are positive for winds blowing from the north and east directions. 184 10 h-rt.. .v * ***., ** # *»• • «*• * • • . * • . % • * * . * * • ^ » » • • • t&frf * *•* * •* • *+ ?rMB$&r**: ****** • * * + + ** * * * * • S^ypSSsA. **. ryv/fih ...... ***** **•• ^ > ^ * * •++ + * * * — • . . , * * * - * . +> **• * . * . .* ? * . * • • -10 h--15 -10 -5 0 5 u wind component (West-East, m/s) 10 15 Figure B.4. Port Hardy hourly wind presented as v (north-south component) versus u (west-east component), u and v are positive for winds blowing from the north and east directions 185 Maximum Annual Wind Speed at 10m (m/s) Figure B.5. Frequency distribution of Port Hardy annual maximum hourly wind speeds for 1953 to 2000. 186 Figure B.6. (a) Hourly wind speed and (b) direction at Island Copper [—] compared to Port Hardy [-] for 1 - 14 January 2000. 187 Figure B.7. (a) u wind component and (b) v wind component at Island Copper [—] compared to Port Hardy [-] for 1-14 January 2000. u and v are positive for winds blowing from the north and east directions. 188 15 •10 CD CU Q . CO - a £ 5 x=a+b.y where a=1.786, b=0.653 and std error =2.106 y=a+b.x where a=1.936, b=Q,35 and std error =1.542 r=0.478 . * + ++ + 5 10 Island Copper wind speed (rn/s) 15 Figure B.8. Correlation between Port Hardy and Island Copper hourly wind speed for 9 December 1999 to 30 November 2000. Winds < 1 m s"1 have been excluded. Lines and equations represent "least-square-fit" regressions, taking Island Copper and Port Hardy wind alternately, as the dependent and independent variables. 189 15 I 1 1 1 1 y=a+b.x where a=0.084, b=0.525 and std error =2.269 x=a+b.y where a=0.611, b=0.801 and std error =2.803 r=0.649 10 h -10 I-_15 I i i i i i -15 -10 -5 0 5 10 15 Port Hardy u component (m/s) Figure B.9. Correlation between Port Hardy and Island Copper hourly u wind component for 9 December 1999 to 30 November 2000. Winds < 1 m s"1 have been excluded. Lines and equations represent "least-square-fit" regressions, taking Island Copper and Port Hardy wind alternately, as the dependent and independent variables. 190 -10 h -15' -15 -10 -5 0 5 Port Hardy v component (m/s) 10 15 Figure B.10. Correlation between Port Hardy and Island Copper hourly v wind component for 9 December 1999 to 30 November 2000. Winds < 1 m s"1 have been excluded. Lines and equations represent "least-square-fit" regressions, taking Island Copper and Port Hardy wind alternately, as the dependent and independent variables. 191 11 | 1 1 1 1 1 1 r 0 50 100 150 200 250 300 350 Direction (°) Figure B. l l . Probable wind direction at Island Copper [—] compared to Port Hardy [--] for 9 December 1999 to 30 November 2000. Winds < 1 m s"1 have been excluded. 192 (a) 351 1 1 1 1 1 1 1 1 1 3 0 - I 2 5 - 1 > s Ii to ii § 2 0 -Q h ~& ! I i Frequency (cycles/day) Figure B.12. (a) Wind spectrum (6 blocks of 60 days, Vi block overlap, Planning window, hourly data) for Island Copper [—] and Port Hardy [--] for 9 December 1999 to 30 November 2000. (b) Wind coherence [—] and phase [--]. 193 • East Quadrant X South Quadrant A West Quadrant Port Hardy wind (m/s) Figure B.13. Comparison of Island Copper and Port Hardy extreme wind speeds (greater than three standard deviations from the mean) for 9 December 1999 to 30 November 2000. 194 APPENDIX C Conductivity Temperature Depth (CTD) Measurements C.l Introduction Conductivity-temperature-depth profiling has been carried out on 18 occasions: July, August, September and November 1997; April, August and December 1998; March, June, September and December 1999; February, March, June, July, August and November 2000; and March 2001. The primary conductivity-temperature-depth instrument (CTD) was an Ocean Sensors, Inc., Model OS200 (Figure F.10). Measurements were supplemented on occasion with Sea-Bird Electronics SBE19 and SBE25 to verify the performance the OS200 and to measure additional parameters such as dissolved oxygen, transmissivity and florescence. This Appendix details the OS200 performance and overviews the calibration procedures. C.2 OS200 Performance CTDs measure conductivity, temperature and depth (pressure) from which salinity and density can be calculated. Salinity was calculated using the Practical Salinity Scale 1978 (PSS78) from conductivity, temperature and pressure (Lewis and Perkin 1978, Perkin and Lewis 1980). The OS200 was operated in an averaging mode and programmed to average 99 raw measurements every -0.5 s to provide stable and accurate measurements. Corrections to calibrated salinity were made using field bottle samples measured with a Guildline Model 8400 Autosal. The resulting salinity accuracy was 0.03 and resolution was 0.001. An additional salinity correction to the PSS78 salinity was required to account for non-seawater sourced ions and is described in Chapter 2. The temperature accuracy and resolution were 0.02 °C and 0.001 °C, respectively. The uncertainty of depth measurements was 0.23 m, 195 which was the standard error of OS200 measurements (corrected for temperature and depth effects) compared to simultaneous measurements made with the SBE19 and thermistors. The depth resolution of 0.15 m was achieved in interface areas by slowing the CTD descent rate. C.3 OS200 Operation OS200 profiles were primarily taken at the West Station buoy - located near the centre and deepest part of the lake (Table C.l , Figure C.l). The OS200 was lowered and raised through the water column by cable and hydraulic winch powered by a petrol generator. The OS200 was attached by nylon ties, electrical tape and shackle to the cable. A thirty-five pound lead ball weight was attached to the cable approximately one metre below the OS200 to ensure the CTD fell vertically through the water column. Depth was calculated by subtracting the atmospheric pressure from the recorded pressure, and dividing by the integrated local density and gravity. Atmospheric pressure was determined by holding the CTD above the water or at a known depth (i.e. 0.1 m). CTD profiles were also performed at the East Station and over the North Injection System and South Injection System diffusers (Table C.l , Figure C.l). These measurements showed the lake was horizontally homogeneous. For this reason, only profiles at the west station are reported in this thesis. In addition, CTD measurements with the instrument held at the interface depth (known as a thermocline hold) were made to identify internal waves. These measurements were not successful and were supplanted by thermistors (Appendix D). The CTD data excluded measurements taken before the instrument was submerged and during the up-cast. Calibration corrections for conductivity, temperature and pressure were applied. Bad data points were removed with a median filter, which replaced outliers with the current median value. Outliers were defined as values exceeding three standard deviations from the current median value, which was the median of a nine data point window. The additional salinity correction for non-seawater sourced ions (Chapter 2) was then added. 196 C.4 OS200 Calibration The OS200 was calibrated before and after each trip by comparing OS200 measurements to more accurate laboratory instruments to determine correction equations. Wilton (1998) calibrated salinity for trips 1-5 (Table C.l), from OS200 and Autosal calculated salinity. However, the salinity calibration was found to be temperature dependent. So calibrations for trips 1-5 were revised into conductivity calibrations and subsequently calibrations done for conductivity instead of salinity. Calibration of the measurable conductivity and temperature, followed by the calculation of salinity, was more rigorous than calibration of salinity, which is a secondary parameter that is dependent on conductivity and temperature. Trips 1-10 (Table C.l) calibrations for salinity and conductivity, were performed in a dewar using dilutions of Burrard Inlet water by comparing OS200 to Autosal measurements in a method described by Wilton (1998). For subsequent trips, calibration of conductivity was performed concurrently with temperature in a temperature controlled and mechanically stirred water bath. A range of conductivity was achieved by using -12 and -24 salinity water and changing the temperature. Laboratory measurements of temperature were made with a Guildline 9540 Platinum Resistance thermometer (accuracy 0.01°C) and a Dymec Model 2801A Quartz thermometer, whereas conductivity was measured with the Guideline Model 8400 Autosal (accuracy 0.005 in S). Conductivity calibrations varied considerably - by 0.2 mS/cm at 30 mS/cm over three years. Of more concern was that pre- and post calibration were often different. The difference in salinity between the OS200 and in situ bottle samples measured with the Autosal increased with layer depth and had an average error of 0.02 and average scatter of ±0.03. This level of accuracy was marginal considering the salinity of the lower layer was changing at -0.05 every four-months. So the OS200 salinity was further improved by calibrating against the Autosal salinity measured from in situ bottle samples. In situ bottle samples were taken for each pit lake layer for most trips. The OS200 was also calibrated for depth (pressure). Depth was corrected as a function of temperature based on experiments in the water-bath, warm and cold 3 m deep 197 swimming pools and in the pit lake, where depth was simultaneous measured with the OS200 and tape-measure. The depth correction for depth varies by 1 m over the 20 °C range of water temperatures expected in the pit lake. Furthermore, depth was corrected as a function of depth from comparisons between OS200 pit lake profiles and thermistor temperatures (at a known depth) and SBE19 profiles. The depth correction for depth was 2.14 %, which was 4.28 m at 200 m depth. The standard error of OS200 measurements (corrected for temperature and depth effects) compared to simultaneous measurements made with the SBE19 and thermistors was 0.22 and 0.23 for the upper and lower interfaces, respectively. Table C.l . Summary of CTD profiles using OS200 unless noted otherwise. Trip Date Cast Name Location Type Other Operator 1 23-Jul-97 WC1 West Station Deep Cast MW 1 23-Jul-97 WH1 West Station Thermocline Hold MW 1 23-Jul-97 NCI North Injection System Deep Cast MW 1 23-Jul-97 NH1 North Injection System Thermocline Hold MW 1 23-Jul-97 EC1 East Station Deep Cast MW 1 23-Jul-97 EH1 East Station Thermocline Hold MW 1 23-M-97 EC2 East Station 10m Very Fast Drops MW 1 23-Jul-97 SCI South Injection System Deep Cast MW 1 23-Jul-97 SHI South Injection System Thermocline Hold MW 1 24-Jul-97 WC3 West Station Deep Cast MW 1 24-Jul-97 WH3 West Station Thermocline Hold MW 1 24-Jul-97 FWC3 west of West Station Deep Cast MW 1 24-Jul-97 VWC3 far west of West Station Deep Cast MW 1 24-M-97 NC3 North Injection System Deep Cast MW 1 24-M-97 ND3 North Injection System Thermocline Drift MW 1 24-M-97 ND4 North Injection System 12m Up Down Drift MW 1 24-Jul-97 ND5 North Injection System 12m Up Down Drift MW 1 24-Jul-97 SC3 South Injection System Deep Cast MW 2 28-Aug-97 WC11 West Station Deep Cast DO MW 2 28-Aug-97 E C U East Station Deep Cast DO MW 2 28-Aug-97 NC11 North Injection System Deep Cast MW 2 28-Aug-97 sen South Injection System Deep Cast MW 2 29-Aug-97 EH11 East Station Thermocline Hold MW 2 29-Aug-97 EH12 East Station Thermocline Hold MW 2 29-Aug-97 EC 12 East Station Deep Cast MW 2 29-Aug-97 NH11 North Injection System Thermocline Hold MW 2 29-Aug-97 NH12 North Injection System Thermocline Hold MW 3 24-Sep-97 WC21 West Station Deep Cast MW 3 24-Sep-97 WC22 West Station Deep Cast MW 3 24-Sep-97 EC21 East Station Deep Cast DO MW 3 24-Sep-97 SC21 South Injection System Deep Cast DO MW 3 24-Sep-97 SC22 South Injection System 20m Very Slow Cast DO MW 4 23-Nov-97 WC31 West Station Deep Cast MW 4 23-Nov-97 WH31 West Station Thermocline Hold MW 4 23-Nov-97 SC31 South Injection System Deep Cast MW 4 23-Nov-97 SH31 South Injection System Thermocline Hold MW 5 3-Apr-98 WC41 West Station Deep Cast MW 5 3-Apr-98 WC42 West Station Deep Cast MW 6 14-Aug-98 WC51 West Station Deep Cast DO*l TF 6 14-Aug-98 WC52 West Station Deep Cast TF 6 14-Aug-98 WC53 West Station Deep Cast TF 7 3-Dec-98 WC61 West Station Deep Cast TF 7 3-Dec-98 WC62 West Station Deep Cast TF 7 3-Dec-98 NH61 North Injection System Thermocline Hold TF 7 3-Dec-98 NH62 North Injection System Thermocline Hold TF 8 3-Mar-99 WC71 West Station Deep Cast DO*l TF 8 4-Mar-99 WC72 West Station Deep Cast DO*l TF 8 4-Mar-99 WC73 West Station Deep Cast TF 8 4-Mar-99 WC74 West Station Shallow Cast TF 8 4-Mar-99 SY71 Rescan - Time 3 Transit 1 Yoyo Transit TF 8 4-Mar-99 SY72 Rescan - Time 3 Transit 3 Yoyo Transit TF 199 Table C.l . Summary of CTD profiles using OS200 unless noted otherwise (cont.). Trip Date Cast Name Location Type Other Operator 8 4-Mar-99 SY73 Rescan - Time 4 Transit 1 Yoyo Transit TF 8 4-Mar-99 SY74 Rescan - Time 4 Transit 3 Yoyo Transit TF 9 9-Jun-99 EC81 East Station Deep Cast TF 9 9-Jun-99 WC81 West Station Deep Cast TF 9 9-Jun-99 WC82 West Station Deep Cast (Fast) TF 9 9-Jun-99 SC81 South Injection System Deep Cast TF 9 9-Jun-99 SH82 South Injection System Thermocline Hold TF 10 22-Sep-99 WC91 West Station Deep Cast TF 10 22-Sep-99 WC92 West Station Deep Cast TF 10 22-Sep-99 EC91 East Station Deep Cast TF 10 22-Sep-99 EC92 East Station Deep Cast TF 11 6-Dec-99 WC101 West Station Deep Cast TF 11 6-Dec-99 SBWC101 West Station Deep Cast SBE25 DJ 11 6-Dec-99 WC102 West Station Deep Cast TF 11 6-Dec-99 SBWC102 West Station Deep Cast SBE25 DJ 11 6-Dec-99 SC101 South Injection System Deep Cast TF 11 6-Dec-99 SBSC101 South Injection System Deep Cast SBE25 DJ 11 6-Dec-99 SH101 South Injection System Thermocline Hold TF 11 6-Dec-99 SBSH101 South Injection System Thermocline Hold SBE25 DJ 12 31-Jan-00 WC111 West Station Deep Cast TF 13 13-Mar-00 WC121 West Station Deep Cast TF 13 13-Mar-OO WC122 West Station Deep Cast TF 14 5-Jun-00 WC131 West Station Deep Cast TF 14 5-Jun-00 WC132 West Station Deep Cast TF 14 5-Jun-00 SC131 South Injection System Deep Cast TF 14 5-Jun-00 SCI 32 South Injection System Deep Cast TF 15 5-Jul-00 WC141 West Station 68m Cast TF 15 5-Jul-00 0705ICW1 West Station 68m Cast Idronaut BB 15 6-Jul-00 WC142 West Station Deep Cast TF 15 5-Jul-00 0705icw2 West Station 68m Cast Idronaut BB 15 5-Jul-00 NW141 North Pit Wall Temp TC Thermocline Hold TF 15 5-Jul-00 0705icnl North Injection System 73m Cast Idronaut BB 15 6-Jul-00 0706icpl South ARD POND 0.4m Cast Idronaut BB 16 29-Aug-00 WC151 West Station Deep Cast TF 16 30-Aug-00 SC151 South Injection System Deep Cast TF 17 29-Nov-OO WC161 West Station Deep Cast TF 17 29-Nov-OO WC162 West Station Deep Cast TF 18 15-Mar-01 WC171 West Station Deep Cast TF 18 15-Mar-01 SBWC171 West Station Deep Cast SBE19 TF 18 15-Mar-01 WC172 South Injection System Deep Cast TF 18 15-Mar-01 SBWC172 South Injection System Deep Cast SBE19 TF 18 15-Mar-01 SC172 South Injection System Deep Cast TF 18 15-Mar-01 SBSC172 South Injection System Deep Cast SBE19 TF MW - Mike Wilton, TF - Tim Fisher, BB - Betram Boehrer, DJ - Dave Jones, DO - dissolved oxygen, Idronaut - Idronaut Ocean Seven 316 (C, T, P, DO, pH, Eh, Turbidity, fluorescence), SBE25 - Sea-bird Electronics Inc. SBE25 (C, T, P, DO), SBE19 - Sea-bird Electronics Inc. SBE19 (C, T, P, fluorescence, transmissivity), Rescan -transects of lake (Rescan 1999b), *1 Unsuccessful. 200 Figure C.l . Island Copper pit lake bathymetry showing locations of CTD casts at West Station, East Station, South Injection System (SIS) diffuser and North Injection System (NIS) diffuser. Horizontal datum is WGS-84 and depths are referenced to the ultimate water level (UWL), which is 2.44 m above mean sea level (adapted from Muggli et al. 2000). 201 APPENDIX D Thermistor Chain Mooring Design and Thermistor Specification and Calibration D . l Introduction The field investigation was intensified beyond the water quality and periodic CTD measurements by deploying two thermistor chains and a meteorological station from December 1999 to November 2000. There were two mooring sites, known as the central mooring thermistor chain (CMTC) and south injection system thermistor chain (SISTC). This appendix provides details on the thermistor chain mooring design and thermistor layout, specifications and calibrations. The background and motivation for thermistor chains are given in the research proposal Fisher and Lawrence (1999) Island Copper pit lake physical processes - proposal for the installation of thermistor chains. In summary, thermistor chains were used to continuously measure changes to the pit lake structure to provide a better understanding of physical process than could be gained from periodic CTD profiles. Concurrent measurements of meteorological conditions allow for the environmental forcing to be quantified. The purpose of the deployment was to improve the understanding of the pit lake in several key areas: upper interface equilibrium depth; anoxia of the intermediate layer; intermediate layer plume dynamics and lower interface mixing dynamics. The thermistor measurements also provide baseline data for the development of a predictive model for the pit lake. 202 D.2 Thermistor Chain Design The CMTC was located in the deepest part of the lake, which was in near the middle of the lake (Figure 3.1) and close to the West Station buoy. The purpose of the CMTC was to measure the evolving thermal stratification of the water column and thermistors were spread from the surface to 250 m. A two point mooring system was designed to position the CMTC near the deepest part of the lake. The CMTC was fixed in place with surface lines to the existing BHP's West Station buoy and to a new mooring buoy and anchor system. The design drawing for the mooring is shown in Figure D.l (note there were minor changes to dimensions and materials from the design drawings during construction The top 20 m of the thermistor chain was 1/2' Tenex rope to which thermistors were permanently attached. The remaining thermistor chain was 3/16' wire rope (impregnated and coated with plastic to give VA" overall) and thermistors were attached via a mounting block or zap straps and electrical tape. The bottom of the thermistor chain was weighted to keep the thermistor chain vertical. Swivels were used at the both the top and bottom of the thermistor chain to prevent twisting. The SISTC was located over the south injection system (SIS) diffuser. The purpose of this thermistor chain was to assess the impact on the south injection system (SIS) of the upper halocline and identify internal waves not seen at the central mooring because of its position close to the centre of the lake. The SISTC was 100 m in length. The position of the SIS diffuser was determined with a Differential Global Positioning System by pulling up on the BHP's existing SIS buoy, which was tied to the end of the diffuser, and sighting along the line of the outfall pipe visible onshore. Fixing the SISTC position over the SIS diffuser was critical so a three point mooring system was designed. The three point mooring system was provided by surface lines to: (1) BHP's existing SIS buoy to the north, which was pulled out of the way by connecting it to the BHP's North Injection System buoy; (2) a temporary mooring buoy and anchor system to the west; and (3) a ground anchor on the shore to the south-east. The SISTC construction was identical to the CMTC. The design drawings for the mooring scheme are shown in Figures D.2, D.3 and D.4 (note there were minor changes to dimensions and materials from the design drawings during construction). 203 D.3 Thermistor Chain Deployment and Recovery Both thermistor chains were recovered by means of a retrieval line attached to the bottom of the thermistor chain. The recovery operation started at the new mooring buoy to which the other end of the retrieval line was attached. By pulling on the retrieval line the bottom of the CMTC was brought to the surface first. Thermistors on the wire rope were removed from the bottom up and the retrieval line and wire rope were allowed back into the lake. At the top of the thermistor chain the 20 m tenex rope with thermistors attached was removed and the wire rope attached to the CMTC buoy. The mooring system, wire rope and retrieval line were left in place while the thermistors (both loose and on the 20 m rope) were removed to a safe dry location. The downloading, checking, battery replacement and reprogramming of the thermistors was a two day process. The redeployment followed the same operation but in reverse (Figure F.12). The thermistor chains were recovered and deployed bottom up to overcome the perceived difficulty of lifting the heavy thermistor chain from the top which was heavily instrumented. Two problems were encountered with deployment and recovery at the CMTC. First, the anchor system on the temporary mooring buoy dragged several times under windy conditions when the boat was tied up to it as this added extra wind loading. The anchor did not provide a strong hold because the depth of water reduced the scope (angle of mooring line) so that when pulled it tended to lift rather than drag as designed. The dragging was worsened by the direction of pull, which was off pit wall benches, so that the anchor was pulled into deeper water further reducing the scope. The second problem was that the lower thermistors became tangled twice, which lifted instruments above their designated depth. Tangling was cause by deploying the instruments too quickly from a position too close to the CMTC buoy. The SISTC system was trouble free because of the shallower depth and greater rigidity of the three point mooring. D.4 Thermistor Instruments The different thermistors used had different specifications (Table D.l) which were suitable for different locations on the thermistor chains. The Onset Stowaways had slow response and low accuracy and were used in the upper layer where large but slow 204 temperature changes were expected. The Brancker TR 1000s were accurate with fast response and were used in the upper and lower interfaces to measure the rapid temperature fluctuations caused by internal waves. TSKA Wadars were accurate but with slow response and were used in areas like the intermediate layer where small and slow temperature changes were expected. Brancker XP200s where used for coarse temperature measurements and pressure (depth) checks. Up to 51 thermistors where used for the first half of the deployment year (Figure F.ll), dropping to 33 in the second half of the deployment as instruments borrowed from University of Victoria were returned (Table D.2). Thermistors turnover also occurred due to repair needs. The type and availability of instruments as well as changing conditions in the pit lake necessitated modifications to the thermistor layout. For example, the rising upper interface position and fluctuating surface water level required that instruments monitoring the upper interface be moved up and down for different deployments. Table D.3 and D.4 summarize the thermistor deployments including instrument serial number, type, depth, sampling rate and miscellaneous information. Note that a Stowaway thermistor was deployed in the SIS sedimentation pond to measure the temperature of the SIS injection system discharge. In addition, a XL200 temperature/pressure instrument was placed in the meteorological station logger box December 1999 - June 2000 to measure atmospheric pressure. D.5 Thermistor Calibration Thermistors were calibrated before and after their deployment to the Island Copper pit lake. Calibration was performed in a temperature controlled and mechanically stirred water bath. Thermistor measurements were compared to temperature readings from a Guildline 9540 Platinum Resistance thermometer (accuracy 0.01 °C) and a Dymec Model 2801A Quartz thermometer. Details of the calibrations are given in Table D.5. Calibration equations to correct thermistors were consistent between pre- and post calibrations, with the exceptions noted in Table D.5. Depths were marked on the thermistor chain before deployment. Depths were corrected in the post-processing to allow for the submergence of the thermistor chain buoys, 205 which was estimated before deployment but found to be different in the field. The SISTC depths were corrected by -0.11 m and the CMTC depths were corrected by -0.07 m. 206 Table D.l. Specifications of Thermistors # Instrument Data Memory (# samples) Accuracy Resolution Time Constant 10 Brancker TR-1000 T 256 k, sc 0.02 °C 0.005 °C 30-40 s 5+12 UVic Brancker TR-1000F T 256 k, sc 0.02 °C 0.005 °C 1-2 s 4 TSKA Wadars T 64 k, sc 0.02 °C 0.005 °C 3 min 14 Stowaways T 8 x 64 k, sc 6 x 8 k, sc 0.2 °C 0.2 °C 10 min 3 + 5 UVic Brancker XL-200 (1 x700 m, 2x50 m, UVic 5 x 100 m) T P 14 k, sc 0.1 °C 0.25 % of max press 0.01 °C 10 min T - temperature, P - pressure, sc - self contained submersible' UVic - University of Victoria instruments used December 1999 to June 2000. Table D.2. Total Instrument Summary (working). Instrument Dec 6ch 1999 - 03rd Feb 2000 - 16th March 2000 8th June 2000 - 1 Sept 2000 -31st Jan 2000 13 March 2000 - 5th June 2000 29 Aug 2000 29 Nov 2000 Stowaway 15 15 15 15 14 Wadar 4 4 4 4 4 XL200 8 8 8 2 2 TR1000 19 24 24 12 15' Total 46 51 51 33 35 207 &0| C •c o o £ 8 ° Q Z 208 c o o C3 s 3 c I c ' 3 U o H C •c o o 03 C U cn 03 H E & 5 Q E l o CN I cn o c E l o -a vo cn cn CN S o 0 0 o o o s o o o s 8 CN cn o oo Tt-oo on o| s vn v o VO >/-> o o o 2 CN m o • a vo cn vo c l o o o CN X CN - * V-) o cn r-CN o o o iS oo o o o s o o o s l o CN J IX VO o O o o s o o o VO 00 VO -3 o o CN X vo m o CN 209 JS 00 CB cj Q H U c o ca £ o o E o E 00 -a <D c tD 6 o —1 o CN CN ro oo •a c 3 O kn ca -a "ob c C3 e CN CN CN ro oo oo 14 o i= o 8 «> _o "u X) CA .g CU • a -a o £ w CO ca 00 oo C > CN 8 ^ ® X tU > c3 c o u E E oo CN ro oo •a c 3 o -a Ji "ob c ctl Q U tu > J4 o tu J S o tu a, vo uo O * o o CN M U tU J3 U T 3 "5b c ca r ® CN CN ® cn oo oo * fc oo c £ E E e E E £ o o ro VO VO o O o O o o <N CN • i CN l l CN CN CN l l CN l l CN l l II _ U II _ U II II tU tU II tU II tU II tU "ob "5b 00 00 00 00 00 s a c a B a c cs ct) ct) ca ca ca ca kH kn kn kH tU tu tU <u <u tu tU <C <£ 4=i <C •a ca ct) ct) ca ca ca ca J3 J3 J3 J3 J3 X x Q. cl a, D. 0 . &. Q. <u tU tU tU tu tu tu T3 T3 T3 T3 T3 - a - a T3 T3 T3 T3 T3 -a -a 4) U U •B 3 3 ta c pa a ta a ta a ta a ta a ta a k" kH t n V J SH t n tn CA C/5 03 tu tU D 0) 4) ca s CN * E E E E E E E o oo —H uo o CN CN CN ro I/O CN CN CN CN CN CN CN ® ® ® ® ® ® ® ro ro U0 VO OO oo ro UO U0 uo uo oo .—i —H —H O •> oo oo 00 00 > > * * * * * x ca tu J 3 ca X CD 00 a cd J 3 O ca X ca x o ca c o X) T 3 tu ro CN •<t UO * o o CN —1 X tu o o ca 4—) 3 X) 3 O - a ca tu T3 tu ta ON uo CO [> [> on c o o tu fc o cj c o tU tu ca tu 'o tu c tU E 3 c ca ca 4 c o ca tu o a. T3 <u > X o c o ta k-x 13 U uo * c .o tj tU fc o CJ c o tu k-tu o C3 k-3 o tj ca -a tu 'o o p. c CU E 3 B ca kH tu ca tu =3 Ji 00 ca j= E O ca o O - a > x o B o ca U vo * CJ tu tU tu o-ca •<d x x g - B 'tn T3 tU O . D . ca tu 4 tu c a o CJ tu J 3 tU e B o cj 00 c X T3 tj tU tn >> ca s ca O o B * C3 T3 tj tU fc o o B cj tU X tu > ca tu E tu 00 00 o 4 o 00 uo tU kH <o 4 & 3 O H E o cj T3 B ca j= o ta 4 4 p. 2 tu o o o CN 4) B 3 210 oo oo oo oo oo oo oo oo oo 211 8 vo O 00 o c. -4—> o ,o E-i 212 Table D.5. Thermistor Calibration Summary Description SN Type Pre-cal Post-cal Comment * StowAway 73151 Sl-Stow,32k cl4 c20 pre- and post calibrations match 15 73153 S2- Stow,32k cl4 c20 pre- and post calibrations match 73154 S3-Stow,32k cl4 c20 pre- and post calibrations match 73155 S4-Stow,32k cl4 c20 pre- and post calibrations match 73156 S5-Stow,32k cl4 c20 pre- and post calibrations match 73157 S6-Stow,32k cl4 c20 pre- and post calibrations match 73159 S7-Stow,32k cl4 - *4, no post cal 4880 S8-Stow,32k cl4 c20 pre- and post calibrations match 1282 S9-Stow,8k c!4 c20 pre- and post calibrations match 3610 S10-Stow,8k cl4 c20 pre- and post calibrations match 0058 Sll-Stow,8k cl4 c20 pre- and post calibrations match 73097 S12-Stow,32k cl4 c20 pre- and post calibrations match 1284 S14-Stow,8k cl4 c20 pre- and post calibrations match 1299 S15-Stow,8k cl4 c20 pre- and post calibrations match 1281 S16-Stow,32k cl4 c20 pre- and post calibrations match Hobopro 277272 Hobopro 72 cl4 c20 pre- and post calibrations match TR1000 6775 TR1000 cl4,cl5 c20,c21 pre- and post calibrations match *6 15 6776 TR1000 cl4 c20 pre- and post calibrations match *6 6777 TR1000 cl4 c20 pre- and post calibrations match *6 6778 TR1000 cl4 c20,c21 pre- and post calibrations match *6 6779 TR1000 cl4 c20 pre- and post calibrations match *6 6780 TR1000 cl6,cl7 ,cl8 c20,c21 Brancker repaired prior to c 16 Brancker repaired prior to cl8 cl8 to postcal -> shift of 0.12°C post cals match *5 6781 TR1000 cl4,cl5 ,cl8 c20,c21 Brancker repaired prior to cl8 cl8 to postcal -> shift of 0.22°C post cals match *5 6782 TR1000 cl4 c20 pre- and post calibrations match *6 6783 TR1000 cl6,cl7 c20 Brancker repaired prior to cl6 pre- and post calibrations match *6 6850 TR1000 cl4 c20 pre- and post calibrations match *6 8834 TR1000F cl6,cl7 c20 new from Brancker prior to cl6 pre- and post calibrations match *6 8837 TR1000F cl6,cl7 c20,c21 new from Brancker prior to cl6 pre- to postcal -> shift/slope, max diff=0.03°C still within specification of 0.05°C precals match, postcals match *5 8838 TR1000F cl4,cl5 ,cl8 c20,c21 Brancker repaired prior to cl8 cl8 to postcal -> shift/slope, max diff=0.01°C still within specification of 0.05°C postcals match *5 8839 TR1000F cl6,cl7 c20,c21 pre- and post calibrations match *6 8840 TR1000F cl6,cl7 c20,c21 pre- and post calibrations match *6 213 Table D.5. Thermistor Calibration Summary (cont.) Description SN Type Pre-cal Post-cal Comment * TR1000 (Uvic) 12 8144 TR1000 cl4 cl8 pre- and post calibrations match 8145 TR1000 cl4 cl8 pre- and post calibrations match but steep slope, max diff = 0.02°C 8146 TR1000 cl4 cl8 pre- and post calibrations match 8147 TR1000 cl4 cl8 pre- and post calibrations match 8148 TR1000 cl4 cl8 pre- and post calibrations match 8149 TR1000 cl4 cl8 pre- and post calibrations match 8150 TR1000 cl4 cl8 pre- and post calibrations match 8151 TR1000 cl4 cl8 pre- and post calibrations match 8152 TR1000 cl4 cl8 pre- and post calibrations match 8153 TR1000 cl4 cl8 pre- and post calibrations match 8155 TR1000 cl4 cl8 pre- and post calibrations match 8157 TR1000 cl4 cl8 pre- and post calibrations match XL200 3 5423 XL200 cl4 - no post cal *3 5424 XL200 cl4 c20 pre- and post calibrations match 7056 XL200 cl4 c20 pre- and post calibrations match XL200 (Uvic) 5 8136 XL200 cl4 c20 pre- and post calibrations match 8137 XL200 cl4 c20 pre- and post calibrations match 8138 XL200 cl4 cl8 pre- and post calibrations match 8139 XL200 cl4 cl8 pre- and post calibrations match 8140 XL200 cl4 cl8 pre- and post calibrations match Wadars 4 W178 Wadar cl5 c20 pre- and post match, used combined cal curve *9 W181 Wadar cl5 c20 pre- and post match, used combined cal curve *9 W182 Wadar cl5 c20 pre- and post match, used combined cal curve *9 W183 Wadar cl5 c20 pre- and post match, used combined cal curve *9 *3 XL200 #5423 failed by corrosion. *4 S7-73159 initially failed to readout, but was successfully readout with a battery change back in the lab. *5 Calibration shift observed in post calibration, which was less than instrument specified accuracy. Use pre-calibration correction. *6 Calibration shift observed in post calibration which exceeds instrument specified accuracy. Use pre-calibration correction. *9 Pre- and post calibrations combined to give greater data spread for calibration curve. 214 215 • IE BUDYS c/) 5 °« 1 OUTFALL DATING L TING SIS SYSTEt^ >-• 3 m </) - NEW FL CDNNEC z L D Z LJ L J CY ON <r 1 — • z > (J> _ i OO • 0 > >-i — (/> OJ • CY i—i Z X Q_ z r~i <c tt l_J CQ PER CTII CH NG AS ISHER VEN" CL L J CY i—< LJ ISHER • • ~> a > z _ J ISHER z C J Z >— <r <c <: L_ ' 1 c/> L Y _ l u r o a z X C L co r— OJ <L 1— _ J OO Z 3 L J • X i— 216 Ld 2 1 [ i | > | • JZ <r C/) <r C/) z <c • Ld C J <L~ Ld U (/) ( X c_> PQ Ld X GO 0> ON 0> Ld PQ Ld > • 0 1 O J 217 218 APPENDIX E Meteorological Measurements E.l Introduction A meteorological station was deployed concurrently with the thermistor chains to measure the weather - the architect of the lake thermal stratification. This Appendix provides details on the meteorological station deployment, instrumentation and calibration. E.2 Meteorology Station Setup The Environmental Fluid Mechanics group's "Terra" meteorological station was operated at the Island Copper Mine from December 1999 to November 2000 (Figure F.9). It is referred to in the thesis as the Island Copper meteorological station, or in reference to meteorological conditions as Island Copper. The Island Copper meteorological station was located onshore near the south injection system at 5605850 N 608400 E (50°36' 127°29') and an elevation of ~ 8 m above mean sea level (Figure 3.1). Chapter 4 discusses the rationale for the location of the meteorological station on land rather than on the lake itself. Measuring the meteorological conditions on the lake surface is generally preferable. However, the variability of wind across the lake due to sheltering and blocking from the pit walls nullified the advantages of deploying the single meteorological station on a raft on the lake surface. It was felt that locating the meteorological station on the lake did not guarantee results that were significantly more representative of the average lake wind field. The on land station was at a similar elevation to the lake water surface and has similar exposure to winds. Importantly, an onshore meteorological station ensured less risk to equipment and greater data security. 219 The meteorological station was located on the southern shore as it is flat and open to Rupert Inlet, whereas the west, north and east sides of the pit lake are surrounded by 60-72 m high pit walls. The specific location was selected due to ease of access and because there was an opening in the alder trees. These trees were planted as part of the reclamation over the beach dump and have grown to ~ 2 m tall. The weather conditions at the meteorological station site were considered similar to conditions in the middle of the pit lake, as both had comparable exposure to the dominant westerly and easterly winds. E.3 Meteorology Station Instrumentation The Island Copper meteorological station was equipped to measure rain, wind speed, wind direction, air temperature and relative humidity, and solar radiation. An R.M Young propeller vane anemometer was used for the wind speed and direction. The anemometer was mounted 4.4 m above the ground. Vaisala temperature and relative humidity (RH) sensors were mounted in a R.M. Young radiation shield. Short wave solar radiation was measured with a Li-cor pyranometer. Measurements from these instruments were recorded with a Terra Sciences data logger at 15 min intervals. Temperature, RH, and solar radiation were spot measurements, whereas wind speed and direction were 15 min averages. Rain was measured with a tipping bucket rain gauge every 0.26 mm and recorded by an Onset Event logger. These were the primary measurements and unless stated otherwise the data discussed in the thesis was from these instruments. Data were not recorded for the time it took to download and restart the data loggers, which was 1-2 hours for the Terra Sciences data logger. The only instrumentation problem occurred with the tipping bucket rain gauge, which failed to record rain for the September to November 2000 deployment. The failure resulted from the tipping bucket not tipping, which was possibly caused by the Onset Event logger wires restricting the motion of the tipping bucket. Secondary measurements were taken as backup. Temperature and RH were measured with an Onset Hobopro at 15 min intervals. Pressure and temperature were measured for the first six months with a XL200 (see Appendix D) placed in the data logger box. The 220 meteorological station measurements, instrumentation and periods of record are summarized in Table E. 1. E.4 Meteorology Station Calibration All meteorological station sensors were calibrated before and after deployment to Island Copper. The R.M Young anemometer specified wind direction accuracy is ± 3°. Correction to the anemometer wind direction was important as the calibration correction was up to 25° at large angles (approaching 360°). Pre- and post-calibrations were constant. The wind direction should be accurate to ± 5°, including the uncertainty of setting the vane zero point. Wind speed was checked against a fixed speed (3 m s"1) motor. Pre- and post deployment wind speeds were the same as that recorded for the instrument after factory calibration. The R.M Young anemometer specified wind speed accuracy is ± 0.3 m s"1 with a threshold of 1 m s"1. Vaisala temperature and RH sensors have specified accuracies of 0.2 °C and 3 %, respectively. The Vaisala temperature required a correction - 2.37 °C but this correction was constant for pre- and post calibrations. The Vaisala RH pre- and post calibrations agreed to within ± 5 % and this better reflects the accuracy. The Li-cor pyranometer specified accuracy is ± 5 % and the pre- and post calibration were within ±3 %. E.5 Port Hardy Airport Meteorological Station Environment Canada operate a permanent weather station at Port Hardy Airport (Station ID 1026270, location 50°41' 127°22' and elevation of 22 m), which is 12 km north-east of Island Copper and on the east coast of Vancouver Island (Figure B.l). Temperature, relative humidity, pressure, precipitation was obtained for 1996 - 2001 and windspeed and direction for 1953 - 2001. These data are detailed in Appendix B. Cloud cover was estimated from daily noon observations of weather conditions at Port Hardy Airport 221 Table E. l . Meteorological station instrumentation and deployment periods. Parameter Instrument 9Dec99-31Jan00 31Jan00-13Mar00 13Mar00-5June00 5June00-09Aug00 29Aug00 30Nov00 Temperature Vaisala X X X X X RH Vaisala X X X X X short-wave solar Li-cor X X X X X wind speed R.M Young X X X X X wind direction R.M Young X X X X X rainfall tipping bucket X X X X X Temperature Hobopro X X X X X RH Hobopro X X X X f Temperature XL200 X X X - -Pressure XL200 X X X - -x = instrument used and data recovered, f = instrument failure and no data, - no instrument. 222 APPENDIX F Photographs Figure F.l . Oblique aerial of Island Copper Mine and Rupert Inlet view from east looking west (August 1993). Figure F .2. Flooding the Island Copper Mine pit with seawater from Rupert Inlet (background) (June 1996) April 9, 1996 Figure F.3. Island Copper Mine pit post closure view looking west (9 April 1996). Figure F.4. Island Copper Mine pit during flooding view looking west (17 July 1996). Figure F.5. Island Copper Mine pit at end of flooding view looking west (29 Julyl996). Figure F.6. North injection system outfall including de-aerator. Figure F.8. South injection system sedimentation pond and outfall intake. 228 

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

Comment

Related Items