Science, Faculty of
Earth, Ocean and Atmospheric Sciences, Department of
DSpace
UBCV
Burns, Sean Raymond Patrick
Master of Applied Science - MASc
University of British Columbia
Abandoned exploration boreholes are commonly found around mine sites in a fractured crystalline rock environment. If the abandoned boreholes have not been properly decommissioned they have the potential to create connections through the rock fractures and influence ground water flow and contaminant transport. A fully three-dimensional discrete fracture model is used to investigate the impact of abandoned boreholes on contaminant transport from a waste-rock pile overlying a fractured rock mass. Dissolved contaminants travel through the fractured rock mass under the influence of a subhorizontal regional hydraulic gradient towards a downstream compliance boundary. A number of different fracture geometries are investigated to gain an understanding of the field situations in which abandoned boreholes can be expected to have an impact. The effect of fracture density, transmissivity contrasts, and borehole diameter and location are studied. The simulation results show that vertical abandoned boreholes are most likely to have an impact when large, sub-horizontal, high-transmissivity features are present in the network. Low fracture density, aperture variability, relatively high horizontal transmissivity, and the presence of major features in the fracture network all lead to abandoned boreholes having a greater overall influence. If an abandoned borehole is transversely offset from the central flow line passing through the source zone the contaminant plume can migrate towards the borehole in a direction not predicted by the average regional hydraulic gradient. In field-scale fracture networks smaller borehole diameters leads to shorter breakthrough times and higher contaminant concentrations at the downstream boundary due to the interplay between the fracture network and borehole void space. The presence of abandoned boreholes can be expected to have important implications in the design of monitoring networks to detect ground water contamination when these fracture network and abandoned borehole properties exist.
7492989 bytes
application/pdf
CONTAMINANT TRANSPORT THROUGH A B A N D O N E D BOREHOLES IN FRACTURED R O C K by SEAN R A Y M O N D PATRICK BURNS B . A S c , University of British Columbia, 1997 A THESIS SUBMITTED IN PARTIAL F U L F I L L M E N T OF THE REQUIREMENTS FOR THE DEGREE OF M A S T E R OF APPLIED SCIENCE in THE F A C U L T Y OF G R A D U A T E STUDIES Department of Earth and Ocean Science We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH C O L U M B I A June 2000 © Sean Raymond Patrick Burns, 2000 In p r e s e n t i n g t h i s t h e s i s i n p a r t i a l f u l f i l m e n t o f t h e r e q u i r e m e n t s f o r an a d v a n c e d degree at the U n i v e r s i t y o f B r i t i s h C o l u m b i a , I agree t h a t t h e L i b r a r y s h a l l make i t f r e e l y a v a i l a b l e f o r r e f e r e n c e and s t u d y . I f u r t h e r agree t h a t p e r m i s s i o n f o r e x t e n s i v e c o p y i n g o f t h i s t h e s i s f o r s c h o l a r l y p u r p o s e s may be g r a n t e d by t h e head o f my depar tment o r by h i s o r h e r r e p r e s e n t a t i v e s . I t i s u n d e r s t o o d t h a t c o p y i n g o r p u b l i c a t i o n o f t h i s t h e s i s f o r f i n a n c i a l g a i n s h a l l not be a l l o w e d w i t h o u t my w r i t t e n p e r m i s s i o n . Department The U n i v e r s i t y o f B r i t i s h C o l u m b i a V a n c o u v e r , Canada Abstract Abandoned exploration boreholes are commonly found around mine sites in a fractured crystalline rock environment. If the abandoned boreholes have not been properly decommissioned they have the potential to create connections through the rock fractures and influence ground water flow and contaminant transport. A fully three-dimensional discrete fracture model is used to investigate the impact of abandoned boreholes on contaminant transport from a waste-rock pile overlying a fractured rock mass. Dissolved contaminants travel through the fractured rock mass under the influence of a sub-horizontal regional hydraulic gradient towards a downstream compliance boundary. A number of different fracture geometries are investigated to gain an understanding of the field situations in which abandoned boreholes can be expected to have an impact. The effect of fracture density, transmissivity contrasts, and borehole diameter and location are studied. The simulation results show that vertical abandoned boreholes are most likely to have an impact when large, sub-horizontal, high-transmissivity features are present in the network. Low fracture density, aperture variability, relatively high horizontal transmissivity, and the presence of major features in the fracture network all lead to abandoned boreholes having a greater overall influence. If an abandoned borehole is transversely offset from the central flow line passing through the source zone the contaminant plume can migrate towards the borehole in a direction not predicted by the average regional hydraulic gradient. In field-scale fracture networks smaller borehole diameters leads to shorter breakthrough times and higher contaminant concentrations at the downstream boundary due to the interplay between the fracture network and borehole void space. The presence of abandoned boreholes can be expected to have important implications in the design of monitoring networks to detect ground water contamination when these fracture network and abandoned borehole properties exist. ii Table of Contents Abstract ii Table of Contents iii List of Tables v List of Figures vi Acknowledgements ix 1.0 Introduction 1 2.0 Review 6 3.0 Methods 13 3.1 Fracture Network Modeling 13 3.2 Overview of the FRACMAN/MAFIC suite 13 3.2.1 FRACWORKS 14 3.2.2 MESHMONSTER 17 3.2.3 EDMESH 18 3.2.4 MAFIC 19 3.3 Solute Transport 20 3.4 Implementation of Abandoned Boreholes 22 3.4.1 Borehole as a Discrete Fracture 23 3.4.2 Borehole-Fracture Intersections 27 3.4.3 Test of Flow through Borehole 28 3.4.4 Test of Particle Tracking through Borehole 29 3.4.5 Effect of Borehole Diameter 31 3.5 Conceptual Model Description 33 4.0 Results 48 4.1 Base Case Fracture Network 48 4.1.1 Fracture Generation Region 50 4.1.2 Base Case Borehole 50 4.1.3 Grid Discretization 52 4.1.4 Case 1 - Base Case Results 54 4.2 Influence of Background Fracture Network 60 4.2.1 Case 2 - Higher Horizontal Transmissivity 60 4.2.2 Case 3 -Higher Vertical Transmissivity 64 4.3 Influence of Major Features 66 4.3.1 Case 4 - Major Feature Transmissivity = lxl0'3 m2/s 67 4.3.2 Case 5 - Threshold Effect 69 4.3.3 Case 6 - Lower Feature Terminates atx = 30 m 70 iii 4.4 Influence of Borehole Location and Diameter 72 4.4.1 Case 7 - Borehole atx = 0 m 72 4.4.2 Case 8 - Borehole atx = -35 m 73 4.4.3 Case 9 - Borehole aty = 15 m 75 4.4.4 Case 10- Influence of Borehole Diameter 76 4.5 Case 11 - Variable Network with No Major Features 77 5.0 Discussion 112 5.1 Results Summary 112 5.2 Fracture Network Structure 116 5.3 Implications for Monitoring Network Design 119 5.4 Proper Borehole Decommissioning 124 6.0 Conclusions 126 References 128 iv List o f Tables Table 3-1: Parameters for analytical problem 32 Table 4-1: Generation parameters for base case fracture network 48 Table 4-2: Flow and transport output parameters for all simulations 55 Table 4-3: Generation parameters for Case 11 fracture network 78 Table 4-4: Borehole parameters for Case 11 80 v List of Figures Figure 1-1: Contaminants leaching from a waste rock pile into a fractured rock aquifer 5 Figure 2-1: Aquifer cross-contamination due to abandoned well in multiaquifer system 12 Figure 3-1: Definition of orientation convention used in F R A C M A N 36 Figure 3-2: Orthogonal view of the base case fracture network 37 Figure 3-3: Parameters for incorporation of a circular borehole as a rectangular fracture 38 Figure 3-4: Discretization algorithm around the intersection of a borehole with a fracture 39 Figure 3-5: Test of flow through a 10 cm diameter borehole 40 Figure 3-6: Fracture network designed to test particle-tracking through the borehole 41 Figure 3-7: Particle breakthrough at downstream boundary for test fracture network 42 Figure 3-8: Pathway of a single particle traveling through the test fracture network 43 Figure 3-9: Analytical model used to test the effect of borehole size on flow rate and residence times 44 Figure 3-10: Flow behaviour of a large void space within a thin conduit 45 Figure 3-11: Hypothetical mine site used to develop the conceptual model 46 Figure 3-12: Cross-section view through the conceptual model showing boundaries and dimensions 47 Figure 4-1: Cross-section of base case fracture network showing fracture traces for background fractures and major features 83 Figure 4-2: Orthogonal view of modeling domain showing flow and transport boundaries 84 Figure 4-3: Diagram of modeling domain and fracture generation region 85 Figure 4-4: Finite element grid discretization 86 Figure 4-5: Case 1 - Particle breakthrough at downstream exit boundary 87 Figure 4-6: Case 1 - Lateral breakthrough distribution of particles across exit boundary 88 Figure 4-7: Pathways of selected particles after exiting the borehole 89 Figure 4-8: Snapshots of particles traveling through the fracture network when no borehole is present 90 Figure 4-9: Case 2 - Particle breakthrough at downstream boundary 91 Figure 4-10: Case 2 - Lateral breakthrough distribution of particles 92 Figure 4-11: Case 3 - Particle breakthrough at downstream boundary 93 Figure 4-12: Case 3 - Lateral breakthrough distribution of particles 94 Figure 4-13: Case 4 - Particle breakthrough at downstream boundary 95 Figure 4-14: Case 4 - Lateral breakthrough distribution of particles 96 Figure 4-15: Case 5 - Particle breakthrough at downstream boundary 97 Figure 4-16: Case 6 - Particle breakthrough at downstream boundary 98 Figure 4-17: Case 6 - Lateral breakthrough distribution of particles 99 Figure 4-18: Case 7 - Particle breakthrough at downstream boundary 100 Figure 4-19: Case 7 - Lateral breakthrough distribution of particles 101 Figure 4-20: Case 8 - Particle breakthrough at downstream boundary 102 Figure 4-21: Case 8 - Lateral breakthrough distribution of particles 103 Figure 4-22: Case 9 - Particle breakthrough at downstream boundary 104 Figure 4-23: Case 9 - Lateral breakthrough distribution of particles 105 Figure 4-24: Case 10 - Influence of borehole diameter on downstream breakthrough curves 106 Figure 4-25: Case 10 - Influence of borehole diameter on total system flow and median particle residence times 107 Figure 4-26: Case 11 - Particle breakthrough at downstream boundary for a random fracture network with no major features 108 vii Figure 4-27: Isometric views of borehole locations for Case 14 109 Figure 4-28: Case 11 - Particle breakthrough at downstream boundary for an entirely random fracture network 110 Figure 4-29: Particle breakthrough at downstream boundary for a constant aperture fracture network 111 viii A c k n o w l e d g e m e n t s I would like to thank my supervisor, Leslie Smith, for his support, advice, guidance and patience during the preparation of this thesis. I am grateful for the support and advice of Petros Gaganis, Bob Parney, and the other graduate students in the U B C hydrogeology group who have given me advice over the past two years. Golder Associates Ltd. developed the discrete fracture modeling code used in this study, and special thanks go to the F R A C M A N group in Seattle, including Tom Doe, Bi l l Dershowitz, Paul La Pointe and Glori Lee, who made the source code available, without which I could not have made the necessary modifications. I would like to give a very special thank you to my parents, Gail and Ray Burns, who have been incredibly supportive over my university career, and who have always encouraged me to follow my dreams. Claire Bradford also deserves special mention for her emotional support and editing prowess. Funding from the National Science and Engineering Research Council provided financial support for this research. ix 1.0 Introduction Protecting and maintaining the quality of natural water supplies is necessary for the well being of all life on earth, including humanity. As the population of the world expands the demand for material goods will increase, and new sources of raw minerals will need to be uncovered and developed. In the past, poor and ill-advised mining practices have resulted in catastrophic pollution to water supplies and long-term environmental damage. While mining practices and accountability have improved significantly in recent years, environmental damage continues to occur at active sites and there is still room for improvement. Ground water is the most vulnerable fresh water source to long-term pollution from mining waste, as residence times in the subsurface can be on the order of tens to hundreds of years. Once a ground water resource has been contaminated it is technically difficult and very expensive to remediate, i f it can be done at all. It is more cost-effective to engage in practices that prevent ground water contamination from occurring in the first place, than to attempt to clean up the problem afterwards. If proper care is taken a balance can be met between meeting the resource requirements of society while sustaining the quality of the ground water for the long term. Mining typically requires the extraction of a small amount of valuable ore minerals from large volumes of sub-economical gangue material. Open-pit mining can result in waste rock piles and tailings dams that need to be managed for decades or even centuries. Waste rock is the material removed to gain access to the ore body, and usually contains 1 metallic minerals such as pyrite, which have no economical value. Figure 1-1 depicts a situation in which a waste rock pile is situated in an elevated area overlying a fractured rock mass. Over long time periods exposure to rainfall and atmospheric oxygen can result in acidic water leaching out of the waste rock piles and mobilizing toxic heavy metals. This process is commonly referred to as acid rock drainage (ARD). If no impervious barrier is present to impede migration the metals will travel vertically through the unsaturated zone and reach the water table, contaminating the underlying ground water. Under the Environmental Assessment Act in British Columbia a potential mine site must be rigorously studied to assess any potential for ground water contamination (RSBC 1996). A compliance boundary for a site is established to delineate the maximum extent at which some degradation to the local ground water quality is deemed acceptable by the regulatory agency. A network of compliance monitoring wells is installed at the compliance boundary and ground water is regularly sampled. If contamination above a threshold value is detected at the compliance boundary then the site owner will face financial penalties on top of the cost of containing and remediating the contamination. The processes governing A R D are still an active area of research, and there is a great deal of uncertainty when predicting contaminant fluxes to the subsurface. Once contaminants have reached the water table there is the potential for off-site migration in the direction of ground water flow. Numerical models are useful to help predict contaminant fluxes from potential sources at a site and assess the risk of exceeding the threshold concentration at 2 the compliance boundary. Ground water modeling for a fractured rock terrain is more complicated than in porous media due to the heterogeneous nature of a fractured rock mass. Flow is often concentrated in a small subset of the total number of fractures and contaminant transport occurs along complicated, tortuous pathways. Gathering hydrogeological characterization data in crystalline rock is very expensive, so the amount of information available to the modeler is usually quite limited. The higher uncertainty in characterizing fractured rock leads to additional engineering effort and expense. One source of uncertainty at a mine site overlying fractured rock is the presence of any abandoned exploration boreholes around the ore body. Before a mine can be developed the ore body must be delineated, so a developed mine site may have tens to hundreds of exploration boreholes. If previous phases of exploration at a site failed to detect economic mineralization then all drilling records for a site could be lost. If the boreholes were not properly grouted before being abandoned they can act as conduits between stratigraphic layers, or form preferential pathways in the fracture network. The potential for cross-contamination of layered aquifers that are perforated by boreholes is well understood. It is hypothesized that in a fractured rock the boreholes could connect fractures creating a preferential pathway and allowing contamination to travel faster to the compliance boundary, or contaminating a previously protected zone. Connections can be formed that cannot be anticipated based solely on knowledge of the statistical fracture geometry and the contaminant plume may move in a direction not anticipated based on conventional porous medium approximations. 3 This thesis investigates the potential impacts that these abandoned boreholes could have at a mine site underlain by a fractured rock mass. The conceptual model assumes that contaminated leachate from a waste rock pile has reached the water table. A three-dimensional discrete fracture code is used to model the flow and contaminant transport through the fracture network and to assess the impacts the boreholes have on the system. The objective of this thesis is to gain an understanding of the types of subsurface conditions where abandoned exploration boreholes have the potential to exacerbate a contamination problem, or reduce the effectiveness of a monitoring network design. The first two chapters present the motivation for the work and a review of related work on fracture network modeling and the effects of abandoned wells and boreholes. In the third chapter the modeling methodology is presented along with the modifications that were necessary for the purposes of this study. The fourth chapter presents the results of the set of simulations used to demonstrate the effects of boreholes in different fracture networks. The final two chapters explain the implications of the simulation results and summarize the conclusions that can be drawn from this study. 4 5 2.0 Review Numerical models for simulating flow and transport in porous media have been available for several decades (e.g. Freeze, 1971; McDonald and Ffarbaugh, 1988; Zheng, 1990). Recently, the potential for siting high-level radioactive waste in deep crystalline rock has prompted interest in modeling flow and transport through fractures (e.g. Dverstorp et al., 1992; Chan et al., 1993; Bodvarsson et al., 1997). Fracture modeling has been applied on a variety of length scales from a single fracture to field-scale networks. The simplest model for fracture flow is that of two parallel plates of constant aperture (e.g. Snow, 1965; Witherspoon et al., 1980). Several studies have been conducted on flow and transport through a single fracture with various complicating factors including spatial variability in the aperture (Brown, 1987; Raven et al., 1988; Moreno et al., 1988), diffusion into the porous rock matrix (e.g. Sudicky and Frind, 1982; Neretnieks, 1980) and surface sorption (e.g. Burkholder, 1976; Freeze and Cherry, 1979; Wels and Smith, 1994). Laboratory experiments on flow through a single fracture have suggested that the cubic law is not valid for most naturally occurring rough-walled fractures (e.g. Raven et al., 1988) but most discrete fracture models make this assumption due to limitations on the complexity that can be modeled over large domains. Field scale models for flow and transport through fracture networks fall under one of three broad categories: (National Research Council, 1996) (1) equivalent continuum models which model the fracture network using average properties (e.g. Neuman, 1987; Carrera et al., 1990), (2) discrete network simulation models which account for each fracture individually (e.g. Smith et al., 1985; Dershowitz et al , 1995), and (3) hybrid 6 techniques which apply some combination of the previous two approaches (e.g. Schwartz and Smith, 1988; La Pointe et al., 1995). Equivalent continuum models do not attempt to model each fracture pathway individually, but use an average representation of the hydraulic conductivity field. The main advantage of these models is that they can be less computationally intensive than discrete fracture models so a larger domain can be modeled. This approach is most appropriate when modeling a rock with many fracture connections or significant matrix permeability. When the network is heterogeneous with large contrasts in hydraulic conductivity, equivalent continuum models require a fine grid discretization and can become as computationally intensive as discrete fracture methods (Naff et al., 1998a,b). In sparse fracture networks the individual fracture connections and tortuous pathways may need to be incorporated to realistically model solute transport. Discrete fracture network models attempt to include every important fracture in the rock mass as an individual feature in the model. The models can be either two-dimensional where the fractures are represented by interconnected line elements, or three-dimensional where the fractures are planar features. There are three-dimensional models available, based on finite difference (e.g. Therrien and Sudicky, 1996) or finite element (e.g. Dershowitz et al., 1995) discretization schemes, that allow for secondary flow through the porous matrix. Nordqvist et al. (1992) presented a three-dimensional model that could simulate aperture variability within individual fractures. The more complex the model, the smaller the size of the domain that can be modeled, due to computational 7 limitations. The main advantage of discrete network models is that they include the physical structure of the fracture network and can account for fracture connections and preferential flow paths that may be of primary importance, especially in sparse networks. The disadvantage to these models is that they are more computationally intensive, and require extensive site characterization data to be any more meaningful than an equivalent continuum approximation. Hybrid models have evolved from an attempt to combine the benefits of both equivalent continuum and discrete fracture models. There are a variety of approaches used to combine the physical accuracy of discrete networks with the computational advantages of a statistical continuum. One approach has been to analyze the fracture network and only include the major conducting conduits (e.g. La Pointe et a l , 1995). This approach is valuable only i f it can be assumed that the smaller fractures are not important to the flow and transport. Another method uses small-scale discrete fracture networks to generate statistics for flow and transport that are subsequently applied to a field-scale continuum model (Schwartz and Smith, 1988; Parney and Smith, 1995). This modeling approach is a subject of ongoing research to try and determine appropriate links between small and large-scale simulations. One common purpose for modeling subsurface flow from a potential contamination source is to aid in the design of monitoring networks. A significant amount of research has been done on optimizing the design of monitoring networks in porous media. Massmann and Freeze (1987) present a framework for the design of monitoring networks 8 which incorporates uncertainty in the selection of a best network alternative. Meyer et al. (1994) provide a comprehensive method for optimizing monitoring network design by minimizing network cost, maximizing the probability of detection, and minimizing the size of the plume at the time of detection. Storck et al. (1997) expanded this work to a fully three-dimensional analysis for porous media. Jardine et al. (1996) expand upon the method of Massmann and Freeze (1987) to provide a decision-analysis framework for the design of monitoring networks in fractured rock, using a two-dimensional discrete fracture model to evaluate the best alternative from a number of alternate monitoring strategies. A comprehensive presentation of monitoring network design in a three-dimensional fractured rock mass has yet to be provided. The potential environmental hazards from abandoned drinking water wells are widely known. Most of these wells contain a steel casing but through time the casing can corrode exposing the contacting sediments to the well bore water. If the well crosses a low permeability unit, such as a clay aquitard, there is the potential for cross contamination into previously protected aquifers. Figure 2-1 shows a typical situation where the potential for cross contamination exists. Pumping from the lower confined aquifer can cause contaminants to migrate downwards from the upper unconfined aquifer through the abandoned well. The ground water resource in the lower aquifer is polluted and the potential exists for contaminants to reach the water supply well. Many jurisdictions have legislation requiring landowners to properly grout and decommission abandoned wells on their property. If the landowner does not comply they can be held liable for any contamination problems brought about from their wells. 9 Numerous studies have been done on the potential for abandoned wells or unplugged boreholes to lead to subsurface contamination of ground water. The potential for water supply wells that connect previously isolated water-bearing zones to facilitate contaminant migration has been well documented (e.g. Gass et al., 1977; Lacombe et al., 1995). Drilling for oil and gas has left a legacy of abandoned wells and exploration boreholes that present a similar risk (e.g. Wait and McCollum, 1963; Warner and McConnell, 1993; Langhus and Snider, 1999). Studies into deep-well injection of hazardous waste have assessed the potential for fast contaminant migration to shallow aquifers through leaky boreholes (Javandel et al., 1988; Raven et al., 1990). In 1972 the Atomic Energy Commission was forced to abandon a potential high-level radioactive waste repository site near Lyons, Kansas, when it was discovered that the site was riddled with old oil and gas exploration wells (Zeller, 1973). Some researchers have developed analytical models to describe the steady-state fluid flow across an aquitard that has been penetrated by a circular well (e.g. Avci , 1992; Brikowski, 1993). Transient flow through a borehole has also been studied when hydraulic gradients exist due to injection wells (Javandel et al., 1988; Avci , 1994). Lacombe et al. (1995) developed a numerical model to study flow and advective-dispersive contaminant transport through a borehole that connected an upper unconfined aquifer with a lower confined aquifer. They showed that the abandoned borehole could cause significant cross-contamination when hydraulic gradients existed across the borehole due to pumping in the lower aquifer. 10 A l l of the above studies deal with the environmental impacts of abandoned boreholes in a porous medium or an equivalent continuum. No attempt has yet been made to assess the impact of abandoned exploration boreholes in a fractured rock environment using a discrete fracture model. The heterogeneous nature of contaminant transport through fractured rock means that equivalent continuum models are not sufficient to properly address the influence of leaky boreholes. The present study uses a discrete fracture model to investigate the impact of abandoned exploration boreholes on contaminant transport from a surface waste rock pile through a fractured rock mass. 11 12 3.0 Methods 3.1 Fracture Network Modeling In this study a three-dimensional discrete fracture model was used to model flow and transport through the fracture network. It was necessary to take account of the fracture scale heterogeneity to properly model the influence of an abandoned borehole in a fractured rock mass. The influence of individual fractures and connected pathways through the rock matrix are an important physical process that needed to be included in the chosen model. For this reason, using an averaging method such as the statistical continuum method would not be appropriate. The F R A C M A N / M A F I C suite of programs developed by Golder Associates (Dershowitz et al. 1995) was chosen as it provides nearly all of the required tools, from fracture network generation to flow and contaminant transport modeling. The main modification required was the addition of a routine to incorporate abandoned boreholes. F R A C M A N / M A F I C provides the advantage of allowing complex fracture geometry to be included in the analysis. This could allow for a wide variety of fracture networks to be investigated by following the method described in this thesis. 3.2 Overview of the FRACMAN/MAFIC suite Four applications from the F R A C M A N suite were used to model the discrete fracture network and solve the flow and transport problem: F R A C W O R K S , MESHMONSTER, EDMESH, and MAFIC. F R A C W O R K S is used to generate the three-dimensional fracture networks, MESHMONSTER is used to create a basic finite-element mesh, E D M E S H is used to make appropriate refinements to the mesh, and M A F I C is used to solve flow and contaminant transport through the fracture network. 3.2.1 FRACWORKS F R A C W O R K S allows a wide variety of discrete features to be generated from deterministic or stochastic descriptions. The base case fracture geometry for the simulations presented in the following chapter was generated using a combination of deterministic fractures to model two major features and two stochastic fracture sets for the background network. The first stage of generating a three-dimensional fracture set is to choose a generation model. The random networks used in this study were generated using the Poisson Rectangle model. This is a simplified version of the Enhanced Baecher model (Dershowitz et al., 1989), with fracture dimensions specified using length and width instead of an effective radius. The Poisson model assumes that the fracture centers are randomly distributed in space. Once the fracture center has been chosen, the fracture geometry is determined by specifying the dimensions and orientation of the fracture. Fracture mechanics suggest that the general shape of a fracture in homogeneous rock will be elliptical (Baecher et al. 1977). The base case geometry assumed square planar fractures for simplicity, although polygonal approximations to elliptical fractures are incorporated in the final simulations with minimal additional computational effort. 14 The fracture density is determined by specifying the number of fractures to generate, and the size of the generation region. In the base case the background network consisted of two sets of 750 fractures, generated in the following region: -62.5 m < x < 62.5 m, -37.5 m < y < 37.5 m, -25 m < z < 25 m. The generation region described above was chosen to be larger than the modeling region to avoid edge effects in the fracture network, as explained in Chapter 4. This results in 1500 fractures in a volume of 148000 m 3 , or 0.003 fractures/m3. The fracture density can be specified using a number of alternate density intensity measures. According to Dershowitz et al. (1995) the preferred measure for fracture intensity is the areal intensity, P32 P 32 = A f / V t (1) where: Af = Total area of fractures, and V t = Total volume. For the base case, the length and width of background fractures are both a constant 10 m, so the areal intensity P 3 2 = [(10 m x 10 m) x 1500] / 148000 m 3 = 1.01 m"1. It is recognised that an individual fracture length of 10 m is larger than is typical for a natural 15 fractured rock mass, but large fractures were necessary to form connected networks due to computational limits in the large fracture volume modeled. The fracture orientations are specified by distributions for the trend and plunge of the fracture pole. The fracture pole is defined as a line normal to the fracture plane. Figure 3-1 defines the trend as the clockwise angle between the negative x direction and the projection of the fracture pole onto the x-y plane. The plunge is the angle between the x-y plane and the pole. For both fracture sets in the base case a constant orientation was adopted for simplicity. When variation in the fracture orientation was incorporated in later simulations a bivariate-normal distribution was used. In the base case fracture network, the fractures in Set 1 were horizontal with a constant pole trend of 180° and a plunge of 90°. Set 2 were vertical and oriented parallel to the y-axis with a trend of 90° and a plunge of 0°. The transmissivity of fractures was specified directly in the fracture generation model, and the corresponding aperture calculated according to the cubic law for flow through parallel plates: (Snow, 1965) Q = _pgb^.dh ( 2 ) 12/i dl where Q is the volumetric flow per unit thickness p is the density of fluid (in this case water) g is the gravitational constant 16 b is the hydraulic aperture of the fracture JJ. is the viscosity of the fluid dh/dl is the hydraulic gradient. The term pgb3/12 p is the transmissivity, T, of the fracture. For water at constant density and viscosity the transmissivity is proportional to the cube of the aperture. In this model 3 2 the density of water used is 1000 kg/m , the gravitational constant is 9.81 m/s , and the 3 2 absolute viscosity of water is 1x10" Ns/m , resulting in the following relationship: T = kxb3 (3) where k\ is a constant equal to 8.175xl05 (ms)"1. Along with the two stochastic sets described above, two large horizontal fractures were generated to provide a fast pathway when an abandoned borehole was present. The large features were generated using a deterministic rectangular model in F R A C W O R K S . The fracture centers were specified at Xo(x, y, z) = (0 m, 0 m, 10 m) and (0 m, 0 m, -8 m). Both fractures had a length of 100 m and a width of 50 m so that they spanned the entire modeling domain. Figure 3-2 shows an orthogonal view of the base case fracture network generated using F R A C W O R K S . 3.2.2 MESHMONSTER This program generates a triangular finite-element mesh from the fracture network data files created by F R A C W O R K S . The algorithm first finds the intersections between all of 17 the fractures, then creates triangular elements to fill in the area of the fracture planes. Triangular elements are then further subdivided until all elements are below a maximum element area specified by the user. The aspect ratio, being the ratio of the longest element side to the shortest side, is also checked to eliminate long, thin fractures, which can cause solution difficulties in the numerical model. The version of the program used is not capable of automatically generating tetrahedral elements for the rock matrix, so i f matrix flow is desired these must be specified manually. A l l simulations in this thesis assumed that the matrix is impermeable. For many types of geologic environments, such as fractured sandstone, this is not a realistic assumption, but for the purposes of assessing the influence of abandoned boreholes matrix diffusion was not considered to be a key physical process. 3.2.3 EDMESH The meshing algorithm used in MESHMONSTER is complicated and imperfect, so the resulting mesh must be edited to remove problems that will interfere with the calculation of a flow solution (Parney, 1999 p. 142). E D M E S H automatically takes the mesh output from MESHMONSTER and removes or merges nodes that are likely to cause problems. Nodes that are not connected to any boundary through the fracture network (i.e. isolated fractures) are removed. E D M E S H can also be used to refine the finite-element mesh in regions where, for example, higher hydraulic gradients exist. E D M E S H was used in this thesis to refine the mesh at fracture intersections, and around the location of abandoned boreholes. The specific details of the mesh generating and editing algorithms are proprietary to Golder Associates. 18 3.2.4 MAFIC MAFIC provides the numerical solution to flow and solute transport through the discrete fracture networks generated with F R A C M A N . The details of the flow and transport solution can be found in Miller et al. (1995). The diffusivity equation describing flow can be written as (Bear, 1972): d_ dx, f dP dz^ — + PS — Sx, dx dP = p(a+(j)f5)— + q dt (4) where: xi = coordinate directions P = fluid density = fluid viscosity = permeability (absolute) P = fluid pressure g = gravitational constant z = vertical direction (upward) a = pore compressibility = porosity P = fluid compressibility q = source/sink term t = time 19 The approximation of an incompressible fluid is applied to fracture flow in two dimensions to simplify equation 4 to a simple volume-conservation equation: S — -TV2h = q (5) dt where: S = fracture storativity h = hydraulic head T = fracture transmissivity q = source/sink term t = time V 2 = two-dimensional Laplace operator MAFIC uses a Galerkin finite-element solution scheme to approximate the solution for this equation. A n Incomplete Cholesky Conjugate Gradient (ICCG) algorithm is used to solve the flow problem (Meijerink and van der Vorst, 1977). Once the flow solution is found and the hydraulic head at each node is known, the groundwater velocity in each element can be calculated for use in the particle-tracking algorithm. In this thesis only steady-state flow is modeled, but the analysis could easily be extended to include transient flow with an associated increase in the required computational effort. 3.3 Solute Transport MAFIC simulates transport through the finite-element mesh with a particle-tracking algorithm based on similar approaches in two dimensions (Smith and Schwartz, 1984; Hull et al., 1987). The longitudinal distance traveled by an individual particle in any 20 given time-step contains a deterministic convective component and a stochastic dispersive component. A second transverse dispersive component, orthogonal to the first dispersive component, is included to determine the total travel distance and direction. If a particle reaches an element edge during a time-step it is placed in the adjacent element and allowed to continue traveling through the fracture network. A particle is removed i f it reaches a sink node, or becomes stuck in the system. For the purpose of this study, i f a particle becomes stuck at any time during its journey through the fracture network it is removed from the calculation of the total number of particles released. When the breakthrough of particles at the downstream boundary is compared to the initial number of particles released, any particles that become stuck are removed from the equation. This was done because the behaviour of stuck particles was not thought to be indicative of any real physical phenomenon, but rather a result of idiosyncrasies in the particle tracking algorithm. Dershowitz et al. (1991) indicate the problems presented by applying the particle tracking algorithm to a finite-element grid. Prediction of nodal flux and head distributions do not depend on accurate small-scale (1-10 cm) heads and fluxes. However, particle tracking solute transport modeling requires that small-scale heads and fluxes be accurate since the local particle motion is controlled by head differences within single elements. The Galerkin approximation applied in the MAFIC code does not ensure a locally consistent flux field, but only global conservation of mass and a consistent head field. Dershowitz 21 et al (1991) indicated the necessity to refine meshes and condition the particle tracking algorithm to develop accurate transport solutions. In this thesis fracture networks that resulted in large amounts of stuck particles were eliminated from the study to avoid introducing bias through the error in the particle-tracking model. Fracture density and individual fracture properties were chosen from a range that produced acceptable model behaviour and minimized the amount of particles becoming stuck while traversing the fracture network. This limited the scope of the fracture networks that could be simulated, but ensured that the results presented in Chapter 4 were not influenced by this problem. 3.4 Implementation of Abandoned Boreholes To accurately model both flow and transport through a borehole in a fracture network, it is necessary to include the borehole as a discrete feature. M A F I C includes a technique for incorporating a well or borehole as a series of nodes with a specified group flux. If the group flux is set to zero (no pumping or injection) then any inflow to the borehole through a node is compensated by an equivalent outflow through the other nodes. This technique works well for the flow system but does not allow particle tracking through the borehole. Any particle reaching a group-flux node is removed from the transport calculations and is assumed to exit at the borehole. Particles cannot enter at one fracture intersection along the borehole and exit at a different intersection. 22 For this reason it was decided to include the borehole as an equivalent discrete fracture, which would allow for flow and transport without any modifications to the numerical model. Instead of modifying the MAFIC code, which could result in the introduction of programming errors in the code, a separate program was written which takes the finite element mesh output from E D M E S H and adds a borehole as a discrete fracture. 3.4.1 Borehole as a Discrete Fracture The parameters needed to specify a single rectangular fracture are aperture, transmissivity, storativity, width, length, and location. The aperture and the transmissivity are essentially only a single parameter as they are directly related in the numerical model by the cubic law. The numerical code checks to make sure that the aperture and transmissivity are consistent. The desired location and length of the borehole govern the location and length of the corresponding discrete fracture. The storativity is not important as all simulations in this study incorporated a steady-state flow field. We are left with two parameters, aperture and width, which must be mathematically related to the diameter of the borehole to ensure an accurate implementation. The scenario investigated assumes a crystalline fractured rock, so the borehole was assumed to be an open conduit with reasonably smooth sides. If the borehole is assumed to be filled with geologic debris, its effective area and aperture would be less than used here. Figure 3-3 shows the parameters needed to simulate a borehole with an equivalent 23 discrete fracture. The appropriate aperture, b, and fracture width, w, are found by matching the volume and flow rate of a borehole with diameter, D. The volume of a fracture of width, w, aperture, b, and length, 1 is simply: Vf = wbl (6) The volume of a borehole of diameter, D, and length, 1 is: K - ^ - l (7) Matching the two equations and solving for the fracture width gives: w = —— (8) 4b Daugherty et al. (1985) gives the equation for the mean velocity in a circular conduit with smooth walls: -=P8_D**L (9) 32ju dl 24 thus the flow rate through the circular conduit is: 128 \x dl Equation 9 is based on the assumption that flow is laminar through the conduit and that the void space is fully saturated. It is practically impossible for turbulent flow to exist when Reynolds number, R, is less than 2000, because any turbulence that is set up will be damped out by viscous friction (Daugherty et al. 1985). Reynolds number is a dimensionless number with the following form: R = ^ P (11) Rearranging and substituting the values for density and viscosity used in equation 2 gives: (12) D where diameter is in meters and velocity is in meters per second. For a typical borehole diameter of 10cm (0.1m) the critical velocity is 0.02 m/s. The highest velocities modeled 25 in this study were 0.0273 m/s and were possibly in the turbulent range with a Reynolds number of 2730. However, the exact critical R is really indeterminate and Daugherty et al (1985) suggest that its value is normally around 4000. Care must be taken when applying this model to check that flow through the borehole is not likely to be turbulent, which would invalidate the underlying conceptual model. Setting the equation for flow through a fracture (equation 2) equal to the equation for flow through a circular conduit (equation 10) gives: (13) Solving for the aperture, b, gives: b3 = 3K D* (14) 32 w Combining equations 8 and 14 gives: (15) Tt D«1.28Z) (16) w 4 26 So we have a direct relationship that provides the width and aperture of a discrete fracture to correspond to a borehole of specified diameter, D . This relationship ensures consistency in the void space volume, volumetric flow rate, velocity, and contaminant residence times. A small error may be introduced around the borehole entrance by the fact that the shape of the conduit is rectangular instead of circular. However, due to the fine discretization around the borehole, this error is likely to be insignificant compared with other systematic error intrinsic in the finite element model approximation. 3.4.2 Borehole — Fracture Intersections The borehole is added to the finite element model after the mesh is generated and refined. The code for adding the boreholes requires the following parameters for each borehole: Xo(x,y,z), the borehole start location, Xi(x.y.z), the borehole end location, and D , the borehole diameter. The program is capable of including any number of boreholes with any orientation. The boreholes can extend all the way across the domain or terminate at any point, but must be linear. The program takes the input parameters and calculates the three-dimensional equation for the line along which the borehole lies. Each fracture within the modeling domain is checked to see if it intersects with the borehole. If a fracture and a borehole intersect, the 27 fracture's elements are checked to find the element center that lies closest to the intersection point. The intersection point is then relocated to the center of this element. This eliminates problems that occur when a borehole intersects a fracture close to the boundary of two elements. Specifying a reasonable mesh refinement around the borehole location reduces any errors caused by adjusting the true intersection point, and ensures a good mesh design by the borehole-meshing algorithm. For every intersection found the borehole program refines the discretization within the intersected element. Positioning the intersection at the nearest element center reduces the likelihood of having elements with high aspect ratios, which can cause problems in the numerical solver. Figure 3-4 demonstrates the algorithm used to discretize a fracture-borehole intersection. The borehole is discretized as a fracture and the intersection nodes are calculated. The length of the intersection is the fracture width, w. The discretization algorithm ensures a consistent set of nodes and triangular elements are input into the flow solver. The number of elements used to discretize the borehole is governed by the maximum borehole element aspect ratio, which can be specified in the program. 3.43 Test of Flow through Borehole A simple scenario with an analytical solution is used to test the implementation of the borehole in the finite element model. The scenario is shown in the bottom-right corner of Figure 3-5. A single borehole is connected on each end by a constant head boundary. The hydraulic gradient across the borehole is varied for ten runs of the numerical model and the circles show the resulting flow rates. The solid line shows the analytical solution 28 for flow through a filled circular conduit. As expected, the numerical model corresponds almost exactly to the analytical solution, except for some insignificant numerical error. The value of the solution when Reynolds number is equal to 2000 is also plotted. Gradients above this value of ~ 0.00065 m /100 m = 6.5x10"6 could possibly result in turbulent flow through the 10 cm diameter borehole. 3.4.4 Test of Particle Tracking through Borehole Another simple scenario was set up to test particle transport through the borehole. A simple fracture network was generated, and although no analytical solution exists for comparison, the network is simple enough to yield intuitive behaviour. The fracture network for this test scenario and the borehole location are shown in Figure 3-6. The left side of the uppermost fracture intersects a constant-head boundary with hi = 100 m, while the right side of the lowermost fracture intersects a boundary with h 2 = 0 m. A l l other boundaries are impermeable. The borehole is a vertical line located at the center of the domain and intersects all three horizontal fractures. A vertical fracture on the right hand side connects the upper and middle fractures and another vertical fracture on the left-hand side connects the middle and lower features. The contaminant source zone is on the upper-left constant-head boundary. The horizontal square fractures are 20 meters on each side. The borehole is 10 centimeters in diameter. When the borehole is not present contaminants must follow a tortuous path through the fracture network. The pathway travels to the right across the domain in the upper fracture, back to the left in the middle fracture, then back towards the exit boundary at the 29 right-hand side of the lower fracture. When the borehole is present it creates a direct connection from the upper fracture to the lower fracture and particles only have to travel once across the domain to reach the downstream boundary. It is expected that the borehole will increase ground water flow and significantly reduce particle residence times in the network. The upper and lower fractures have a transmissivity of lx l0" 4 m2/s 6 2 while all of the other fractures have a transmissivity of 1x10" m /s. Figure 3-7 shows the normalized breakthrough of particles at the downstream boundary, both with and without the borehole. It can be seen from this plot that the borehole reduces particle residence times significantly. The peak of the breakthrough curve is also higher with the borehole, which corresponds to higher contaminant concentrations. The median travel time for particles is 7150 s without the borehole and 131s with the borehole. The velocity in the borehole is 0.299 m/s and its length between the fractures is 10 m, so the borehole residence time is close to 33 s, or about VA of the median residence time in the system. Figure 3-8 shows the pathway a single particle takes through the fracture network when the borehole is present. It can be seen from this plot that particles are indeed traveling through the borehole as expected. When the borehole is present it creates a shortcut through the network, as well as increasing the overall conductivity. The flow rate through the system increases from 4.93xl0"5 m3/s to 2.37xl0"3 m3/s with the addition of the borehole. The flow pattern in the lower fracture is analogous to that of an injection 30 well in a confined aquifer under a regional gradient. Particle paths are radial away from the borehole but bend in the direction of the downstream boundary. 3.4.5 Effect of Borehole Size It was found that i f a larger diameter borehole is used in the above scenario, the overall hydraulic conductivity of the fracture network is not increased significantly, as the thin fractures account for most of the resistance in the network. However, larger boreholes increase the volume of the void space in the fracture network, and therefore increase residence times. Intuitively one might think that a larger borehole connecting the fracture network will increase both the flow rate and the particle velocities, but the latter is not always the case. The aperture of the high transmissivity fractures in the above example is calculated using Equation 3 to be 5.0xl0"4 m. Likewise the aperture of the lower-transmissivity fractures is l . l x l 0"4 m. Using these values, the total volume of the fracture void space is calculated as 0.49 m . The volume of the 10 cm diameter borehole is 0.016 m , which is small in comparison. The volume of a 25 cm diameter borehole would be 0.08 m which is starting to become significant in comparison to the fracture volume. Smaller boreholes which increase the hydraulic conductivity of the network by creating a new, fast pathway without significantly increasing the system volume will have the greatest impact in decreasing contaminant residence times. 31 Figure 3-9 shows a scenario that is used to demonstrate how a large void-space, such as a borehole, can behave when connected to constant head boundaries through relatively thin fractures. A fracture consisting of three parts has a total length of L. The outer, thin sections both have an aperture of bo and a length of lo. The middle section has a larger aperture of bi and length of \\ = L - 21o- The length of the middle section is varied but the total length remains constant. The analytical solution for this problem was solved and the flow rate and hydraulic residence time calculated for different lengths of the middle section, l i . Figure 3-10 shows a plot of the residence time and flow rate as the length of the middle, large aperture fracture is varied. The parameters for this plot are given in Table 3-1. Description Parameter Value Units Thin fracture aperture bo lxlO" 4 m Thick fracture aperture bi lxlO" 1 m Upstream head Ho 100 m Downstream head Hi 0 m Total Length L 100 m Table 3-1. Parameters for analytical problem It can be seen from Figure 3-10 that the total flow rate through the system increases exponentially as the relative length of the thick portion is increased. The residence time 32 increases with the thick fracture length up to a maximum, then decreases again. On the left-hand side of this plot, when the thick fracture is small, the thin fracture dominates the flow. In this scenario the hydraulic conductivity of the total system is governed by the thin fracture, and the thick fracture acts as a 'reservoir', slowing down flow. On the right-hand side of the plot the thick fracture is dominating the flow, as the thin low-conductivity portion of the system is small. In this scenario the thick fracture acts as a 'pipeline', and increasing its length will significantly increase the conductivity (reducing the( length of the resistive thin fracture) and residence times will decrease. In most of the scenarios investigated in this study the behaviour more closely reflects the left side of Figure 3-10. The borehole is not connected directly to the constant head boundaries, so flow and transport must take place through the lower-conductivity fractures. Increasing the size of the borehole does not increase the hydraulic conductivity of the system significantly as the system is dominated by the background fractures. Of course, i f the borehole were directly connected to the constant-head boundaries the behaviour would more closely resemble the right-hand side of Figure 3-10. 3.5 Conceptual Model Description The simulations presented in the following section all use the same basic conceptual model. A schematic view of a hypothetical mine site is shown in Figure 3-11. A waste rock pile sits above a flow divide at the top of a hill and acts as a source of contaminants to the ground water system. The underlying fractured rock unit is 25 m thick at the left side of the domain, and 20 m thick at the right side of the domain. The simulation 33 domain is 100 m long by 50 m wide. The total volume of the simulation domain is 112500 m 3 . A layer of low-permeability surficial deposits overlies the fractured rock unit downstream of the recharge area. Figure 3-12 is a cross section through the center of the conceptual model that shows the boundaries and the source zone. The domain is surrounded by impermeable boundaries except for the two constant head boundaries on the top-left and right sides of the domain. The constant head boundary at the top-left acts as the source of water for the flow system, and corresponds to recharge at the top of the hill. The boundary on the left is impermeable because it is assumed to be on a groundwater divide. The constant head boundary on the right hand side of the domain acts as the exit boundary for flow and contaminant transport, and is meant to represent a vertical equipotential. A contaminant source zone is situated at the upstream constant head boundary where the waste rock pile contacts the fractured-rock mass. The size of the source zone is representative of a 'hot-spot' in the waste rock pile where a higher than normal flux of contaminants are released, due to preferential pathways through the pile. Contaminants released at the source zone will percolate down through the open fractures and then move horizontally across the fractured crystalline unit to the downstream boundary. The lower boundary corresponds to impermeable bedrock, while the boundary on the top right along the slope is due to low-permeability overburden. 34 It is recognized that the conceptual model presented here is somewhat idealized compared to a real field situation. This conceptual model was chosen to include the key processes that were thought to be important, while minimizing the influence of other, secondary effects. The low-permeability boundary that covers most of the upper surface of the fractured rock unit was necessary to confine the aquifer. When a specified head boundary was used to represent an unconfined aquifer the majority of the particles exited at the water table. The downstream boundary could represent a compliance surface where contaminant concentrations are of the most concern. If the boreholes were allowed to intersect a specified head boundary at the upper surface then they would no longer be acting simply as connections through the network, but as a source of recharge to the fracture network. In this case no particles would travel into the borehole, as the head in the borehole would be higher than that in the surrounding fractures. In many field situations the top of the borehole would be capped or caved in preventing water from entering at the surface so it is realistic to assume that the boreholes do not intersect a constant head boundary. 35 Up Figure 3-1: Definition of orientation convention used in F R A C M A N 36 B 100% of fractures displayed X -s V -z -25% of fractures displayed V -Figure 3-2: Isometric view of the base case fracture network 37 Borehole Fracture D i a m e t e r , D W i d t h , w IN I N A p e r t u r e , b Figure 3-3: Parameters for incorporation of a circular borehole as a rectangular fracture 38 Figure 3-4: Discretization algorithm around the intersection of a borehole with a fracture 39 Figure 3-6: Fracture network designed to test particle-tracking through the borehole 41 o a; o _c CD O m o m CO o + LU o co o + LU O 00 V CO o + LU O CO O + LU o CO CO o + LU O CO o + LU O CO o + LU O CO cu E o c <D -t-> s-t a3 T 3 O co o + LU o Csi CO O + LU O o O CM O CO O CD d CM CO O CO o o CN O O O + LU O » S bxi 42 Figure 3-9: Schematic of analytical model used to test the effect of borehole size on flow rate and residence times 44 ( s / c i u ) a j e y w v o y 47 4.0 Results 4.1 Base Case Fracture Network The base case fracture network consists of two orthogonal fracture sets and two large horizontal fractures. The large fractures represent major features of high conductivity in the horizontal plane, such as sedimentary bedding planes or large fissures caused by erosional unloading. The smaller fractures are randomly sampled from two sets of statistical parameters, and represent the background fracture network. Figure 4-1 shows the fracture traces on a cross-section through the center of the modeling domain. The major features are connected through the background fracture network, but no single background fracture is large enough to directly connect the two major features. Therefore, all flow and transport must take place along some indirect pathway through the network. A single, transmissive fracture that connected the two major features would create a preferential pathway through the network, and produce similar results to including a borehole in the network. The generation parameters for all of the fractures in the base case network are provided in Table 4-1. Parameter Fracture set one Fracture set two Generation Region 125 m x 65 m x 40 m 125 m x 65 m x 40 m Number of Fractures 750 750 Fracture Model Poisson Rectangle Poisson Rectangle Generation Mode Centers Centers Truncation Mode Off Off Number of Sides 4 4 Pole (trace, plunge) 0°, 90° 0°, 0° Pole distribution Constant Constant Fracture dimension l O m x 10m l O m x 10m Transmissivity lxl0" 6 m 2 /s lxl0" 6 m 2 / s Aperture l x l 0 " 4 m l x l 0 " 4 m Table 4-1: Generation parameters for base case fracture network 48 The statistics for both sets of background fractures are identical, except for orientation. The base case network is meant to represent the simplest possible situation, therefore all fractures in a given set have the same orientation, aperture, and size. The center of each fracture is randomly positioned within the fracture generation region. It is possible to sample any of these parameters from probability distributions, but constant values are used in the base case simulation for simplicity and to ease comparisons with later simulations. In the base case scenario, both of the large fractures have a transmissivity of l x l 0"4 m2/s and extend laterally across the entire modeling domain. The transmissivity of the major features is two orders of magnitude higher than that of the background fractures which, combined with their large areal extent, causes them to dominate the flow system. The uppermost feature located at z = 10 m is bounded on all sides by impermeable boundaries, while the lower feature located at z = -8 m intersects a constant head boundary on the right-hand side. The contaminant source zone is located at the top of the domain on the left-hand side as shown in Figure 4-2. For contaminants to reach the exit boundary on the right-hand side they must travel at least partly through the background fracture network to reach the major features. Most of the particles representing the contaminant release will travel partly along the major features and partly along the background fractures to reach the exit boundary, although a small percentage will bypass the major features and travel only in the background network. 4 9 4.1.1 Fracture Generation Region The generation region for the background fractures is larger than the modeling domain to avoid 'boundary effects' within the modeling domain. The generation algorithm randomly positions a fracture center within the generation region. Because fracture centers will not be positioned outside the boundaries, a certain amount of fractures that overlap the boundaries will be missing, resulting in a sparser fracture network near the generation region boundaries. If the generation region is larger than the modeling by the Vi width of the largest fracture generated this error will be eliminated. For this reason the fractures were generated in a 125 x 65 x 40 meter region for the 100 x 50 x 25 meter modeling domain as shown in Figure 4-3. 4.1.2 Base Case Borehole For all of the different fracture networks investigated in this study, the flow and transport model is first solved without any boreholes present, and then the same model is solved with a borehole in a specified location. The results from the two simulations are then compared to determine the effects of the borehole on the fracture network flow and contaminant transport. The vertical borehole for the base case fracture network is located at x = -25 m and y = 0 m. As can be seen in Figure 4-1, the borehole intersects the two major conductive zones at z = 8 m and -10 m respectively and any background fractures in between. Flow can enter or leave the borehole at any intersection, depending on the local hydraulic head distribution. In the base case simulation the flow down the borehole varies only by a 50 small percentage between the major features, representing the limited hydraulic interaction of the borehole with the background network. If the flow through the borehole decreases by 1% after a fracture intersection then there is a 1% chance that particles will leave the borehole at that intersection and enter the background fracture network. There is a 99% chance they will continue traveling down the borehole. The borehole diameter for this simulation was 10 cm. The mathematical algorithm allows for boreholes of any size and orientation to be included in the finite-element model, as long as the finite element grid is designed to ensure an appropriate refinement around any borehole-fracture intersections. As explained in Chapter 3, i f a borehole's diameter is large it may not act as a fast-pathway, but rather can increase contaminant residence times. The large radius borehole will contain a significant volume of water compared to the rest of the fracture network, and will act like a reservoir or holding-tank in the system. This will be the case if the volume of void space in the borehole is large in comparison to the volume of void space in the rest of the fracture network, and the borehole is only connected to the constant head boundaries through the thin background fractures. While a larger borehole radius will raise the overall hydraulic conductivity of the fracture network and therefore the total flow through the fracture network, the volume of water will be much larger and the residence times will increase significantly. Small diameter boreholes may form previously unavailable pathways, while not significantly increasing the volume of void space in the system, thereby decreasing contaminant residence times and increasing 51 maximum concentrations at the downstream boundary. The effect of using different borehole diameters in the base case fracture network is investigated later in this chapter. Figure 4-1 shows that there will be numerous possible paths that will lead from the source to the exit boundary. However, no individual fracture is long enough to connect the major features, and particles must travel along both of the background fracture sets. It should be noted that a two-dimensional trace plot such as this will always appear less connected than the corresponding three-dimensional network as there is another dimension for connections to be made in. This base case fracture network is well connected with many possible paths for particles to travel along to reach the exit boundary. 4.1.3 Grid Discretization Figure 4-4 shows the finite-element grid discretization of the lower large high-conductivity fracture. Note the further grid refinement around the borehole location and the fracture intersections. This refinement was done to reduce model errors where the hydraulic gradients were likely to be large. The grid refinement shown was chosen as the point at which further grid refinement did not lead to significant changes in the flow or particle transport results. This level of discretization resulted in approximately 20000 nodes and 35000 elements for the base-case scenario. Some subsequent simulations with more complicated fracture networks required close to 100000 elements. Available computer memory allowed simulations of up to 300000 elements, but this level of discretization was not found to be necessary for the fracture networks considered here. 52 An increase in the fracture densities over those used here would be approaching the limits of a Pentium II, 450 M H z machine with 196 M B of R A M . A single base-case simulation including fracture generation, discretizing, editing and solving the finite element model takes approximately 15 minutes on the above machine. Some of the more complex simulations took several hours to complete. It was found that approximately 1000 particles were needed to provide a reliable breakthrough curve at the exit boundary for a single realization of the fracture geometry. Certain simulations required significantly more than 1000 particles to be injected at the source due to some particles becoming stuck in the system. Using a constant concentration at the source zone for all simulations would result in a higher number of particles entering the system for fracture networks with high flow rates. For the purposes of this study the amount of mass entering the system for different fracture networks is not as important as the behaviour of the mass once it has entered the network. For this reason the number of particles traveling through the system was kept constant rather than the source concentration, so concentration values shown on breakthrough curves should be interpreted in a relative sense. The longitudinal and transverse dispersivities were set at 1.0 m and 0.707 m respectively. Test simulations found that most of the plume spreading was due to tortuous pathways in the fracture network and not due to the dispersivity. Setting the dispersivities several orders of magnitude lower yielded similar results. 53 4.1.4 Case 1 - Base Case Results Table 4-2 shows a summary of the results of all simulations described in this chapter, including the base case scenario. Q s y s is the total system flow, representing the flow entering through the upstream constant head boundary and leaving through the downstream boundary. V b 0 r is the downward ground water velocity in the borehole, and Qbor is the volumetric flow rate through the borehole. The parameter % Q b 0 r is equal to Qbor / Qsys x 100% and represents the percentage of the total system flow that travels through the borehole. The T parameters are as defined in Chapter 3. It can be seen from Table 4-2 that the total flow through the fracture network increases by a factor of 3.3 from 3.07x10"5 m3/s to 1.02xl0"4 m3/s when the borehole is included, while the mean contaminant breakthrough time decreases by a factor of 7.8. When the borehole is present in the base case scenario it controls the flow system, channeling 87% of the network flow and 98.5% of the particles released at the source. Figure 4-5 shows the breakthrough of particles at the exit boundary for the base case simulation with and without a borehole. Both plots are normalized by dividing the number of particles exiting during a specific time period by the total number of particles travelling through the network. For instance, the spike in the breakthrough curve for the case with a borehole means that 9% of the particles exited during that particular time period. Each time period is 2000 s long. The curve is smoothed using a moving average of the three closest time periods. 54 T3 ii u © o a a Q CS S 3 Vi u u C N C N 0 S. 55 The effects of adding a borehole in this scenario can be quite clearly seen on this plot. The simulation with the borehole has a much higher spike of particles exiting in a short time period, which represents a higher concentration of contaminant than the case without the borehole. Furthermore, contaminants reach the downstream boundary earlier in time for the case with the borehole. In fact, most of the particles have already exited by t = l x l 0 5 s in the case with the borehole, whereas in the case without the borehole the first few particles are just beginning to reach the downstream boundary at this time. The parameters used to quantify the shape of the breakthrough curve are Tio, T5o and T9o. Tio is the time when 10% of the particles have exited at the downstream end of the modeling domain, and represents the nose of the breakthrough curve. T50 is the time when half of the particles have exited and is identical to the median of the particle breakthrough times. T90 is the time when 90% of the particles have exited and represents the tail of the breakthrough curve. The difference between T90 and Tio quantifies the spread of the breakthrough curve. A lower spread would correspond to a more tightly contained plume, and higher contaminant concentrations at the downstream boundary, as is seen in the base case scenario with the borehole. The parameters T25, T 7 5 , and T95 are also included to give a complete description of the breakthrough curve. Table 4-2 gives the values of the parameters described above for all simulations, including the base case. Figure 4-6 is a plot of the exit location of each particle against its exit time. The vast majority of the particles reach the downstream boundary while travelling along the lower major conduit, and exit at the same longitudinal and vertical location (x = 50m, z = -10 m). For this reason the plot shows the transverse, or y-component of the exit location vs. 56 time. The time axis is logarithmic. This plot shows the temporal and transverse spread of the particles at breakthrough, and where along the y-axis particles show up first. The most apparent difference between the two scenarios is the order-of-magnitude longer time for particles to travel through the fracture network when no borehole is present. Another interesting feature that could not be seen on the breakthrough plots is the structure of the contaminant plume when a borehole is present. It can be seen from this plot that the borehole controls the flow system in the large fracture, and that the flow lines in the fracture are similar to those that occur in a confined aquifer with an injection well. Particles close to the center of the domain tend to exit first as they follow the most direct pathway from the borehole to the exit boundary. Particles that exit further from the center of the domain take longer to reach the boundary, as they have to travel along a longer arc, even backwards for part of the journey. The particles spread laterally almost across the entire domain. The impermeable boundaries at either side of the domain interfere with the flow system and limit the spread of particles. These boundaries represent undisturbed flow lines in the regional flow system, and while it would be desirable to have the boundaries as far away as possible, computational requirements limit the size of the domain. One interesting feature to note for the case with the borehole is that the temporal spread of the plume is even less than depicted on the breakthrough curve. The breakthrough curve for any given lateral location will encompass a smaller time spread than when the breakthrough is averaged across the entire y-axis. This means that the breakthrough 57 curve is underestimating the maximum concentrations through averaging along the lateral direction, and concentrations for the case with the borehole will be even higher than already shown. Figure 4-7 shows the pathway a single particle takes from the source zone to the exit boundary. After the particle exits the borehole it travels backwards along the x-axis to the upstream impermeable boundary. As it approaches the boundary the ground water velocity carries it outwards in the negative y direction. Near the lateral boundaries the velocity is in the positive x direction which carries the particle towards the exit boundary again. The particle exists at y = -24.1 m. Particles exiting towards the outskirts of the domain must travel along a much longer path resulting in a later exit-time for both the first and last particles to reach the boundary. Particles that leave the borehole on the downstream side travel directly to the downstream boundary under a strong gradient, and exit the system at the earliest time. When a particle leaves an element at any intersection, the probability of it entering any other intersecting element is based on the proportion of flow to that element. Stream-tube routing is not incorporated in the M A F I C flow and transport code. Some clustering can be observed in the lateral exit-location of particles for the case without the borehole. This is due to heterogeneity in the discrete fracture network caused by the random distribution of fractures in the domain. Although all fractures have the same equivalent radius and transmissivity, areas with greater than average fracture densities will result in higher local permeability and preferential flow paths. 58 Figure 4-8 shows a frame-by-frame cross-section of particle transport through the fracture network when no borehole is present. It can be seen from this plot that there are two major pathways through the network along with many secondary ones. One concentrated packet of particles travels almost straight down through the background network until it meets the lower major feature. At this point it moves horizontally towards the exit boundary. Another concentrated packet of particles travels vertically until it reaches the upper major feature, and then horizontally across the feature until the point where the feature intersects the upper impermeable boundary. At this point the particles again travel downwards through the background fracture network to the lower major feature and out through the exit boundary. The effect of the two major pathways can be seen on Figure 4-5. The breakthrough curve contains two peaks corresponding to the two distinct paths. When the borehole is present the particles are concentrated along a single preferential pathway through the fracture network. As mentioned above, 98.5% of the particles travel along this pathway through the borehole in the base case network. This will result in less of the fracture network being exposed to contaminants and will greatly reduce the fracture surface area contacted by the plume. Although sorption is not included in this analysis, the implications of this behaviour on the mobility of sorbing contaminants can be inferred. Reducing the available fracture surface area will reduce the amount of sorption, increasing the mobility of sorbing contaminants and reducing residence times further when the borehole is present. 59 The following set of simulations investigates the effect of changing the fracture network transmissivities, while keeping the geometrical properties the same. To facilitate comparison with the base case scenario the same random-number seed was used to generate the fracture network, resulting in the exact same locations for the fractures. 4.2 Influence of Background Fracture Network 4.2.1 Case 2 — Higher horizontal transmissivity The next simulation investigates the effect of anisotropy in the background fracture network. In this case the transmissivity of the horizontal fracture set was raised by two orders of magnitude to l x l 0"4 m2/s (corresponding to an aperture of 5xl0" 4 m), while the vertical set was left at 1x10" m /s. The transmissivity of the two large fractures was kept at lx l0" 4 m2/s. This alteration to the base case scenario represents an introduction of anisotropy into the system. A typical field scenario that corresponds to this fracture network is horizontal fractures along consolidated sedimentary bedding planes, connected by tighter vertical fractures. Other than the above modifications the parameters for this simulation are the same as provided in Table 4-1. This alteration to the background fracture network raised the overall permeability of the geologic unit and causes the flow through the system to increase by 64% over the base case, to 5.04x10"5 m3/s when no borehole is present. When the borehole is present at x = -25 m the flow through the system is LlOxlO" 4 m3/s. The total flow through the system is higher both with the borehole and without compared to the base case. However, the 60 relative impact of adding the borehole in this scenario is reduced. In the base case the borehole increased the total system flow by a factor of 3.3, whereas in this scenario the borehole increased the system flow by a factor of 2.2. The proportion of the total system flow that channeled through the borehole was also less in this case at 77%. It is expected that a vertical borehole would have a more significant effect where the horizontal fracture transmissivities are higher than the vertical. A highly conductive vertical conduit should be more effective in this situation because with low vertical fracture transmissivities the competing preferential flow paths through the background network are at a disadvantage. Any connected flow-path through the background network must include some vertical low-transmissivity fractures. The reason that the borehole has less of a relative effect on the system flow than the base case is not the anisotropy but the overall increase in background network conductivity. The significance of the borehole relies upon the fact that it connects two relatively high transmissivity fractures and creates a preferential pathway that did not exist before. The greater the contrast between the background transmissivity and the major features the greater an impact the borehole can be expected to have. This relationship is investigated later in this chapter. Figure 4-9 shows the breakthrough of particles at the downstream boundary for Case 2. The breakthrough time parameters are provided in Table 4-2. Generally, the behaviour of these curves closely resembles the breakthrough curves for the base case. The borehole has the effect of allowing contaminants to reach the downstream boundary sooner and at 61 much higher concentrations. The breakthrough curve for the simulation without the borehole shows that the particle breakthrough is dispersed over a much longer time frame, but with lower concentrations. More specifically, although the borehole was shown to have less of an effect on the system flow in this simulation vs. the base case, the borehole has a greater effect on the contaminant breakthrough as measured by the median residence time, T50. T50 decreases by a factor of 8.5 in Case 2 with the inclusion of the borehole, but only by a factor of 7.8 in the base case. This result is more significant when taken in conjunction with the borehole's decreased effect on the flow in Case 2. The reason for this increased borehole effect on transport in Case 2 is the longer residence times for Case 2 without the borehole. It can be seen from Table 4-2 that Q s y s and the breakthrough times are similar for Case 2 and the base case when the borehole is present. This is expected, because the higher transmissivity horizontal fractures will not have a large effect when the borehole is present. The majority of the system flow will take place through the borehole, the major features, and through the vertical background fractures from the source zone to the upper major feature. 97.5% of the particles released from the source zone travel along this path through the borehole. However, except for the initial breakthrough, the breakthrough times are longer for Case 2 when no borehole is present. This is due to the overall volume increase caused by the larger aperture horizontal fractures being introduced. Increasing the fracture transmissivity by two orders of magnitude results in an increase in fracture volume of 4.64 according to the cubic law (Tocb ). The total system conductivity increased by a factor of 1.6 as 62 evidenced by the increase in the flow rate in Case 2, which was not enough to account for the higher volume of the flow paths, hence the residence times increased. The fact that Tio and T25 are smaller for Case 2 indicates that the fastest pathways are faster in this scenario, although the majority of the pathways are slower. The lateral breakthrough distribution for Case 2 is shown in Figure 4-10. The lateral pattern of breakthrough for the case with the borehole is once again similar to the base case, with the average particle exit time increasing with distance from the center of the domain. This is due to the borehole being located in the center of the domain. In later sections the borehole is located off-center and the effect on the lateral breakthrough distribution is investigated. Without the borehole the lateral breakthrough distribution is not structured, but rather is spread fairly evenly across the domain and across an extended time period. Increasing the transmissivity of the horizontal background fracture set had the effect of increasing the overall conductivity of the system. The borehole had a smaller effect on the system flow, but a larger effect on contaminant transport than in the base case. This is due to the complex interplay between fracture aperture, network conductivity, preferential pathways and contaminant residence times. The next section expands upon the concept of anisotropy in the background fracture network. 63 4.2.2 Case 3 - Higher Vertical Transmissivity In this case the transmissivity of the vertical fracture set was raised by two orders of magnitude over the base case to lx l0" 4 m2/s, while the horizontal set was left at lxlO" 6 2 * 4 2 m /s. The transmissivity of the two large fractures was maintained at 1x10" m /s. Increasing the transmissivity of the vertical fracture set had a larger effect on the flow system than increasing the horizontal transmissivities in Case 2. Without the borehole, the total flow through the system increased by a factor of 4.8 to 1.47x10"4 m3/s. The reason for this is that in the base case the vertical background fracture set represents the longest constriction in the major flow-path. Most of the horizontal distance from the source zone to the downstream constant-head boundary can be traversed along one of the major conductive features. To reach the lower feature and the exit boundary, however, water must travel along the vertical background fractures from the source zone to the upper feature and again between the upper and lower features. These vertical fractures act as the bottleneck in the base case, and increasing their transmissivities by two orders of magnitude opens this bottleneck significantly. When the vertical borehole was included in the system the total flow was 2.24x10"4m3/s, or 2.2 times the base case scenario with the borehole. The overall effect of the borehole on the flow rate in Case 3 was therefore an increase by a factor of 1.5 (relative to Case 3 with no borehole), which means the borehole had a less significant effect here than in Case 2 or the base case. The percentage of the flow through the borehole was also less in this case at 55%. The reason for this is the overall increase in the hydraulic conductivity 64 of the background fracture network, which makes the preferential pathway created by the borehole less significant. It can be seen by Figure 4-11 though, that the borehole is still causing a large spike in the particle breakthrough curve early in time. In fact, in Table 4-2 it can be seen that the early breakthrough with the borehole is even faster than that of Case 1. However, the effect of the borehole on contaminant transport, measured by the change in the median particle residence time, is only a factor of 3.8. This is due to faster particle transport for Case 3 without the borehole compared to the base case without the borehole. The proportion of particles travelling through the borehole is reduced to 89%, because in this situation there are several conductive alternate pathways. This effects the tail of the breakthrough curve for the case with the borehole, as particles traveling through the background network pathways still take longer than those traveling through the borehole. Figure 4-12 shows the lateral breakthrough distribution for Case 3. On this plot it is evident that the borehole is having a less significant impact on contaminant transport than in the base case. For the case with the borehole the lateral distribution is less structured due to the larger percentage of particles that bypass the borehole and travel through the background fracture network. There is an overlap of particle exit times for the simulations with and without the borehole that was not seen in the base case simulations. Even the particles that do travel through the borehole exhibit a greater spread in their residence times in Case 3. 65 From the previous two simulations a number of observations can be made. It is obvious that a borehole will have a much more significant effect on the flow system when there is less choice of alternative fast-pathways for the flow and transport. For a vertical borehole, the greatest impact will occur when the background network has a lower vertical hydraulic conductivity, as the borehole will create a fast vertical pathway where none existed previously. As was observed in Case 2, increasing the transmissivity of a fracture set does not necessarily result in faster contaminant transport, although it will result in a greater hydraulic conductivity and hence a greater flow rate. The important consideration is whether any restrictions along the preferential pathways are increased. When a borehole was present, increasing the vertical transmissivities resulted in earlier breakthrough and higher concentrations. The vertical fractures were the bottleneck in the preferential pathway along the section between the source area and the upper major feature. Increasing the horizontal transmissivities did not have much of an effect when the borehole was present because the main pathway through the system does not extensively sample the horizontal background set. 4.3 Influence of Major Features The next set of simulations investigates the effect of changing the properties of the two main features located at z = 10 m and z = -8 m for the same background network considered in the base case. The first two cases look at the effect of increasing and decreasing the transmissivity of the major fractures, to see how the influence of the borehole depends on these pathways. The third case investigates the effect of a high contrast between the major features and the background by raising the transmissivity of 66 the major features while lowering that of the background fractures. The fourth case looks at the influence of the abandoned borehole on the base case scenario when the lower major feature does not extend all the way to the downstream constant-head boundary. 3 2 4.3.1 Case 4 - Major feature transmissivity = 1x10'm fs In Case 4 the fracture network is identical to the base case except that the transmissivity of the two major features was increased to lxlO" 3 m2/s. Increasing the transmissivity by one order of magnitude corresponds to increasing the fracture aperture, and volume, by a factor of 2.15. The major features are only connected through the low transmissivity background fractures, so shorter contaminant residence times will only be expected if the increase in flow rate is large enough to counteract the increased volume of void space. Table 4-2 shows the flow rate through the system has increased by 7% over the base case to 3.29xl0"5 m3/s when no borehole is present. When the borehole is present the system flow has increased by 214% to 2.18xl0"4 m3/s. Increasing the transmissivity of the major fractures has a much larger effect when a borehole is present to connect them together. In this case there is only a short length of flow that must take place through the background fracture network, and the increased hydraulic conductivity can counteract the restriction in the flow-path caused by the background network. This behaviour results in the borehole increasing the flow rate through the system by a factor of 6.6 in Case 4, in comparison to a factor of 3.3 for the base case. 67 The contaminant breakthrough curves for Case 4 are given in Figure 4-13. The breakthrough curve for the case with the borehole reaches a higher peak and occurs earlier in time than in the base case. The time difference can be seen more clearly looking at the Tio and T5o values in Table 4-2. The breakthrough curve without the borehole occurs later in time than in the base case, leading to an overall greater contrast between the borehole and no-borehole simulations for Case 4. For the no-borehole case contaminant residence times are longer because the increased hydraulic conductivity from the more transmissive major features does not overcome the larger volume of the void space. The background fractures connecting the major features constrict the flow. When a borehole is present the conductivity is much greater due to the new pathway it forms between the major features. This causes contaminant travel-times to be shorter than in the base case. The borehole controls both the flow and transport through the system when it is present, channeling 99% of the system flow and 98.6% of the particles. The lateral breakthrough distribution for Case 4 given in Figure 4-14 shows the extent that the borehole controls the shape of the plume at breakthrough. The breakthrough pattern is very structured and reflects the flow-paths coming out of the borehole in the lower major fracture. Because the borehole contains nearly all of the flow in the system, there is limited interference in the pattern due to localized flow from background fractures into the major fracture. 68 4.3.2 Case 5 - Threshold effect In this case the transmissivity of the major features was lowered from the base case to the threshold value where they just begin to have an effect on the system. When the major fracture transmissivity was set at the same value as the background fractures there was no flow through the borehole and all particles traveled down one of the alternate pathways to the exit boundary. The transmissivity was subsequently increased until the borehole began to conduct flow. The threshold value for thus fracture geometry was T = 1.86xl0"6 m2/s. At this point the total flow through the system was found to be 8.38x10" m Is without the borehole, and 8.63x10"6 m3/s with the borehole; a difference of 3%. The flow through the borehole accounted for 13% of the total system flow. It can be seen from Figure 4-15 that there is no appreciable difference in the breakthrough curve when the borehole is included in the system. No particles travel through the borehole because it does not create a more preferable pathway to those already existing in the background fracture network. Even though the major features that the borehole connects extend across the entire domain and have a transmissivity that is 86% higher than the rest of the fractures, the network is already well-connected and the pathway created by the borehole is not transmissive enough to channel the particles. It is interesting to note that, when no borehole is present, the contaminant residence times were longer than the base case both when the major-features' transmissivity was raised 69 in Case 4 and lowered in Case 5. This implies that, for this realization of the fracture network geometry, there is a single transmissivity value that will maximise the flow through the fracture network. The flow rate increases with increasing major feature transmissivity. This is the same behaviour that was demonstrated in the simple analytical scenario presented in Chapter 3. 4.3.3 Case 6 - Lower feature terminates atx = 30 m In this case the large feature at z = -8 m no longer extends across the entire modeling domain but extends longitudinally only from x = -50 m to x = 30 m. This lower feature no longer intersects the constant head boundary at x = 50 m so that now both major features are hydraulically connected to the boundaries only through the background network. The base case background fracture network is used in this simulation, and the transmissivity of the major features is lxl0" 4 m 3 /s . There are now restrictions in the flow-path at both ends when the borehole is present which would be expected to reduce its effect on the flow and transport. If the upper major feature was connected to the upstream boundary and the lower major feature was connected to the downstream boundary then adding the borehole would have the largest effect, creating a direct high T pathway all the way from the source to the sink. The greater the proportion of the flow path that is along low T fractures, the less effective the borehole can be in increasing the overall conductivity, and the more the volume of the borehole is going to act to lengthen residence times. 70 The flow and transport parameters from this simulation are given in Table 4-2. The flow through the system is less than the base case both with and without the borehole. The relative effect of the borehole on the system flow is also less, increasing the flow by 13% compared with 330% for the base case. The borehole still controls 77% of the flow in this system. The borehole makes less of an impact in Case 6 than in the base case, but it still controls the fracture flow when it is present. The breakthrough of contaminants occurs much later in time when compared to the base case. T50 for Case 6 is 9.22x105 s compared with 2.48x105 s for the base case. As can be seen in Figure 4-16, the borehole still has the effect of shortening the particle residence times in the fracture network, but the breakthrough is more dispersed in time and the concentrations are lower. The borehole had the overall effect of accelerating the median breakthrough by a factor of 1.7 in Case 6, compared with a factor of 7.8 in the base case. This case is representative of any situation in which two or more relatively large features are connected by a borehole within the fracture network. When the features do not form part of any high conductivity pathway to the source or sink the borehole has some effect on contaminant transport, but it is substantially diminished. The plot of the lateral breakthrough distribution for Case 6 is given in Figure 4-17. The particles have to find their way to the exit boundary through the background network along discrete fractures, making the lateral distribution of particle exit locations more clustered. In the base case the locations were spread across the entire width of the domain, making the plume easier to detect. In this case clustering leaves gaps in the 71 plume, which could allow contaminants to travel past a discrete monitoring location undetected. The particles exit in similar locations for the case with the borehole and without, but the borehole causes the particles to exit earlier in time. 4.4 Influence of Borehole Diameter and Location The next group of cases considers the effect of changing the borehole properties for the base case scenario. The effect of both the borehole's location and its diameter on the flow and transport through the fracture network are investigated. Because the fracture network used in this section is identical to the base case, the flow and transport properties without the borehole are the same as the base case. For this reason only the results with the borehole are compared to the base case scenario when the borehole was located at x = - 25 m, y = 0 m with a diameter of 10 cm. 4.4.1 Case 7 - Borehole atx = 0m In the base case scenario the borehole was located at x = -25 m which is only 5 meters downstream from the upper constant-head boundary and 15 meters downstream of the contaminant source zone. In this case the borehole was moved 25 meters further downstream from the flow and contaminant source zones to x = 0 m. The results of the simulation are summarized in Table 4-2. Positioning the borehole further downstream had very little effect on the groundwater flow through the fracture network. The total flow through the system and the flow through the borehole are essentially the same as the base case. 72 Figure 4-18 shows the breakthrough of contaminants at the downstream boundary for Case 7 plotted against the base case results. It can be seen from this plot that there is no significant difference in the temporal breakthrough across the downstream boundary when the borehole is located farther away from the source zone. The geometry of the system means that there is no change in the length or conductivity of the main pathway when the borehole is moved downstream along the centerline. There is a slightly lower peak concentration with the borehole at x = 0 m and a longer tail to the breakthrough curve. Even though the proportion of flow through the borehole is the same as the base case, a higher percentage of particles bypass the borehole and travel through the background network when it is located further away. In this case 93.9% of the particles traveled through the borehole compared with 98.5% in the base case. Figure 4-19 shows the lateral breakthrough distribution for Case 7. The main difference between this plot and the base case is that the particles are spread out more through time when the borehole is present. The front of the plume reaches the downstream boundary at the same time, but the tail is more spread out along the entire width of the domain. 4.4.2 Case 8 - Borehole atx = -35 m This case investigates the effect of moving the borehole 10m closer to the contaminant source zone to x = -35 m. The constant head boundary along the top of the domain extends from x = -50 to -30 m, so in this case the top of the borehole is actually intersecting the constant-head boundary. Based on the results of moving the borehole further away in the previous case, simply moving the borehole closer would only be 73 expected to raise peak concentrations somewhat, and not effect breakthrough times. However, because the high conductivity conduit now intersects a constant head boundary, the flow system and contaminant transport are significantly altered. In this simulation 92% of the groundwater flow enters the domain through the borehole. The borehole connects the upper constant head boundary to the lower major feature, which in turn is connected to the downstream constant head boundary. This creates a high conductivity pathway through the system, significantly increasing the total system flow in this case. The system flow with the borehole is 2.26xl0"4 m3/s for this scenario compared to 1.02xl0"4 m3/s for the base case. In previous examples the borehole acted as a connection between two fractures, so that flow entered the borehole at one fracture intersection and traveled under a constant gradient to another fracture intersection. In the present scenario, the borehole is connected to a constant-head boundary, and acts as a source of groundwater to all of the surrounding intersecting fractures. The hydraulic head in the borehole is higher than that in the surrounding fracture elements, therefore flow is out of the borehole at all fracture intersections and no particles can enter the borehole. If the borehole intersected the contaminant source zone then there would be a large influx of fast moving contaminants traveling along the major flow path, and reaching the downstream boundary at high concentrations and at an early time. This would represent a worse case scenario but is unlikely to occur in practice. Most contaminant source areas, such as a tailings pond at a 74 mine site, would be lined with some form of an impermeable barrier, and any known boreholes in this area would be properly sealed from groundwater flow. The contaminant breakthrough curve for Case 8 is provided in Figure 4-20. The borehole in this case actually slows down contaminant transport, as contaminants are pushed away from the fast pathway due to the higher hydraulic head in the borehole. Particles must take a longer and more tortuous route through the background network after they are forced backwards and to the lateral edges of the domain. The lateral breakthrough distribution in Figure 4-21 shows evidence of this. When the borehole is present the majority of particles are still located close to the domain edge, even though the downstream boundary is 85 meters from the borehole location. 4.4.3 Case 9 - Borehole aty = 15 m In this case the effect of locating the borehole off of the x-axis is investigated. The parameters for this simulation are the same as the base case, except for the borehole location which is now at x = -25 m y = 15 m. The borehole still extends vertically across the entire modeling domain. Figure 4-22 shows the breakthrough of particles in this scenario averaged across the entire outflow boundary. Comparing this curve to the breakthrough with the borehole at y = 0 m shows that the temporal breakthrough is similar for both cases. The main difference caused by the off-center borehole location is shown in the lateral breakthrough distribution plot, Figure 4-23. It can be seen in this plot that the borehole's position has 75 caused the plume to deflect laterally towards the positive y direction. The earliest breakthrough no longer occurs at x = 0 m but approximately at x = 15 m, the lateral position of the borehole. The impermeable boundary located at y = 25 m is interfering with the flow system and preventing particles from moving further away from the x-axis, and preventing the true effect of the borehole location being demonstrated. It is clear from this plot that the borehole controls the contaminant plume exit location, and in this scenario an effective monitoring network design would need to take into account the influence of the borehole. 96.4% of the particles travel through the borehole, which is less than the base case scenario, but still enough to influence the downstream structure of the plume. The behaviour described above has implications for the design of detection monitoring networks, which are discussed in Chapter 5. 4.4.4 Case 10 - Influence of Borehole Diameter The next set of simulations test the influence of the borehole diameter on the contaminant transport through the base case network. The borehole is located at x = -25 m, as in the base case. The diameter was set at 2 cm, 5 cm, 10 cm (base case) and 20 cm, and the flow and transport modeled. A summary of the results of the individual simulations is given in Table 4-2. The contaminant breakthrough curves for the different borehole diameters are given in Figure 4-24. It can be seen from this plot that the larger the borehole, the longer contaminants take to reach the downstream boundary. The peak concentration varies between the different runs, but is likely related more to the random variations between 76 runs due to the random dispersion process, than to any significant physical phenomenon. The peak concentration and the tail of the breakthrough curve are related to the proportion of particles that travel through the borehole, whereas the nose of the breakthrough curve depends more on the borehole properties. Figure 4-25 shows the total system flow rate and the mean contaminant breakthrough time plotted against the borehole diameter. The flow rate varies linearly with the borehole diameter over this range, which means there is a linear increase in the total hydraulic conductivity of the system. The interesting behaviour shown on this plot is that the median contaminant residence time also increases with borehole diameter. Intuitively one would expect that a more conductive system would result in shorter residence times. However, while the residence time is proportional to the flow rate, it is inversely proportional to the conducting volume. For the geometry considered here, the increased volume of a larger borehole outweighs the increased conductivity, resulting in longer residence times. In effect, when there are constrictions in the main flow path a larger borehole will slow down contaminant transport and act more like a holding tank than a fast conduit. 4.5 Case 11 - Variable Network with No Major Features In previous studies of the effect of abandoned boreholes in porous media the borehole connected two aquifers across a relatively impermeable aquitard (Lacombe et al., 1995; Avci , 1994). A hydraulic gradient was established across the aquitard due to either a pumping or injection well or natural hydraulic head differences between the two aquifers. 77 Contaminants t raveled f r om the contaminated aquifer, across the borehole , and into the prev ious ly uncontaminated aqui fer due to the hydrau l ic gradient. In the present study, the in f luence o f the borehole on contaminant transport was not cont ro l led b y an extract ion w e l l , but rather was based on the connect ions and preferent ial pathways created in the network under a un i f o rm regional gradient. So far, the cases invest igated i n this study have incorporated large determinist ic features i n the fracture geometry that create a s ingle preferent ial pathway w h e n a ver t ica l borehole is present. The present case demonstrates that abandoned boreholes can in f luence f l o w and contaminant transport i n an ent irely stat ist ical ly generated fracture network. The rock mass is stat ist ical ly homogeneous, but w i l l have loca l areas o f h igher fracture density and permeab i l i t y due to random fracture c luster ing and t ransmiss iv i ty d ist r ibut ion. The stat ist ical parameters used to generate the fracture network for Case 11 are g iven in Tab le 4 -3 . Parameter Fracture set one Fracture set two Generat ion R e g i o n 125 m x 65 m x 40 m 125 m x 65 m x 40 m N u m b e r o f Fractures 300 500 Fracture M o d e l Enhance Baecher Enhanced Baecher Generat ion M o d e Centers Centers Truncat ion M o d e O f f O f f N u m b e r o f S ides 6 6 Po le (trace, p lunge) 0 , 9 0 0 , - 1 0 P o l e d is t r ibut ion B ivar ia te N o r m a l B ivar ia te N o r m a l D is t : K l , K 2 , K 3 10, 1 0 , 0 10, 1 0 , 0 S ize D is t r ibu t ion L o g n o r m a l L o g n o r m a l M e a n Rad ius 8 m 5 m Std. D e v . Rad ius 5 1 Aspec t Ra t io 1 1 Termina t ion % 0 0 Transmiss iv i t y D is t . L o g n o r m a l L o g n o r m a l M e a n Transmiss iv i t y l x l 0 " 5 m 2 / s l x l 0 " 6 m 2 / s Std. dev. T ransmiss iv i t y 5 x l 0 " 5 5 x l 0 " 7 Table 4-3: Generat ion parameters for Case 11 fracture network 78 There are two main fracture sets, both of which have a bivariate-normal distribution for orientation and a lognormal distribution for equivalent radius and transmissivity. Set 1 is a sub-vertical set consisting of more, smaller fractures with lower transmissivity. Set 2 is a sub-horizontal set consisting of less, larger fractures with higher transmissivity. Figure 4-26 shows three cross-sections through the fracture network at y = 10 m, 0 m, and -10 m respectively. *» The flow and transport through the network was solved for the fracture network alone and for the case when a total of eight boreholes were present in the domain. The locations and sizes of the boreholes are given in Table 4-4. The boreholes have various sizes, locations and orientations representing a range of possible exploration holes. None of the boreholes penetrate the upper constant head boundary or the source zone. In previous cases a single borehole had influential behaviour because the major features extended across the entire domain. As was demonstrated in Case 7 the position of the borehole was not as important as the fact that it connected the two major features somewhere within the domain. In the present case the fractures are randomly located within the domain, so the location of the main conducting features are not know a priori, i f they exist. The influence of a single borehole in this situation would change significantly with its position depending on whether or not it happened to form a dominant hydraulic pathway. For this reason a number of boreholes were used in this case to sample a wider distribution of the fracture network. Figure 4-27 shows a plot of the borehole distribution and orientations in an isometric view of the modeling domain. 79 Number X 0 (x , y, z) (m) x i ( x > y>z) ( m ) Orientation (Tr, PI) Diameter (cm) 1 -25, 5, 50 -25, 5, -50 0°, 90° 5.0 2 -15, 10,50 -15, 10, -50 0°, 90° 4.0 3 15,-5, 50 15,-5,-50 0°, 90° 5.0 4 0, 0, 50 0, 0, -50 0°, 90° 6.0 5 -29, 4, 20 -10, -5,-20 205°, 62° 5.0 6 -15,-10, 50 -5,-10, -50 180°, 87° 5.0 7 35,0, 50 35,0, -50 0°, 90° 4.0 8 0,-15, 50 10,-10, -50 153°, 84° 7.0 able 4-4: Borehole parameters for Case 11 Figure 4-28 shows the particle breakthrough curve at the downstream boundary. It is obvious from this plot that the boreholes are having an influence, even when the discrete major features of the previous simulations are not present. The boreholes have increased the flow rate through the network from 1.16x10" m /s to 1.38x10" m/s. The median residence time has decreased from 4.84xl0 5 s to 3.53xl0 5 s. The maximum concentration at the downstream boundary is higher for this fracture network when the boreholes are present, however in previous cases the influence of the borehole was much more pronounced. In previous simulations the borehole created a single preferential pathway through the fracture network that dominated the flow. In the present case, however, there are a number of alternate pathways through the eight boreholes and the 80 background fracture network, which creates a dispersive effect that lowers the overall peak concentration. The variability of both fracture size and transmissivity creates a fracture network with preferential pathways and flow channeling. The larger, more transmissive horizontal fractures provide the potential for the sub-vertical boreholes to make preferable connections. When the same simulation was solved using only a single vertical borehole at x = -25 m, y = 0 m there was no noticeable change from the case without the borehole. It is not the number of boreholes in the domain that is important in determining the flow and transport behaviour, but rather i f any individual borehole forms a significant preferential pathway. The large number of boreholes increases the probability of a connection being made. A second situation was investigated to test whether the fracture sizes and locations were more important than aperture variability in creating preferential pathways with the boreholes present. A n identical realization of the fracture geometry was used, but the fracture transmissivity was set as a constant. To ease comparison, the transmissivity value was adjusted through trial and error until the flow rate through the network without the boreholes was the same as for the case with a variable transmissivity. The constant 6 2 value of transmissivity to achieve this was found to be 2.85x10" m /s. The network was solved both without any abandoned boreholes and with the same borehole locations as given in Table 4-4. In this case the presence of the boreholes increased the flow rate from 1.16x10" m/s to 1.70x10" m/s. Figure 4-29 shows the breakthrough curves at the 8 1 downstream boundary both with and without the abandoned boreholes. It is clear from this plot that the boreholes are no longer having any effect on the transport of contaminants through the fracture network. The breakthrough curves are essentially identical except for small random fluctuations. For the same fracture geometry in the network, the abandoned boreholes only had an influence when aperture variability was included. It should be noted that the fracture network in these simulations is well connected. Abandoned boreholes are expected to have an important influence in sparse fracture networks based on the fracture geometry alone as they can create connections through the network where non existed before. In a relatively dense fracture network aperture variability is important because the network is already well-connected. This simulation demonstrates the flexibility of the modeling approach, and the fact that abandoned boreholes can have an influence in random fracture networks when the right conditions exist. It would be possible to gather probability data on the impact of one or several boreholes for a given fracture geometry by performing a Monte Carlo simulation using different realizations of the fracture statistics and systematically varying the borehole location. This information would be useful when assessing the potential impact of abandoned boreholes at a real world site where the fracture geometry is known, within some uncertainty, but the borehole locations are not. Such an analysis would require a significant computational effort and therefore a large project budget to be feasible, but would be necessary to analyse the impact of abandoned boreholes in a decision analysis framework. 82 83 on CD 'S c CD R e l <u g c R CL> T3 cn on CD on O s-< j i CD & .2 s-i O e l o o '5b ca e l o '3 S-l CL> e l CD CD CD CD a on o ccj S - l a o '5b CD l-i 60 _ej O e l o '5b CD e l o ta t-l CD e l CD 60 CD r>-l T3 3 cS a o - d &o e l •»—i "5 T3 O B <+-< o 60 ccj CD I. 3 85 00 o X u s-l on eg o o 1 on U CD O 00 O o o ID o o co o o o o sapi^JBd # pazjjeujJON S-9 87 88 x - a x i s Figure 4-7: Pathway of a particle through a borehole in the base case fracture network 89 ; j i i ' J | | | ]_ _j ! , = £ o I t [ i * I j j } [ ? ! i I ' s p 5 4 • ; * « a — : ; : j ] -% i ! i 1 ;J : . i ! i i j* j " 4 ^ fa |^ : ; i i ! : i :; .. * 1 | ;! 1 .! j i • % ' = 2 8 8 i i c CD Vi CD UH P H Vi CD "o CD t-i o x> o c a CD I M s-c o CD c CD iri c3 CD CD Vi cd ,£> CD •B Xi 60 o M l .s u 1 OH O on -*-» O Xi Vi m 90 i Tt CD S-3 91 92 93 o + o o o o o o o c ! CO CM T - T - CM CO 1 -i i i (Ul) UOtteOO - ! |BJ8}B-| U) 94 CD O + LU O i -O + LU O LO o + LU o CO 1 05 O d co o ci o CD o LO o o co o CM o o ci LO o + LU O LO o + LU O CD i i g o c LO H LO o + LU O LO o + LU O LO O + LU O CM LO O + LU O O O + LU O O Sd|3jued # pazjieuuoN o O •a O i t-l Vi U u S bl) 95 96 97 98 o + LU 99 100 101 102 103 104 $ • t • cu o cu SZ cu o 1 si o £ CQ o o • oo z • o • % '.•44;. V •It • • • 1 o o o c o + LU o CD O + LU O jo E H ( I U ) U011BOO-! iBjajB-) 1 o £3 O is 'H T3 o 5 <u i-1—1 OS <u oo U rs • « IM 3 105 106 (s/eui) ajey M O U o o o o o o o o O J u J l i j l J J L L J L L J L L I L L I I s - C D i n f C O C M T - O O O O O O O O O T _ to L O O L O O L O O L O O L O O a> (soooi-) 3UJ | i q6nojm>|Eeja UBipai/y M to 107 108 Figure 4-27: Isometric views of borehole locations for Case 14 109 110 I l l 5.0 Discussion 5.1 Results Summary In the previous section the effect of boreholes on flow and transport through a fracture network with a variety of parameters was studied. The effect of various processes was presented including the structure of the fracture network, the relative importance of large individual features, the influence of the location and diameter of the abandoned boreholes, and aperture variability. The majority of the cases where the borehole had an important influence had a similar structure, consisting of a background fracture network of statistically generated fractures, and two discrete features with relatively high transmissivity which extend across the entire horizontal plane. The final case also showed that abandoned boreholes could have a significant effect on flow and transport in a fracture network with homogeneous fracture statistics throughout. The first case studied, the base case, consisted of two sets of relatively low transmissivity fractures and two high transmissivity major features. In the base case network a single borehole at x = - 25 m was found to have a significant effect on both the flow and contaminant transport. The single borehole increased the total flow through the fracture network, reduced contaminant residence times and increased contaminant concentrations at the downstream boundary. The borehole also influenced the shape, spread and exit location of the contaminant plume resulting in a more structured and less disperse plume reaching the downstream boundary through the lower major feature. 112 The next set of cases investigated the effect of changing the parameters of the background fracture networks. It was found that the borehole had a much greater effect on contaminant transport when the transmissivity of the vertical fracture set was low and the transmissivity of the horizontal fracture set was high. When the situation was reversed and the vertical transmissivity was high the borehole was much less significant. In this case the path created by the vertical borehole was competing with a large number of conductive pathways created by the high-transmissivity vertical fractures. Increasing the transmissivity of the major features while the parameters of the background fracture network were fixed resulted in the borehole having a greater effect on flow and transport. When the transmissivity of the major features was increased from lx l0" 4 m2/s to lx l0" 3 m2/s, adding the borehole increased the flow by twice the amount of the base case. The ratio of contaminant residence time without the borehole to residence time with the borehole was also twice as large as the base case. When the transmissivity of the major features was reduced to being only slightly larger than that of the background fractures, the total flow rate was only increased by 3% with the addition of the borehole. The contaminant residence time and particle distribution were not significantly changed as no particles traveled through the borehole. It was found that the contrast in transmissivity between the background fractures and major features was the important factor in determining the significance of the borehole. 3 2 When the major feature transmissivity was left at 1x10" m /s as above, but the 8 2 transmissivity of the background fractures was further lowered to 1x10" m /s the effect 113 of the borehole was further increased. This was due to the much lower flow rate and longer residence times when the borehole was not present. When the borehole was included the pathway it created was nearly as conductive as before, because the background fractures only make up a small portion of the path length. When neither of the major features directly intersected a constant head boundary the influence of the borehole on the contaminant breakthrough curve was much less pronounced. Contaminant residence times are still reduced by the borehole, but less so than in the base case. The breakthrough concentrations are no longer significantly higher with the borehole. The flow rate is also still effected, but much less significantly than in the base case. The borehole has the greatest effect when it joins up a highly transmissive flow path. The more restrictions in the flow path from low transmissivity fractures, the less of an effect the borehole will have. Adjusting the longitudinal position of the borehole showed that the distance of the borehole from the contaminant source was not important. The length and total conductivity of the major pathway was the same unless the borehole was located close enough to be directly connected to either a constant head boundary or the contaminant source zone. When the borehole was connected to the upper constant head boundary the flow system was significantly altered and contaminant transport followed a completely different path. Flow rates ended up being much higher, but residence times were longer due to the flow out of the borehole causing particles to initially travel upstream and 114 around the source area. If a borehole were connected directly to the contaminant source zone a major flux of contaminants through the system would be expected. Situating a borehole away from the source zone centerline resulted in the contaminant plume shifting transversely towards the borehole location. This result shows that abandoned boreholes in unknown locations could result in maximum contaminant concentrations at a location that would not be predicted by looking at the average groundwater flow direction alone. Defining a hydraulic gradient is always a difficult exercise in a complex heterogeneous flow system such as a fractured aquifer, but abandoned boreholes may add to the complexity by channeling contamination away from the flow direction predicted based on the regional water table configuration. For instance, in Case 8 when the borehole was located at y = 10 m the borehole channeled 96.4% of the released particles, but only 89% of the flow. The effect of the borehole diameter depends on how directly the borehole is connected to the constant head source and exit boundaries. If the borehole is well connected through high-transmissivity fractures then a larger borehole will result in a great increase in the flow through the system and rapid migration of contaminants to the exit boundary. The more likely scenario in a fractured rock environment is that an abandoned borehole will be indirectly connected to the flow boundaries through a series of fractures, including low-transmissivity restrictions in the flow path. In this case a larger borehole diameter will have a small effect on raising the overall conductivity and will result in longer transport times as contaminants reside in the large void space of the bore. If a borehole is 115 directly connected at the upper end to contamination at the ground surface in a recharge zone then the borehole has the potential to cause significant contamination of any underlying intersected fracture zones. In this case the borehole could serve to decrease travel times to a compliance boundary and increase the ultimate concentrations. Larger boreholes would be expected to worsen the contamination problem. The influence of the borehole was also examined for an entirely stochastic fracture network with the same boundary conditions as the previous simulations. In this case the influence of the borehole was not as pronounced but similar behaviour was observed. The base case was designed to isolate and investigate the key interesting processes within the limitations of the modeling approach. Therefore, the effect of the borehole is somewhat exaggerated over what would be expected for most field situations. However, it was shown in the final simulation that the key processes are still present in a completely random fracture network and that abandoned boreholes can have a significant influence in more realistic scenarios. This simulation also demonstrated the importance of aperture variability in forming preferential pathways with the abandoned boreholes. 5.2 Fracture Network Structure The set of simulations presented in the previous section show the results of adding a single borehole in a number of different fracture networks. Most of the fracture networks investigated have a similar structure, with a low transmissivity background network and two horizontal high transmissivity features. It was found that for a single borehole to have a significant probability of influencing the flow and transport it was necessary to 116 have the major features present in the fracture network. Part of the cause of this phenomenon can be attributed to the geometric properties of the domain studied. The total vertical extent of the fractured-rock unit is only 25 meters at its thickest point. The horizontal distance from the source zone to the exit boundary is 100 meters. Vertical boreholes in a relatively dense, homogeneous network will not have a large impact as the vertical distance traveled is significantly less important that the horizontal distance. The near-horizontal orientation of the regional hydraulic gradient is the most important factor in determining this behaviour. In a situation where the hydraulic gradient was in a near-vertical alignment, vertical boreholes would be expected to have an important influence on the transport of contaminants in a downward direction. Previous studies performed on the effect of boreholes or wells in porous media have concentrated on the case where an upper unconfined aquifer is separated from a lower confined aquifer by a low permeability aquitard, (e.g. Avci , 1992; Lacombe et al, 1995). The borehole is situated vertically so that it crosses the aquitard and creates a preferable hydraulic connection across the two aquifers. There is a gradient in the hydraulic head across the two aquifers caused by either pumping or injection so that groundwater will flow from the contaminated to the uncontaminated formation. The effect of a vertical borehole in a homogeneous porous medium has not been investigated, as in such a case a borehole would not be expected to have an important impact. This is especially true when the dominant groundwater direction is horizontal and vertical gradients are negligible. 117 Similarly, in fractured rock a vertical borehole would not be expected to significantly influence contaminant transport, except in the case where it creates a connection between two permeable zones (i.e. horizontal fractures) across a lower permeability zone (i.e. sparsely fractured bedrock). Certain statistically homogeneous fracture networks can be treated as porous media approximations i f a representative elementary volume can be defined (Schwartz and Smith, 1988). It is expected that the same fracture networks that can be modeled entirely as equivalent porous medium approximations will not be significantly influenced by the presence of abandoned boreholes. Conversely, some fracture networks contain fractures on a variety of length scales for which a representative elementary volume cannot be defined. Barton (1995) describes fracture networks that are fractal in nature. They have the property that the larger the volume of a fracture network that is sampled, the larger the maximum fracture size is likely to be. This type of fracture network description results in many small fractures and a few larger ones. The larger fractures and the connections between them are likely to govern the flow system. A network of this type would have a much higher likelihood of being influenced by a vertical borehole creating a discrete highly transmissive pathway. If two or more large fractures are intersected by the borehole a new dominant pathway could be created causing an increase in the permeability of the entire network, and possibly a rerouting and acceleration of contaminant transport through the system. The base case scenario is also analogous to the case where two highly fractured zones are vertically separated by a zone of relatively low permeability. More dispersion would take 118 place along the tortuous pathways in a fracture zone relative to the large single fracture used in the model, but a similar temporal and spatial behaviour would be expected. Without any boreholes present, flow and transport would take place primarily in the fracture zones. The sparsely fractured zone separating the two fracture zones would limit contaminant transport into the lower aquifer. In the presence of one or more vertical boreholes however, a new conductive pathway would be created allowing contaminants to reach the previously uncontaminated zone. Furthermore, contaminants may show up in an area where they were not previously expected to be, which has implications for the design of detection monitoring networks. Reducing the hydraulic head in the lower fractured reservoir through pumping would likely cause contaminants to travel vertically through the borehole, contaminating the lower aquifer and eventually reaching the pumping well. 5.3 Implications for Monitoring Network Design A typical problem faced by consulting engineers is choosing the best location for a limited number of monitoring wells to provide plume detection before contaminants reach a compliance boundary. If contaminants are shown to reach a compliance boundary at concentrations above a predetermined level a site owner may face large penalties on top of the cost of containment and remediation. However, the high cost of drilling and completing detection wells in hard rock will limit the number of wells that can be included in any monitoring network design, often to only a few wells. For this reason it is important for the designer to be able to predict the most likely plume location, the plume size and the maximum concentrations at the detection well. The results of this 119 study have shown that abandoned boreholes could have a significant impact on all three of these parameters. The ideal plume for the purposes of detection has a large footprint at high concentrations so that it is easy to intersect and detect contaminated water. This conflicts with the objectives of containment and remediation because a larger plume will be much more expensive to remove from the subsurface. There is a tradeoff therefore between locating detection wells close to the contaminant source where plumes will be difficult to detect and relatively inexpensive to remediate or locating the detection wells close to the compliance boundary where less wells will be needed but the cost of remediation will be higher. There is also the possibility that once a plume reaches a monitoring network close to the compliance boundary it will be below the threshold concentration and remedial action will not be required. The best way to handle the conflicting objectives for monitoring network design is through a decision analysis approach. Massmann et al. (1991) provide a technique for applying decision analysis to the multi-criterion problem of monitoring network design. Jardine et al. (1996) apply this technique for siting monitoring networks in fractured rock aquifers. A decision analysis for monitoring network design could be extended to three dimensions following the method used in this study, but would require the generation of a large number of particles to provide useful concentration values in the finite-element network. In this study the particle tracking algorithm required as much as 15 minutes to 1 2 0 execute for 1000 particles, so increasing the number of particles significantly would result in large simulation times for the decision analysis approach. In the base case scenario the abandoned borehole causes the contaminant plume to reach the downstream boundary faster and at higher concentrations. The plume was also less spread out temporally and initial breakthrough took place along the centerline of the flow system, with initial breakthrough time increasing as the distance from the centerline increased. Some of these effects could be considered positive and some negative from the frame of reference of the owner of a site with a potential contamination problem. A positive influence would be any that reduces the total cost, including cost of monitoring, remediation and any fines levied for compliance failure. The decrease in residence time would reduce the cost of a monitoring network design simply because groundwater sampling and testing would occur over a shorter time period before initial detection occurs at a detection monitoring well. However, the total cost of such a situation will be high, as the concentrated plume will need to be contained and remediated before it reaches the compliance boundary. According to standard economic discounting principles the earlier in time a cost occurs, the greater its impact. The present value of any future cost can be calculated by the following equation: FV PV = -?—- (17) (1 + 0 " 121 where FV is the future value in real terms, PV is the present value in today's dollars, n is the number of compound periods and i is the discount rate. For instance, i f a compliance boundary is breached 100 years from now, and the fine is one million dollars at that time, then it would be possible to cover that cost by investing $73 today assuming it is possible to earn an interest rate of 10% over this time period. Obviously this represents an extreme scenario and the choice of the discount rate is subject to debate. A more conservative estimate of the discount rate might be 3%, in which case the present cost of failure 100 years from now would be $52000. Either way it is detrimental to the site owner when the abandoned borehole decreases residence times. The higher concentrations present in the plume with the abandoned borehole can be an advantage because they facilitate detection of the plume and may prevent non-compliance. However, i f the contaminant concentrations were very low without the borehole then they may have been below the compliance limit. In this case the borehole would be severely detrimental as without it the entire cost of remediation could be avoided. The smaller temporal spread of the plume can be both a benefit and a detriment, as there is a tradeoff in monitoring and remediation costs. A smaller plume means that samples must be taken more frequently to ensure detection. If the plume passes through the detection monitoring point between two sampling periods then the site owner could be faced with the high cost of non-compliance. The influence of the borehole means that more frequent samples must be taken and therefore raises monitoring costs. However, 122 once the plume has been detected, i f containment and remediation must take place then a smaller plume will result in substantially lower cleanup costs. The increased cost of sampling must be weighed against the decreased cost of remediation to determine whether or not the abandoned boreholes have increased costs. It is necessary to first determine probabilities of detection and failure through stochastic simulation to make this calculation. When a borehole is located off to the side of the main flow direction from the source zone it was shown that the plume can be displaced transversely at the exit boundary. If the abandoned borehole locations are not known then a more complex and expensive monitoring network will be required. Usually the optimal location of a monitoring well will be along the centerline of the flow direction, where the maximum plume concentrations are most likely to occur. When abandoned boreholes are scattered throughout a site then the initial breakthrough of the plume could occur in line with a borehole. A number of monitoring wells would need to be located along the transverse direction to ensure detection of the front of the plume. Overall it is clear that abandoned boreholes at unknown locations.will have a negative impact on the design of a monitoring network. The boreholes increase the complexity and uncertainty of the flow and transport system and therefore increase the cost of ensuring compliance. If the geometric properties of the fracture network make it prone to being influenced by abandoned boreholes, then it may be beneficial to properly abandon all drill holes at a site to prevent future costs and design difficulties. It may also be cost-123 effective to locate and seal existing boreholes at potentially contaminated sites, as reducing uncertainty and complexity will result in lower design and engineering costs. 5.4 Proper Borehole Decommissioning Boreholes drilled into any subsurface environment with the potential for groundwater contamination should be properly decommissioned. The easiest way to decommission the type of small, uncased boreholes that are typically used for subsurface exploration is to fill the entire borehole with an approved sealing material. The approved method for choosing fill materials for decommissioning boreholes in British Columbia is given in A S T M standard D5299-99. A 2" diameter borehole only requires 0.02 cubic feet of fdling material per foot of length. One bag of bentonite will be enough to fill 35 feet of well. This is a relatively cheap option and will save substantial costs later in trying to locate undocumented wells. It is desirable to properly decommission the boreholes at the time of drilling and core extraction, thereby removing one factor of uncertainty from future environmental studies at the site. A more difficult decision needs to be made at a site where abandoned boreholes may already exist in unspecified locations. If the underlying geological environment is such that abandoned boreholes are not likely to cause concern then it may be more cost effective to do nothing at the site and assume the small amount of risk associated with the boreholes. However, i f the underlying geology is prone to being affected by the abandoned boreholes, as discussed in the previous section, then it may be cost-effective to attempt to detect and plug the boreholes at the site. 124 It will likely be difficult to attempt to locate small exploration boreholes once the drilling records have been lost. The first stage in attempting to detect the boreholes is to investigate historical data at the site. Interviewing employees who have worked at the site, especially those who may have been involved with exploration can give important leads. Old exploration maps and pre-construction reports may delineate the locations of exploration holes. Total magnetic surveys have proven useful in detecting abandoned wells in the subsurface, but only when a steel casing exists. This technique will not be useful for exploration boreholes that are typically uncased. Williams and Hecker (1998) had success using ground-penetrating radar to locate old oil and gas exploration boreholes. Geophysical methods may be useful for detecting abandoned boreholes but the associated cost could be prohibitive. It may be necessary to simply evaluate the increased risk at the site from the abandoned boreholes and incorporate this risk into the facility and detection monitoring design. 125 6.0 Conclusions The present modeling study shows the possible impacts of abandoned boreholes for a specific range of fracture geometry and borehole properties. Abandoned boreholes were found to have the greatest influence when they connect large, high-transmissivity fractures and form a preferential pathway through the fracture network. In this situation the presence of abandoned boreholes had the following impacts: • Earlier arrival of contaminants at a downstream compliance boundary • Higher peak contaminant concentrations at the boundary • More compact plume with a smaller breakthrough period and shorter tail • Less surface area contact in the fracture network therefore earlier arrival for sorbing contaminants • Transverse displacement of the plume towards the borehole location • Plume transport in a direction other than that which would be predicted based on the regional hydraulic gradient • Increased cost of detection monitoring network design to account for increased uncertainty when borehole locations are unknown. Improperly abandoned vertical boreholes are expected to have an impact in a horizontally oriented flow domain if there are large horizontal fractures present which dominate the flow system. The boreholes will only have an effect if the network is such that the borehole creates a shortcut between two permeable pathways, or creates a new pathway. A multi-scale or fractal network could be studied stochastically to analyze the 126 probabilistic effects of abandoned boreholes on the flow and transport. The boreholes would likely have a significant effect on those realizations where the larger fractures happen to be intersected by the borehole. For other realizations where the boreholes do not intersect more than one major feature, and essentially are a dead-end, there would be no noticeable effect. Stochastic simulations require significantly more computational time and effort, so deterministic simulations were used in this study. A representative set of realizations was presented which demonstrates the expected range of behaviour but does not provide probabilities for the output parameters based on a given set of fracture network statistics. The results of flow and transport modeling demonstrate that there is a potential for abandoned boreholes to form connections in a fracture network and create new fast pathways for solute transport. Sub-vertical boreholes are most likely to have an effect when the major conducting fractures are sub-horizontal. This is especially true when the major flow direction is in the horizontal plane. Perhaps equally important is the fact that in many fracture networks abandoned boreholes are not expected to have any impact. In dense, uniform networks with a significant amount of vertical fractures a vertical borehole will not significantly increase the conductivity. The potential for contamination problems due to abandoned boreholes must be assessed on a site-specific basis. It is hoped that the information provided in this thesis will help guide the risk assessment process at any site in fractured rock. 127 References A S T M D5299-99, Standard guide for decommissioning of ground water wells, vadose zone monitoring devices, boreholes, and other devices for environmental activities, American Society for Testing and Materials, West Conshohocken, PA. Avci, C. B., Flow occurrence between confined aquifers through improperly plugged boreholes, J. Hydrol, 139, 97-114, 1992. Avci, C. B., Evaluation of flow leakage through abandoned wells and boreholes, Water Resour. Res., 30(9), 2565-2578, 1994. Baecher, G. B., N . A . Lanney, and H. Ff. Einstein, Statistical description of rock properties and sampling, Proc. of 18th U.S. Sym. On Rock Mechanics, American Institute of Mining Engineers, 5C1-8, 1977. Barton, C.C., Fractal analysis of the scaling and spatial clustering of fractures in rock: in Barton, C.C. and La Pointe, P.R., eds., Fractals in the Earth Sciences: Plenum Press, New York, 141-177, 1995. Bear, J., Dynamics of Fluids in Porous Media, American Elsevier Publishing Co., New York, N Y , 1972. Brikowski, T., Flow between aquifers through filled cylindrical conduits: Analytical solution and application to underground nuclear testing sites, J. Hydrol., 146, 115-130, 1993. Brown, R. Fluid flow through rock joints: the effect of surface roughness, J.G.R., 92(B2), 1337-1347, 1987. Bodvarsson, G. S., T. M . Bandurraga, and Y . S. Wu (Eds.), The Site-Scale Unsaturated Zone Model of Yucca Mountain, Nevada, for the Viability Assessment, Rep. L B N L -40376, Lawrence Berkeley Natl. Lab. Berkeley, Calif, 1997. Burkholder, H. C , Methods and data for predicting nuclide migration in geologic media, paper presented at International Symposium on Management of Wastes from the L W R Fuel Cycle, Denver, Colo., 1976. Carrera, J., J. Heredia, S. Vomvoris, and P. Hufschmied, Modeling of flow on a small fractured monzonitic gneiss block, pp. 115-167 in Hydrogeology of Low Permeability Environments, International Association of Hydrogeologists, Hydrogeology: Selected Papers, vol. 2, S. P. Neuman and I. Neretnieks, eds. Hannover: Heise, 1990. Chan, T., P. Griffith, B. W. Nakka, and K. R. Khair, Finite-element modelling of geomechanical and hydraulic responses to the Room 209 heading extension excavation 128 response. Experiment II: Post-excavation analysis of experimental results, Atomic Energy of Canada Limited., Report no. 9566, -4 July, 1993. Daugherty, R. L. , J. B. Franzini, and E. J. Finnemore, Fluid Mechanics with Engineering Applications, eighth edition, McGraw-Hill Inc., 1985. Dershowitz, W. S., J. Geier, and K. Lee, Field validation of conceptual models for fracture geometry, submitted to Rock Mechanics and Rock Engineering, 1989 Dershowitz, W. S., P. Wallman, J. E. Geier, and G. Lee, Preliminary discrete fracture network modeling of tracer migration experiment at the SCV site, SKB Technical Report 91-23, Swedish nuclear fuel and waste management co., 1991. Dershowitz, W. S., G. Lee, J. Geier, T. Foxford, P. LaPointe, and A . Thomas, F R A C M A N interactive discrete feature data analysis, geometric modeling, and exploration simulation, User Documentation, Version 2.5, Golder Associates, 1995. Dverstorp, B. , J. Andersson, and W. Nordqvist, Discrete fracture network interpretation of field tracer migration in sparsely fractured rock, Water Resour. Res., 28(9), 2327-2343, 1992. Freeze, R. A. , Three-dimensional, transient, saturated-unsaturatediflow in a groundwater basin, Water Resour. Res., 7(2), 347-366, 1971. Freeze, R.A., and J.A. Cherry, Groundwater, Englewood Cliffs, N.J. , Prentice-Hall, 1979. Gass, T. E., J. E. Lehr, and H. W. Heiss Jr., Impact of abandoned jwells on ground water, Rep. EPA-600/3-77-095, U.S. Environ. Prot. Agency, Washington, D. C , 1977. Hull, L. C , J. D. Miller, and T. M . Clemo, Laboratory and simulation studies of solute transport in fracture networks, Water Resour. Res., 23(8), 1505-1513, 1987 Jardine K. , L. Smith, T. Clemo, Monitoring networks in fractured rocks: A decision analysis approach, Ground Water, 34(3), 504-518, 1996 Javandel I., C. F. Tsang, P. A . Witherspoon, D. Morganwalp, Hydrologic detection of abandoned wells near proposed injection wells for hazardous waste disposal, Water Resour. Res., 24(2), 261-270, 1988. Lacombe, S., E. A . Sudicky, S. K. Frape, and A . J. A . Unger, Influence of leaky boreholes on cross-formational groundwater flow and contaminant transport, Water Resour. Res., 31(8), 1871-1882, 1995. Langhus, B., and A . Snider, ZOEI, a computer model to calculate the zone of endangering influence of class II injection wells, AAPG Bulletin, 83(7), 1999. 129 La Pointe, P. R., S. Follin, J. Hermanson, and P. Wallman, A method for identifying and analyzing conductive pathways in fractured rock masses for performance assessment transport models, SKB Technical Report 95-28, Swedish nuclear fuel and waste management co., 1995. Maloszewski, P., and A . Zuber, Tracer Experiments in Fractured Rocks: Matrix Diffusion and the Validity of Models, Water Resour. Res., 29(8), 2723-2735, 1993. Massmann, J., and R. A . Freeze, Groundwater contamination from waste management sites: The interaction between risk-based engineering design and regulatory policy. 1. Methodology, Water Resour. Res., 23(2), 351-367, 1987. Massmann, J, R. A . Freeze, J. L. Smith, T. Sperling, and B. R. James, Hydrogeological decision analysis: 2. Applications to groundwater contamination, Ground Water, 29(4), 536-548, 1991. McDonald, M . G., and A . W. Harbaugh, A modular three-dimensional finite-difference ground-water flow model: U.S. Geological Survey Techniques of Water Resources Investigations, Book 6, Chapter A l , 1988. Meijerink, J. A . , and H. A. van der Vorst, A n iterative solution method for linear systems of which the coefficient matrix is a symmetric M-matrix, Math. Of Comp., 31(137), 148-168, 1977. Meyer, P. D., A . J. Valocchi, and J. W. Eheart, Monitoring network design to provide initial detection of groundwater contamination, Water Resour. Res., 30(9), 2647-2659, 1994. Miller, I., G. Lee, W. Dershowitz, and G. Sharp, MAFIC Matrix/fracture Interaction Code with Solute Transport, User Documentation, Version pi.5, Golder Associates, 1995. Moreno, L. , Y .W. Tsang, C.F. Tsang, V.F. Hale, and I, Neretnieks, Flow and tracer transport in a single fracture: a stochastic model and its relation to some field observations, Water Resour. Res., 24(12), 2033-2048, 1988. Naff, R. L. , D. F. Haley, and E. A . Sudicky, High-resolution Monte Carlo simulation of flow and conservative transport in heterogeneous porous media, 1. Methodology and flow results, Water Resour. Res., 34(4), 663-667, 1998a. Naff, R. L. , D. F. Haley, and E. A . Sudicky, High-resolution Monte Carlo simulation of flow and conservative transport in heterogeneous porous media, 2. Transport Results, Water Resour. Res., 34(4), 669-697, 1998b. 130 Neretnieks, I., Diffusion in the rock matrix: A n important factor in radionuclide retardation?, J.G.R., 85,4379,1980. Neuman, S., Stochastic continuum representation of fractured rock permeability as an alternative to the R E V and fracture network concepts, Rock Mechanics: Proc. Of the 28th U.S. Symposium, 533-561, Balkema, Rotterdam, 1987. Nordqvist, A . W., T. W. Tsang, C. F. Tsang, B. Dverstorp, and J. Andersson, A variable aperture fracture network model for flow and transport in fractured rock, Water Resour. Res., 28(6), 1704-1714, 1992. Parney, R. W., Statistical continuum modeling of mass transport through fractured media, in two and three dimensions, Ph.D. Thesis, University of British Columbia, 1999. Parney, R., and L. Smith, Fluid velocity and path length in fractured media, Geophysical Research Letters, 22(11), 1437-1440, 1995. Raven, K .G . , K.S. Novakowski, and P.A. Lapcevic, Interpretation of field tests of a single fracture using a transient solute storage model, Water Resour. Res., 24(12), 2019-2032, 1988. Raven K. G., D. W. Lafleur, and R. A . Sweezey, Monitoring well into abandoned deep-well disposal formations at Sarnia, Ontario, Canadian Geotechnical Journal, 27(1), 105-118, 1990. Schwartz, F. W., and L. Smith, A continuum approach for modeling mass transport in fractured media, Water Resour. Res., 24(8), 1360-1372, 1988. Smith, L. and F. W. Schwartz, An analysis of the influence of fracture geometry on mass transport in fractured media, Water Resour. Res., 20(9), 1241-1252, 1984. Smith, L. , C. W. Mase, F. W. Schwartz, and D. Chorely, A numerical model for transport in networks of planar fractures, Sym. On Hydrogeology of Rocks of Low Permeability, 666-675, Tucson, Ariz., International Association of Hydrogeologists, 1985. Snow, D. T., A parallel plate model of fractured permeable media, Ph.D. Thesis, University of California, 1965. Storck, P., J. W. Eheart, and A . J. Valocchi, A method for the optimal location of monitoring wells for detection of groundwater contamination in three-dimensional heterogeneous aquifers, Water Resour. Res., 33(9), 2081-2088, 1997. Sudicky, E.A., and E.O. Frind, Contaminant transport in fractured porous media: Analytical solutions for a system of parallel fractures, Water Resour. Res., 18(6), 1634-1642, 1982. 131 Therrien, R., and E. A . Sudicky, Three-dimensional analysis of variably-saturated flow and transport in discretely-fractured porous media: Model development and illustrative examples, Jour. Contam. Hydrol., 23(1-2), 1-44, 1996. Wait, R. L., and M . J. McCollum, Contamination of fresh water aquifers through an unplugged oil-test well in Glynn County, Georgia, Ga. Geol. Surv. Mineral Newsletter, 16(3-4), 74-80. 1963. Warner, D. L. , and C. L. McConnell, Assessment of environmental implications of abandoned oil and gas wells, Journal of Petroleum Technology, 45(9), 874-880, 1993. Wels, C , and L. Smith, Retardation of sorbing solutes in fractured media, Water Resour. Res., 30(9), 2547-2563, 1994. Williams, E. G., and B. Hecker, Use of ground penetrating radar to locate old oil & gas wells; a case study, AAPG Bulletin, 82(9), 1998. Witherspoon, P.A., J.S.Y. Wang, K . Iwai, and J.E. Gale, Validity of cubic law for fluid flow in a deformable rock fracture, Water Resour. Res., 16(6), 1016-1024, 1980. Zeller, E. J., The disposal of nuclear waste, California Geology, 24(4), 79-87, 1973. Zheng, C , MT3D, A modular three-dimensional transport model for simulation of advection, dispersion, and chemical reactions of contaminants in groundwater systems, Report to the Kerr Environmental Research Laboratory, US Environmental Protection Agency, Ada, OK., 1990. Zimmerman, R., and I. W. Yeo, Fluid flow in rock fractures: cubic law, lubrication equation, and stokes equation, paper presented at International Symposium on Dynamics of Fluids in Fractured Rocks: Concepts and Recent Advances, Ernest Orlando Lawrence Berkeley National Laboratory, 1999. 132
Thesis/Dissertation
eng
Geological Engineering
Vancouver : University of British Columbia Library
For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use.
Graduate
Contaminant transport through abandoned boreholes in fractured rock
Text