Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Monitoring contamination in surface water and groundwater in a catchment with an unconfined porous aquifer… Chen, Yaming 2009

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

Item Metadata

Download

Media
24-ubc_2009_fall_chen_yaming.pdf [ 4.29MB ]
Metadata
JSON: 24-1.0052863.json
JSON-LD: 24-1.0052863-ld.json
RDF/XML (Pretty): 24-1.0052863-rdf.xml
RDF/JSON: 24-1.0052863-rdf.json
Turtle: 24-1.0052863-turtle.txt
N-Triples: 24-1.0052863-rdf-ntriples.txt
Original Record: 24-1.0052863-source.json
Full Text
24-1.0052863-fulltext.txt
Citation
24-1.0052863.ris

Full Text

MONITORING CONTAMINATION IN SURFACE WATER AND GROUNDWATER IN A CATCHMENT WITH AN UNCONFINED POROUS AQUIFER OVERLYING FRACTURED BEDROCK  by Yaming Chen B.Sc., Lanzhou University, 1986 M.Sc., Lanzhou University, 1989  A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in The Faculty of Graduate Studies (Geological Science) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) May 2009 © Yaming Chen, 2009  ABSTRACT  The factors influencing the design of an integrated surface water and groundwater contaminant monitoring network are examined for a system where a permeable unconfined aquifer overlies fractured bedrock in a headwater stream catchment. The unconfined porous aquifer is modeled as a homogeneous and isotropic porous medium with deterministic properties. The fractured bedrock is modeled as a lower permeability rock mass, but with a stochastic fracture network composed of three orthogonal fracture sets. The surface domain consists of a V-shaped overland flow zone and a linear stream channel with constant width. A continuous contaminant source zone is situated close to the land surface. The analysis is based on a steady state representation of the groundwater and surface water flow system and transport of a non-reactive solute in both the surface and subsurface domains. The probability of plume detection is defined in terms of the likelihood of detecting contamination in a performance monitoring network prior to its detection at a downstream compliance boundary. The ratio of the permeability of the unconfined porous aquifer to the bulk permeability of the underlying fractured bedrock is a key factor influencing detection probabilities. A higher permeability unconfined aquifer reduces detection probabilities in the stream, and increases detection probabilities in the subsurface. Fracture network connectivity and aperture size have a significant influence on detection probabilities in the subsurface, especially in the bedrock. Dilution of contamination in the headwater stream by inflowing non-contaminated water has a strong influence on detection probabilities. Due to the effect of dilution, sampling both stream water and streambed groundwater will enhance detection probabilities in the stream. The contaminant detection threshold is a key factor in influencing the stream water and groundwater monitoring. A lower detection threshold will increase detection probabilities in the surface water and groundwater systems, especially in the stream water. The results in this study demonstrates the importance of integrating stream water and groundwater samples in the design of a performance monitoring network and documents the factors that determine contaminant detection probabilities in the surface stream, the unconfined aquifer, and the underlying fractured bedrock.  ii  TABLE OF CONTENTS  ABSTRACT .............................................................................................................................ii TABLE OF CONTENTS ........................................................................................................iii LIST OF TABLES.................................................................................................................... v LIST OF FIGURES ................................................................................................................. vi ACKNOWLEDGEMENTS.................................................................................................... xv 1  INTRODUCTION ............................................................................................................. 1 1.1 Motivation and review of previous research ............................................................. 1 1.2 Objectives and contributions ..................................................................................... 7 1.3 Thesis outline............................................................................................................. 8  2  METHODOLOGY TO EVALUTE GROUNDWATER MONITORING NETWORKS ................................................................................................................... 12 2.1 Introduction.............................................................................................................. 12 2.2 Model domain for the groundwater system ............................................................. 13 2.3 Boundary conditions for groundwater flow and solute transport ............................ 14 2.3.1 Boundary conditions for flow....................................................................... 14 2.3.2 Boundary conditions for solute transport ..................................................... 14 2.4 Description of the numerical model: FRAC3DVS .................................................. 15 2.5 Approach for evaluation of groundwater monitoring network................................ 17 2.6 Scenarios and model parameters ............................................................................. 19 2.6.1 Base Case and the parameters ...................................................................... 19 2.6.2 Other scenarios and the parameters .............................................................. 22 2.7 Discretization of model domain and Monte Carlo simulations ............................... 23 2.7.1 Discretization of model domain ................................................................... 23 2.7.2 Number of Monte Carlo realizations ............................................................ 26 2.8 Summary.................................................................................................................. 27  3  FACTORS INFLUENCING GROUNDWATER MONITORING ................................ 41 3.1 3.2 3.3 3.4 3.5  Introduction.............................................................................................................. 41 Statistical properties of all the GW scenarios.......................................................... 42 Groundwater flow, transport and detection probabilities in the Base Case............. 44 Influence of hydraulic conductivity of unconfined porous aquifer ......................... 47 Influence of fractured bedrock................................................................................. 49 3.5.1 Influence of fractured bedrock on detection in unconfined porous aquifer.. 49 3.5.2 Influence of fracture network connectivity .................................................. 51 3.5.3 Influence of fracture apertures...................................................................... 54 3.5.4 Influence of rock matrix permeability .......................................................... 57 3.6 Influence of contaminant detection threshold.......................................................... 58 3.7 Influence of source zone location ............................................................................ 60 iii  3.8 Summary.................................................................................................................. 61 4  METHODOLOGY TO EVALUATE INTEGRATED SURFACE WATER AND GROUNDWATER MONITORING NETWORKS .............................................. 92 4.1 Introduction.............................................................................................................. 92 4.2 Model domain of the SW-GW system..................................................................... 93 4.3 Boundary conditions for coupled SW-GW flow and transport ............................... 98 4.3.1 Boundary conditions for coupled SW-GW flow .......................................... 98 4.3.2 Boundary conditions for coupled SW-GW transport ................................... 99 4.4 Description of numerical model: HydroGeoSphere ................................................ 99 4.5 Approach to evaluate integrated SW-GW monitoring network ............................ 102 4.6 Scenarios and model parameters ........................................................................... 105 4.6.1 Base Case and the parameters .................................................................... 105 4.6.2 Other scenarios and the parameters ............................................................ 109 4.7 Summary................................................................................................................ 110  5  FACTORS INFLUENCING INTEGRATED SURFACE WATER AND GROUNDWATER MONITORING ............................................................................. 120 5.1 Introduction............................................................................................................ 120 5.2 Summary of statistical properties of the coupled SW-GW scenarios.................... 121 5.3 SW-GW flow, transport and detection probabilities in the Base Case.................. 123 5.3.1 SW-GW flow and transport in the first realization..................................... 123 5.3.2 Statistics of flow and transport in all realizations....................................... 131 5.3.3 Detection probabilities for the Base Case monitoring network.................. 133 5.4 Influence of hydraulic conductivity of unconfined porous aquifer ....................... 134 5.5 Influence of fracture apertures in the bedrock....................................................... 139 5.5.1 Influence of aperture of horizontal fracture set .......................................... 139 5.5.2 Influence of apertures of vertical fracture sets ........................................... 142 5.6 Influence of dispersivities in stream water ............................................................ 145 5.7 Influence of detection threshold ............................................................................ 148 5.8 Influence of source zone location .......................................................................... 151 5.9 Summary................................................................................................................ 154  6 CONCLUSIONS ............................................................................................................ 199 6.1 Conclusions: groundwater monitoring .................................................................. 199 6.2 Conclusions: integrated surface water and groundwater monitoring .................... 201 6.3 Suggestions for future research ............................................................................. 206 REFERENCES ..................................................................................................................... 208 APPENDICES ...................................................................................................................... 218 Appendix A. Appendix B. Appendix C. Appendix D. Appendix E. Appendix F. Appendix G.  Results of flow and solute transport in Scenario SW-GW-1 .................. 218 Results of flow and solute transport in Scenario SW-GW-2a................. 220 Results of flow and solute transport in Scenario SW-GW-2b ................ 223 Results of solute transport in Scenario SW-GW-3a................................ 226 Results of solute transport in Scenario SW-GW-3b ............................... 227 Results of solute transport in Scenario SW-GW-4 ................................. 228 Results of solute transport in Scenario SW-GW-5 ................................. 229 iv  LIST OF TABLES Table 2.1 Parameters for unconfined porous aquifer in the Base Case .................................. 29 Table 2.2 Distribution parameters for fracture network in the Base Case.............................. 29 Table 2.3 Dispersivities for transport simulation in fractures in the Base Case..................... 29 Table 2.4 Parameters for simulations in rock matrix in the Base Case .................................. 30 Table 2.5 Investigated scenarios and features for groundwater simulations .......................... 30 Table 3.1 Statistics of the fracture networks in each scenario................................................ 65 Table 3.2 Bulk hydraulic conductivities of the bedrock in each scenario .............................. 65 Table 3.3 Statistics of detection time and location at the compliance boundary in each scenario................................................................................................................... 66 Table 4.1 Parameter values for simulation of surface water ................................................ 113 Table 4.2 Investigated scenarios and features for coupled SW-GW simulations................. 113 Table 5.1 Bulk hydraulic conductivity of the bedrock in each scenario .............................. 158 Table 5.2 Stream water flow in the first realization of each scenario .................................. 158 Table 5.3 Properties of stream water flow of all realizations in each scenario .................... 159 Table 5.4 Monte Carlo realizations, solute concentrations in stream water, detection time and location at compliance boundary in each scenario ................................ 159  v  LIST OF FIGURES  Figure 1.1 Conceptual diagrams of the catchment and the performance monitoring network .............................................................................................................. 11 Figure 2.1 Model domain for groundwater simulation in the system with an unconfined porous aquifer overlying fractured bedrock. Dimension (m). ........................... 31 Figure 2.2 Diagrams of upstream inflow boundary and downstream outflow boundary in the model domain .......................................................................................... 31 Figure 2.3 Locations of source zone and water table in the domain. A Three dimensions, B cross-section at Y=80m. Dimension (m). ...................................................... 32 Figure 2.4 Locations of performance monitoring wells in groundwater monitoring network .............................................................................................................. 32 Figure 2.5 Histograms of random data in truncated distributions for fracture locations, lengths and apertures in the first realization of the Base Case. A, B and C Data for fracture centre locations in uniform distributions, D and E data for fracture lengths in exponential distributions, and F data for fracture apertures in a lognormal distribution................................................................. 33 Figure 2.6 Fracture network in the first realization of the Base Case for simulation of groundwater flow and transport. A Three dimensions, B and C crosssections at Y=100m and X=250m. Aperture (m), Dimension (m). ................... 34 Figure 2.7 Simulation of the small fractured domain for test of grid spacing. A Fractured domain, B spatial discretization, C simulated hydraulic head (m), and D source zone location. Dimension (m). ............................................................... 35 Figure 2.8 Nodal concentrations of a monitoring well in the small fractured domain with different element sizes ....................................................................................... 35 Figure 2.9 Nodal concentrations of monitoring well #5 in the first realization of the Base Case with different element sizes in the bedrock domain ................................. 36 Figure 2.10 Nodal concentrations of monitoring well #7 in the first realization of the Base Case with different vertical element sizes in the unconfined aquifer ....... 36 Figure 2.11 Nodal concentrations of monitoring well #143 in the first realization of the Base Case with different vertical element sizes in the unconfined aquifer ....... 37 Figure 2.12 Nodal concentrations of monitoring well #122 in the first realization of the Base Case with different time weightings (TW) ............................................... 37 Figure 2.13 Detection time at compliance boundary vs. number of Monte Carlo realizations in the Base Case. “Minimum”, “Maximum” and “Average” represents the earliest, the latest and the arithmetic average detection time, for each set of additional realizations. ............................................................... 38 Figure 2.14 Detection probabilities of monitoring wells in the unconfined aquifer vs. number of Monte Carlo realizations in the Base Case. A Cross-section at X=350m, B cross-section at Y=100m. .............................................................. 39 Figure 2.15 Detection probabilities of monitoring wells in the bedrock vs. number of Monte Carlo realizations in the Base Case. A Cross-section at X=100m, B cross-section at Y=100m. .................................................................................. 40 Figure 3.1 Hydraulic head and saturation profile in first realization of the Base Case. A Head (m) in 3D domain, B head (m) in cross-section at Y=80m, C vi  saturation in 3D domain, and D saturation in cross-sections at X=60m and Y=80m (crossing source zone). Dimension (m)................................................ 67 Figure 3.2 Saturation profile in the depth at the location of source zone centre in first realization of the Base Case............................................................................... 67 Figure 3.3 Plume in groundwater (concentrations above detection threshold) at the time of detection at the compliance boundary in the first realization of the Base Case. A Three dimensions, B cross-section at Y=80m (source zone). Concentration (g/L), Dimension (m). ................................................................ 68 Figure 3.4 Detection probabilities of monitoring wells in the Base Case. A Unconfined porous aquifer, B fractured bedrock. ................................................................. 69 Figure 3.5 Relationship between the bedrock hydraulic conductivity and the detection time at the compliance boundary in 200 realizations of the Base Case. A Detection time vs Kx, B detection time vs Ky, and C detection time vs Kz. ....... 70 Figure 3.6 Plume in groundwater (concentrations above detection threshold) at the time of detection at the compliance boundary in the first realization of Scenario GW-1 (higher Kuc). A Three dimensions, B cross-section at Y=80m (source zone). Concentration (g/L), Dimension (m). ..................................................... 71 Figure 3.7 Detection probabilities of monitoring wells in Scenario GW-1 (higher Kuc). A Unconfined porous aquifer, B fractured bedrock. ......................................... 72 Figure 3.8 Outer boundaries of the plumes (defined by the detection threshold) in the fractured bedrock and the overlying unconfined porous aquifer, for the first & tenth realizations of the Base Case at the time of detection at the compliance boundary. A The first realization, B the tenth realization.............. 73 Figure 3.9 Comparison of the outer boundaries of the plumes (defined by the detection threshold) in the unconfined porous aquifer between the case with the underlain fractured bedrock and the case with the underlain equivalent porous medium of the bedrock, for the first & tenth realizations of the Base Case at the time of detection at the compliance boundary. A The first realization, B the tenth realization. .................................................................... 74 Figure 3.10 Contours of hydraulic head (m) in the unconfined porous aquifer (crosssection at X=120m) and in the underlain fractured bedrock (cross-section at X=103m) for the first realization of the Base Case. Dimension (m)................. 75 Figure 3.11 Contours of hydraulic head (m) in the unconfined porous aquifer (crosssection at X=120m) and the underlain equivalent porous medium of the fractured bedrock (cross-section at X=103m) for the first realization of the Base Case. Dimension (m). ............................................................................... 75 Figure 3.12 Plume in groundwater (concentrations above detection threshold) at the time of detection at the compliance boundary in the first realization of Scenario GW-2a (lower horizontal fracture density). A Three dimensions, B crosssection at Y=80m (source zone). Concentration (g/L), Dimension (m). ........... 76 Figure 3.13 Detection probabilities of monitoring wells in Scenario GW-2a (lower horizontal fracture density). A Unconfined porous aquifer, B fractured bedrock. ............................................................................................................. 77 Figure 3.14 Plume in groundwater (concentrations above detection threshold) at the time of detection at the compliance boundary in the first realization of Scenario GW-2b (lower vertical fracture densities). A Three dimensions, B crosssection at Y=80m (source zone). Concentration (g/L), Dimension (m). ........... 78 vii  Figure 3.15 Detection probabilities of monitoring wells in Scenario GW-2b (lower vertical fracture densities). A Unconfined porous aquifer, B fractured bedrock. ............................................................................................................. 79 Figure 3.16 Plume in groundwater (concentrations above detection threshold) at the time of detection at the compliance boundary in the first realization of in Scenario GW-3a (smaller horizontal fracture apertures). A Three dimensions, B cross-section at Y=80m (source zone). Concentration (g/L), Dimension (m). .................................................................................................. 80 Figure 3.17 Detection probabilities of monitoring wells in Scenario GW-3a (smaller horizontal fracture apertures). A Unconfined porous aquifer, B fractured bedrock. ............................................................................................................. 81 Figure 3.18 Plume in groundwater (concentrations above detection threshold) at the time of detection at the compliance boundary in the first realization of Scenario GW-3b (smaller vertical fracture apertures). A Three dimensions, B crosssection at Y=80m (source zone). Concentration (g/L), Dimension (m). ........... 82 Figure 3.19 Detection probabilities of monitoring wells in Scenario GW-3b (smaller vertical fracture apertures). A Unconfined porous aquifer, B fractured bedrock. ............................................................................................................. 83 Figure 3.20 Plume in groundwater (concentrations above detection threshold) at the time of detection at the compliance boundary in the first realization of Scenario GW-4 (higher rock matrix permeability). A Three dimensions, B crosssection at Y=80m (source zone). Concentration (g/L), Dimension (m). ........... 84 Figure 3.21 Detection probabilities of monitoring wells in Scenario GW-4 (higher rock matrix permeability). A Unconfined porous aquifer, B fractured bedrock. ...... 85 Figure 3.22 Plume in groundwater (concentrations above detection threshold) at the time of detection at the compliance boundary in the first realization of Scenario GW-5a (higher detection threshold: 1% of source concentration). A Three dimensions, B cross-section at Y=80m (source zone). Concentration (g/L), Dimension (m). .................................................................................................. 86 Figure 3.23 Detection probabilities of monitoring wells in Scenario GW-5a (higher detection threshold: 1% of source concentration). A Unconfined porous aquifer, B fractured bedrock. ............................................................................. 87 Figure 3.24 Plume in groundwater (concentrations above detection threshold) at the time of detection at the compliance boundary in the first realization of Scenario GW-5b (lower detection threshold: 0.01% of source concentration). A Three dimensions, B cross-section at Y=80m (source zone). Concentration (g/L), Dimension (m)......................................................................................... 88 Figure 3.25 Detection probabilities of monitoring wells in Scenario GW-5b (lower detection threshold: 0.01% of source concentration). A Unconfined porous aquifer, B fractured bedrock. ............................................................................. 89 Figure 3.26 Plume in groundwater (concentrations above detection threshold) at the time of detection at the compliance boundary in the first realization of Scenario GW-6 (source zone on bedrock). A Three dimensions, B cross-section at Y=80m (source zone). Concentration (g/L), Dimension (m). ........................... 90 Figure 3.27 Detection probabilities of monitoring wells in Scenario GW-6 (source zone on bedrock). A Unconfined porous aquifer, B fractured bedrock. .................... 91  viii  Figure 4.1 Illustration of a small headwater catchment with stream flow (photo from the Bertrand Creek Watershed, Langley, BC, Canada, taken by Cindy Starzyk) . 114 Figure 4.2 Types of stream water flow and streamlines in open channel (plane view)....... 114 Figure 4.3 Mechanisms affecting solute concentrations in stream water (point source)..... 114 Figure 4.4 Mechanisms affecting solute concentrations in stream water (continuous source).............................................................................................................. 115 Figure 4.5 Surface domain (overland zone and stream channel) of the catchment for coupled SW-GW simulations. Dimension (m). ......................................... 115 Figure 4.6 Subsurface domain of the catchment with an unconfined porous aquifer overlying fractured bedrock for coupled SW-GW simulations. Dimension (m).................................................................................................................... 116 Figure 4.7 Illustration of the methods for discretization of surface and subsurface domains (Therrien et al., 2005) ....................................................................... 116 Figure 4.8 Illustration of the approaches for coupling of surface and subsurface domains (Therrien et al., 2005)...................................................................................... 117 Figure 4.9 Layout of the surface water observation points at the outlet of stream channel . 117 Figure 4.10 Layout of the surface water observation points along the centerline of stream channel ................................................................................................. 118 Figure 4.11 Locations of the performance groundwater monitoring wells for the integrated SW-GW monitoring........................................................................ 118 Figure 4.12 Fracture network in three dimensions in the first realization of the Base Case for simulation of coupled SW-GW flow and transport. Aperture (m), Dimension (m). ................................................................................................ 119 Figure 5.1 Evolution of hydraulic heads in the stream water and in the streambed groundwater at observation point (at the node at X=100m, Y=100m) in the channel for the first realization of the Base Case ............................................ 160 Figure 5.2 Evolution of stream flux at observation point (at the node at X=100m, Y=100m) in the channel for the first realization of the Base Case.................. 160 Figure 5.3 Total discharge along the stream channel at the steady state for the first realization of the Base Case............................................................................. 161 Figure 5.4 Water depth along the stream channel at the steady state for the first realization of the Base Case............................................................................. 161 Figure 5.5 Flow velocity along the stream channel at the steady state for the first realization of the Base Case............................................................................. 162 Figure 5.6 Flow velocity profiles at stream cross-sections at the steady state for the first realization of the Base Case. A Stream with the bedrock in subsurface, B stream with the equivalent porous medium. .................................................... 163 Figure 5.7 Coupled steady-state surface water flow in the first realization of the Base Case. A Hydraulic head (m) and flow velocity vectors; B water depth (m); C, D and E surface flow velocity Vx, Vy and Vz (m/d); and F flux exchange with the subsurface (m/d). Dimension (m). ..................................................... 164 Figure 5.8 Coupled steady-state groundwater flow in the first realization of the Base Case. A Hydraulic head (m) in the 3D domain, B saturation profile, C head (m) in cross-section at Z=120m (in the unconfined aquifer), D head (m) in cross-section at Z=110m (in the bedrock), E head (m) and flow velocity vectors in cross-section at Y=100m (beneath the streambed), and F head  ix  (m) and flow velocity vectors in cross-sections at X=60m, X=250m and X=500m. Dimension (m)................................................................................. 165 Figure 5.9 Evolution of plumes (concentrations above detection threshold) in surface water and groundwater before detection of the contaminant at the compliance boundary in the first realization of the Base Case. A 1500 days, B 2000 days, C 2500 days, and D 3184 days (detection at the stream outlet). Concentration (g/L), Dimension (m)................................................... 166 Figure 5.10 Plume in the SW-GW system (concentrations above detection threshold) at time of detection at compliance boundary in the first realization of the Base Case. A plume in stream water (2D); B plume in groundwater (crosssection Y=70m, source zone); C solute entering stream water from groundwater through streambed (not through overland pathways), crosssections at X=50m, 70m, 80m, 100m, 150m, 200m, 250m; and D solute entering stream water through streambed, cross-section Y=100m (domain centerline). Concentration (g/L), Dimension (m). ........................................... 167 Figure 5.11 Concentrations along the stream centreline (in stream water and streambed) at the time of detection at the compliance boundary in the first realization of the Base Case............................................................................................... 168 Figure 5.12 Concentration profiles in stream water at the time of detection at the compliance boundary in the first realization of the Base Case. Crosssections at A X=100m, B X=200m, C X=300m, D X=400m, and E X=500m. .......................................................................................................... 169 Figure 5.13 Statistics of water depth, discharge and flow velocity at stream outlet in different number of realizations for the Base Case. A Water depth, B total discharge, and C mean velocity (Vx). “Minimum”, “Maximum” and “Average” represents the smallest, the largest and the arithmetic average values of each parameter, for each set of additional realizations. ................... 170 Figure 5.14 Statistics of the maximum solute concentration in stream at the time of detection at the compliance boundary, and detection time at compliance boundary in different number of realizations for the Base Case. A Maximum solute concentration, B detection time. “Minimum”, “Maximum” and “Average” represents the smallest, the largest and the arithmetic average values of each parameter, for each set of additional realizations....................................................................................................... 171 Figure 5.15 Detection probabilities at the stream outlet in the Base Case ........................... 172 Figure 5.16 Detection probabilities of groundwater monitoring wells in the Base Case. A Unconfined porous aquifer, B fractured bedrock. ........................................... 172 Figure 5.17 Coupled steady-state SW-GW flow in the first realization of Scenario SWGW-1 (higher Kuc). A Surface water depth (m), B surface flow velocity Vx (m/d), C flux exchange (m/d) between surface and subsurface, D saturation profile in the subsurface, E hydraulic head (m) in the subsurface domain, and F groundwater hydraulic head (m) and flow velocity vectors in crosssection at Y=100m (beneath the streambed). Dimension (m). ........................ 173 Figure 5.18 Comparison of stream water flow parameters for the first realizations of Scenario SW-GW-1 (higher Kuc) and the Base Case. A Water depth, B total discharge, and C mean velocity (Vx)................................................................ 174  x  Figure 5.19 Plume in the SW-GW system (concentrations above detection threshold) at time of detection at compliance boundary in the first realization of Scenario SW-GW-1 (higher Kuc). A Plume in stream water (2D), B plumes in stream water and groundwater (3D), C plume in groundwater (cross-section Y=70m, source zone), and D solute entering stream water through streambed (cross-section Y=100m, domain centerline). Concentration (g/L), Dimension (m)....................................................................................... 175 Figure 5.20 Concentrations along the stream centreline (in stream water and streambed) at the time of detection at the compliance boundary in the first realization of Scenario SW-GW-1 (higher Kuc)................................................................. 176 Figure 5.21 Concentration profiles in stream water at the time of detection at the compliance boundary in the first realization of Scenario SW-GW-1 (higher Kuc). Cross-sections at A X=337m, B X=350m, C X=400m, D X=450m, and E X=500m................................................................................................. 177 Figure 5.22 Detection probabilities at the stream outlet in Scenario SW-GW-1 (higher Kuc)................................................................................................................... 178 Figure 5.23 Detection probabilities of groundwater monitoring wells in Scenario SWGW-1 (higher Kuc). A Unconfined porous aquifer, B fractured bedrock. ....... 178 Figure 5.24 Comparison of stream water flow parameters for the first realizations of Scenario SW-GW-2a (larger horizontal fracture apertures) and the Base Case. A Water depth, B total discharge, and C mean velocity (Vx)................. 179 Figure 5.25 Plume in the SW-GW system (concentrations above detection threshold) at time of detection at compliance boundary in the first realization of Scenario SW-GW-2a (larger horizontal fracture apertures). A Plume in stream water (2D), B plumes in stream water and groundwater (3D), C plume in groundwater (cross-section Y=70m, source zone), and D solute entering stream water through streambed (cross-section Y=100m, domain centerline). Concentration (g/L), Dimension (m). ........................................... 180 Figure 5.26 Concentrations along the stream centreline (in stream water and streambed) at the time of detection at the compliance boundary in the first realization of Scenario SW-GW-2a (larger horizontal fracture apertures)........................ 181 Figure 5.27 Concentration profiles in stream water at the time of detection at the compliance boundary in the first realization of Scenario SW-GW-2a (larger horizontal fracture apertures). Cross-sections at A X=160m, B X=200m, C X=300m, D X=400m, and E X=500m. ........................................................... 182 Figure 5.28 Detection probabilities at the stream outlet in Scenario SW-GW-2a (larger horizontal fracture apertures)........................................................................... 183 Figure 5.29 Detection probabilities in groundwater monitoring wells in Scenario SWGW-2a (larger horizontal fracture apertures). A Unconfined porous aquifer, B fractured bedrock. ........................................................................................ 183 Figure 5.30 Comparison of stream water flow parameters for the first realizations of Scenario SW-GW-2b (larger vertical fracture apertures) and the Base Case. A Water depth, B total discharge, and C mean velocity (Vx). ......................... 184 Figure 5.31 Plume in the SW-GW system (concentrations above detection threshold) at time of detection at compliance boundary in the first realization of Scenario SW-GW-2b (larger vertical fracture apertures). A Plume in stream water (2D), B plumes in stream water and groundwater (3D), C plume in xi  groundwater (cross-section Y=70m, source zone), and D solute entering stream water through streambed (cross-section Y=100m, domain centerline). Concentration (g/L), Dimension (m). ........................................... 185 Figure 5.32 Concentrations along the stream centreline (in stream water and streambed) at the time of detection at the compliance boundary in the first realization of Scenario SW-GW-2b (larger vertical fracture apertures)............................ 186 Figure 5.33 Concentration profiles in stream water at the time of detection at the compliance boundary in the first realization of Scenario SW-GW-2b (larger vertical fracture apertures). Cross-sections at A X=100m, B X=200m, C X=300m, D X=400m, and E X=500m. ........................................................... 187 Figure 5.34 Detection probabilities in groundwater monitoring wells in Scenario SWGW-2b (larger vertical fracture apertures). A Unconfined porous aquifer, B fractured bedrock. ............................................................................................ 188 Figure 5.35 Comparison of the concentrations profiles in stream water at the time of detection at the compliance boundary in the first realization of Scenario SW-GW-3a (smaller αl in stream water) and the Base Case. Cross-sections at A X=100m, B X=200m, C X=300m, D X=400m, and E X=500m............. 189 Figure 5.36 Comparison of the concentrations profiles in stream water at the time of detection at the compliance boundary in the first realization of Scenario SW-GW-3b (smaller αt in stream water) and the Base Case. Cross-sections at A X=100m, B X=200m, C X=300m, D X=400m, and E X=500m............. 190 Figure 5.37 Detection probabilities at the stream outlet in Scenario SW-GW-3b (smaller αt in stream water) ........................................................................................... 191 Figure 5.38 Plume in the SW-GW system (concentrations above detection threshold) at time of detection at compliance boundary in the first realization of Scenario SW-GW-4 (lower detection threshold 1.0×10-4g/L). A Plume in stream water (2D), B plumes in stream water and groundwater (3D), C plume in groundwater (cross-section Y=70m, source zone), and D solute entering stream water through streambed (cross-section Y=100m, domain centerline). Concentration (g/L), Dimension (m). ........................................... 191 Figure 5.39 Concentrations along the stream centreline (in stream water and streambed) at the time of detection at the compliance boundary in the first realization of Scenario SW-GW-4 (lower detection threshold 1.0×10-4g/L) .................... 192 Figure 5.40 Concentration profiles in stream water at the time of detection at the compliance boundary in the first realization of Scenario SW-GW-4 (lower detection threshold 1.0×10-4g/L). Cross-sections at A X=100m, B X=200m, C X=300m, D X=400m, and E X=500m......................................................... 193 Figure 5.41 Detection probabilities at the stream outlet in Scenario SW-GW-4 (lower detection threshold 1.0×10-4g/L) ..................................................................... 194 Figure 5.42 Detection probabilities in groundwater monitoring wells in Scenario SWGW-4 (lower detection threshold 1.0×10-4g/L). A Unconfined porous aquifer, B fractured bedrock. ........................................................................... 194 Figure 5.43 Plume in the SW-GW system (concentrations above detection threshold) at time of detection at compliance boundary in the first realization of Scenario SW-GW-5 (source zone on bedrock). A Plume in stream water (2D), B plumes in stream water and groundwater (3D), C plume in groundwater (cross-section Y=70m, source zone), and D solute entering stream water xii  Figure Figure  Figure Figure  Figure  Figure  Figure  Figure  Figure  through streambed (cross-section Y=100m, domain centerline). Concentration (g/L), Dimension (m). .............................................................. 195 5.44 Comparison of locations for solute entering stream water in Base Case and Scenario SW-GW-5 (source zone on bedrock) at the time of detection at the compliance boundary................................................................................. 196 5.45 Concentration profiles in stream water at the time of detection at the compliance boundary in the first realization of Scenario SW-GW-5 (source zone on bedrock). Cross-sections at A X=100m, B X=200m, C X=300m, D X=400m, and E X=500m................................................................................. 197 5.46 Detection probabilities in groundwater monitoring wells in Scenario SWGW-5 (source zone on bedrock). A Unconfined porous aquifer, B fractured bedrock. ........................................................................................................... 198 A.1 Statistics of water depth, discharge and flow velocity at stream outlet in different number of realizations for Scenario SW-GW-1 (higher Kuc). A Water depth, B total discharge, and C mean velocity (Vx). “Minimum”, “Maximum” and “Average” represents the smallest, the largest and the arithmetic average values of each parameter, for each set of additional realizations....................................................................................................... 218 A.2 Statistics of the maximum solute concentration in stream at the time of detection at the compliance boundary, and detection time at compliance boundary in different number of realizations for Scenario SW-GW-1 (higher Kuc). A Maximum solute concentration, B detection time. “Minimum”, “Maximum” and “Average” represents the smallest, the largest and the arithmetic average values of each parameter, for each set of additional realizations...................................................................................... 219 B.1 Coupled steady-state SW-GW flow in the first realization of Scenario SWGW-2a (larger horizontal fracture apertures). A Surface water depth (m), B surface flow velocity Vx (m/d), C flux exchange (m/d) between surface and subsurface, D saturation profile in the subsurface, E hydraulic head (m) in the subsurface domain, and F groundwater hydraulic head (m) and flow velocity vectors in cross-section at Y=100m (beneath the streambed). Dimension (m). ................................................................................................ 220 B.2 Statistics of water depth, discharge and flow velocity at stream outlet in different number of realizations for Scenario SW-GW-2a (larger horizontal fracture apertures). A Water depth, B total discharge, and C mean velocity (Vx). “Minimum”, “Maximum” and “Average” represents the smallest, the largest and the arithmetic average values of each parameter, for each set of additional realizations...................................................................................... 221 B.3 Statistics of the maximum solute concentration in stream at the time of detection at the compliance boundary, and detection time at compliance boundary in different number of realizations for Scenario SW-GW-2a (larger horizontal fracture apertures). A Maximum solute concentration, B detection time. “Minimum”, “Maximum” and “Average” represents the smallest, the largest and the arithmetic average values of each parameter, for each set of additional realizations. ............................................................. 222 C.1 Coupled steady-state SW-GW flow in the first realization of Scenario SWGW-2b (larger vertical fracture apertures). A Surface water depth (m), B xiii  Figure C.2  Figure C.3  Figure D.1  Figure E.1  Figure F.1  Figure G.1  surface flow velocity Vx (m/d), C flux exchange (m/d) between surface and subsurface, D saturation profile in the subsurface, E hydraulic head (m) in the subsurface domain, and F groundwater hydraulic head (m) and flow velocity vectors in cross-section at Y=100m (beneath the streambed). Dimension (m). ................................................................................................ 223 Statistics of water depth, discharge and flow velocity at stream outlet in different number of realizations for Scenario SW-GW-2b (larger vertical fracture apertures). A Water depth, B total discharge, and C mean velocity (Vx). “Minimum”, “Maximum” and “Average” represents the smallest, the largest and the arithmetic average values of each parameter, for each set of additional realizations...................................................................................... 224 Statistics of the maximum solute concentration in stream at the time of detection at the compliance boundary, and detection time at compliance boundary in different number of realizations for Scenario SW-GW-2b (larger vertical fracture apertures). A Maximum solute concentration, B detection time. “Minimum”, “Maximum” and “Average” represents the smallest, the largest and the arithmetic average values of each parameter, for each set of additional realizations. ............................................................. 225 Statistics of the maximum solute concentration in stream at the time of detection at the compliance boundary, and detection time at compliance boundary in different number of realizations for Scenario SW-GW-3a (smaller αl in stream water). A Maximum solute concentration, B detection time. “Minimum”, “Maximum” and “Average” represents the smallest, the largest and the arithmetic average values of each parameter, for each set of additional realizations...................................................................................... 226 Statistics of the maximum solute concentration in stream at the time of detection at the compliance boundary, and detection time at compliance boundary in different number of realizations for Scenario SW-GW-3b (smaller αt in stream water). A Maximum solute concentration, B detection time. “Minimum”, “Maximum” and “Average” represents the smallest, the largest and the arithmetic average values of each parameter, for each set of additional realizations...................................................................................... 227 Statistics of the maximum solute concentration in stream at the time of detection at the compliance boundary, and detection time at compliance boundary in different number of realizations for Scenario SW-GW-4 (lower detection threshold 1.0×10-4g/L). A Maximum solute concentration, B detection time. “Minimum”, “Maximum” and “Average” represents the smallest, the largest and the arithmetic average values of each parameter, for each set of additional realizations. ............................................................. 228 Statistics of the maximum solute concentration in stream at the time of detection at the compliance boundary, and detection time at compliance boundary in different number of realizations for Scenario SW-GW-5 (source zone on bedrock). A Maximum solute concentration, B detection time. “Minimum”, “Maximum” and “Average” represents the smallest, the largest and the arithmetic average values of each parameter, for each set of additional realizations...................................................................................... 229  xiv  ACKNOWLEDGEMENTS  First, I would like to express the deep appreciation from my heart to my co-supervisors: Dr. Leslie Smith and Dr. Roger Beckie for giving me the great opportunity, the generous support, the valuable advice, the love and the patience over the whole period of my PhD study at UBC. I also would like sincerely to thank my Supervising Committee member, Dr. Ulrich Mayer for his kind assistance in implementing my research project, and Dr. Marc Bustin and Dr. Richard Chase in my Candidacy Exam Committee, Dr. Robert Millar, Dr. Erik Eberhardt, Dr. Michael Ward and Dr. Kent S. Novakowski in my final Doctoral Examination Committee. Then, I would like to say thanks to the kind donors of UBC Lorntzsen and MacKay Scholarship, George E Winkler Memorial Scholarship, Lisle and Sheila Jory Award in Geological Engineering, which I won throughout my course of study here. In addition, the research was funded by the grants from the Natural Sciences and Engineering Research Council of Canada and the Canadian Water Network of the Networks of Centres of Excellence of Canada, as well as other grants from Dr. Leslie Smith and Dr. Roger Beckie, and it was implemented with the world-class computer resources of the Westgrid of Canada. I appreciate the great help from my research collaborators Dr. Rob McLaren and Dr. Young-Jin Park at the University of Waterloo in providing the numerical codes: FRAC3DVS and HydroGeoSphere used in my research, and the assistance from the Westgrid’s staff Dr. Roman Baranowski in running the programs. Also I am grateful for Dr. Jean Cho in providing the data used in my research and my UBC Hydrogroup colleagues Rich Amos, Cindy Starzyk, Sergi Molin and EOS colleagues Sukhi Hundal, Karim Damani, David Jones and Albert Cui for their help throughout my study. Finally, I want to say thanks to everyone in my life who has been giving me the love, care and support in achieving my goal of success. xv  1 INTRODUCTION  1.1 Motivation and review of previous research  Monitoring of leachate migration from landfills or other waste disposal facilities is carried out to minimize the risk of contaminants entering surface water and groundwater systems. The design of sampling frameworks for surface water and groundwater has been studied for a number of years, especially as it relates to facilities located above permeable sediments. However, guidance for design of monitoring networks in fractured bedrock or in an integrated porous and fractured system is quite limited. This study investigates the factors that influence detection of the contamination using a groundwater monitoring network and an integrated surface water and groundwater monitoring network for a system where a permeable, unconfined porous aquifer is located above fractured bedrock in a headwater stream catchment. Monitoring for contamination in both surface water and the groundwater system is considered. In this section, previous research is briefly reviewed for each of the areas relevant to this thesis: (1) groundwater flow and solute transport in fractured media, (2) design of contaminant monitoring networks for groundwater systems in granular porous media, (3) design of contaminant monitoring networks for fractured media, and (4) surface water monitoring systems. Review papers describing the evolution of concepts and models used to predict groundwater flow and solute transport in fractures and fractured rock systems include Smith and Schwartz (1993), Chiles and de Marsily (1993), Adler and Thovert (1999), Berkowitz (2002), Bodin et al. (2003a, 2003b), Neuman (2005) and MacQuarrie and Mayer (2005). Of 1  particular importance to this study has been the significant progress in the numerical analysis of three-dimensional, variably saturated groundwater flow and solute transport in discretelyfractured porous media (Therrien and Sudicky, 1996; Therrien et al., 2003, Therrien et al., 2005). An early comprehensive review of the literature on the groundwater monitoring network designs was carried out by Loaiciga et al. (1992). A valuable summary of methodologies for groundwater monitoring network design and associated guidance for implementation was presented by the Task Committee on the State of the Art in Long-term Groundwater Monitoring Design (ASCE, 2003). This report noted significant progress has been made in the design of contaminant monitoring networks for groundwater systems in porous media by using statistical and numerical models. Massmann and Freeze (1987a, 1987b) presented a framework based on a risk-cost-benefit analysis for landfill management and assessment of alternative groundwater monitoring strategies. They formulated the problem to maximize the net present value of a stream of costs (facility construction and operation), risks (probability of failure) and the benefits (revenue for service), and integrated monitoring into the framework by calculating the probabilities of plume detection and the benefits of early detection, given a particular monitoring network. Their framework allows for the identification of the best design among a number of different monitoring well locations. The application they presented considered only advective solute transport. Meyer and Brill (1988) adopted the framework of Massmann and Freeze (1987a) and developed a framework for optimal groundwater monitoring network design in a porous medium by linking simulation and optimization models. They used the objectives of maximization of the probability of detection of plumes and minimization of the number of monitoring wells, to examine trade-offs between contaminant plume detection and the 2  number of wells to be installed. Meyer et al. (1994) extended this optimization analysis to examine the initial detection of groundwater contamination at a waste disposal facility by optimizing the multiple conflicting objectives including maximization of the probability of detecting a contaminant leak, minimization of the number of monitoring wells, and minimization of the expected area of contamination at the time of detection. They solved the problem under conditions of uncertainty in both the leak source location and the spatial distribution of the saturated hydraulic conductivity in a two-dimensional porous medium. Their analysis ignored the effect of dispersion. Storck et al. (1997) extended the work of Meyer et al. (1994) to a three-dimensional porous medium, but also did not consider the effect of dispersion. Hudak and Loaiciga (1993) carried out an optimal groundwater monitoring network design to facilitate early detection near the boundaries of a contaminant source in a multilayered aquifer system. Hudak et al. (1995) developed an optimal groundwater monitoring network design for a regional-scale aquifer, where the zone of potential contamination exceeded a distance of 1km from the contaminant source. Yenigül et al. (2005) carried out a study to assess the reliability of groundwater monitoring systems in a two-dimensional porous medium. They concluded that lateral dispersivity of the medium has a significant influence on the reliability of the monitoring system, that the number and the location of the groundwater monitoring wells are dependent on the heterogeneity of the medium and size of the contaminant leak, and that the reliability of the common practice of adopting three downgradient monitoring wells is inadequate for detection of groundwater contamination due to landfills. Further research into groundwater monitoring network designs, using methods not relevant to the study here, can be found in Morisawa and Inoue (1991), Dougherty and  3  Marryott (1991), Marryott et al. (1993), James and Gorelick (1994), Cieniawski et al. (1995), Reed et al. (2000) and Hudak (2002). Despite the fact that the conceptual framework for groundwater monitoring network design has been significantly advanced, the previous studies on groundwater monitoring network design were limited to fully-saturated isotropic and homogeneous porous media. None of these studies dealt with the problem in a three-dimensional variably-saturated aquifer system, and none of them considered the effects of dispersion and diffusion in the design of monitoring networks in three-dimensional media because of the huge computation effort. Few studies have examined monitoring network design in fractured bedrock. Jardine et al. (1996) adopted the decision framework developed by Massmann and Freeze (1987a) to examine relationships between monitoring network design, probabilities of plume detection, and the characteristics of the fracture system. They concluded that locating monitoring zones in boreholes at the locations where fractures had the largest apertures, or in zones with higher fracture density, yielded higher probabilities of plume detection. Their analysis was limited to two-dimensions and did not consider the influence of flow and solute transport in the rock matrix, or the effect of dispersion and diffusion in solute transport. Most of the research on surface water monitoring has focused on the issues related to sampling frequency and sampling network design. Dixon and Chiswell (1996) reviewed surface water quality monitoring program design publications from 1970 to early 1995. They pointed out that sampling frequency is of great importance in designing a monitoring program, and it must be tied to river flow (event-based sampling for an event flow, regular sampling for a steady flow), and that the spatial distribution of sampling sites within a catchment should be determined on the basis of the monitoring program goals. Claessen 4  (1997) compared the similarities and differences in monitoring surface water and groundwater including the layout of monitoring networks and monitoring frequencies, and concluded that monitoring systems for surface water and groundwater were seldom integrated into one system. Kristensen and Bøgestrand (1996) and Alexander et al. (1998) summarized surface water quality monitoring practice in Europe and the US, respectively. They identified that the sampling frequency differs substantially depending on the purpose of the monitoring programs, and it can be annually, quarterly, monthly or even weekly. Studies on the issues related to surface water and groundwater interactions have been carried out by hydrologists and hydrogeologists for a few decades. Winter (1999) and Sophocleous (2002) summarized the effects of geology on the interaction of surface water and groundwater, and they concluded that the local geological conditions control the streamaquifer recharge-discharge processes through the surface water beds. Oxtobee and Novakowski (2002) carried out a field study of groundwater/surface water interaction to examine the exchange processes between fractured bedrock aquifer and a bedrock stream, and they discovered that groundwater discharge in a bedrock stream environment occurs primarily through discrete point sources associated with open fractures, as compared to more diffuse, or continuous seepage zones often observed in a porous media environment. Conant et al. (2004) presented a case study of a groundwater plume discharging to a river, and concluded that sampling of water in streambed deposits is necessary to accurately characterize plumes discharging to the river, and sampling of surface water alone is insufficient due to the rapid dilution of contaminants in the river. Stephenson et al. (2006) implemented a study of hydraulic characterization of injection/extraction well system for steam enhanced remediation of contamination in a limestone having low matrix porosity and a sparse distribution of fractures, and identified that success of the remedial method 5  dependant on the fracture network for delivery of remedial fluid or extraction of contaminants. Chapman et al. (2007) carried out a field investigation to examine the impact of surface water and groundwater interaction on contaminant transport, and they suggested that surface water and groundwater concentration measurements together provide greater confidence in identifying and quantifying transport processes, rather than groundwater measurements alone. Significant progress has been made recently in developing numerical models such as HydroGeoSphere (Therrien et al., 2005) and in applying these types of numerical codes to simulate fully-coupled surface-subsurface flow and solute transport in porous media and fractured media for characterizing and quantifying groundwater contributions to stream flow, surface-subsurface exchange fluxes, as well as rainfall-runoff event responses at various scales (VanderKwaak and Loague, 2001; Bauer et al., 2006; Kollet and Maxwell, 2006). The most recent attempts to simulate flow and solute transport in coupled surface water and groundwater systems were carried out by Jones et al. (2006), Jones et al. (2008), Li et al. (2008) and Sudicky et al. (2008) in watershed-scale unconfined porous aquifers. The advance of the numerical techniques provides the tool for implementation of this study, and the findings in these numerical studies help to better understand the issues related to water and solute exchanges between the surface and the subsurface. However, none of the previous research simulated coupled surface-subsurface flow and solute transport in a system with both porous media and fractured media. No studies have examined factors influencing groundwater monitoring or an integrated groundwater and surface water monitoring in a catchment with an unconfined porous aquifer overlying fractured bedrock.  6  For detection of contaminant leakage from a landfill or a hazardous waste disposal site, a monitoring network may be designed to simply monitor contamination in surface water rather than groundwater, based on the perception of faster contaminant transport in surface water, or the monitoring network may be designed to monitor contamination in a shallow unconfined aquifer rather than the deeper fractured bedrock based on the lower cost of drilling wells. However, there is a lack of a scientific basis to justify the monitoring network design and the strategies for monitoring contamination in surface water and groundwater, or in an upper unconfined porous aquifer and underlying fractured bedrock.  1.2 Objectives and contributions  Figure 1.1 is a conceptual diagram illustrating the problem examined in this thesis for monitoring the contamination in a headwater stream catchment where an unconfined porous aquifer overlies fractured bedrock. The downstream limit of the domain represents a regulatory compliance boundary. A contaminant source is present in the upper part of the catchment. A performance monitoring network is to be designed to detect contamination prior to detection at the compliance boundary. Monitoring sites can be positioned in both the surface stream and groundwater wells, and the groundwater wells can be completed in either the upper porous aquifer or the fractured bedrock. There is lack of guidance for practitioners to follow in practice to monitor the contamination in such a system. The objectives of this research are to investigate what factors are important and how they interact to influence the design of monitoring networks to detect contamination in both surface water and groundwater. Factors influencing design of surface water and groundwater monitoring networks that are examined include longitudinal and transverse dispersivities in 7  stream water, hydraulic conductivity of unconfined porous aquifer, stochastic properties of fracture networks (fracture connectivity and apertures), rock matrix permeability, contaminant detection threshold, and source zone location. The contributions of this study are: (1) This study is the first to evaluate an integrated surface water and groundwater monitoring network in a catchment with an unconfined porous aquifer overlying fractured bedrock, and to examine how the factors of both surface water and groundwater systems interact in influencing the integrated surface water and groundwater monitoring; (2) This study is the first to fully account for advective and dispersive/diffusive effects in simulation of contaminant transport in both surface water and groundwater and in both porous aquifer and fractured bedrock for design of an integrated surface water and groundwater monitoring network. All previous studies were limited to groundwater monitoring network designs in two or three-dimensional porous media, typically ignoring the effect of dispersion, or were limited to the investigation of factors influencing groundwater monitoring in fractured media in two-dimensions, also ignoring the effect of dispersion. (3) The study helps to better understand coupled surface water and groundwater flow and transport patterns in a three-dimensional variably-saturated porous and fractured system, and provides guidelines for design of integrated monitoring networks for detection of contamination in surface water and groundwater.  1.3 Thesis outline  The thesis covers two topics:  8  •  Evaluation of groundwater (GW) monitoring networks, and the factors that influence groundwater monitoring in a catchment with an unconfined porous aquifer overlying fractured bedrock, without integrated surface water;  •  Evaluation of integrated surface water and groundwater (SW-GW) monitoring networks, and the factors influencing integrated surface water and groundwater monitoring in a headwater stream catchment with an unconfined porous aquifer overlying fractured bedrock.  The thesis is composed of six chapters: Chapter 1 provides the motivation, a brief review of previous research, the research objectives and contributions, as well as the outline of the thesis; Chapter 2 describes the methodology for evaluation of a groundwater (GW) monitoring network, the numerical model and parameters used for simulations of groundwater flow and transport without coupling to surface water; Chapter 3 presents the results of simulations of multiple realizations for each of 10 scenarios of steady-state groundwater flow and non-reactive contaminant transport without coupling to surface water, and discusses how the factors interact to influence monitoring of contamination in groundwater systems; Chapter 4 describes the methodology for evaluation of an integrated surface water and groundwater (SW-GW) monitoring network, the numerical model and the parameters used for simulations of coupled surface water and groundwater (SW-GW) flow and transport; Chapter 5 presents the results of simulations of multiple realizations for each of 8 scenarios of coupled steady-state surface water and groundwater (SW-GW) flow and non-reactive contaminant transport, and discusses how the factors 9  in surface water and groundwater systems interact to influence the integrated SW-GW monitoring; Chapter 6 summarizes the conclusions regarding groundwater monitoring and integrated surface water and groundwater monitoring, discusses the implications of the results for hydrogeological practice, as well as providing suggestions for the future research.  10  Performance Monitoring Network  Source  Unconfined Porous Aquifer  Vertical Exaggeration: 8×  Fractured Bedrock  Location of groundwater monitoring wells Location of surface water monitoring stations  Regulatory Compliance Boundary  Figure 1.1 Conceptual diagrams of the catchment and the performance monitoring network  11  2 METHODOLOGY TO EVALUTE GROUNDWATER MONITORING NETWORKS  2.1 Introduction  The evaluation of a performance monitoring network in a groundwater system includes an estimation of the likelihood for any given monitoring well location that the contaminant concentration will exceed a detection threshold prior to the plume reaching a regulatory compliance boundary. This chapter presents the methodologies used to examine the likelihood of contaminant detection in a catchment where an unconfined porous aquifer is located above fractured bedrock. Monte Carlo simulation techniques are used to estimate plume detection probabilities as a function of the location of monitoring wells downgradient from a contaminant source zone. Section 2.2 describes the model domain of the groundwater system; Section 2.3 specifies the boundary conditions for groundwater flow and solute transport; Section 2.4 gives a brief description of the numerical model used; Section 2.5 outlines the approach for evaluation of groundwater monitoring network; Section 2.6 provides the details of the scenarios investigated and the model parameters used in the scenarios; Section 2.7 is an evaluation of the model domain discretization and number of Monte Carlo realizations by means of test simulations; and Section 2.8 summarizes the methodologies for evaluation of the groundwater monitoring network.  12  2.2 Model domain for the groundwater system  The model domain used for simulation of groundwater flow and solute transport is shown in Figure 2.1. The model domain is 500m long and 200m wide with a V-shaped land surface. The elevation is 100m at the base of the domain, while the land surface elevations at the upstream inflow boundary (X=0m) and at the downstream outflow boundary (X=500m) are shown in Figure 2.2. The top surface is symmetric about the centerline of the domain with a slope of 3m/500m in the X direction (parallel to the valley axis) and 1m/100m in the Y direction (across the valley axis). The subsurface is composed of an unconfined homogeneous and isotropic porous aquifer overlying fractured bedrock. The bedrock is 15m thick (from elevation at 100m to 115m), and the unconfined porous aquifer ranges in thickness from 6m to 10m. The groundwater flow system is at steady-state. The fractured bedrock contains a fracture network with randomly-distributed fractures embedded in a homogeneous and isotropic rock matrix with low (but finite) permeability. The fracture network is composed of three fracture sets (one horizontal set and two vertical sets) with fixed orientations parallel to the XY, XZ and YZ planes. The three sets of fractures are assumed to be spatially uncorrelated (no joints are persistent). All fractures are represented by idealized two-dimensional parallel plates with a rectangular form. No fault zones are included that could extend across the entire domain from the upstream to the downstream boundary. The procedures and statistical distributions used to generate the fracture sets are described in detail in Section 2.6. The size of the domain and the thickness of the bedrock and unconfined aquifer to be modeled are chosen based on the availability of the computer resources for handling the 13  simulations of the three-dimensional catchment-scale variably-saturated groundwater (and later coupled surface water and groundwater in Chapter 4 and 5) in the discretely-fractured bedrock with thousands of fractures and about three million elements.  2.3 Boundary conditions for groundwater flow and solute transport  2.3.1 Boundary conditions for flow  The boundary conditions for groundwater flow are: (1) Constant hydraulic heads at the upstream inflow boundary (X=0m) and the downstream outflow boundary (X=500m) specified with the values shown in Figure 2.2 that vary linearly from a high at the valley sides and a low at the valley centre. In this way, both horizontal and vertical hydraulic gradients exist on the water table at each of the boundaries. (2) No flow boundaries representing groundwater divides at the side boundaries (Y=0m and Y=200m) and an impermeable boundary at the base of the model domain (Z=100m). (3) A very small recharge (1.0×10-4m/yr) applied on the ground surface. Recharge at this rate improved the stability of the numerical solution, without affecting groundwater flow and solute transport in any significant way.  2.3.2 Boundary conditions for solute transport  The boundary conditions for groundwater transport are: 14  (1) A non-reactive contaminant source defined as a zone with a constant concentration of 1.0g/L and a size of 20m×20m×2m, extending in the unconfined porous aquifer from X=50m to X=70m, Y=60m to Y=80m (see Plot A in Figure 2.3) and from Z=121.5m to Z=123.5m in the region just above the water table (see Plot B in Figure 2.3). The size of the source zone is chosen relative to the size of the model domain and also the space of the fractures beneath the unconfined porous aquifer. (2) Zero concentration at the inflow boundary plane at X=0m. (3) The 3rd type boundary at the outflow boundary at X=500m, and no mass flux at all other boundaries. (4) Zero initial (background) concentration.  2.4 Description of the numerical model: FRAC3DVS  The numerical code FRAC3DVS was chosen to simulate steady-state groundwater (GW) flow and non-reactive solute transport in the porous and fractured system shown in Figure 2.1. It is a robust numerical model for simulation of three-dimensional, variably-saturated groundwater flow and solute transport in porous and discretely-fractured media, and has been applied successfully in solving many groundwater flow and transport problems. The details of the code including the governing equations can be found in Therrien and Sudicky (1996) and Therrien et al. (2003). Here is only a brief description about how the code works. In FRAC3DVS, variably-saturated flow in porous media (the unconfined aquifer and the rock matrix) is described with Richards equation, and variably-saturated flow in fractures is described with an analogy of Richards equation. Solute transport in both the porous medium and discrete fractures is described by the classical advection-dispersion equation. In this 15  study, the unconfined aquifer is variably saturated, and the fractured bedrock is fullysaturated with the flow boundary conditions specified (see Figure 3.1 in Chapter 3). The unconfined porous aquifer and the lower permeable rock matrix are represented by three-dimensional blocks, and the fractures are represented by two-dimensional rectangular parallel plates. The horizontal and vertical fractures are incorporated in the grid by superimposing two-dimensional face elements onto the three-dimensional grid of block elements. In order to fully couple the fractures with the porous matrix, faces and blocks share common nodes along the fracture walls, and the model assumes that along the surface comprising a fracture-matrix interface, hydraulic head and concentration in the matrix along the surfaces are equal to the hydraulic head and concentration in the adjoining fractures. The model takes into account advective transport, mechanical dispersion and molecular diffusion in both porous media and fractured media. In solving the flow and transport equations with FRAC3DVS, both finite element and finite difference approaches are available. Because of lower CPU memory requirements and faster execution, the finite difference option was used in the simulations reported here, in combination with an adaptive time-stepping control with a maximum time-step of 90 days. Therrien and Sudicky (1996) performed a variety of groundwater simulations in discretelyfractured porous media with this numerical program, and concluded that the finite-element and finite-difference approaches yield identical results. Test simulations using the first realization of the Base Case in this study also showed there was no difference in flow and transport solutions using the two approaches.  16  2.5 Approach for evaluation of groundwater monitoring network  The groundwater monitoring network (Figure 2.4) consists of 171 potential well locations for performance monitoring. The monitoring locations are placed in 9 rows between the source zone and the compliance boundary (from X=100m to X=500m), with 19 wells in each row (from Y=10m to Y=190m). The distance is 50m between neighbouring rows in the X direction, and 10m between neighbouring wells in each row in the Y direction. The smaller well spacing is selected in the Y direction in order to get better estimation of detection of the contaminant plume in that direction considering the possible effects of lateral advective transport in the fractured bedrock and lateral dispersion in the unconfined aquifer. To evaluate the groundwater monitoring network, a detection threshold is defined in terms of a concentration value above the background concentration, and expressed relative to the source concentration, such as 0.1% of the source concentration. It is applied at both the compliance boundary and the performance monitoring wells to examine if the contaminant is detected at the compliance boundary and also at a monitoring well. It is assumed that all the performance monitoring wells extend to the basal boundary with full vertical resolution in both the unconfined aquifer and the fractured bedrock. Monitoring points in each of the wells correspond to the nodes in the numerical grid at the well location. Each monitoring point (at each node) in a monitoring well is assumed to represent a monitoring port in a multi-level sampler (sampling at each port or node). This approach neglects the additional uncertainty in plume detection related to the vertical position of a monitoring zone relative to the location of the plume. In addition, within FRAC3DVS, it is assumed that the monitoring wells do not affect the flow field or influence solute transport. The nodal concentrations are compared directly with the detection threshold 17  (not weighted based on discharge in the wells), and no accord is taken of effects of vertical mixing and borehole dilution in the wells. The monitoring network is evaluated by assuming that the regulatory compliance boundary is located at X=500m, corresponding to the outflow boundary of the simulation domain. Detection of contaminant at the compliance boundary is assumed to occur if the concentration(s) at any of the nodes on this boundary exceeds the detection threshold. The transport simulation is terminated once this happens to permit identification of which unit (the unconfined aquifer or the bedrock) the contaminant is most likely to be detected in, before the contaminant reaches the compliance boundary. Detection of the contamination in a performance monitoring well is assumed to occur if the concentration value at any node representing a monitoring port exceeds the detection threshold concentration during the same timestep that contaminant is detected at the compliance boundary. This time is called the “detection time” hereafter. This detection criterion is appropriate for the continuous source considered here, but may be inappropriate for an instantaneous or a finite duration source, in which zones close to the release could be flushed by the time the plume reaches the compliance boundary. Detection probabilities are computed for each monitoring well location as the proportion of realizations in the Monte Carlo simulation where the detection is recorded at that location. The detection probabilities are evaluated separately for the unconfined aquifer and the fractured bedrock. It is important to emphasize that the detection probabilities in the performance network are not defined in an absolute sense, but rather as the probability of detecting the plume at the same time the contamination is detected at the compliance boundary. For example, if a performance monitoring well at location “A”, as shown in Figure 2.4, has a detection probability of 0.7 in the bedrock, this means that there is a 70% probability of detecting the 18  plume in the bedrock at the location “A” at the same time the contaminant is detected anywhere on the compliance boundary. The result also implies that there is a 30% probability that monitoring at the location “A” will fail to detect contamination at the time it is detected at the compliance boundary. If the detection probability in a performance monitoring well is 0.7 in the bedrock and 0.1 in the unconfined aquifer, this indicates that the contaminant is more likely to be detected in the bedrock at this location as the plume reaches the compliance boundary. If the detection probability at the location “A” is 0.9 in the bedrock and 0.9 in the unconfined aquifer, this indicates the contaminant is very likely to be detected by both the bedrock monitoring points, and by those in the unconfined aquifer.  2.6 Scenarios and model parameters  2.6.1 Base Case and the parameters  The parameters describing the unconfined porous aquifer (Table 2.1), the fracture network (Table 2.2), dispersion within the fracture planes (Table 2.3), and flow and transport in the rock matrix (Table 2.4) are representative of a site with a moderately permeable unconfined aquifer overlying low permeability bedrock containing a relatively sparse network of open fractures. In all the simulations, the contaminant is assumed to be nonreactive, and the free-water diffusion coefficient is 1×10-9m2/s. Distinct values of tortuosity (τ) are applied for diffusion in the unconfined aquifer and the fractured bedrock. The van Genuchten model (van Genuchten, 1980) is used to characterize the moisture retention curve in the unsaturated zone of the unconfined aquifer. The parameters listed in Table 2.1 are representative of a sandy silt. 19  For each fracture set in the bedrock, fracture center locations are generated from a uniform distribution, fracture lengths are generated from an exponential distribution, and fracture aperture from a lognormal distribution by using the parameters in Table 2.2. These distributions are uncorrelated and independent, which neglects the impact of the correlation of fractures on reducing connectivity of the fracture network and permeability of the bedrock (de Dreuzy et al., 2004). The fracture length and fracture aperture distributions are truncated with defined lower and upper bounds. The lower bound on the fracture length distribution is chosen to be consistent with the feasible grid discretization used in the numerical solution, and described below. The fracture statistics follow the guidance that can be obtained from field data presented by Raven (1986) and Bonnet et al. (2001). A domain extending 30m in the X and Y axes and 15m in the Z axis beyond the boundaries of the flow domain was used for generation of the fracture network. In this way, fractures whose centers lie outside the domain but extend into the domain are included, thus eliminating boundary effects. Figure 2.5 shows the histograms of fracture centre locations, lengths and apertures in the first realization of the Base Case. Plot A, B and C show the data in uniform distributions for the centre locations of fractures (parallel to the XY plane) in the X, Y and Z axes; Plot D and E show the data for the lengths of fractures (parallel to the XY and XZ planes), and the data is in exponential distributions although the plots do not have this appearance because of the truncation bounds selected for the fracture lengths; Plot F shows the data in lognormal distribution for the apertures of fractures (parallel to the XY plane). The data for other sets of fractures is not shown here because of the similarity in the distributions. The fracture network in the first realization of the Base Case incorporates a total of 3897 fractures (1554 fractures parallel to the XY plane, 1158 fractures parallel to the XZ plane and 1185 fractures parallel to the YZ plane). The network is shown in Figure 2.6. While the three 20  dimensions shown in Plot A of Figure 2.6 might give the impression of a dense network of open fractures; the cross-sections shown in Plot B and C of Figure 2.6 and the scan-line density in Table 2.2 indicate a relatively sparse network of open fractures. The saturated hydraulic conductivity for a single fracture is calculated with the “Cubic Law” (Snow 1969):  Kf = ρg(2b)2/(12µ)  (2.1)  where Kf is the saturated hydraulic conductivity of a single fracture (m/s), 2b is the fracture aperture (m), ρ is the fresh water density (1000kg/m3), µ is the dynamic viscosity of water [1.124×10-3kg/(m⋅s)], and g is the acceleration due to gravity (9.8m/s2). The bulk hydraulic conductivity (Kx, Ky and Kz) of the bedrock in the X, Y, Z directions is calculated from the simulation results using Darcy’s law:  Ki = Qi /[Ai(dh/dxi)]  (2.2)  where Ki is the bulk hydraulic conductivity in i direction (i = x in the x-direction, i = y in the y-direction, and i = z in the z-direction), Qi is the volume of water (m3/s) flowing through the cross-section area Ai (m2) under the specified hydraulic gradient dh/dxi (equal to 1/500) across the fractured bedrock domain in the i direction. The detection threshold for the Base Case is predetermined to be 1.0×10-3g/L, which is 0.1% of source zone concentration (1.0g/L). This threshold value is chosen to represent a contaminant concentration attributable to the source zone and also with the consideration of  21  the practical detection capability of many common contaminants relative to their background concentrations.  2.6.2 Other scenarios and the parameters  Table 2.5 outlines 9 scenarios used to examine the factors influencing the probabilities of plume detection for the groundwater monitoring in comparison to the Base Case. These factors include the hydraulic properties of the unconfined aquifer, the fracture network and rock matrix, the contaminant detection threshold, and the source zone location. All the scenarios use the same domain and the same parameters as in the Base Case, except for the key parameters to be investigated. Scenario GW-1 is designed to investigate the influence of hydraulic conductivity (Kuc) of the unconfined porous aquifer on detection probabilities by increasing Kuc one order of magnitude from the Base Case. Scenario GW-2a and Scenario GW-2b are designed to examine the influence of fracture network connectivity on detection probabilities due to the reduction of the densities by a factor of 2 for the set of horizontal fractures and for the sets of vertical fractures, respectively. Scenario GW-3a and Scenario GW-3b are designed to examine the influence of the aperture of the horizontal and the vertical fracture sets, respectively, by reducing the mean aperture one order of magnitude from the Base Case but keeping the lower and upper bounds unchanged in the distribution. Scenario GW-4 is designed to investigate the influence of matrix permeability (Km) of the bedrock on detection probabilities with Km being increased by 2 orders of magnitude from the Base Case. Scenario GW-5a and Scenario GW-5b are designed to investigate the influence of detection threshold concentrations on detection probabilities by increasing and reducing the 22  detection threshold by one order of magnitude from the Base Case. Scenario GW-6 is designed to examine the influence of the location of the contaminant source zone on detection probabilities with the source zone sitting on the bedrock.  2.7 Discretization of model domain and Monte Carlo simulations  2.7.1 Discretization of model domain  The fractured bedrock is discretized into block elements with the sizes ∆X=∆Y=∆Z=1m. The unconfined aquifer is discretized into block elements with sizes ∆X=∆Y=1m and ∆Z=0.6m-1m to accommodate the variable ground surface. The domain contains 2.5 million elements with 2,618,226 nodes, which is the upper limit that the available state-of-art Westgrid computer resources in Canada can handle. Using this discretization, the bedrock unit is represented by 15 layers, and the unconfined aquifer is divided vertically into 10 layers. A small-domain simulation was carried out to test the suitability of a grid spacing of 1.0m in the bedrock for the Base Case analysis. As shown in the plots of Figure 2.7, the domain, 10m by 10m by 10m, incorporated 8 fractures, each with an aperture of 26.4µm. The rock matrix was assigned with a hydraulic conductivity of 5.8×10-10m/s. All other material properties were equal to those assigned in the Base Case. A constant hydraulic gradient of 0.3 was specified between the inflow boundary (X=0m) and the outflow boundary (X=10m) for a steady-state flow simulation in a fully-saturated domain. The source zone (2m by 4m by 2m) with the specified constant concentration of 1.0g/L was located at X=0m-2.0m, Y=4.0m8.0m and Z=4.0m-6.0m. The simulation time was 10 years. In a sequence of simulations, the 23  domain was discretized into element sizes of ∆X=∆Y=∆Z=2m, 1m, 0.5m and 0.25m, respectively. Figure 2.8 shows the nodal concentrations versus element size in a monitoring well located at X=6m, Y=6m, Z=0-10m. The results demonstrate that the discretization with element size of 2m is inappropriate for the fractured domain while nodal concentrations have converged when element sizes 1m or smaller are used. The differences between nodal concentrations are less than 0.5% in the fractures and 2% in the rock matrix for those computed with the element sizes of 1m and 0.5m. A transport simulation was also carried out for the fractured bedrock domain using the first realization of the Base Case with element sizes of ∆X=∆Y=∆Z=1m and ∆X=∆Y=∆Z=0.5m. Figure 2.9 shows that differences are small between the nodal concentrations computed with the element sizes of 1.0m and 0.5m, which demonstrates that the element size of ∆X=∆Y=∆Z=1m is sufficient for spatial discretization of the bedrock. This is the finest spatial discretization to accommodate as many fractures as possible in the bedrock unit and maintain the numerical accuracy in the catchment-scale aquifer system with the Westgrid computation capacity. Certainly, it would be great to use even finer grids when more powerful computation resource is available in the future for better characterization of solute transport in both the fractures and rock matrix. A recent study carried out by Weatherill et al. (2008) for simulation of solute transport in a single fracture with the code HydroGeoSphere showed that the fine spatial discretization at the scale of the fracture aperture is ideal to characterize the low matrix diffusion in bedrock. Furthermore, to test the suitability of the spatial discretization in the unconfined aquifer, the simulations were carried out in the first realization of the Base Case with the element sizes of ∆X=∆Y=1m but different ∆Z by dividing the unconfined aquifer vertically into 5, 6, 24  10, 12, 13, and 14 layers (the upper limit that the Westgrid network computer system used for this research could accommodate). The simulations were implemented by using the constant hydraulic heads 121.50m and 117.35m, respectively, for the inflow boundary and outflow boundary (the averages of the hydraulic heads in Figure 2.2). The results in Figure 2.10 and 2.11 show that differences in the nodal concentrations are very small in the saturated zone of the unconfined aquifer and in the bedrock, as well as in the unsaturated zone of the unconfined aquifer when a vertical element size of the ∆Z smaller than 1m is used, which demonstrates that the spatial discretization of the unconfined aquifer with element sizes of ∆X=∆Y=1m and ∆Z=0.6m-1m is appropriate. Also in Figure 2.10, large differences in the concentration calculations are observed only at the nodes in the unsaturated zone of the unconfined aquifer when the vertical element size of 1m-2m is used, but the differences beneath the water table (2.6m below the ground surface at the well location) are very small when using the vertical element size of 1m-2m or smaller than 1m. This result indicates that the vertical discretization with element size smaller than 2m is sufficient to obtain reliable results in contaminant transport and detection in the saturated zone of the unconfined aquifer. It should be mentioned that with the values of the dispersivities in Table 2.1, 2.3 and 2.4 and the spatial discretization in the unconfined aquifer (∆X=∆Y=1m and ∆Z=0.6m-1m) and in the bedrock (∆X=∆Y=∆Z=1m), the Peclet number constraints are slightly violated in the transverse directions in the unconfined aquifer as well as the fractures and the matrix in the bedrock. However, the results in Figure 2.8, 2.9, 2.10 and 2.11 demonstrated that with the small element sizes used, the violation of the Peclet number requirement does not affect the detection results in the system. Finally, by using the same spatial discretization of the model domain, the simulations of solute transport were also carried out for the first realization of the Base Case to test the 25  influence of different time-stepping procedures on the detection results. Figure 2.12 shows that the differences between the nodal concentrations in monitoring wells are very small when using a time weighting factor of 0.51, and a fully-implicit weighting factor of 1.0. Therefore, the time-weighting of 0.51 is selected for transport simulations in all the scenarios, considering that the fully-implicit weighting may generate more numerical dispersion.  2.7.2 Number of Monte Carlo realizations  The Monte Carlo approach is used to characterize the contaminant detection probability in the system attributable to uncertainty in defining the fracture networks in the bedrock, while the unconfined aquifer is treated deterministically. The approach involves generating a unique fracture network in the bedrock for each realization of a scenario, and then simulating groundwater flow and solute transport in the unconfined aquifer and the bedrock to obtain a realization of the contaminant plume for calculation of detection probabilities at the monitoring well locations in each of the units. To identify the minimum number of realizations required for reliable statistics in evaluation of the groundwater monitoring network, detection time at the compliance boundary and detection probabilities at the performance monitoring wells were calculated in a sequence of simulations with up to 500 realizations for the Base Case. Figure 2.13 shows the relationship between the number of realizations and the detection time at the compliance boundary, including the earliest (minimum) detection time, the latest (maximum) detection time, and the arithmetic average detection time, for each set of additional realizations. The result suggests that 200 realizations are sufficient to characterize 26  the variability in the detection time. Figure 2.14 shows the detection probabilities at the monitoring wells in the unconfined aquifer, for selected cross-sections at X=350m (Plot A) and Y=100m (Plot B). Figure 2.15 shows the detection probabilities at the monitoring wells in the bedrock, for selected cross-sections at X=100m (Plot A) and Y=100m (Plot B). The results indicate that the detection probabilities in both the unconfined aquifer and the bedrock converge when the sequence of simulations reach 200 realizations. Therefore, 200 realizations are used for simulations of flow and transport in all the scenarios reported here.  2.8 Summary  This chapter has described the methodologies for evaluating the groundwater monitoring network. They include description of the problem geometry and the numerical techniques used for simulation of steady-state, variably-saturated groundwater flow and non-reactive contaminant transport in a catchment with an unconfined homogeneous and isotropic porous aquifer overlying sparsely fractured bedrock. They also include the techniques for estimation of the probabilities of plume detection with the performance monitoring wells in the unconfined aquifer and the fractured bedrock. To investigate the factors influencing detection probabilities in this type of system, 10 scenarios are defined here, which differ from each other in values of key parameters. The factors investigated include hydraulic conductivity of the unconfined aquifer, fracture network connectivity, fracture apertures, rock matrix permeability, detection threshold, and source zone location. The results and implications on how these factors interact to influence the detection probabilities are presented in the next chapter. The following three key points are emphasized: 27  (1) The properties of the unconfined aquifer are treated in a deterministic fashion. Uncertainty in prediction of contaminant transport originates from the stochastic representation of the fracture networks in the bedrock, and it is characterized with Monte Carlo simulations of fracture networks within each of the 10 scenarios. (2) Detection is defined in terms of the probability of observing a contaminant concentration above the detection threshold in a performance monitoring well at the same time contamination is detected at the compliance boundary. The probability can also be interpreted in terms of the likelihood that at a given well location, the contaminant will reach the compliance boundary prior to detection in the performance well monitoring network. (3) Advective transport, hydrodynamic dispersion and molecular diffusion are accounted in the evaluation of detection probabilities in a three-dimensional porous medium overlying a fractured medium. This approach provides a more realistic estimation of contaminant plume geometry and detection. The effects of dispersion and diffusion were ignored in the earlier studies of monitoring network design. (4) Simulations are most reliable with saturated groundwater flow; values for unsaturated flow may be affected by mesh density.  28  Table 2.1 Parameters for unconfined porous aquifer in the Base Case Effective Saturated Hydraulic Porosity  Conductivity Kuc (m/s)  van Genuchten  Dispersivities, Tortuosity  Parameters  & Bulk Density  α = 1.1/m 0.38  β = 6.08  -6  2.5×10  αl = 2.0m αth = 0.2m αtv = 0.02m  γ = 0.84 Swr = 0.17  τ = 0.5 ρb = 1,690kg/m3  Table 2.2 Distribution parameters for fracture network in the Base Case  Fracture  Orientation  Set  Volumetric Density  Length (m)  Aperture (µm)  min, max, Lamda λ  min, max, mean, stdev  (volume of fractures m3 in per m3 of bedrock)  1  // XY Plane  2, 40, 0.01  20, 800, 250, 150  1.2×10-4  2  // XZ Plane  2, 15, 0.02  20, 800, 250, 150  1.2×10-5  3  // YZ Plane  2, 15, 0.02  20, 800, 250, 150  1.2×10-5  Scan-line Density in Average 6 fractures at a vertical scan-line of the depth of 15m 5 fractures at a horizontal scanline of per 100m 5 fractures at a horizontal scanline of per 100m  * “min”, “max” represents the truncation value of the distributions, “stdev” represents “standard deviation”.  Table 2.3 Dispersivities for transport simulation in fractures in the Base Case Parameter  Value (m)  Longitudinal dispersivity (αl)  2.0  Transverse dispersivity (αt)  0.2  29  Table 2.4 Parameters for simulations in rock matrix in the Base Case Effective Saturated Hydraulic Porosity  Conductivity Km (m/s)  Dispersivities, Tortuosity & Bulk Density αl = 0.5m αth = 0.05m  0.05  -10  5.8×10  αtv = 0.005m τ = 0.2 ρb = 2,000kg/m3  Table 2.5 Investigated scenarios and features for groundwater simulations Scenario  Intent  Features and Parameters  Base Case  Base Case  See Section 2.6.1, detection threshold 1.0×10-3g/L  GW-1  Hydraulic  Conductivity  Hydraulic conductivity of unconfined porous aquifer (Kuc)  of Unconfined Aquifer  increased by 1 order of magnitude  GW-2a  Fracture Connectivity  Density for horizontal fracture set reduced by factor 2  GW-2b  Fracture Connectivity  Densities for vertical fracture sets reduced by factor 2  GW-3a  Fracture Aperture  Mean of horizontal fracture apertures reduced by factor 10  GW-3b  Fracture Aperture  Mean of vertical fracture apertures reduced by factor 10  GW-4  Rock Matrix Permeability  Matrix permeability (Km) increased by 2 orders of magnitude  GW-5a  Detection Threshold  Increased to 1% of source concentration  GW-5b  Detection Threshold  Reduced to 0.01% of source concentration  GW-6  Source Zone Location  Contaminant source zone on top of the bedrock  30  Unconfined Porous Aquifer  Fractured Bedrock  Figure 2.1 Model domain for groundwater simulation in the system with an unconfined porous aquifer overlying fractured bedrock. Dimension (m).  125m  125m 124m  121.9m  122m  122m  121.9m 121.1m  121m 118.9m  Unconfined Porous Aquifer 115m  115m  115m  Fractured Bedrock  100m  115m  Fractured Bedrock  100m Upstream Inflow Boundary (at X=0m)  118.9m 118.3m Unconfined Porous Aquifer  100m  100m  Downstream Outflow Boundary (at X=500)  Figure 2.2 Diagrams of upstream inflow boundary and downstream outflow boundary in the model domain  31  A  Source Zone Concentration 1.0g/l  B  Source Zone Concentration 1.0g/l  Figure 2.3 Locations of source zone and water table in the domain. A Three dimensions, B cross-section at Y=80m. Dimension (m).  A  Location of Monitoring Well  Compliance Boundary  Figure 2.4 Locations of performance monitoring wells in groundwater monitoring network 32  B  Data of Centre Locations at the X Axis for Fractures in the XY Plane  Data of Centre Locations at the Y Axis for Fractures in the XY Plane  300  300  250  250  200  200  # of data  # of data  A  150 100  150 100 50  50 0  0  -30.0  44.3  118.6  192.9  267.2  341.5  415.8  490.1  -30.0  2.2  Location at the X axis (m)  D  Data of Centre Locations at the Z Axis for Fractures in the XY Plane  300  150  250  120  200  90 60  131.2  163.5  195.7  34.0  39.3  676.7  785.5  150 100  30  50  0  0  80.0 85.9 91.9 97.8 103.7 109.7 115.6 121.5 127.5 133.4  2.0  7.3  12.7  Location at the Z axis (m)  18.0  23.3  28.6  Fracture length (m)  F  Data of Lengths for Fractures in the XZ Plane  Data of Apertures for Fractures in the XY Plane  800  300 250 200  # of data  # of data  99.0  Data of Lengths for Fractures in the XY Plane  180  E  66.7  Location at the Y axis (m)  # of data  # of data  C  34.5  150 100 50  600 400 200 0  0 2.0  3.8  5.6  7.5  9.3  11.1  Fracture Length (m)  12.9  14.8  24.0  132.8  241.6  350.4  459.1  567.9  Fracture aperture (µm)  Figure 2.5 Histograms of random data in truncated distributions for fracture locations, lengths and apertures in the first realization of the Base Case. A, B and C Data for fracture centre locations in uniform distributions, D and E data for fracture lengths in exponential distributions, and F data for fracture apertures in a lognormal distribution.  33  A  B  C  Figure 2.6 Fracture network in the first realization of the Base Case for simulation of groundwater flow and transport. A Three dimensions, B and C cross-sections at Y=100m and X=250m. Aperture (m), Dimension (m).  34  B  A  ∆X=∆Y=∆Z =1m  Fractures and matrix  C  D  Source Zone Concentration (1.0 g/l)  Figure 2.7 Simulation of the small fractured domain for test of grid spacing. A Fractured domain, B spatial discretization, C simulated hydraulic head (m), and D source zone location. Dimension (m).  Nodal Concentrations vs Element Size in Small Fractured Domain (Monitoring Well at X=6m,Y=6m, Z=0-10m) 10  Z (m)  8 6 4 2 0 0.0  0.2  0.4  0.6  0.8  1.0  Concentration (g/l) ∆X=∆Y=∆Z=2.0m ∆X=∆Y=∆Z=0.5m  ∆X=∆Y=∆Z=1.0m ∆X=∆Y=∆Z=0.25m  Figure 2.8 Nodal concentrations of a monitoring well in the small fractured domain with different element sizes 35  Nodal Concentration vs Element Size in the Bedrock (Monitoring Well 5 at X=100m, Y=50m, Z=100m-115m) 115  Z (m)  112 109 106 103 100 0.00  0.02  0.04  0.06  0.08  0.10  0.12  0.14  Concentration (g/l) ∆X=∆Y=∆Z=0.5m  ∆X=∆Y=∆Z=1m  Figure 2.9 Nodal concentrations of monitoring well #5 in the first realization of the Base Case with different element sizes in the bedrock domain  Nodal Concentration vs Vertical Element Size in the Unconfined Aquifer (Monitoring Well 7 at X=100m, Y=70m, Z=100m-123.7m) 125  Z (m)  120 115 110 105 100 0.0  0.1  0.2  0.3 ∆Z=1.74m ∆Z=0.79m  0.4 0.5 0.6 Concentration (g/l) ∆Z=1.45m ∆Z=0.70m  0.7  0.8  0.9  1.0  ∆Z=0.87m ∆Z=0.62m  Figure 2.10 Nodal concentrations of monitoring well #7 in the first realization of the Base Case with different vertical element sizes in the unconfined aquifer  36  Nodal Concentration vs Vertical Element Size in the Unconfined Aquifer (Monitoring Well 143 at X=450m, Y=100m, Z=100m-121.3m) 125  Z (m)  120 115 110 105 100 0.0E+00  2.0E-04  4.0E-04  6.0E-04  8.0E-04  1.0E-03  1.2E-03  Concentration (g/l) ∆Z=1.26m ∆Z=0.53m  ∆Z=1.05m ∆Z=0.48m  ∆Z=0.63m ∆Z=0.45m  Figure 2.11 Nodal concentrations of monitoring well #143 in the first realization of the Base Case with different vertical element sizes in the unconfined aquifer  Nodal Concentration vs Time Weighting (Monitoring Well 122 at X=400m, Y=80m, Z=100m-121.8m) 124 120 Z (m)  116 112 108 104 100 0.0E+00  4.0E-04  8.0E-04 Concentration (g/l) TW=1.0  1.2E-03  1.6E-03  TW=0.51  Figure 2.12 Nodal concentrations of monitoring well #122 in the first realization of the Base Case with different time weightings (TW)  37  Detection Time at Compliance Boundary vs Number of Realizations 120  Time (yrs)  100 80 60 40 20 0 0  50  100  150  200 Minimum  250 300 Realizations Maximum  350  400  450  500  Average  Figure 2.13 Detection time at compliance boundary vs. number of Monte Carlo realizations in the Base Case. “Minimum”, “Maximum” and “Average” represents the earliest, the latest and the arithmetic average detection time, for each set of additional realizations.  38  A  Wells in Cross-Section at X=350m 0.6  Probability  Y=60m Y=70m  0.4  Y=80m Y=90m  0.2  Y=100m Y=110m  0.0 0  50  100  150  200  250  300  350  400  450  500  Realizations  B  Wells in Cross-Section at Y=100m 1.0  Probability  0.8  X=100m  0.6  X=200m  0.4  X=300m X=400m  0.2 0.0 0  50  100  150  200  250  300  350  400  450  500  Realizations  Figure 2.14 Detection probabilities of monitoring wells in the unconfined aquifer vs. number of Monte Carlo realizations in the Base Case. A Cross-section at X=350m, B cross-section at Y=100m.  39  A  Wells in Cross-Section at X=100m  Probability  0.8 Y=60m Y=70m Y=80m Y=90m Y=100m Y=110m  0.6 0.4 0.2 0.0 0  50  100  150  200  250  300  350  400  450  500  Realizations  B  Wells in Cross-Section at Y=100m  Probability  1.0 0.8  X=100m  0.6  X=200m  0.4  X=300m X=400m  0.2  X=500m  0.0 0  50  100  150  200  250  300  350  400  450  500  Realizations  Figure 2.15 Detection probabilities of monitoring wells in the bedrock vs. number of Monte Carlo realizations in the Base Case. A Cross-section at X=100m, B cross-section at Y=100m.  40  3 FACTORS INFLUENCING GROUNDWATER MONITORING  3.1 Introduction  This chapter presents the analysis of monitoring of groundwater (GW) contamination in the system with the unconfined porous aquifer overlying fractured bedrock, where there is no complementary monitoring of contamination in surface water. Factors examined include the hydraulic conductivity of the unconfined aquifer, the fracture network connectivity (a function of density, length, orientation and correlation), the aperture of the horizontal and the vertical fracture sets, the hydraulic conductivity of the rock matrix, the detection threshold concentration and the location of the source zone. The influence of these factors on detection probabilities of the performance monitoring network is evaluated on the basis of detection of the contamination in both the unconfined aquifer and the bedrock at the time that the contaminant reaches the compliance boundary with a concentration above the detection threshold. In this chapter, Section 3.2 presents a statistical summary of the generated fracture networks, the bulk hydraulic conductivities of the bedrock, as well as the time and location for detection of the contaminant at the compliance boundary in all the scenarios; Section 3.3 presents the discussion of the results of groundwater flow and transport & detection probabilities in the Base Case; Section 3.4 presents the discussion of the influence of hydraulic conductivity of the unconfined porous aquifer; Section 3.5 presents the discussion of the influence of fracture connectivity, fracture apertures and matrix permeability in the bedrock, as well as the influence of the fractured bedrock on the contaminant transport and detection in the overlying unconfined porous aquifer; Section 3.6 presents the discussion of 41  the influence of detection threshold; Section 3.7 presents the discussion of the influence of source zone location; Section 3.8 summarizes how these factors interact to influence the groundwater monitoring.  3.2 Statistical properties of all the GW scenarios  Given the distribution properties for fracture orientation, length and aperture in Table 2.2, 200 fracture networks are generated. Each realization reflects a relatively sparse network of open fractures in the bedrock. The average and the coefficient of variation of the number of fractures in each fracture set are given in Table 3.1 for each of the scenarios. All the scenarios use the same set of fracture network realizations as those in the Base Case except for Scenario GW-2a and Scenario GW-2b where the volumetric densities of the horizontal and vertical fracture sets are reduced. With the generated fracture networks and the given matrix permeability, the saturated bulk hydraulic conductivities of the fractured bedrock are calculated for all scenarios using the method described in Section 2.6.1. The statistical results are given in Table 3.2. Compared with the Base Case, the average bulk hydraulic conductivity of the bedrock decreases when the fracture densities (Scenario GW-2a and Scenario GW-2b) and the fracture apertures (Scenario GW-3a and Scenario GW-3b) are reduced, but the variability (the coefficient of variation) is increased in these cases. Reduction of the mean of the fracture apertures results in a significant decrease in the average of the bedrock permeability, but a dramatic increase in the variability (coefficient of variation of as high as 50%) due to the fact that the lower and upper bounds of the distribution of the apertures are unchanged, and the locations of those fractures with apertures close to the upper bound vary from realization to 42  realization. The results also show that increase of the rock matrix permeability by two orders of magnitude (Scenario GW-4) only leads to a small increase in the average bulk hydraulic conductivity of the bedrock, but a small decrease of the variability. Among the fracture properties, fracture aperture, especially the horizontal fracture aperture, has the most significant impact on the bulk hydraulic conductivity of the bedrock. The averages of Kx and Ky are over 2 orders of magnitude smaller when the mean aperture of the horizontal fracture set is reduced by a factor of 10 in Scenario GW-3a. The changes of the hydraulic conductivities of the bedrock will inevitably cause impact on groundwater flow and solute transport in the catchment system, and also the detection probabilities of the performance monitoring network. Table 3.3 gives the statistics of the time that the contaminant is detected at the compliance boundary (detection time) in each of the scenarios, and also the location where the contaminant is detected at the compliance boundary (detection location). The average of the detection time at the compliance boundary for the Base Case is 58.3 years, and detection at the compliance boundary occurs only in fractured bedrock in this case. Recall that detection time is defined in terms of the first arrival of the contaminant at the compliance boundary with a concentration above the detection threshold. Significant changes in the detection time and location can be observed in Table 3.3 for other scenarios which investigate different key parameters. These results are discussed in the following sections of this chapter.  43  3.3 Groundwater flow, transport and detection probabilities in the Base Case  The results of the groundwater flow simulation in the first realization of the Base Case are shown in Figure 3.1 and Figure 3.2. Plot A and B of Figure 3.1 show the hydraulic head contours in the 3D domain and in the 2D vertical cross-section at Y=80m (crossing the source zone), indicating that the fracture network creates a quite irregular groundwater flow field in the fractured bedrock at the local scale. The hydraulic head contours in the unconfined aquifer are also non-uniformly distributed in the dominant flow direction. Plot C and D of Figure 3.1 show the saturation profiles in the 3D domain and in the 2D vertical cross-sections at X=60m and Y=80m (crossing the source zone), indicating that the bedrock is fully saturated, and the unconfined aquifer is partially saturated near the ground surface. Figure 3.2 shows the saturation profile across the full depth of the system at the location of the source zone centre (at X=60m and Y=70m). The saturation at this location is 44% on the top, 48% at the elevation 123.5m, 96% at the elevation 122.6m, and 100% under the elevation 121.7m. The results indicate that the source zone contamination (between elevation 121.5m and 123.5m) is situated on the water table and extends into the partially saturated zone, which means solute transport will occur in both the saturated and unsaturated zones. The results from simulations of 200 realizations in the Base Case show that the average inflow rate at the upstream boundary is 12.81m3/d, and the outflow rate at the downstream boundary is 12.86m3/d. The average recharge rate across the water table is approximately 0.03m3/d. With the specified flow boundary conditions, there is no significant change in the position of the water table from realization to realization, but there are significant changes in the local flow field in the bedrock due to the variation of the fracture networks.  44  Figure 3.3 shows the plume (concentrations above the detection threshold) in the first realization of the Base Case at the time that the contaminant is detected at the compliance boundary. It is a representation of the snapshot of the plume at the time where the threshold of 0.1% of the source zone concentration is reached at the compliance boundary. Plot A shows the 3D envelop around the perimeter of the plume, and Plot B shows the plume in the vertical cross-section at Y=80m where the source zone is located. It should be pointed out that these plots are shown with the vertical exaggeration 10× (same as the plots in all other scenarios), and only the concentrations above the detection threshold. Solute mass originating from the source zone is present outside this envelop, but at concentrations below the detection limit. The results indicate that the contaminant transport occurs in both the unsaturated and saturated zones of the unconfined aquifer, and also in the fractured bedrock (the blobs shown in the cross-section representing the plume in the vertical fractures parallel to the XZ plane). It is clear that the concentrations are lower than the detection threshold in most part of the rock matrix. The contaminant first reaches the compliance boundary in the fractured bedrock with the concentration(s) above the detection threshold. Because of the fracture network, the geometry of the plume is extremely irregular in the bedrock, and also irregular in the unconfined aquifer. Figure 3.4 shows the detection probabilities of the performance groundwater monitoring network in the system at the time the contaminant is detected at the compliance boundary. They are computed separately for the unconfined aquifer (Plot A) and the bedrock (Plot B). The detection probabilities show how likely the contaminant will be detected at this time at a potential performance monitoring well location in each unit. In other words, they indicate where a monitoring well should be located when a monitoring network is designed to detect the contaminant before it reaches the compliance boundary. As shown in the plots, detection 45  probabilities of over 90% occur in those wells located downstream of the source zone between Y=50m and Y=90m in both the unconfined aquifer and the fractured bedrock; the probabilities of detecting the contaminant outside this area are lower. The area encompassed by the detection probability contours is about 30m wider in the bedrock than that in the unconfined aquifer, indicating that because of the fracture networks, the contaminant is likely to have more lateral spreading and to be detected in a wider area in the bedrock. In addition, it can be noticed from Figure 3.4 that detection probabilities near the source zone are higher in the unconfined aquifer and lower in the bedrock, because it takes a certain distance downstream from the source for plume to enter the bedrock (about 50m-80m in this case). Moreover, both the plume in the first realization (Figure 3.3) and the detection probability contours (Figure 3.4) indicate that detection of the contaminant at the compliance boundary is likely to occur in the bedrock, rather than in the unconfined aquifer in the Base Case. The saturated hydraulic conductivity (Kuc) of the unconfined aquifer in this case (Table 2.1) is of the same magnitude as the average of the bulk hydraulic conductivity (Kx, Ky, Kz) of the bedrock (Table 3.2). The ratio Kuc/Kx, Kuc/Ky and Kuc/Kz is approximately equal to 1.0, 0.9 and 3.6, respectively. In addition, results show that the average of the detection time at the compliance boundary is 58.3 years in the Base Case, and the variability of the detection time at the compliance boundary (a coefficient of variation of 24.0%, Table 3.3) is relatively large compared to that of the bulk hydraulic conductivity of the bedrock (a coefficient of variation of 5.1% for Kx, 5.2% for Ky, 4.2% for Kz, Table 3.2). No correlation is observed between the detection time at the compliance boundary and the bulk hydraulic conductivity (Kx, Ky and Kz) of the bedrock in the 200 realizations considered (see Figure 3.5). The results demonstrate that change of the fracture networks from realization to realization may not 46  cause a significant change in the permeability of the bedrock, but have a significant impact on contaminant transport and detection. The results indicate that the fracture networks in the bedrock play the key role in controlling the detection of the contaminant at the compliance boundary in this case. In addition, it should emphasized again that although the surface topography used in the groundwater monitoring analysis is in symmetric V-shape, no interactions between groundwater and surface water is considered. The simulations with a uniform topography in the Y-direction also showed the topography has no effects on the detection probabilities in the groundwater system.  3.4 Influence of hydraulic conductivity of unconfined porous aquifer  The influence of the hydraulic conductivity of the unconfined porous aquifer on plume detection is investigated in Scenario GW-1 using the same fracture networks in the bedrock as those in the Base Case, but with a hydraulic conductivity one order of magnitude higher in the unconfined aquifer. The saturated hydraulic conductivity (Kuc) of the unconfined aquifer is increased to 2.5×10-5m/s. As a result, the ratio of the hydraulic conductivity of the unconfined aquifer to that of the bedrock Kuc/Kx, Kuc/Ky and Kuc/Kz is approximately 10, 9 and 36, respectively. Figure 3.6 shows the plume (concentrations above the detection threshold) in the first realization at the time the contaminant is detected at the compliance boundary, in which Plot A shows the 3D envelop of the plume surface, and Plot B shows the plume in the crosssection at Y=80m where the source zone is located. Figure 3.7 shows the detection probabilities at the performance monitoring well locations in the unconfined aquifer (Plot A) 47  and in the bedrock (Plot B). Both the plume geometry and the detection probability contours indicate that the contaminant is first detected in the unconfined aquifer at the compliance boundary. In addition, as shown in Figure 3.7, the detection probabilities of the performance monitoring network measured in terms of the area encompassed by the detection contours is substantially higher in the unconfined aquifer than the bedrock. The results indicate that contaminant transport is dominantly controlled by the unconfined aquifer with higher permeability, and that the detection of the contaminant occurs mainly in this aquifer at the time the contaminant reaches the compliance boundary. Consequently, the average of detection time at the compliance boundary is significantly reduced to 10.3 years due to the increase of permeability of the unconfined aquifer in this case, and it is much shorter than the average of detection time 58.3 years in the Base Case (Table 3.3). The variation in the detection time is also reduced significantly (a coefficient of variation of 1.9%, Table 3.3), in comparison to that in the Base Case (a coefficient of variation of 24.0%). Overall, the results demonstrate that when the permeability of the unconfined aquifer is one order of magnitude higher than that of the bedrock, detection of the contamination is more likely to occur in the unconfined aquifer than in the bedrock, and the time for detection of the contaminant at the compliance boundary is shorter. The results imply that in a case where the permeability of the unconfined aquifer is high in comparison to the permeability of the bedrock, monitoring the contamination in the unconfined aquifer may be sufficient to detect the contaminant before it reaches the compliance boundary, and monitoring in the bedrock may not be needed in such a case. Finally, it was observed that in general, when ratio of the hydraulic conductivity of the unconfined aquifer to that of the bedrock (Kuc/Kx) is equal and greater than 10, detection at 48  the compliance boundary occurs in the unconfined aquifer. When the ratio is less than 1, detection at the compliance boundary occurs in the bedrock. When the ratio is between 1 and 10, detection at the compliance boundary can be either in the unconfined aquifer or in the bedrock, or in both.  3.5 Influence of fractured bedrock  3.5.1 Influence of fractured bedrock on detection in unconfined porous aquifer  A comparison of detection probabilities in Figure 3.4 and 3.7 suggests there is a greater degree of transverse migration in the unconfined aquifer in the Base Case than in the case of Scenario GW-1 where the fracture networks are identical but the unconfined aquifer is assigned a value of hydraulic conductivity one order of magnitude higher. The area encompassed by the detection contours in the unconfined aquifer is about 20m wider in the Base Case than that in Scenario GW-1. To examine this issue further, simulations were carried out for the first 10 realizations of Base Case by replacing the fractured bedrock with a hydraulically equivalent homogeneous and anisotropic porous medium, which has the corresponding Kx, Ky and Kz of fractured bedrock in each of the realizations. Figure 3.8 shows the outer boundaries of the plume in the bedrock and the overlying unconfined porous aquifer (defined by the detection threshold of 0.1% source concentration), at the time the contaminant is detected at the compliance boundary, for the first realization (Plot A) and the tenth realization (Plot B) of the Base Case. Figure 3.9 presents a comparison of the contaminant plume in the unconfined porous aquifer for the first realization (Plot A)  49  and the tenth realization (Plot B) of the Base Case, with the case where the bedrock is replaced by the equivalent porous medium. Figure 3.8 illustrates that the fracture networks in the bedrock produce the irregular transport patterns in both the bedrock and the overlying unconfined porous aquifer. The transport pattern in the unconfined porous aquifer is not the symmetric ellipse shape commonly predicted in homogeneous and isotropic porous media, but has a tendency to follow the transport pattern in the bedrock instead. Figure 3.9 demonstrates that the particular structure of the fracture network has an influence on the plume geometry in the overlying unconfined aquifer, tending to increase transverse spreading in comparison to the case where the fracture bedrock is replaced with the equivalent porous medium. This linkage will also act to influence the detection probabilities in the unconfined aquifer. Examination of the other eight realizations also shows such a phenomenon. This phenomenon can be explained by a comparison of the hydraulic head contours in Figure 3.10 and 3.11 for the first realization of the Base Case. Figure 3.10 shows that the fracture network creates an irregular groundwater flow field in the bedrock, and a nonuniform flow field in the overlying homogeneous and isotropic porous aquifer (also see Plot A and B in Figure 3.1). As a result, more transverse migration occurs, generating a wider and irregularly-shaped plume in the unconfined porous aquifer. A similar behaviour was also observed by Therrien and Sudicky (1996) in simulation of solute transport in a porous aquifer under a sparsely-fractured aquitard. In addition, it can be noticed that the head contours in the unconfined aquifer tend to resemble the patterns of head contours in the bedrock. This explains why the plume in the unconfined porous aquifer has a tendency to follow the transport pattern in the bedrock. By contrast, none of the effects above exists in the flow field for the overlying porous aquifer when the bedrock is replaced with an equivalent 50  homogeneous and anisotropic porous medium, as shown in Figure 3.11. A similar comparison can also be made with the results from simulation of the first realization of Scenario GW-1 where the hydraulic conductivity of the unconfined porous aquifer is one order of magnitude higher than that in the Base Case. This comparison indicated that the influence of the fracture network on flow and transport in the overlying unconfined aquifer becomes less when the permeability of the unconfined aquifer is higher.  3.5.2 Influence of fracture network connectivity  The influence of fracture network connectivity on detection probabilities of the groundwater monitoring network is investigated by reducing by half the volumetric density of the horizontal fracture set (Scenario GW-2a) and the vertical fracture sets (Scenario GW2b). The scan-line density (Table 2.2) and the average number of fractures (Table 3.1) are both reduced by almost one half for the horizontal fracture set and the vertical fracture sets, respectively, in comparison to the Base Case. The fractures that remain in Scenario GW-2a and GW-2b are the same as in the Base Case. As described in Berkowitz et al. (2000) and Wealthall et al. (2001), fracture network connectivity is a function of fracture density, length, orientation, as well as spatial correlation among fractures. Since the three sets of fractures in the bedrock are assumed to be spatially uncorrelated and independent from each other, with the fixed fracture lengths and orientations, the reduction of fracture densities will inevitably reduce connectivity of the fracture networks and also the permeability of the bedrock. As shown in Table 3.2, in comparison to the Base Case, the bulk hydraulic conductivity of the bedrock is reduced approximately by a factor of 3 for Kx and Ky, and 1.4 for Kz on 51  average in Scenario GW-2a, and a factor of 2 for Kx, Ky and Kz on average in Scenario GW2b, respectively. The variability in bulk hydraulic conductivity of the bedrock, however, is increased due to the reduction of the fracture connectivity. The coefficient of variation is 10.5%, 9.2% and 5.7% for Kx, Ky and Kz in Scenario GW-2a; 7.0%, 6.4% and 6.7% for Kx, Ky and Kz in Scenario GW-2b, higher than those in the Base Case. As a result, the ratios of the hydraulic conductivities of the unconfined aquifer and the bedrock (Kuc/Kx) in these scenarios are between 1 and 10. As the fracture densities and fracture connectivity are reduced, the tortuosity of the connected fracture networks will be increased, causing more tortuous migration of the contaminant in the bedrock. Consequently, as shown in Table 3.3, the average of the detection time at the compliance boundary is increased (by 23.2% in Scenario GW-2a and 18.0% in Scenario GW-2b), in comparison to the Base Case. The variability of detection time (the coefficient of variation) is also slightly increased (by 1.4% in Scenario GW-2a and 3.0% in Scenario GW-2b). The lower fracture densities and connectivity result in longer detection time on average at the compliance boundary, with larger variability. The plume (concentrations above the detection threshold) at the time that the contaminant is detected at the compliance boundary is shown in Figure 3.12 for the first realization in Scenario GW-2a and in Figure 3.14 for the 1st realization in Scenario GW-2b, respectively. The detection probabilities of the performance monitoring network at the time that the contaminant is detected at the compliance boundary are shown in Figure 3.13 for Scenario GW-2a and in Figure 3.15 for Scenario GW-2b, respectively. The results show that in comparison to the Base Case, the reduction of the fracture densities and the permeability of the bedrock has a significant effect on contaminant transport and detection in the bedrock and also in the overlying unconfined aquifer. 52  The reduction of fracture densities, especially for the set of horizontal fractures, results in more lateral contaminant transport in the bedrock (Plot A of Figure 3.12) but lower continuity of the plume as defined by the detection threshold (Plot B of Figure 3.12), indicating that a lower fracture connectivity will cause more uncertainty in detection of the contaminant in the bedrock. The reduction of the permeability of the bedrock with the lower fracture densities, especially for the sets of the vertical fractures, results in greater contaminant transport in the overlying unconfined aquifer (the plots of Figure 3.14), indicating that there will be greater probability in detection of the contaminant in the overlying unconfined aquifer when the fracture connectivity in the bedrock is lower. The plots in Figure 3.13 and 3.15 show that the area encompassed by the detection contours in the bedrock (Plot B) for the two cases is about 20m wider than in the Base Case. This also suggests that the contaminant is likely to have more lateral spreading and to be detected in a wider area in the bedrock, due to the increase in detection time at the compliance boundary in these cases with lower fracture densities. In addition, the plots show that due to reduction of the bedrock permeability, the detection probabilities in the overlying unconfined aquifer (Plot A) are slightly increased in both cases. The detection probabilities (Figure 3.13 and 3.15) together with the plumes (Figure 3.12 and 3.14) indicate that detection of the contaminant at the compliance boundary can occur in both the bedrock and the unconfined aquifer in both scenarios, and monitoring the contamination in both units is needed to maximize the likelihood of detecting the contaminant before it reaches the compliance boundary. The results reinforce the conclusion that detection at the compliance boundary can occur in either the unconfined aquifer or the bedrock or both, when the ratio of hydraulic conductivities of the unconfined aquifer to the bedrock (Kuc/Kx) is between 1 and 10. Finally, it should be pointed out again that the 53  probability of detection calculated with the assumption of perfect vertical resolution yields the maximum detection probability in the sparsely-fractured bedrock.  3.5.3 Influence of fracture apertures  The influence of fracture apertures on plume detection is investigated with scenarios where the mean aperture is reduced by a factor of 10 for the set of horizontal fractures (Scenario GW-3a) and for the sets of vertical fractures (Scenario GW-3b), but the lower and upper bounds of the aperture distribution keep unchanged. The reduction of the mean aperture will lower the permeability of the bedrock, and influence the contrast in the bulk hydraulic conductivity of the bedrock relative to the unconfined aquifer. The results in Table 3.2 show that in comparison to the Base Case, the average bulk hydraulic conductivity of the bedrock is reduced by over two orders of magnitude for Kx and Ky, and by a factor of 3.5 for Kz in Scenario GW-3a, due to the reduction of the mean of the horizontal fracture apertures. The average bulk hydraulic conductivity of the bedrock is reduced by about one order of magnitude for Kx and Ky, and a factor of 35 for Kz in Scenario GW-3b, due to the reduction of the mean of the vertical fracture apertures. However, the coefficient of variation of the bulk hydraulic conductivity of the bedrock is increased to 23.2%, 20.0% and 10.3% for Kx, Ky and Kz in Scenario GW-3a, and to 15.0%, 12.6% and 38.5% for Kx, Ky and Kz in Scenario GW-3b, respectively. This is due to the fact that the lower and upper bounds of the aperture distribution are unchanged, and the locations of those fractures with apertures close to the upper bound vary from realization to realization. The results demonstrate that the permeability of the bedrock is mainly controlled by the apertures especially for the set of the horizontal fractures. Although the fracture networks are 54  the same, the reduction of fracture apertures results in a dramatic decrease in the average of the bedrock permeability in both cases. Consequently, the ratio of the hydraulic conductivity of the unconfined aquifer to that of the bedrock (Kuc/Kx) is increased to be greater than 10 in both scenarios. Figure 3.16 and 3.18 show the plume (concentrations above the detection threshold) in three-dimensional envelop (Plot A) and in cross-section (Plot B) at the time the contaminant is detected at the compliance boundary in the first realization for Scenario GW-3a and Scenario GW-3b, respectively. The geometries of the contaminant plumes indicate that apertures for both the horizontal fracture set and the vertical fracture sets are important in controlling the contaminant transport patterns in the bedrock. The plots in Figure 3.16 show that the horizontal fracture apertures control how far the contaminant will migrate horizontally in the bedrock, and determine if first detection will occur in the bedrock at the compliance boundary. The plots in Figure 3.18 show that the vertical fracture apertures control how deep the contaminant will migrate vertically from the unconfined aquifer into the bedrock, and also have a significant impact on detection of the contaminant at the compliance boundary. The results indicate that due to the reduction of the fracture apertures, the contaminant transport is dominantly in the unconfined aquifer in both cases, and detection of the contaminant at the compliance boundary also occurs in this aquifer. As shown in Table 3.3, the reduction of the fracture apertures results in a significant decrease in the variability of detection time at the compliance boundary (a coefficient of variation of 1.6% for Scenario GW-3a, 1.8% for Scenario GW-3b), in comparison to the coefficient of variation of 24.0% for the Base Case. The smaller variation observed in the detection time is similar to that in Scenario GW-1 (with a more permeable unconfined aquifer). The results again reinforce the point that when the ratio of the hydraulic 55  conductivity of the unconfined aquifer to that of the bedrock (Kuc/Kx) exceeds 10, the unconfined aquifer plays the dominant role in contaminant transport and detection in the groundwater system. Furthermore, the results in Table 3.3 show that the average of the detection time at the compliance boundary is significantly increased (109.7 years for Scenario GW-3a, 107.1 years for Scenario GW-3b), in comparison to both the Base Case (the detection time 58.3 years in average) and Scenario GW-1 (the detection time 10.3 years in average). This increase reflects the longer travel time in the unconfined aquifer, relative to the bedrock pathway in the Base Case. The results demonstrate again that the permeability of the unconfined aquifer is a key factor in determining the time for detection of the contaminant at the compliance boundary. Figure 3.17 and 3.19 show the detection probabilities in the unconfined aquifer (Plot A) and the bedrock (Plot B) at the time that the contaminant is detected at the compliance boundary for Scenario GW-3a and Scenario GW-3b, respectively. Because of the reduction of the fracture apertures, the detection of the contaminant at the compliance boundary occurs in the unconfined aquifer, not in the bedrock for both scenarios. The detection probabilities and the areas encompassed by the probability contours in the unconfined aquifer in these scenarios are similar to those in Scenario GW-1. The results suggest that in the case where the unconfined aquifer plays a dominant role in controlling contaminant transport when the ratio of the hydraulic conductivity of the unconfined aquifer to that of the bedrock (Kuc/Kx) exceeds 10, the fractured bedrock becomes much less significant in influencing detection at the compliance boundary as well as the detection probabilities in the overlying unconfined porous aquifer. In such cases, the necessity in monitoring for contamination in the bedrock will be lower, and monitoring the contamination in the unconfined aquifer may be sufficient to detect the contaminant before it reaches the compliance boundary. 56  3.5.4 Influence of rock matrix permeability  The influence of rock matrix permeability on plume detection is investigated in Scenario GW-4 using the same fracture networks but with the matrix hydraulic conductivity (Km) two orders of magnitude higher than in the Base Case. As shown in Table 3.2, in comparison to the Base Case, the average of the bulk hydraulic conductivity of the bedrock is increased by 19.2% for Kx, 18.5% for Ky, and 14.3% for Kz, due to the increase of the matrix permeability. The variability (the coefficient of variation) is slightly reduced by 0.9% for Kx, 1.0% for Ky, and 0.4% for Kz. The ratio of the hydraulic conductivity of the unconfined aquifer to that of the bedrock is less than 1 (Kuc/Kx= 0.8) in this scenario. Figure 3.20 shows the plume (concentrations above the detection threshold) in the first realization of Scenario GW-4. The 3D envelop of the plume geometry (Plot A) has no significant difference from that in the Base Case. However, the cross-section (Plot B) shows that the plume has higher continuity in the bedrock, indicating that the probability in detection of the contaminant in the bedrock will increase when the permeability of the rock matrix is higher. Figure 3.21 shows the detection probabilities of the monitoring network in the unconfined aquifer (Plot A) and the bedrock (Plot B) in this scenario. Similar to the Base Case, the contaminant is first detected at the compliance boundary in the bedrock. The results support the point that when the ratio of hydraulic conductivity of the unconfined aquifer to that of the bedrock is less than 1, detection at the compliance boundary occurs in the bedrock. The results also demonstrate that for the case studied, the rock matrix permeability has a relatively small influence on the plume detection in comparison with the fracture 57  permeability. In other words, the calculation of detection is not sensitive to a two order of magnitude increase in matrix permeability in this case. In addition, as shown in Table 3.3, the average of detection time at the compliance boundary in this scenario is increased by 7.5%, but the coefficient of variation in detection time is reduced by 11.7%, in comparison to the Base Case. The results indicate that in the case where the rock matrix permeability is higher, more contaminant is transferred by advection and diffusion from the fracture walls into the matrix blocks, leading to a later arriving of the contaminant at the compliance boundary with a concentration above the detection threshold. As more solute mass enters into the rock matrix, the probability of detecting the contaminant in the bedrock will increase.  3.6 Influence of contaminant detection threshold  The influence of the contaminant detection threshold on groundwater monitoring is investigated in Scenario GW-5a and Scenario GW-5b using different detection threshold concentrations applied both at the compliance boundary and in the performance monitoring wells. The detection threshold is 1%, 0.1% and 0.01% of the source zone concentration in Scenario GW-5a, in Base Case, and in Scenario GW-5b, respectively. The fracture network realizations used in the simulations of these scenarios are identical to those in the Base Case. Figure 3.22 and 3.24 show the 3D geometry (Plot A) and the 2D cross-section (Plot B) of the plume (concentrations above the detection threshold) in the first realization at the time the contaminant is detected at the compliance boundary in Scenario GW-5a (with the higher detection threshold of 1%) and in Scenario GW-5b (with the lower detection threshold of 0.01%), respectively. The comparison of the plots in these scenarios with those in the Base 58  Case demonstrates that the plume geometry in the unconfined aquifer and the bedrock is very sensitive to the detection threshold. In the case with the detection threshold equal to 1% of the source concentration (Scenario GW-5a), the plume in the unconfined aquifer appears much larger than that in the bedrock, and the contaminant is first detected at the compliance boundary in the unconfined aquifer (the plots in Figure 3.22). When the detection threshold is 0.1% of the source concentration (the Base Case, the plots in Figure 3.3) or 0.01% of the source concentration (Scenario GW-5b, the plots in Figure 3.24), the size of the plume (defined by the detection threshold) in the bedrock is much larger than that in the unconfined aquifer, and the contaminant is first detected at the compliance boundary in the bedrock. The results in Table 3.3 show that for the case with detection thresholds of 1%, 0.1% and 0.01% of source zone concentrations, the average of detection time at the compliance boundary is 115.7 years, 58.3 years and 35.2 years, respectively, and the coefficient of variation is 8.8%, 24.0% and 30.1%, respectively. The higher the detection threshold used, the longer would be the expected detection time at the compliance boundary, and the smaller would be the variability in the detection time at the compliance boundary. Figure 3.23 and 3.25 show the detection probabilities of the monitoring wells in the unconfined aquifer (Plot A) and in the bedrock (Plot B) for Scenario GW-5a and Scenario GW-5b, respectively. The results show that the areas encompassed by the detection probabilities contours in both the unconfined aquifer and the bedrock are closely related to the detection thresholds. The results indicate that a higher detection threshold leads to lower probabilities in detection of the contaminant in the bedrock but higher probabilities in detection of the contaminant in the unconfined aquifer.  59  3.7 Influence of source zone location  The impact of the source zone location on groundwater monitoring is investigated in Scenario GW-6, which is designed with all the same parameters used in Scenario GW-1, but with a different source zone location in the Z axis. In Scenario GW-6, the source zone is at Z=115m to Z=117m, located at the sediment – bedrock interface. In Scenario GW-1, the source was at Z=121.5m to Z=123.5m, located on the water table in the unconfined aquifer. It is already known from Section 3.4 that when the source zone is located on the water table in Scenario GW-1 (with the hydraulic conductivity of 2.5×10-5m/s in the unconfined aquifer, one order of magnitude higher than the bulk hydraulic conductivity of the bedrock), the contaminant plume (concentrations above the detection threshold) is predominantly in the unconfined aquifer (the plots in Figure 3.6). Detection of the contaminant at the compliance boundary occurs in the unconfined aquifer, and the detection probabilities in this aquifer are much higher than those in the bedrock at the time the contaminant is detected at the compliance boundary (the plots in Figure 3.7). However, the plume (concentrations above the detection threshold) in the first realization (Figure 3.26) and the detection probability contours (Figure 3.27) in Scenario GW-6 illustrate that when the source zone is located at the sediment – bedrock interface, the contaminant is first detected at the compliance boundary in the bedrock, rather than in the unconfined aquifer, despite the fact that the ratio of hydraulic conductivities of the unconfined aquifer and the bedrock is over 10 in this case. The size of the plume and the area encompassed by the detection probability contours in the bedrock are much larger than those in the unconfined aquifer at the time the contaminant is detected at the compliance boundary. That is opposite to what was observed in Scenario GW-1. The results confirm the expected result 60  that the closer the source zone is to the bedrock, the more likely it is that the contaminant will be detected in the bedrock. Moreover, as shown in Table 3.3, when the source zone is located at the sediment – bedrock interface, the variability of the detection time at the compliance boundary is significantly increased (a coefficient of variation of 12.3%) in this case, in comparison to Scenario GW-1 (a coefficient of variation of 1.9%), despite that the increase in the average detection time at the compliance boundary is small.  3.8 Summary  Using the methodology presented in Chapter 2, the probabilities of detection of the contaminant in performance wells of a groundwater monitoring network prior to detection at the compliance boundary were evaluated in a catchment where an unconfined aquifer overlies sparsely fractured bedrock. Factors influencing plume detection in the groundwater system were discussed based on the results from Monte Carlo simulations of steady-state groundwater flow and solute transport in 10 scenarios. The results demonstrate that the unconfined aquifer and the fractured bedrock should be viewed as an integrated hydrologic system. The probability of plume detection in the system is predominantly influenced by the ratio of the hydraulic conductivity of the unconfined aquifer to that of the fractured bedrock, coupled with the contaminant detection threshold and source zone location. With the source zone on the water table, the unconfined aquifer 8.94m thick at the source location, and the detection threshold of 0.1% of the source zone concentration, it was observed that when the ratio of the hydraulic conductivity of the unconfined aquifer to that of 61  the bedrock (Kuc/Kx) exceeds 10, detection at the compliance boundary occurs in the unconfined aquifer. When this ratio is less than 1, detection at the compliance boundary occurs in the bedrock. When the ratios are between 1 and 10, detection at the compliance boundary can be either in the unconfined aquifer or in the bedrock, or in both. The hydraulic conductivity of the unconfined aquifer has a significant impact on contaminant transport and detection probabilities. When the unconfined aquifer is more permeable (with a hydraulic conductivity one order of magnitude higher than that of the bedrock), monitoring the contamination in the unconfined aquifer may be sufficient to reliably detect the contaminant before it reaches the compliance boundary, provided that the wells in the unconfined aquifer are properly located (downstream of the source zone), and that the heterogeneity in the unconfined aquifer is not a significant factor. Among the fracture properties, fracture network connectivity and fracture aperture have significant impacts on detection probabilities of the groundwater monitoring network, and the rock matrix permeability has a minor effect. Lower fracture network connectivity with smaller densities tends to increase the tortuous migration of contaminant in the bedrock. Consequently, the detection time at the compliance boundary increases (with larger variability), and the contaminant is likely to be detected in a wider area in the bedrock prior to detection at the compliance boundary. The influence of the bedrock is dependent upon the fracture apertures. The horizontal fracture apertures control how far the contaminant will migrate in the bedrock, and determine if detection will occur in the bedrock at the compliance boundary. The vertical fracture apertures control how deep the contaminant will enter from the unconfined aquifer into the bedrock, and also have a significant impact on likelihood to detect the contaminant in the bedrock.  62  The fractured networks also have a significant impact on groundwater flow, transport and detection of the contaminant in the overlying unconfined porous aquifer. The fracture networks produce the irregular flow patterns and plume geometries, and increase transverse spreading in both the units. The transport and detection patterns in the unconfined aquifer tend to follow those in the bedrock. The results suggest the fracture networks not only contribute to the uncertainty in detection of the contamination in the bedrock but also in the unconfined aquifer. The results also imply that in designing a groundwater monitoring network in such a system, the locations of monitoring wells in the unconfined aquifer and the bedrock should be selected with consideration of the fracture networks in the bedrock. The influence of fractured bedrock becomes less significant when the unconfined aquifer has a high permeability. The plume geometry and the detection probabilities in both the unconfined aquifer and the bedrock are sensitive to the detection threshold. The detection time and location at the compliance boundary are also dependent upon the detection threshold. A higher detection threshold leads to lower probabilities in detection of the contaminant in the bedrock, and increases the detection time at the compliance boundary. Finally, the source zone location has a significant influence on the performance of the groundwater monitoring network and the detection location at the compliance boundary. The closer the source zone is to the bedrock, the more likely is that detection of the contaminant will occur in the bedrock. For monitoring near the source zone, probability of detection depends upon downgradient distance from source before the contaminant reaches the bedrock with concentrations above the detection threshold. Overall, the results indicate that design of the groundwater monitoring network requires an understanding of the large-scale hydraulic conductivities of both the unconfined aquifer 63  and the fractured bedrock. It should not be assumed a priori that if the network of open fractures in the bedrock appears sparse, monitoring the bedrock pathways may not be important.  64  Table 3.1 Statistics of the fracture networks in each scenario  Scenario  # of Fractures // XY Plane  # of Fractures // XZ Plane  # of Fractures // YZ Plane  Average  Coefficient of Variation (%)  Average  Coefficient of Variation (%)  Average  Coefficient of Variation (%)  1581  3  1117  3  1135  3  GW-2a  793  5  1117  3  1135  3  GW-2b  1581  3  561  5  569  5  Base Case GW-1 GW-3a GW-3b GW-4 GW-5a GW-5b GW-6  Table 3.2 Bulk hydraulic conductivities of the bedrock in each scenario  Scenario  Kx (×10-6m/s)  Ky (×10-6m/s)  Kz (×10-6m/s)  Average  Coefficient of Variation (%)  Average  Coefficient of Variation (%)  Average  Coefficient of Variation (%)  2.6  5.1  2.7  5.2  0.7  4.2  GW-2a  0.9  10.5  1.0  9.2  0.5  5.7  GW-2b  1.6  7.0  1.7  6.4  0.3  6.7  GW-3a  0.01  23.2  0.02  20.0  0.2  10.3  GW-3b  0.2  15.0  0.3  12.6  0.02  38.5  GW-4  3.1  4.2  3.2  4.2  0.8  3.8  Base Case GW-1 GW-5a GW-5b GW-6  65  Table 3.3 Statistics of detection time and location at the compliance boundary in each scenario Detection Time  Detection Location  Base Case  Average (yrs) 58.3  Coefficient of Variation (%) 24.0  Unconfined Porous Aquifer No  Fractured Bedrock Yes  GW-1  10.3  1.9  Yes  No  GW-2a  71.8  25.4  Yes  Yes  GW-2b  68.8  27.0  Yes  Yes  GW-3a  109.7  1.6  Yes  No  GW-3b  107.1  1.8  Yes  No  GW-4  62.7  21.2  No  Yes  GW-5a  115.7  8.8  Yes  Yes  GW-5b  35.2  30.1  No  Yes  GW-6  11.4  12.3  No  Yes  Scenario  * “Yes”, “No” represents detection of contaminant occurs in the unit or not.  66  A  B  C  D  Figure 3.1 Hydraulic head and saturation profile in first realization of the Base Case. A Head (m) in 3D domain, B head (m) in cross-section at Y=80m, C saturation in 3D domain, and D saturation in cross-sections at X=60m and Y=80m (crossing source zone). Dimension (m).  Saturation vs Depth at the Location of Source Zone Centre (X=60m, Y=70m, Z=100m-123.94m) 125  Z (m)  120 115 110 105 100 0.0  0.2  0.4  0.6  0.8  1.0  Saturation  Figure 3.2 Saturation profile in the depth at the location of source zone centre in first realization of the Base Case  67  A  B  Figure 3.3 Plume in groundwater (concentrations above detection threshold) at the time of detection at the compliance boundary in the first realization of the Base Case. A Three dimensions, B cross-section at Y=80m (source zone). Concentration (g/L), Dimension (m).  68  A  B  Figure 3.4 Detection probabilities of monitoring wells in the Base Case. A Unconfined porous aquifer, B fractured bedrock.  69  A  Detection Time at Compliance Boundary vs Bedrock Hydraulic Conductivity 120  Time (yrs)  100 80 60 40 20 0 2.0E-06  2.2E-06  2.4E-06  2.6E-06  2.8E-06  3.0E-06  2.8E-06  3.0E-06  3.2E-06  Kx (m/s)  B 120  Time (yrs)  100 80 60 40 20 0 2.2E-06  2.4E-06  2.6E-06 Ky (m/s)  C 120  Time (yrs)  100 80 60 40 20 0 6.0E-07  6.5E-07  7.0E-07  7.5E-07  8.0E-07  8.5E-07  9.0E-07  Kz (m/s)  Figure 3.5 Relationship between the bedrock hydraulic conductivity and the detection time at the compliance boundary in 200 realizations of the Base Case. A Detection time vs Kx, B detection time vs Ky, and C detection time vs Kz.  70  A  B  Figure 3.6 Plume in groundwater (concentrations above detection threshold) at the time of detection at the compliance boundary in the first realization of Scenario GW-1 (higher Kuc). A Three dimensions, B cross-section at Y=80m (source zone). Concentration (g/L), Dimension (m).  71  A  B  Figure 3.7 Detection probabilities of monitoring wells in Scenario GW-1 (higher Kuc). A Unconfined porous aquifer, B fractured bedrock.  72  A  Plumes in the Unconfined Porous Aquifer and the Fractured Bedrock (the 1st Realization) 200  Y (m)  160 120  Source Zone  80 40 0 0  50  100  150  200  250 X (m)  300  350  400  450  500  unconfined porous aquifer fractured bedrock  B  Plumes in the Unconfined Porous Aquifer and the Fractured Bedrock (the 10th Realization) 200 Y (m)  160 120  Source Zone  80 40 0 0  50  100  150  200  250 X (m)  300  350  400  450  500  unconfined porous aquifer fractured bedrock  Figure 3.8 Outer boundaries of the plumes (defined by the detection threshold) in the fractured bedrock and the overlying unconfined porous aquifer, for the first & tenth realizations of the Base Case at the time of detection at the compliance boundary. A The first realization, B the tenth realization.  73  A  Plume in the Unconfined Porous Aquifer (the 1st Realization) 200  Y (m)  160 120  Source Zone  80 40 0 0  50  100  150  200  250  300  350  400  450  500  400  450  500  X (m) with underlain fractured bedrock with underlain equivalent porous medium  B  Plume in the Unconfined Porous Aquifer (the 10th Realization) 200 Y (m)  160 120  Source Zone  80 40 0 0  50  100  150  200  250 X (m)  300  350  with underlain fractured bedrock with underlain equivalent porous medium  Figure 3.9 Comparison of the outer boundaries of the plumes (defined by the detection threshold) in the unconfined porous aquifer between the case with the underlain fractured bedrock and the case with the underlain equivalent porous medium of the bedrock, for the first & tenth realizations of the Base Case at the time of detection at the compliance boundary. A The first realization, B the tenth realization.  74  Figure 3.10 Contours of hydraulic head (m) in the unconfined porous aquifer (cross-section at X=120m) and in the underlain fractured bedrock (cross-section at X=103m) for the first realization of the Base Case. Dimension (m).  Figure 3.11 Contours of hydraulic head (m) in the unconfined porous aquifer (cross-section at X=120m) and the underlain equivalent porous medium of the fractured bedrock (crosssection at X=103m) for the first realization of the Base Case. Dimension (m).  75  A  B  Figure 3.12 Plume in groundwater (concentrations above detection threshold) at the time of detection at the compliance boundary in the first realization of Scenario GW-2a (lower horizontal fracture density). A Three dimensions, B cross-section at Y=80m (source zone). Concentration (g/L), Dimension (m).  76  A  B  Figure 3.13 Detection probabilities of monitoring wells in Scenario GW-2a (lower horizontal fracture density). A Unconfined porous aquifer, B fractured bedrock.  77  A  B  Figure 3.14 Plume in groundwater (concentrations above detection threshold) at the time of detection at the compliance boundary in the first realization of Scenario GW-2b (lower vertical fracture densities). A Three dimensions, B cross-section at Y=80m (source zone). Concentration (g/L), Dimension (m).  78  A  B  Figure 3.15 Detection probabilities of monitoring wells in Scenario GW-2b (lower vertical fracture densities). A Unconfined porous aquifer, B fractured bedrock.  79  A  B  Figure 3.16 Plume in groundwater (concentrations above detection threshold) at the time of detection at the compliance boundary in the first realization of in Scenario GW-3a (smaller horizontal fracture apertures). A Three dimensions, B cross-section at Y=80m (source zone). Concentration (g/L), Dimension (m).  80  A  B  Figure 3.17 Detection probabilities of monitoring wells in Scenario GW-3a (smaller horizontal fracture apertures). A Unconfined porous aquifer, B fractured bedrock.  81  A  B  Figure 3.18 Plume in groundwater (concentrations above detection threshold) at the time of detection at the compliance boundary in the first realization of Scenario GW-3b (smaller vertical fracture apertures). A Three dimensions, B cross-section at Y=80m (source zone). Concentration (g/L), Dimension (m).  82  A  B  Figure 3.19 Detection probabilities of monitoring wells in Scenario GW-3b (smaller vertical fracture apertures). A Unconfined porous aquifer, B fractured bedrock.  83  A  B  Figure 3.20 Plume in groundwater (concentrations above detection threshold) at the time of detection at the compliance boundary in the first realization of Scenario GW-4 (higher rock matrix permeability). A Three dimensions, B cross-section at Y=80m (source zone). Concentration (g/L), Dimension (m).  84  A  B  Figure 3.21 Detection probabilities of monitoring wells in Scenario GW-4 (higher rock matrix permeability). A Unconfined porous aquifer, B fractured bedrock.  85  A  B  Figure 3.22 Plume in groundwater (concentrations above detection threshold) at the time of detection at the compliance boundary in the first realization of Scenario GW-5a (higher detection threshold: 1% of source concentration). A Three dimensions, B cross-section at Y=80m (source zone). Concentration (g/L), Dimension (m).  86  A  B  Figure 3.23 Detection probabilities of monitoring wells in Scenario GW-5a (higher detection threshold: 1% of source concentration). A Unconfined porous aquifer, B fractured bedrock.  87  A  B  Figure 3.24 Plume in groundwater (concentrations above detection threshold) at the time of detection at the compliance boundary in the first realization of Scenario GW-5b (lower detection threshold: 0.01% of source concentration). A Three dimensions, B cross-section at Y=80m (source zone). Concentration (g/L), Dimension (m).  88  A  B  Figure 3.25 Detection probabilities of monitoring wells in Scenario GW-5b (lower detection threshold: 0.01% of source concentration). A Unconfined porous aquifer, B fractured bedrock.  89  A  B  Figure 3.26 Plume in groundwater (concentrations above detection threshold) at the time of detection at the compliance boundary in the first realization of Scenario GW-6 (source zone on bedrock). A Three dimensions, B cross-section at Y=80m (source zone). Concentration (g/L), Dimension (m).  90  A  B  Figure 3.27 Detection probabilities of monitoring wells in Scenario GW-6 (source zone on bedrock). A Unconfined porous aquifer, B fractured bedrock.  91  4 METHODOLOGY  TO  EVALUATE  INTEGRATED  SURFACE  WATER AND GROUNDWATER MONITORING NETWORKS  4.1 Introduction  Based on the insight gained into the evaluation of a performance monitoring network in the groundwater system in Chapter 2 and 3, an integrated surface water and groundwater (SW-GW) monitoring network is considered in the next two chapters. Properties of the surface water system and those factors that were identified to have a significant influence on the groundwater monitoring network are analyzed in the context of an integrated SW-GW monitoring network. Aside from the low cost, an advantage of monitoring surface water is that when a stream is fed by groundwater discharge, it is a good, easily accessed integrator for all the complex flow paths in the subsurface. The factors to be investigated include longitudinal and transverse dispersivities in the stream water, the hydraulic conductivity of the unconfined aquifer, the apertures of the horizontal and vertical fractures in the bedrock, the detection threshold, and the location of the source zone of the contamination. This chapter presents the methodologies used to examine the probability of contaminant detection with an integrated SW-GW monitoring network in a headwater stream catchment with an unconfined porous aquifer overlying fractured bedrock (as shown conceptually in Figure 1.1). The evaluation of integrated SW-GW monitoring network involves an assessment of the likelihood for any given surface water monitoring station (observation point) and groundwater performance monitoring well location that the contaminant concentration will exceed a detection threshold prior to the plume reaching a regulatory  92  compliance boundary at the outflow boundary in the surface and subsurface domains. A numerical model is used to simulate coupled SW-GW flow and transport, and Monte Carlo techniques are used to estimate plume detection probabilities at surface water observation points or at groundwater monitoring well locations. In this chapter, Section 4.2 first introduces the types of stream water flow in open channels and the mechanisms affecting solute concentrations in stream water, and then describes the SW-GW system domain to be modeled; Section 4.3 specifies the boundary conditions for simulation of coupled SW-GW flow and transport; Section 4.4 describes the numerical model used for the simulation; Section 4.5 outlines the approach to evaluate the integrated SW-GW monitoring network; Section 4.6 provides the details of the SW-GW scenarios and the model parameters used in the scenarios; Section 4.7 summarizes the methodologies to evaluate the integrated SW-GW monitoring networks.  4.2 Model domain of the SW-GW system  In a shallow and small natural headwater stream, as illustrated conceptually in Figure 4.1, stream flow can be laminar or turbulent depending on flow rates and velocities. Laminar flow in an open channel occurs at low flow rates with small fluid velocities. The flow velocity is usually the highest at or near the stream centreline and diminishes away from the centreline (Figure 4.2A and 4.2B). Turbulent flow in an open channel usually occurs at high flow rates with larger flow velocities and with swirling in the flow, and the highest flow velocity may not occur at or near the stream centre (Figure 4.2C). The type of stream flow in an open channel can be determined with the criterion known as the Reynolds number (Re), which is defined as Re=4ρVd/µ, where V is the mean of stream flow velocities (m/d), d is the 93  average depth of flow (m), ρ is the fresh water density (1000kg/m3), and µ is the dynamic viscosity of water [1.124×10-3kg/(m⋅s)]. In general, the Reynolds number is low for laminar flow (usually less than 500 for open channel flow) and high for turbulent flow (greater than 2000), depending on the factors such as stream channel geometry (Dingman 1994 and Akan, 2006). Stream flow in an open channel can be steady if the flow characteristics such as water depth and velocity profile remain constant with respect to time at a given flow section, or unsteady if the flow characteristics change with time. The stream flow can be uniform if flow characteristics remain the same along the channel, or non-uniform if flow characteristics change along the channel. The slope of the stream water surface in an open channel can be equal to, less than or greater than the slope of the channel bottom. However, in general, the flow in an open channel is non-uniform and turbulent, and it is relatively shallow at the upstream end and deeper at the downstream end, and the depth and velocity vary at different sections along the channel. The mechanisms that affect the concentrations of a non-reactive solute in stream water are illustrated with the diagrams in Figure 4.3 (point source) and Figure 4.4 (continuous source). They mainly include: (1) advection, the contaminant moving downstream in the channel at the mean stream flow velocity; (2) dispersion, the longitudinal and transverse spreading of the contaminant within the stream channel due to variations of flow velocities (longitudinally, laterally and vertically) at different locations (such as faster flow along the stream centreline and slower flow further away from the centreline, faster flow on the surface and slower flow on the stream bottom), and also due to the effect of turbulence in a turbulent flow system; (3) diffusion, migration of the contaminant in the stream due to the concentration gradients at different locations in the channel; (4) dilution, instantaneous lateral 94  and vertical mixing of the solute mass in the stream with clean water from precipitation and overland runoff entering the stream directly on the surface, and also from groundwater discharge through the streambed, as the contaminant moves to the downstream along the channel; (5) exchange, the contaminant in the stream water entering the adjacent porous medium through the streambed, or mixing with the new solute mass entering the lower reach of the stream with groundwater discharge through the streambed due to the incremental increase of flow along the stream channel. Several other mechanisms may also affect the solute concentrations in the stream, such as physical sorption of the solute to the suspended solids and sediments in the channel. Because of surface water and groundwater interactions, the properties influencing flow and solute transport in groundwater such as the streambed permeability will inevitably have an impact on flow and solute transport, as well as detection probabilities of the contamination in surface water. In order to examine how the factors in the surface water and groundwater system will influence integrated SW-GW monitoring, a model domain, as shown in Figure 4.5 (surface domain) and Figure 4.6 (subsurface domain), is designed to represent a small headwater catchment with a small stream on the land surface, as illustrated in Figure 1.1 and Figure 4.1. The two-dimensional surface domain is composed of an overland zone and a stream channel, which is symmetric about the centerline of the domain. The surface gradient in the overland zone is 10m/500m in the X direction (parallel to the stream channel) and 5m/100m in the Y direction (perpendicular to the stream channel). The surface elevation at the upstream boundary (X=0m) varies from 131m at the centre of stream channel to 136m at the edge, and the elevation at the downstream boundary (X=500m) varies from 121m at the centre of stream channel to 126m at the edge. The stream channel is located in the centre of 95  the surface domain with slope gradients of 10m/500m on the streambed (in the X direction) and 0.25m/1m on the stream bank (in the Y direction). The channel is 500m long, 10m wide at top, 8m wide at the flat bottom, and 0.25m deep (also refer to Figure 4.9 and 4.10 in Section 4.5). For simplicity, the channel is represented as a linear feature in the landscape with the constant width along the entire reach. With this assumption, the influence of changes in the channel geometry including the channel width on stream water depth, flow velocity and discharge is neglected. The subsurface domain is 500m long, 200m wide with a V-shaped symmetric surface. The elevation is 100m at the bottom. The system has the same structure as that used for evaluation of the groundwater monitoring network in Section 2.2, except that the slope gradients on the top surface are slightly larger to support overland and stream flow generation. The subsurface is composed of an unconfined homogeneous and isotropic porous aquifer overlying discretely fractured bedrock. The bedrock is 15m thick (from elevation 100m to 115m), and the unconfined aquifer ranges in thickness from 6m to 21m. The fracture networks in the bedrock have the same features as described in Section 2.2 for the groundwater monitoring network. The bedrock is discretized into block elements with the sizes ∆X=∆Y=∆Z=1m, and the unconfined aquifer is discretized into block elements with the sizes ∆X=∆Y=1m, ∆Z=0.6m2.1m to adapt to the variable ground surface. The surface domain is discretized into twodimensional elements for both the overland zone and the stream channel with the sizes ∆X=∆Y=1m, corresponding to the element sizes ∆X=∆Y=1m in the subsurface domain. The surface and subsurface domains are coupled with the dual node approach (see Section 4.4 for a description of this approach).  96  As described in Section 2.7.1, the spatial discretization of the subsurface domain was determined to be sufficient for simulation of groundwater flow and transport in the unconfined porous aquifer and the fractured bedrock. The spatial discretization of the surface domain is sufficient to simulate surface water flow in the small catchment of 0.1km2, as the recharge rate is specified to be uniform on the land surface, and the surficial properties (see Section 4.6.1) of the catchment that control surface water flow is homogeneous and isotropic. In other studies with simulations of surface water flow, di Giammarco et al. (1996) used the element sizes of 100m×100m and 100m×20m to discretize overland (1000m×800m) and stream channel (1000m×20m) in a tilted symmetric V-catchment; VanderKwaak and Loague (2001) used the element sizes of 14m and 7m to discretize overland and stream in a small catchment of 0.1km2; and Sudicky et al. (2008) used the element sizes of 75m and 10m to discretize overland and streams in a watershed of 17km2. Furthermore, Kollet and Maxwell (2006) pointed out that vertical discretization of subsurface domain near the land surface has no significant impact on stream outflow when runoff is produced by excess saturation if recharge rate on the land surface is lower than saturated hydraulic conductivity of the subsurface. In this study, the recharge rate on the land surface is 0.5m/yr (1.6×10-8m/s), which is lower than the saturated hydraulic conductivity of 2.5×10-6m/s of the unconfined aquifer in the Base Case (see Section 4.3.1 and 4.6.1). Using this discretization, the bedrock unit is represented by 15 layers, and the unconfined aquifer is divided vertically into 10 layers. The catchment is discretized into 2.5 million elements with 2,618,226 nodes in the subsurface domain and 100,000 elements with 100,701 nodes in the surface domain. It was a challenge to solve the linear and non-linear equations simultaneously for coupled SW-GW transient flow and solute transport with such a large number of nodes and the available computer resources. 97  4.3 Boundary conditions for coupled SW-GW flow and transport  4.3.1 Boundary conditions for coupled SW-GW flow  The boundary conditions for coupled SW-GW flow simulations are: (1) Specified constant uniform flux of 0.5m/yr on the land surface, representing the net flux rate of the long-term average annual rainfall minus evapotranspiration. (2) Zero-depth gradient of 0.02 in the X direction at the outflow boundary (at X=500m) in the surface domain. (3) Drain nodes at the outflow boundary (at X=500m) throughout the full vertical extent of the subsurface domain. (4) No flow on the bottom (at Z=100m) of the subsurface domain, and also on the side boundaries (at X=0m, Y=0m and Y=200m) of the surface and subsurface domains. The flux specified on the land surface contributes to water flow in both the surface and subsurface domains. The areal recharge or discharge on the top of the subsurface domain is to be determined by the water exchange between the two domains (see more details in Section 4.4). The zero-depth gradient boundary condition is used to simulate water flow out of the surface domain through the outflow boundary, where the water surface slope is specified to be equal to the bed slope. The zero-depth gradient used in this study is 0.02 in the X direction, which means that the water surface slope at the outflow boundary is equal to the slope of the overland zone and the stream channel from the upstream to the downstream.  98  4.3.2 Boundary conditions for coupled SW-GW transport  The boundary condition for coupled SW-GW transport simulations are: (1) Non-reactive contaminant source with continuous constant concentration (1.0g/L) at the location of X=50m to X=70m, Y=60m to Y=80m, Z=121.5m to Z=123.5m in the subsurface. The source zone is between 6.5m and 8.5m above the bedrock, between 7.8m and 9.8m below the ground surface, and 20m left to the stream channel. The source zone is situated in the saturated zone of the unconfined porous aquifer. As mentioned in the previous groundwater monitoring analysis, the size of the source zone is chosen relative to the size of the model domain and also the space of fractures under the source zone in the unconfined porous aquifer. (2) The 3rd type boundary on the land surface, and also at the outflow boundary (at X=500m) in the surface and subsurface domains. (3) Zero mass flux on the bottom (at Z=100m) of the subsurface domain, and also on the side boundaries (at X=0m, Y=0m and Y=200m) of the surface and subsurface domains. (4) Zero initial (background) concentration in the surface and subsurface domains.  4.4 Description of numerical model: HydroGeoSphere  The numerical code HydroGeoSphere is an extension of the model FRAC3DVS with capacity for simulation of fully-integrated surface water (SW) and groundwater (GW) flow and solute transport (Therrien et al., 2005).  99  In HydroGeoSphere, two-dimensional shallow surface water flow is described with the depth-averaged diffusion-wave approximation of the Saint Venant equation. The governing equations for groundwater flow and transport are the same as those for FRAC3DVS (Therrien and Sudicky, 1996; Therrien et al., 2003), three-dimensional variably-saturated subsurface flow is described with Richards equation, and solute transport is governed by the advection-dispersion equation in all domains. Flow and transport in the subsurface is represented as two-dimensional in a fracture, and three-dimensional in a porous medium. The model solves the surface and subsurface flow and solute transport equations simultaneously in an iterative fashion at every time step. The methods for discretization of surface and subsurface flow and transport domains are illustrated in Figure 4.7. The coupled surface and subsurface domains consist of: (1) a single layer of surface nodes, shown by triangles in Figure 4.7, situated on the land surface, (2) layers of subsurface nodes, shown by circles in Figure 4.7, representing the vadose zone and saturated zone. The 2D surface flow grid is draped over the subsurface 3D mesh to maintain the areal topography of the land surface and to ensure the nodes in the surface grid are coincident with those at the top of the subsurface mesh. Discretization of the surface and subsurface flow and transport equations is identical except for the difference in dimensionality. Two approaches can be used in HydroGeoSphere to couple the surface and subsurface flow and transport domains: a common node approach and a dual node approach. Figure 4.8 illustrates the concepts how the two domains are coupled using these two approaches. The common node approach uses a numerical superposition principle whereby the top layer of nodes represents both surface and subsurface domains (the same node numbers for the two domains), with the assumption of the continuity of hydraulic head and solute concentration, 100  and without the need to compute explicitly the water flux and solute exchange between the two domains. The dual node approach uses separate nodes (and different node numbers) to represent the surface domain and the top layer of the subsurface domain, with the assumption that they are separated by a (possibly) thin layer of porous material, and the water flux and solute exchange terms between the two domains are computed explicitly based on the nodal hydraulic head gradients and concentration gradients across this layer. Although there is no difference in the results of flow and transport using the two different methods, the dual node approach is more general because it does not assume continuity of hydraulic head and concentration between the two domains. Therefore, it is adopted to couple the surface and subsurface water flow and solute transport in this study. In HydroGeoSphere, rainfall and evapotranspiration are considered as source and sink terms in the surface water flow equations. Water exchange between surface and subsurface domains is expressed as infiltration (a negative exchange, representing flow from the surface to the subsurface) and exfiltration (a positive exchange, representing flow from the subsurface to the surface). The infiltration is considered as recharge on the top layer of the subsurface domain in the groundwater flow equations. The surface water flow is simulated with the assumptions of depth-averaged flow velocities and a hydrostatic vertical head distribution. The stream flow is considered to occur once the water depth and the flow velocity at the uppermost node(s) in the stream are greater than zero. In the coupled SW-GW transport simulation with HydroGeoSphere, mass in the surface water system comes from the contaminant source zone in the groundwater system. It is determined by solute exchange between the surface and the subsurface. Both advection and dispersion are considered in computation of the solute exchange using the dual node  101  approach. A positive solute exchange represents the contaminant from the groundwater system entering the surface water system. To solve the coupled SW-GW flow and transport equations, both the finite element and finite difference approaches are available in the HydroGeoSphere. However, with the large number of nodes (over 2.7 millions) in the model domain for this study, the finite difference approach was the only option in the code that could be accommodated with the available computer resources. Completion of one realization for the Base Case simulation required approximately 120 hours of CPU time, 4GB of processor memory and 10GB of disk space on the Westgrid computer system in Canada. The computer resources only allows simulation of the coupled SW-GW flow and transport in the given model domain with the limit of three million nodes. More details about the code HydroGeoSphere including the governing equations for flow and solute transport in coupled surface water and groundwater and the numerical approaches for solving the equations can be found in Therrien et al. (2005). Some of the good examples of applying the code to simulate coupled SW-GW flow and solute transport can be found in Jones et al. (2006), Jones et al. (2008), Li et al. (2008) and Sudicky et al. (2008).  4.5 Approach to evaluate integrated SW-GW monitoring network  The integrated SW-GW monitoring network is comprised of surface water (SW) monitoring and groundwater (GW) monitoring. The regulatory compliance boundary for the monitoring is assumed to be located at the outflow boundary (at X=500m) in both the surface and subsurface domains, and the detection threshold specified at the compliance boundary for the surface water monitoring is assumed to be the same as that for the groundwater 102  monitoring. The Monte Carlo approach is used to characterize the uncertainties in the SWGW flow and transport in association with the stochastic fracture networks in the bedrock. Surface water monitoring is carried out at 9 observation points located on the bottom (at the elevation 121m, from Y=96m to Y=104m) of the stream outlet (at X=500m, the compliance boundary) with an interval of 1m (equal to node spacing), as shown in Figure 4.9. Surface water monitoring is also carried out at additional 4 observation points, as shown in Figure 4.10, situated on the bottom of the stream channel at X=100m, 200m, 300m and 400m along the centreline to examine solute concentration distribution and estimate detection probabilities along the stream. As the code HydroGeoSphere simulates two-dimensional areal surface flow only, the concentration in the stream water is depth-averaged. Therefore, the spatial variability of stream water concentration in the vertical direction is neglected at a specific location. In other word, no vertical mixing in the stream is considered at different depths. Groundwater monitoring is carried out with the monitoring network shown in Figure 4.11, consisting of 171 potential locations of performance monitoring wells between the source zone and the compliance boundary. The layout of the monitoring well locations is the same as was presented in Section 2.5. It is assumed that the monitoring wells also have the full vertical resolution in both the unconfined aquifer and the fractured bedrock. The monitoring wells only examine concentrations along the well nodes and do not affect the flow field or influence solute transport. With these assumptions, the uncertainty related to the vertical position of a monitoring zone relative to the location of the plume in each of the two units is ignored. Detection probabilities computed directly by comparison of the simulated concentrations to the detection threshold reflect “the best case” without consideration of dilution effects within the monitoring wells. 103  The factors influencing integrated SW-GW monitoring are evaluated in 8 scenarios (see Section 4.6), by using the steady-state flow field derived from transient SW-GW flow simulations with fixed average annual net recharge rate on the land surface. In this way, contaminant detection is based on long-term monitoring in surface water and groundwater, neglecting the rapid fluctuation of solute concentrations in response to the temporal variations of stream water flow due to individual rainfall events. The coupled SW-GW transport simulations are terminated once contamination is detected anywhere at the compliance boundary, in the surface domain or the subsurface domain, in order to identify where the contaminant is most likely detected in the catchment (in surface water or in groundwater, in the unconfined aquifer or in the bedrock). Detection of contamination at the compliance boundary is considered to occur if the concentration(s) at any of the surface or subsurface nodes on this boundary exceeds the detection threshold, a value defined as the percentage (0.1%) of the source zone concentration. Similarly, detection at a surface water observation point or at a groundwater performance monitoring well location is considered to occur if the concentration(s) in any node at these locations exceeds the detection threshold at the same time the contaminant is detected at the compliance boundary. Detection probabilities are computed for the surface water observation points and the groundwater performance monitoring wells in the unconfined aquifer and the bedrock, as the proportion of realizations in the Monte Carlo simulation where the detection is recorded at that location. Similar to the point discussed in Section 2.5, detection probabilities in the SW-GW monitoring network are not defined in an absolute sense, but rather as the probability of detecting the plume in the surface domain or the subsurface domain at the same time that the contaminant is detected at the compliance boundary. For example, if the surface water 104  observation point “A”, as shown in Figure 4.9, has a detection probability of 0.8, this means there is a 80% probability of detecting the plume in surface water at this point when the contaminant reaches the compliance boundary in either the surface domain or the subsurface domain with a concentration above the detection threshold. It also implies there is a 20% probability of failure in detection of the contamination at this point. The same rule applies to the groundwater monitoring well at the location “B”, as shown in Figure 4.11. If the detection probability at this well location is 0.2 for the bedrock, it means there is a 20% probability of detecting the plume in groundwater (in bedrock) at this location, but a 80% probability of failure in detecting the plume in the bedrock at this location, at the time contaminant is detected at the compliance boundary in either the surface domain or the subsurface domain.  4.6 Scenarios and model parameters  4.6.1 Base Case and the parameters  (1) Parameters for Simulation of Surface Water  The model parameters used in the Base Case for simulation of the surface water component in the coupled SW-GW flow and transport are summarized in Table 4.1. The Manning roughness coefficient is the same in both the X and Y directions, and is a parameter used in solving the governing equation for surface water flow. This parameter is assigned a value of 3.0×10-6 day/m1/3 for both the overland zone and the stream channel in the Base Case, similar to the value (3.5×10-6 day/m1/3) used in the field-scale study of overland flow in 105  a grassy watershed (Abdul, 1985). Test simulations of the coupled SW-GW flow using the first realization of the Base Case in this study showed that the stream water flow velocity and water depth were more sensitive to the roughness in the channel than in the overland zone. The rill storage and obstruction heights represent the amount of water stored in depressions and obstructions on ground surface. Theoretically, no lateral surface water flow will occur unless the rill storage and the obstruction area are filled. As the surface of the model domain in this study is symmetric V-shaped without any depression or obstruction in the microtopography, both the rill storage and obstruction heights are specified to be zero, which is considered as the conservative case for the purpose of monitoring the contamination in surface water. The coupling length is the leakance factor between the surface and subsurface domains when the dual node approach is adopted, and the value 0.1m is given in Table 4.1 by referring to Therrien et al. (2005) in simulation of the field-scale problem of Abdul (1985). Test simulations of the coupled SW-GW flow using the first realization of the Base Case in this study showed the differences in the stream flow velocities and water depths were very small (0.03%) when using the coupling lengths of 1.0m (close to that of 1.35 in Smith and Woolhiser, 1971), 0.1m (from Therrien et al., 2005) and 1×10-4m (from VanderKwaak, 1999), which indicated that the coupling length is not a critical factor in influencing the stream water flow in the small catchment. The longitudinal and transverse dispersivities in Table 4.1 are used to account for the dispersive effect of contaminant transport in surface water due to the spatial variations of flow velocities, and their values are specified by referring to the studies carried out by VanderKwaak (1999), Sudicky et al. (2008) and the test simulation of surface water transport with HydroGeoSphere (personal communication with the code developers). In this study, 106  longitudinal and transverse dispersivities in the stream water are investigated to examine how they influence solute transport and concentration distribution in stream water, and how much impact they have on detection of the contaminant in the stream channel.  (2) Parameters for Simulation of Groundwater  The model parameters used in the Base Case for simulation of the groundwater component of the coupled SW-GW flow and transport in the unconfined homogeneous isotropic porous aquifer and the fractured bedrock are the same as those listed in Table 2.1, Table 2.2, Table 2.3 and Table 2.4 in Section 2.6.1 for the Base Case of groundwater monitoring network evaluation, except that the hydraulic conductivity of the rock matrix is increased by two orders of magnitude to 5.8×10-8m/s for the coupled flow and transport simulations. It has been demonstrated with the results in Section 3.5.4 that increase of the matrix permeability by two orders of magnitude has no significant impact on plume detection for groundwater monitoring. A test simulation of the coupled SW-GW flow using the first realization of the Base Case showed the differences were very small in the results of surface water flow (hydraulic head, stream water depth, stream flow velocity, water exchange between the surface and the subsurface) when using the hydraulic conductivity 5.8×10-8m/s or 5.8×10-10m/s for the rock matrix. However, the computational time required for simulation of the transient SW-GW flow in the case with the hydraulic conductivity of 5.8×10-10m/s for the matrix was about 3-4 times longer and close to the time limit of 240 hours for running a simulation on the Westgrid computer system. The properties of the stochastic fracture networks used for the Monte Carlo simulations of the coupled SW-GW flow and transport in the Base Case are exactly the same as those in 107  the Base Case for evaluation of the groundwater monitoring network (Table 2.2 in Section 2.6.1). Figure 4.12 shows the first realization of the fracture network in three dimensions for the Base Case, which has the same scan-line fracture densities in the cross-sections, as what were shown in Figure 2.6 in Section 2.6.1 for those in the Base Case for the groundwater monitoring network analysis. The bedrock is sparsely-fractured in the Base Case for simulation of coupled SW-GW flow and transport.  (3) Parameters for Control of Flow and Transport Simulations  Using the model domain and the given parameters, the coupled SW-GW transient flow is simulated first to compute the steady-state flow solution in the surface and subsurface domains, followed by the coupled SW-GW solute transport simulations to compute the solute concentrations in surface water and groundwater. The coupled transient SW-GW flow simulation is carried out with the adaptive timestepping procedure, using the head control of 0.1m (the maximum change of hydraulic heads in each timestep) and the timestep control of 90 days (the maximum change of timestep size in each of the iterations) for both the surface and subsurface domains, and the simulation is run for 100 years, which has been found to be sufficiently long for the system to achieve a steady state (see Section 5.3.1 in Chapter 5). The coupled SW-GW solute transport simulation is carried out with adaptive timestepping, using the concentration control of 1.0×10-3g/L (the maximum change of solute concentrations in surface water and groundwater in each timestep). It was found to be sufficient in controlling the size of timestep (smaller than 3 days in the first realization of  108  transport simulation). The transport implicit-explicit time weighting is 0.51, the same as that in Section 2.7.1 for the groundwater monitoring network analysis.  (4) Detection Thresholds for the Integrated SW-GW Monitoring  To compute the detection probabilities of the integrated SW-GW monitoring network, the same detection threshold 1.0×10-3g/L, 0.1% of the source zone concentration (1.0g/L), is assumed for both surface monitoring and groundwater monitoring in the Base Case. Again, this threshold is chosen based on the practical detection capability of many common contaminants relative to the background concentrations. The transport simulation is run until the contaminant reaches the compliance boundary with a concentration above the detection threshold at any node(s) on this boundary in surface water and groundwater. In this way, the likelihood of detecting the contaminant in stream water and groundwater (in the unconfined aquifer and the bedrock) at the time the contaminant is detected at the compliance boundary can be identified.  4.6.2 Other scenarios and the parameters  Table 4.2 outlines 7 scenarios used to examine the key factors influencing the integrated SW-GW monitoring in comparison to the Base Case. These factors are the hydraulic conductivity of the unconfined porous aquifer, the apertures of the horizontal and vertical fractures in the bedrock, the longitudinal and transverse dispersivities in stream water, the detection threshold and the source zone location.  109  All the scenarios in Table 4.2 use the same surface and subsurface model domains and the same parameters as in the Base Case, except for the key parameters to be investigated. Also the same stochastic properties are used to generate the Monte Carlo realizations of fracture networks in the bedrock for the coupled SW-GW flow and transport simulations in all the SW-GW scenarios, except for Scenario SW-GW-2a and 2b, where the horizontal and vertical fracture apertures are increased. The methods for computation of the saturated hydraulic conductivity of a single fracture and that of the bedrock are the same as presented in Section 2.6.1 in Chapter 2.  4.7 Summary  This chapter has described the methodologies for the evaluation of integrated surface water (SW) and groundwater (GW) monitoring network, including the problem geometry and the numerical techniques used for simulation of fully-coupled SW-GW flow and transport in a headwater stream catchment with an unconfined homogeneous and isotropic porous aquifer overlying sparsely fractured bedrock. The techniques for estimation of the detection probabilities of the performance monitoring network in the surface water and groundwater system have also been described. The factors influencing the integrated SW-GW monitoring are evaluated using 8 scenarios, which differ from each other in key parameters to be examined. These factors include the longitudinal and transverse dispersivities in stream water, the hydraulic conductivity of the unconfined aquifer, the fracture apertures in the bedrock, the detection threshold, and the location of the source zone contamination. The results and discussion on how these factors interact to influence the SW-GW monitoring are presented in Chapter 5. 110  The following key points are emphasized: (1) The integrated SW-GW monitoring network focuses on the long-term monitoring of the contamination in surface water and groundwater. It is carried out on the basis of the fully-coupled SW-GW non-reactive contaminant transport simulation with a steady-state flow solution. The physical significance of modeling steady-state flow and transient transport is that the evaluation of the SW-GW monitoring network in this study does not account for the influence of storm event flow causing temporal changes of solute concentrations in stream water. (2) The parameters of the surface water flow system and the properties of the unconfined porous aquifer are treated deterministically. Uncertainty in the contaminant transport in the catchment originates only from the stochastic fracture networks in the bedrock, and it is characterized with Monte Carlo simulations of fracture networks in 8 scenarios. The same set of realizations of the fracture networks are used in all the scenarios in order to examine the influence of the key parameters in the SW-GW monitoring by comparing the results in different scenarios. (3) Detection is defined in terms of the probability of observing a solute concentration above the detection threshold in the surface and subsurface domains at the same time contamination is detected at the compliance boundary. Detection probabilities of the integrated SW-GW monitoring network can be interpreted as the likelihood for any given surface water monitoring station (observation point) and groundwater performance monitoring well location that the contaminant concentration will exceed the detection threshold prior to the plume reaching the regulatory compliance  111  boundary at the outflow boundary at concentrations above the detection threshold in the surface water and groundwater system.  112  Table 4.1 Parameter values for simulation of surface water Parameter Name  Overland  Stream Channel  Manning roughness coefficient (day/m1/3)  3.0×10-6  3.0×10-6  Rill storage height (m)  0.0  0.0  Obstruction height (m)  0.0  0.0  Coupling length (m)  0.1  0.1  Longitudinal dispersivity (m)  1.0  10.0  Transverse dispersivity (m)  0.1  1.0  Table 4.2 Investigated scenarios and features for coupled SW-GW simulations Scenario  Intent  Features and Parameters  Base Case  Base Case  See Section 4.6.1, detection threshold 1.0×10-3g/L  SW-GW-1  Hydraulic Conductivity  Kuc increased by 1 order of magnitude  of Unconfined Aquifer SW-GW-2a  Fracture Aperture  Mean aperture of horizontal fracture set increased by a factor of 1.5 (from 250µm to 375µm)  SW-GW-2b  Fracture Aperture  Mean aperture of vertical fracture sets increased by a factor of 1.5 (from 250µm to 375µm)  SW-GW-3a  Dispersivity in Stream  Longitudinal dispersivity (αl) in stream water reduced by half  SW-GW-3b  Dispersivity in Stream  Transverse dispersivity (αt) in stream water reduced by half  SW-GW-4  Detection Threshold  Reduced by 1 order of magnitude (0.01% of source concentration)  SW-GW-5  Source Location  Contaminant source zone on top of the bedrock  113  Figure 4.1 Illustration of a small headwater catchment with stream flow (photo from the Bertrand Creek Watershed, Langley, BC, Canada, taken by Cindy Starzyk)  Y  V  V  X  Uniform Laminar Flow (A)  Non-Uniform Laminar Flow (B)  Turbulent Flow (C)  Figure 4.2 Types of stream water flow and streamlines in open channel (plan view)  Runoff  Clean Groundwater Discharge  Runoff Precipitation Runoff  Dispersion/Diffusion  Y  Z  Vx Advection  X  Dilution Plume  X  Mixing with Clean Water  Runoff  Clean Groundwater Discharge  (Plane View)  Exchange of Stream Water and Groundwater Runoff  (Cross-Section)  Figure 4.3 Mechanisms affecting solute concentrations in stream water (point source)  114  Runoff  Advection  Y  Vx  New Solute Mass Entering the Stream  Clean Groundwater Discharge  Precipitation  Runoff  Dispersion/Diffusion  Z  Dilution  X  Runoff  Clean Groundwater Discharge  Mixing  Mixing with New Plume  Plume  X Exchange of Stream Water and Groundwater  Runoff  (Plane View)  (Cross-Section)  Figure 4.4 Mechanisms affecting solute concentrations in stream water (continuous source)  Surface Domain  10m 5m  5m  0.25m 1m  8m  1m  95m  Overland Zone  10m  500m  Stream Channel  Figure 4.5 Surface domain (overland zone and stream channel) of the catchment for coupled SW-GW simulations. Dimension (m).  115  Subsurface Domain  Z=136m  Z=131m  Z=136m  Z=126m  Unconfined Porous Aquifer  Z=121m  Z=126m  Fractured Bedrock  Figure 4.6 Subsurface domain of the catchment with an unconfined porous aquifer overlying fractured bedrock for coupled SW-GW simulations. Dimension (m).  Figure 4.7 Illustration of the methods for discretization of surface and subsurface domains (Therrien et al., 2005)  116  Figure 4.8 Illustration of the approaches for coupling of surface and subsurface domains (Therrien et al., 2005)  Rainfall  Overland Flow Y=95m  Y=105m  Overland Flow  10m Z  Stream Flow Y  X  0.25m  A Y=96m  Observation Point  Y=104m  8m  Z=121m  1m  Figure 4.9 Layout of the surface water observation points at the outlet of stream channel  117  10m Z=131.25m  Rainfall  500m  8m Z=131m Z=121.25m  Z  10m  Y X Stream flow Observation Point  8m  1m  Z=121m  Figure 4.10 Layout of the surface water observation points along the centerline of stream channel  B  Location of Monitoring Well  Compliance Boundary  Figure 4.11 Locations of the performance groundwater monitoring wells for the integrated SW-GW monitoring  118  Figure 4.12 Fracture network in three dimensions in the first realization of the Base Case for simulation of coupled SW-GW flow and transport. Aperture (m), Dimension (m).  119  5 FACTORS INFLUENCING INTEGRATED SURFACE WATER AND GROUNDWATER MONITORING  5.1 Introduction  This chapter presents the results of fully-coupled surface water (SW) and groundwater (GW) flow and transport simulations, and detection probabilities at stream water observation points and at groundwater monitoring well locations in the unconfined porous aquifer and the fractured bedrock for each of the 8 scenarios described in Chapter 4. The factors examined include the longitudinal and transverse dispersivities in the stream water, the hydraulic conductivity of the unconfined aquifer, the apertures of the horizontal and vertical fractures, the contaminant detection threshold and the source zone location. Their influence on monitoring is evaluated in terms of the detection location at the compliance boundary, and the detection probabilities in stream water and groundwater at the time contaminant is detected at the compliance boundary. In the scenarios examined, the hydraulic conductivity of the unconfined aquifer has a significant influence on the length of open water in the stream channel in this headwater catchment, which in turn will influence detection probabilities. Section 5.2 presents a summary of the statistical properties of the coupled SW-GW flow and transport simulations, including the bulk hydraulic conductivity of the bedrock; the water depth, total discharge and mean velocity of the stream water flow; the maximum solute concentration in the stream water, and the detection time at the compliance boundary. Section 5.3 presents the flow, transport and detection probabilities in the surface water and groundwater system for the Base Case. Section 5.4 examines the influence of the hydraulic  120  conductivity of unconfined porous aquifer, while Section 5.5 documents the influence of horizontal and vertical fracture apertures in the bedrock. Section 5.6 examines the influence of longitudinal and transverse dispersivities in stream water. A discussion of the influence of contaminant detection threshold is given in Section 5.7, and a discussion of the influence of source zone location is given in Section 5.8. A summary of how these factors interact to influence integrated surface water and groundwater monitoring is given in Section 5.9.  5.2 Summary of statistical properties of the coupled SW-GW scenarios  The statistics of the fracture networks in the Base Case for the coupled surface water and groundwater (SW-GW) simulations are the same as those in the Base Case for the groundwater simulations presented in Chapter 3 (Table 3.1). The fracture networks for the coupled SW-GW simulations in the other scenarios are the same as in the Base Case, except for Scenario SW-GW-2a and 2b, where the mean aperture of the horizontal and vertical fracture sets is increased by a factor of 1.5 in comparison to the Base Case. Table 5.1 gives the statistics of the bulk hydraulic conductivities of the fractured bedrock. Because the same set of fracture networks is used, the values are the same in each scenario except for Scenario SW-GW-2a and 2b with the larger fracture apertures, where the bulk hydraulic conductivities are higher. Table 5.2 presents the characteristics of the surface water flow for the first realization of each scenario, including stream length, water depth, total discharge, the cross-sectional mean of the depth-average flow velocity (Vx) at the stream outlet, and the Reynolds number of the stream flow. Table 5.3 summarizes the statistics of stream water flow among all the realizations in each scenario. Table 5.4 reports the number of Monte Carlo simulations in each scenario, the statistics of the detection time, and the 121  location at which the contaminant is detected at the compliance boundary in the surface and subsurface domains. Also given is the maximum solute concentration in the stream water along the entire reach at the time the contaminant reaches the compliance boundary at a concentration above the detection threshold(s). The cross-sectional mean of the depth-averaged stream flow velocity reported in Table 5.2 and 5.3 is computed with the simulated water depth and total discharge at the hydrograph nodes at the stream outlet, using the formula V=Q/(wd), where V is the mean flow velocity (m/d), Q is the total discharge (m3/d), d is the water depth (m), and w is the stream width (m). Because the average and the variability of the water depths are small in all the scenarios (Table 5.3), the bottom width (8m) of the channel is used to represent the average of the stream width. It is obvious that the depth-averaged stream flow velocity being calculated is dependent upon the channel width, which is constant along the entire reach. As discussed in Chapter 4, all the parameters for flow and transport in the surface water and groundwater system are treated deterministically within each of the scenarios, except for the stochastic fracture networks. Therefore, the observed variability in Table 5.1, 5.3 and 5.4 originates from the fracture networks, which not only control the magnitude and the variation of the bulk hydraulic conductivities of the bedrock, but also cause the variations of flow, solute concentration in stream water and detection time at the compliance boundary between different scenarios and between the individual realizations in each scenario. The main objective of the analysis is to understand how the properties of the fracture bedrock might influence the design of integrated surface water – groundwater (SW-GW) monitoring network. Within each scenario, Table 5.3 indicates the properties of the stream vary little from one realization to the next.  122  The bulk hydraulic conductivities (Kx and Ky) of the bedrock in Table 5.1 are greater than the saturated hydraulic conductivity of the unconfined aquifer 2.5×10-6m/s for all scenarios, hence the ratio of the hydraulic conductivities of the unconfined aquifer to the bedrock (Kuc/Kx) is always less than 1. Reynolds numbers (Re) in Tables 5.2 and 5.3 for all scenarios are close to the upper bound of the laminar flow, with which the threshold is typically 500 for open channel flow (Dingman, 1994 and Akan, 2006). The maximum solute concentration in the stream water in Table 5.4 is about 16 times higher than the assigned detection threshold in the Base Case.  5.3 SW-GW flow, transport and detection probabilities in the Base Case  5.3.1 SW-GW flow and transport in the first realization  (1) Coupled SW-GW flow  Coupled surface water (SW) and groundwater (GW) flow is simulated in a transient mode, but with the specified steady state boundary conditions. The simulation is run until the flow regime in the surface and the subsurface converges to an equilibrium condition. The steady-state solution of the SW and GW flow fields is used for the transport simulations to calculate plume detection probabilities. Figure 5.1 shows the evolution of hydraulic heads in the stream water and in the streambed groundwater at an observation point (at the node at X=100m, Y=100m) in the stream channel. The streambed groundwater refers to the water in the top layer of the nodes in the subsurface domain. The result indicates that a steady state is achieved in about 74 years 123  since the start of the simulation. The actual time at which the equilibrium of the flow fields occurs is not significant to the study. An output time of 100 years is sufficiently long for the flow system to stabilize. The relation between the hydraulic heads in the surface water and the streambed groundwater shown in Figure 5.1 indicates that groundwater discharges into the stream. As shown in Figure 5.2, the stream flow first appears in the channel at this node at the time of about 62 years from the start of the simulation, and reaches steady-state at about 74 years. The total discharge, the water depth and the mean of the depth-averaged flow velocity along the stream channel at the steady state for the first realization are shown in Figure 5.3, 5.4 and 5.5, respectively. The results indicate that the total discharge, water depth and flow velocity increase in a downstream direction, and the stream is relatively shallow with low flow rates along the entire reach of 465m (from X=35m to X=500m) in the small catchment area (0.1km2, 500m×200m). As shown in Table 5.2, the water depth, total discharge and flow velocity are 1.18cm, 107.8m3/d and 1141.9m/d at the stream outlet, respectively. The Reynolds number is 555. Again, it should be mentioned that the influence of any variations in channel geometry on the stream flow is neglected in the computation of these parameters. Plot A in Figure 5.6 shows the steady-state stream flow velocity profiles in several crosssections along the channel (at the locations of X=100m, X=200m, X=300m, X=400m and X=500m) for the 1st realization of the Base Case, while Plot B shows the stream flow velocity profiles in the same set of cross-sections for the case where the fractured bedrock is replaced with a hydraulically equivalent homogeneous and anisotropic porous medium. The velocity profiles shown in Plot A indicate that the stream flow is slightly slower on the source zone side (the right side, looking downstream) than the other side in the cross-sections at X=100m, 200m, 300m and 400m, but it is opposite at the stream outlet cross-section at 124  X=500m. The stream flow velocity distributions are non-symmetric in all the cross-sections and vary in different locations in the case with the fractured bedrock in the subsurface, despite the fact that there is no influence of the channel geometry on the stream flow. The results demonstrate that the stream flow in the channel is non-uniform, because the existence of the fracture network in the bedrock influences the spatial pattern of groundwater discharge entering the stream channel. By contrast, the velocity profiles shown in Plot B of Figure 5.6 indicate that the stream flow velocity distributions are symmetric in all the cross-sections in the case where the fractured bedrock is replaced with the equivalent porous medium, and the stream flow is uniform. The two-dimensional steady-state surface water flow field in the first realization of the Base Case is shown in Figure 5.7, in which Plot A shows the hydraulic head contours together with the flow velocity vectors in the surface domain; Plot B shows the water depths in the surface domain; Plot C show the surface water flow velocity in the dominant flow direction from the upstream to the downstream; and Plot D shows the exchange of water between the surface and the subsurface (a negative value representing flow from the surface to the subsurface). The results indicate that the overland flow is in the direction towards the centerline of the channel (Plot A), the surface flow mainly occurs inside the channel (Plot B, C), and overland runoff occurs only in the area close to the stream banks (Plot D). Surface water recharges the groundwater system in most of the area outside the channel especially in the upstream, and groundwater discharges into the channel, forming a gaining stream along the entire reach (Plot D). The three-dimensional steady-state groundwater flow field in the first realization of the Base Case is shown in Figure 5.8, in which Plot A shows the hydraulic head contours in the 3D subsurface domain; Plot B shows the saturation profile in the domain; Plot C and D show 125  the head contours at the 2D horizontal cross-sections in the unconfined porous aquifer (at Z=120m) and in the bedrock (at Z=110m); Plot E and F show the head contours together with the velocity vectors at the 2D vertical cross-sections in the two units (at Y=100m beneath the streambed, at X=60m across the contaminant source zone, at X=250m in the centre of the domain, and at X=500m on the outflow boundary). The hydraulic head contours shown in the plots of Figure 5.8 demonstrate that due to the existence of the fracture network, the groundwater flow field is irregular in both the bedrock and the overlying unconfined porous aquifer (Plot D, E and F), and the head contours in the unconfined aquifer tends to resemble the patterns of head contours in the bedrock (Plot C and D). Test simulations of coupled SW-GW flow using the first realization of the Base Case showed that such a phenomenon did not exist when the fractured bedrock was replaced with an equivalent porous medium. This is similar to what was observed in the Base Case for the previous groundwater analysis (Figure 3.10 and 3.11 in Chapter 3). The head contours in the plots of Figure 5.8 also show that the subsurface flow is dominantly in the direction from the upstream to the downstream towards the centreline of the domain where stream channel is located (Plot A, C, D and F). Groundwater discharges into the stream channel along the entire reach, and groundwater flow velocity beneath the stream (at the top of the subsurface domain) is generally smaller in the upstream and larger in the downstream (Plot E and F), but the orientation of the flow velocity vector varies from location to location due to the influence of the fracture network. The subsurface domain is unsaturated near the top edges, as shown in the saturation profile in Plot B, but it is fullysaturated in the bedrock and also along the centreline on the top of the unconfined aquifer, where the excess saturation produces the stream flow in the channel.  126  (2) Factors influencing solute concentrations in stream channel  Prior to a discussion of plume detection probabilities, it is important to understand the mechanisms and main factors that influence solute concentrations in the stream channel. Figure 5.9 and 5.10 illustrate solute transport processes in the surface water and groundwater system at a time before the contaminant is detected at the compliance boundary in the first realization of the Base Case. Figure 5.11 shows the solute concentrations in the stream water and in the streambed groundwater along the stream channel centreline in this realization, and Figure 5.12 shows the concentration profiles in the stream water at the cross-sections along the stream channel. Note that the streambed groundwater refers to the groundwater beneath the stream channel in the top layer of the subsurface nodes, and also that all the plots are shown with the vertical exaggeration 6× (same as the plots in all other SW-GW scenarios). Figure 5.9 shows the migration of the plume (with concentrations above the detection threshold of 1.0×10-3g/L) in the surface water and groundwater system, indicating that the contaminant is first carried downward from the source zone by recharge to the unconfined porous aquifer and the bedrock, and then is carried upward with groundwater discharge into the stream water. It takes approximately 1500 days (Plot A) for the contaminant released from the source to reach the stream channel at a concentration above the detection threshold (Note that there is solute entering into the stream water before this time, but at concentrations below the detection threshold). The plume at concentrations above the detection threshold moves for about 60m downstream in the channel in the time period from 1500 days to 2000 days (Plot B), and about 100m further downstream in the channel in the time period from 2000 days to 2500 days (Plot C), until it is detected at day 3184 at the stream outlet compliance boundary with a concentration above the detection threshold (Plot D). It took 127  about 1684 days for the plume to migrate to the stream outlet, a distance of about 450m from the location where the contaminant first entered the stream with a concentration above the detection threshold, despite the fact that the mean flow velocity at the stream outlet is 1141.9m/d (Table 5.2, Figure 5.5). As the stream is gaining flow from groundwater along the entire reach (Plot D of Figure 5.7, Plot E of Figure 5.8), the result here indicates the important effect of dilution in the solute transport process in the stream water, causing the delay in detection of the contaminant at the stream outlet (compliance boundary). Figure 5.10 shows more details of the plume (with concentrations above the detection threshold) in stream water and groundwater at the time the contaminant is detected at the compliance boundary. Plot A shows the plume in stream water, indicating the solute is diluted by over one order of magnitude from the highest concentration of about 1.7×10-2g/L near the head of the stream to the lowest concentration 1.0×10-3g/L (the detection threshold) at the stream outlet. Plot B shows the plume in groundwater in a cross-section at Y=70m (the source zone centre location), indicating the contaminant transport occurs in both the unconfined aquifer and the bedrock. Plot C shows the plume in stream water together with the plume in groundwater in cross-sections at X=50m, X=60m, X=70m, X=90m, X=150m, X=200m and X=250m, which demonstrates that contaminant from the source zone first enters the stream channel through the subsurface flow paths (with groundwater discharge), not through an overland flow path. Plot D shows the plume in stream water together with the plume in groundwater in a cross-section at Y=100m (beneath the streambed), indicating contaminant at concentrations above the detection threshold enters the stream water from the subsurface at two locations along the stream centreline. An enlarged figure of Plot D shows that at both locations (from X=80m to X=120m, from X=140m to X=220m), the solute mass enters the stream water through the streambed on the source zone side (the right side, looking 128  downstream) between the stream bank and the stream centreline. Both Plot C and D show that as the contaminant migrates along the channel, the plume in the stream water is mixing with new solute mass entering the stream channel with groundwater discharge at different locations. Figure 5.11 shows the distribution of solute concentration in the stream water (for the plume shown in Plot A in Figure 5.10) and also the concentration in the streambed along the channel centreline (at the locations of the five SW observation points and also the locations of X=80m and X=180m). The relation between the concentrations in the stream water and the streambed indicates that the solute in the streambed groundwater discharges into the stream water at several locations along the centreline of the stream channel (for example at X=80m and X=180m, highlighted in 2 rectangles in Figure 5.11). The contaminant first enters the stream channel with groundwater discharge at X=80m (near the head of the stream and the source zone). At this location, groundwater in the streambed has a solute concentration of 1.6×10-2g/L, while the solute concentration in the stream water is 8.5×10-3g/L (at the time the contaminant is detected at the compliance boundary). The solute concentration in the stream water reaches a higher value of 1.0×10-2g/L at the stream centre at the downstream location of X=100m. It is diluted by 2.4 times to a concentration of 4.2×10-3g/L at the location X=180m, where the solute concentration in the streambed groundwater is 4.8×10-3g/L (slightly higher than that in the stream water). This result indicates that additional solute mass is entering the stream with groundwater discharge at this location, and mixing with the plume in the stream water. From that location to the stream outlet, the plume is further diluted by more than 4 times in the distance of 320m (from X=180m to X=500) due to the increase of stream flow in the channel (note that groundwater entering the stream channel with solute concentrations less than the detection threshold contributes to the calculation of the stream 129  water concentrations). The result indicates that because of the rapid dilution, sampling streambed groundwater is necessary to maximize the detectability of the contaminants entering the stream, and sampling stream water alone may not be enough. This finding is similar to what Conant et al. (2004) found in their field study of groundwater plume discharging to a river. Figure 5.12 shows the solute concentration distribution in the stream water for the plume shown in Plot A of Figure 5.10, at the cross-sections of X=100m, X=200m, X=300m, X=400m and X=500m along the stream channel. The concentration profile at the crosssection X=100m (Plot A) shows that solute mass enters the stream from the side of the continuous source zone (the right side bank, looking downstream) at the location of Y=96m, with the highest concentration 1.3×10-2g/L in the stream water, but the solute is diluted by a factor of about two when it spreads laterally to the other side of the stream with the lowest concentration of 7.1×10-3g/L at the location of Y=104m (the left side bank). The concentration profile in the cross-section X=200m (Plot B) shows that as the solute is carried downstream by the stream flow, the location of the peak concentration shifts 1m to the centreline (on the source zone side), and the plume with higher concentrations are still on the source zone side of the stream. The concentration profile at the cross-section X=300m (Plot C) shows that the location of the peak concentration shifts to the stream centreline, and the plume with higher concentrations is not on the source zone side of the stream, but on the other side of the stream. The concentration profile at the cross-section at X=400m (Plot D), however, shows that the plume with higher concentrations moves back to the source zone side again at this location because of dilution (by flow augmentation) plus dispersion (by the faster stream flow on the other side of the stream, see Plot A in Figure 5.6). As the plume moves farther away from the source zone to the stream outlet cross-section at X=500m, the 130  solute concentrations are distributed more evenly across the stream (Plot E), because more dilution of the plume occurs at this location with the increased volume of water in the stream channel. Overall, it is shown in Figure 5.12 that contaminant dilution is greater than one order of magnitude (13 times) from the cross-section X=100m to the stream outlet, while the total discharge of the stream increases by a factor of about 18 from 6.0m3/d to 107.8m3/d (Figure 5.3). The concentrations near the source zone and the stream centreline are higher, and diminish further downstream and towards the stream banks. The results demonstrate that the influence of the source zone location on the solute concentration in the stream water is reduced gradually from the upper reach to the lower reach. The influence of the fracture network in the bedrock on the stream water flow is also reflected in the solute transport and the concentration distributions in the stream.  5.3.2 Statistics of flow and transport in all realizations  Figure 5.13 illustrates the statistics of the water depth (Plot A), the total discharge (Plot B) and the mean flow velocity (Plot C) at the stream outlet with the number of realizations used in the Monte Carlo analysis. Figure 5.14 illustrates the statistics of the maximum (the highest) solute concentrations in the stream (Plot A) and the detection time at the compliance boundary (Plot B) with the number of realizations. Note that the minimum, the maximum and the average shown in Figure 5.13 and 5.14 represents the smallest, the largest and the arithmetic average values of each parameter, for each set of additional realizations. The statistics of these parameters all have converged at 100 realizations in the Base Case.  131  Table 5.3 shows that on average, the water depth, the total discharge, the mean velocity and the Reynolds number at the stream outlet is 1.19cm, 109.9m3/d, 1154.6m/d and 565.9, respectively, indicating that the stream is shallow with a small flow rate but with relative high flow velocity, and the flow is the upper bound of laminar in the Base Case. Table 5.4 shows that the detection time at the compliance boundary is 7.9 years on average in the Base Case, and the highest solute concentration in the stream is 1.62×10-2g/L on average (at the time the contaminant is detected at the compliance boundary), which is over one order of magnitude (16 times) higher than the detection threshold concentration (1.0×10-3g/L). The variability in the stream flow and solute concentration, and the detection time at the compliance boundary originates from the variability of the fracture networks from realization to realization, as the unconfined aquifer is homogeneous, isotropic and deterministic, and the recharge on the surface is uniform and deterministic. The coefficient of variation is 0.5%, 1.0% and 0.6%, respectively, for the water depth, total discharge and mean velocity of stream flow (Table 5.3). The coefficient of variation is 22.0% and 9.1%, respectively, for the maximum solute concentration in the stream and the detection time at the compliance boundary (Table 5.4). The results demonstrate that change of the fracture networks from realization to realization causes a small variability in the stream flow, but relatively large variability in the solute concentrations in the stream. This implies that the solute concentrations in the stream are more sensitive to the variability in the fracture networks in the bedrock, in comparison to the volume of the stream flow.  132  5.3.3 Detection probabilities for the Base Case monitoring network  Figure 5.15 shows the detection probabilities in the stream water for the Base Case at the 9 observation points at the stream outlet (compliance boundary). The results indicate that in the Base Case, there is a 100% probability of detecting the plume at the stream outlet (for the assumed value of the detection threshold) when sampling at the stream centreline. The probability of detecting the plume is 50% for a sampling location 1m to the side of the centreline of the plume, and 0% if the sampling location is at the stream bank. The concentration profiles in Figure 5.12 suggest that sampling the stream water in the upper reach on the source zone side and in the lower reach at and near the stream centreline will have a greater chance to detect the contaminant. The shorter is the distance from the stream to the source zone, the greater is the benefit to sample on the source zone side. Figure 5.16 shows the detection probabilities in the groundwater system of the unconfined aquifer (Plot A) and the bedrock (Plot B) for the Base Case. The ratio of saturated hydraulic conductivity Kuc (2.5×10-6m/s) of the unconfined aquifer to the bulk hydraulic conductivity Kx of the bedrock is 0.8 (see Table 5.1). For such a case in the earlier evaluation of the groundwater monitoring network (Section 3.4 of Chapter 3), it was concluded that detection at the compliance boundary occurs in the bedrock when the ratio of the hydraulic conductivity of the unconfined aquifer to that of the fractured bedrock (Kuc/Kx) is less than 1. The results in Figure 5.16, however, show that detection of the contaminant at the compliance boundary occurs in the surface water system. This response indicates that detection of the contaminant in the coupled SW-GW system is more complex than the case if only groundwater is considered, and the factors influencing the solute concentrations in the stream water and the surface water – groundwater interactions play important roles in 133  determination of the probabilities of detecting the contaminant in the catchment. Figure 5.16 shows the area encompassed by the detection probability contours in the bedrock is larger than that in the unconfined aquifer, and the detection probabilities in the bedrock are generally higher than those in the unconfined aquifer at the same locations of the monitoring wells. However, the detection probability contours indicate that, when measured in terms of concentrations above the detection threshold, the plume in the groundwater system does not move as far from the source in both the unconfined aquifer and the bedrock, due to the time of detection in the surface water system. In addition, the results in Figure 5.16 show that due to the presence of the fracture network, the detection probability contours are irregular in both the bedrock and the overlying unconfined aquifer, and the probability contours in the unconfined aquifer tends to resemble those in the bedrock. This phenomenon is similar to what was observed in the earlier analysis of the groundwater monitoring, and it can be explained with the patterns of the irregular flow fields in the two units shown in Figure 5.8.  5.4 Influence of hydraulic conductivity of unconfined porous aquifer  From the earlier analysis of the groundwater monitoring without coupling surface water, it was observed that the hydraulic conductivity of the unconfined porous aquifer is an important factor influencing groundwater monitoring. To evaluate the integrated SW-GW monitoring network, the hydraulic conductivity (Kuc) of the unconfined aquifer is increased by one order of magnitude to 2.5×10-5m/s in Scenario SW-GW-1, such that the hydraulic conductivity (Kuc) of the unconfined aquifer is about one order of magnitude higher than the  134  average of the bulk hydraulic conductivity of the bedrock (8.1 for Kuc/Kx, 7.8 for Kuc/Ky, and 31.3 for Kuc/Kz). The results of the coupled SW-GW steady-state flow in the first realization of this scenario are shown in Figure 5.17. The surface water depth (Plot A), surface flow velocity Vx (Plot B), and the flux exchange between the surface and the subsurface (Plot C) all indicate that in comparison to the Base Case, the length of the headwater stream becomes shorter when the unconfined porous aquifer is more permeable. The increase of the permeability of the unconfined aquifer results in a deeper water table and a larger unsaturated zone especially in the upstream part of the headwater basin (Plot D), indicating that the recharge tends to infiltrate into the deeper zone in the upstream, with less lateral flow towards the channel (Plot E). Consequently, the stream gains flow from groundwater discharge only in the lower reach of the channel (Plot F). Figure 5.18 shows the comparison of the stream water depth, total discharge and mean flow velocity along the stream channel for the first realizations of Scenario SW-GW-1 and the Base Case, indicating that in the case with a higher permeability unconfined aquifer, the stream is shorter (210m long from X=290m to X=500m) with smaller water depth, discharge and flow velocity along the stream channel except at and near the stream outlet. The water depth, the total discharge, the mean velocity and the Reynolds number at the stream outlet are 1.24cm, 119.1m3/d, 1200.6m/d and 613.2, respectively (Table 5.2), which are slightly larger than those in the Base Case. The results demonstrated that in the case with a higher permeability unconfined aquifer, more groundwater discharges into the channel at and near the stream outlet, while less groundwater discharges into the channel in the upper part of the basin.  135  The plume (concentrations above the detection threshold) in the surface water and groundwater system at the time the contaminant is detected at the compliance boundary for the first realization of Scenario SW-GW-1 is shown in Figure 5.19. Plot A shows the plume in stream water, Plot B shows the 3D envelop of the plumes in stream water and groundwater, Plot C shows the plume in the unconfined aquifer and the bedrock in the crosssections at Y=70m (source zone centre location), and Plot D shows the plume in stream water together with the plume beneath the streambed in groundwater in the cross-section at Y=100m (the centreline of the domain). In comparison to the results in Plot D of Figure 5.9 and Plot A, B, D of Figure 5.10 for the first realization of the Base Case, the contaminant migrates farther longitudinally and laterally in both the higher permeability unconfined aquifer and the bedrock before entering the stream water with groundwater discharge in the lower reach of the channel. The plume almost reaches the compliance boundary at concentrations above the detection threshold in both subsurface units when the contaminant is first detected in the stream water at the stream outlet at 15.2 years, an increase of 74.7% compared with the detection time at the compliance boundary (3184 days, approximately 8.7 years) in the first realization of the Base Case, because of a longer travel path in the subsurface due to the higher permeability of the unconfined aquifer. The highest concentration in the stream water for this case is reduced to about 2.6×10-3g/L (slightly higher than the detection threshold concentration 1.0×10-3g/L, Plot A in Figure 5.19), as the headwater of the stream is farther away from the source zone in the catchment. The relation between the solute concentration in the stream water and the streambed groundwater along the channel, as shown in Figure 5.20, indicates that the solute mass enters the stream water with groundwater discharge through the streambed in the section of the stream from X=330m to X=430m. The highest concentration in the streambed groundwater 136  in this section is 6.7×10-3g/L (at the location of X=350m), 3.7 times higher than the concentration of 1.8×10-3g/L in the stream water at the same location. Dilution occurs as the solute enters the stream water from the streambed groundwater. The concentration profiles in the selected cross-sections along the stream channel in Figure 5.21 show that at the time the contaminant is detected at the compliance boundary, concentrations above the detection threshold first appear on the source zone side of the stream at the cross-section X=337m (Plot A). The concentrations are higher at the crosssection X=350m (Plot B), as solute enters the stream on the source zone side in the section between X=337m and X=400m. The highest concentration (2.4×10-3g/L) of the plume occurs in the cross-section of X=400m (Plot C). As more groundwater discharges into the lower reach of the stream (especially at and near the stream outlet), the plume is diluted, and the locations of the highest concentrations gradually shift to the centreline from the cross-section X=450m (Plot D) to the cross-section X=500 at the stream outlet (Plot E). For the purpose of monitoring, the results in Figure 5.20 indicate that sampling both the stream water and the streambed groundwater will enhance the detectability in the stream before the contaminant reaches the compliance boundary, similar to what was found in the Base Case. In addition, the results in Figure 5.21 indicate that the sampling should be taken on the source zone side in the upper reach, but the location should be changed to positions near or at the stream centreline as the plume moves to the lower reach. The statistics of the stream flow parameters, the maximum concentration in stream water and the detection time at the compliance boundary have converged in simulation of 150 realizations of coupled SW-GW flow and transport in Scenario SW-GW-1 (Figure A.1 and A.2 in Appendix A). The smaller variations shown in Table 5.3 and 5.4 for this scenario demonstrate that with the higher permeability unconfined aquifer, the influence of the 137  fracture networks on the stream flow and solute concentrations is reduced. Due to the longer travel path from the source zone to the stream, the maximum solute concentration is reduced to 2.1×10-3g/L on average, compared to the Base Case (16.2×10-3g/L). The average detection time at the compliance boundary is increased by 83.5% from the Base Case (7.9 years). The permeability of the unconfined aquifer has a significant impact on stream flow and solute transport. Figure 5.22 shows the detection probabilities in the stream water at the stream outlet cross-section, while Figure 5.23 shows the detection probabilities in the unconfined aquifer (Plot A) and the bedrock (Plot B) at the time the contaminant is detected at the compliance boundary. Together with the plumes shown in Figure 5.19, the results suggest that when the hydraulic conductivity of the unconfined aquifer is one order of magnitude higher than that of the bedrock, the contaminant can be detected at the compliance boundary, for the basin geometry, in either the stream water or in the groundwater system (in both the unconfined aquifer and the bedrock), instead of only in the stream water in the Base Case. Figure 5.22 shows that the detection probabilities in the stream water are less than 1, indicating that there are some realizations where contaminants are detected at the compliance boundary in the subsurface rather than at the stream outlet, and that sampling the stream water on the source zone side will have a greater probability to detect the contaminant in the stream. In comparison to the Base Case (Figure 5.15), the skew of detection probability distribution shown in Figure 5.22 indicates that the influence of the source zone location on detection of the plume in the stream water is increased when the unconfined aquifer is more permeable. Figure 5.23 shows that the detection probabilities in groundwater are dramatically increased in the higher permeability unconfined aquifer (Plot A) and the bedrock (Plot B). At the time the contaminant is detected at the compliance boundary, there is very high 138  possibility that the contaminant will be detected in both units in those wells directly downgradient from the source zone. The results imply that in the case with a higher permeability unconfined aquifer, monitoring the contamination in the groundwater system will have greater probability of detection than monitoring the contamination in the headwater of the stream. The areas encompassed by the detection probability contours and the high detection probabilities indicate that the uncertainty in detection of the contamination in both subsurface units is reduced significantly, due to the increase of the permeability of the unconfined aquifer and the ratio of the permeabilities of the two units. Figure 5.23 also shows that due to the fracture networks, the detection probability contours are irregular in both the bedrock and the overlying unconfined porous aquifer, and the probability contours in the unconfined aquifer also tends to resemble those in the bedrock, similar to the phenomenon observed in the Base Case. In addition, the high detection probabilities and the area encompassed by the probability contours in the higher permeability unconfined aquifer are similar to that shown in Figure 3.7 for the similar case in Scenario GW-1 for the groundwater analysis (Section 3.4 in Chapter 3), but they are different in the bedrock.  5.5 Influence of fracture apertures in the bedrock  5.5.1 Influence of aperture of horizontal fracture set  In the earlier analysis of a groundwater monitoring network, it was observed that apertures of both horizontal and vertical fractures were important in influencing groundwater monitoring. To evaluate the integrated surface water and groundwater monitoring network, 139  the mean aperture of the horizontal fracture set is increased by a factor of 1.5 (Scenario SWGW-2a). This leads to an increase in the bulk hydraulic conductivity of the bedrock by a factor of 2.5 on average for Kx and Ky, and 1.6 for Kz (Table 5.1). The coefficient of variations increases by 1.0% for Kx, 1.1% for Ky, but 9.8% for Kz, in comparison to the Base Case. The increase of the horizontal apertures causes more variability in the permeability of the bedrock, especially in the vertical direction. The ratio of hydraulic conductivity of the unconfined porous aquifer to that of the bedrock (Kuc/Kx) is reduced to 0.3 from that in the Base Case (Kuc/Kx=0.8). A comparison of steady-state stream flow in the first realization for this scenario with the Base Case in Figure 5.24 shows that with the increase of the fracture apertures and the permeability of the bedrock, the stream water depth, total discharge and mean flow velocity are smaller along the entire reach. They are reduced by 9.3%, 10.5% and 1.3%, respectively, at the stream outlet (Table 5.2). The stream is 408m long (from X=92m to X=500m), 57m shorter than the stream in the Base Case. The results show that while there is less water flow in the stream, more water flows in the subsurface to the downstream outflow boundary, especially through the more permeable bedrock. Although the difference is hardly noticeable in the visualized groundwater flow system in Figure B.1 in Appendix B, the fact is the increase of the apertures of the horizontal fracture set will significantly increase the permeability of the fractures according to the cubic law, which will lead to more and faster preferential groundwater flow in the bedrock (in both horizontal and vertical directions). The contaminant plume (concentrations above the detection threshold) in the surface water and groundwater system of the first realization in this scenario in Figure 5.25 shows exactly the effects of the increased preferential flow in the fractures. Plot A shows the plume in the stream water, Plot B shows the 3D envelop of the plume in the groundwater system, 140  Plot C shows the plume in the unconfined aquifer and the bedrock in the cross-section at Y=70m (surface zone centre location), and Plot D shows the plume in stream water together with the plume beneath the streambed in the cross-section of Y=100m (the centreline of the domain). It can be seen that the contaminant is first detected at the compliance boundary in the bedrock (Plot B), not in the overlying unconfined aquifer or at the stream outlet (Plot A). The fractured bedrock with the larger apertures for the horizontal fracture set has a significant impact on contaminant transport in both the stream water and the groundwater system (Plot C, D). An enlarged figure of Plot D in Figure 5.25 shows that the solute mass enters the stream on the source zone side through the streambed at 9 locations in total along the stream channel between X=160m and X=440m, but the concentrations above the detection threshold are only present in the stream section from X=200m to X=400m. Figure 5.26 shows the solute exchange between the stream water and the streambed along the stream centreline. The highest concentration is 3.7×10-3g/L in the streambed and 1.3×10-3g/L in the stream water (slightly higher than the detection threshold). Both values are much lower than those in the first realization in the Base Case (Figure 5.11 and 5.12) because a greater proportion of the solute migrates to the compliance boundary through the fractured bedrock. Figure 5.27 shows that the solute concentrations in the stream water are higher on the source zone side in all the cross-sections from X=160m to X=400m, but the concentrations above the detection threshold only appear at the cross-sections at X=300m and X=400m. These results together with Plot A of Figure 5.25 indicate that although the contaminant enters the stream, it is not detected at the stream outlet (the compliance boundary) due to the effect of dilution by flow augmentation. The statistics of the stream flow parameters, the maximum concentration in stream water and the detection time at the compliance boundary have converged in simulation of 200 141  realizations of coupled SW-GW flow and transport in Scenario SW-GW-2a (Figure B.2 and B.3 of Appendix B). The results in Table 5.3 and 5.4 show that the increase in the horizontal fracture apertures in the bedrock results in a significant reduction of the stream flow, a significant increase in the variability of the detection time at the compliance boundary, and the variability in solute concentration in the stream water. Figures 5.28 and 5.29 show the detection probabilities of the SW-GW monitoring network in the unconfined aquifer (Plot A) and the bedrock (Plot B) in this scenario. The increase of the average of apertures of the horizontal fracture set by a factor of 1.5 has a significant impact on monitoring the contamination in stream water and groundwater. The highest detection probability at the stream outlet is 60% (Figure 5.28), which indicates that there is a 40% probability of failure to detect the contaminant at this location. However, detection probabilities in the subsurface (Figure 5.29) show that with the larger horizontal fracture apertures in the bedrock, the likelihood to detect the contaminant in the bedrock is much higher than that in the unconfined aquifer, before the contaminant reaches the compliance boundary. Also there is more likelihood to detect the contaminant in the bedrock at the compliance boundary.  5.5.2 Influence of apertures of vertical fracture sets  In Scenario SW-GW-2b, the mean apertures of both sets of vertical fractures are increased by a factor of 1.5, which leads to an increase of the bulk hydraulic conductivity of the bedrock by a factor of 1.3 on average for Kx and Ky, and 2.0 for Kz, in comparison to the Base Case (Table 5.2). The coefficient of variation is increased by 0.4%, 0.5% and 13.6% for Kx, Ky and Kz, respectively. The ratios of hydraulic conductivities (Kuc/Kx, Kuc/Ky and Kuc/Kz) 142  between the unconfined porous aquifer and the bedrock are 0.6, 0.6 and 1.6, respectively. In comparison to the case with larger horizontal fracture apertures in Scenario SW-GW-2a, the increase of the vertical fracture apertures results in a smaller change in the average and the variation of Kx and Ky, but a larger change in the average and the variation of Kz. The comparison of steady-state stream flow in the first realization for this scenario and the Base Case in Figure 5.30 shows that the increase of the vertical fracture apertures also results in less stream flow along the entire reach. As shown in Table 5.2, the stream water depth, total discharge and flow velocity are slightly reduced. The stream is 451m long (from X=49m to X=500m), 14m shorter than in the Base Case. The results indicate that in comparison to the horizontal fracture apertures, the vertical fracture apertures have less impact on the overall surface water and groundwater flow system in the case studied (Figure C.1 of Appendix C). This is because the total volumetric density of the two sets of vertical fractures in the bedrock is 5 times smaller than that of the horizontal fracture set (Table 2.2). The contaminant plumes (concentrations above the detection threshold) in surface water and groundwater systems of the first realization in this scenario are shown in Figure 5.31. Plot A shows the plume in the stream water, Plot B shows the 3D envelop of the plume in the groundwater system, Plot C shows the plume in the unconfined aquifer and the bedrock in the cross-section at Y=70m (surface zone centre location), and Plot D shows the plume in stream water together with the plume beneath the streambed in the cross-section of Y=100m (the centreline of the domain). The contaminant is first detected at the compliance boundary at the stream outlet as in the Base Case (Plot B). The larger vertical fracture apertures in the bedrock result in a plume with lower concentrations in the shorter stream (Plot A), and a larger plume in the bedrock and in the unconfined aquifer (Plot B, C, D). An enlarged figure of Plot D shows the solute mass enters the stream water on the source zone side along the 143  entire upper and middle reach of the stream (between X=80m and X=280m). The highest concentration 1.1×10-2g/L (at X=200m) in the streambed and 3.8×10-3g/L (at X=100m) in the stream water (Figure 5.32), and the solute concentrations in all of the cross-sections along the stream channel (the plots in Figure 5.33) are smaller than those shown in Figure 5.12 for the Base Case, indicating that the larger vertical fractures result in less solute transport in the stream water and more solute transport in the subsurface, especially in the bedrock. The statistics of the stream flow parameters, the maximum concentration in stream water and the detection time at the compliance boundary have converged in simulation of 200 realizations of coupled SW-GW flow and transport in Scenario SW-GW-2b (Figure C.2 and C.3 in Appendix C). The results in Table 5.3 and 5.4 show that in comparison to the Base Case, the increase in the vertical fracture apertures results in only a fairly small decrease of the stream water flow, but an increase of the detection time at the compliance boundary because of lower solute concentrations in the stream water. However, in comparison to the case with larger horizontal fracture apertures, the increase of the vertical fracture apertures leads to much smaller changes in flow and solute transport in the system. As shown in Figure 5.34, the detection probabilities are only slightly higher in the unconfined aquifer (Plot A) and the bedrock (Plot B) due to the increase in the vertical fracture apertures, but the detection probabilities in stream water at the stream outlet are identical as those in Figure 5.15 for the Base Case. Detection at the compliance boundary in all the realizations occurs at the stream outlet, not in the subsurface. In addition, the result in Figure 5.32 indicates that at the time the contaminant is detected at the compliance boundary (the stream outlet), sampling both the stream water and the streambed groundwater would enhance the probability to detect the contaminant. The result in Figure 5.33 indicates sampling the stream water on the source zone side generally increases the probability of 144  detecting the contaminant in the upper reach of the stream, while sampling at and near the stream centreline provide a greater chance to detect the contaminant as the contaminant migrates to the lower reach of the stream.  5.6 Influence of dispersivities in stream water  As discussed at the beginning of Chapter 4, dispersion (longitudinal and transverse spreading of contaminants in stream water) due to local variations of flow velocities is one of the major mechanisms influencing solute concentrations in stream water. Longitudinal and transverse dispersivities are important parameters in characterizing solute transport in stream water (Fischer et al., 1979). From the Base Case analysis in Section 5.3.1, it is known that the stream flow is nonuniform with flow velocities varying across the channel, which will affect solute dispersion in the stream. Scenario SW-GW-3a and Scenario SW-GW-3b investigate the influence of dispersion on solute transport and the concentration distribution in the stream, and the importance of dispersion relative to dilution by flow augmentation along the stream channel. The values of longitudinal dispersivity and transverse dispersivity for solute transport in the stream water are reduced by a factor of two, respectively, in comparison to the Base Case. The contaminant is first detected at the stream outlet compliance boundary in both scenarios, the same as in the Base Case. The plume geometries in the groundwater system for the two cases at the time the contaminant is detected at the compliance boundary are very similar to what was shown in Panel D of Figure 5.9 and Figure 5.10 for the Base Case. The results indicate that the transport processes in the stream water have an insignificant effect on  145  the contaminant transport in the groundwater system, in comparison to the influence of the groundwater system on the stream water solute concentrations. However, the comparison of solute concentrations in the stream water in the first realizations of Scenario SW-GW-3a with the Base Case in Figure 5.35 shows that on average, the concentration values in the cross-sections at X=100m (Plot A), X=200m (Plot B), X=300m (Plot C) and X=400m (Plot D) increase by 6.2%, 3.1%, 1.9% and 1.5%, respectively, in the case with the smaller longitudinal dispersivity. The magnitude of this increase in the concentration in all the cross-sections is relatively small, and it diminishes from the upper reach to the lower reach until there is no difference at the stream outlet crosssection at X=500m (Plot E), due to dilution by flow augmentation along the stream channel. In addition, the results show that the small increase in concentration in each of these crosssections is approximately the same for the source zone side of the stream and the other side of the stream. In comparison to the reduction of solute concentration due to dilution by flow augmentation along the stream channel (over one order of magnitude from upstream to downstream), the small magnitude of increase in solute concentration (with smaller longitudinal dispersivity) demonstrates that dispersion has less of an effect than dilution in influencing solute concentrations in the stream water, especially at the stream outlet where the observation points are located for the compliance monitoring. The comparison of solute concentrations in the stream water in the first realization of Scenario SW-GW-3b with the Base Case is shown in Figure 5.36. The results indicate that there is less lateral solute dispersion across the stream when the transverse dispersivity in the stream water is smaller. The concentration in the cross-section at X=100m (Plot A) increases on the source zone side of the stream where solute mass enters the stream from the continuous source, but decreases on the other side of the stream due to the less lateral 146  spreading of the plume. The concentration in the cross-section at X=200m increases in the central part of the plume (left to the stream centreline from Y=97m to Y=100m) due to the influence of the source zone, but decreases at the edge of the plume near the stream banks. In the meantime, there is no increase in the solute concentration in the centre of the plume (along the stream centreline) in the cross-sections at X=300m (Plot C), X=400m (Plot D) and X=500m (Plot E) due to dilution by flow augmentation along the stream channel, but there is a decrease in solute concentration at the edge of the plume near the stream banks in all three cross-sections. Less lateral dispersion results in a narrower plume along the stream channel. The change in concentration in each of the cross-sections due to the less lateral dispersion is relatively small, in comparison to the large magnitude of reduction in concentration due to dilution by flow augmentation. The statistics of the maximum concentration in stream water and the detection time at the compliance boundary have converged in simulation of 100 realizations of coupled SW-GW flow and transport in both Scenario SW-GW-3a (Figure D.1 in Appendix D) and SW-GW-3b (Figure E.1 in Appendix E). The results in Table 5.4 show that in comparison to the Base Case, with the smaller longitudinal dispersivity in the stream water in Scenario SW-GW-3a, the average detection time at the stream outlet is longer (by 36.5 days, an increase of 1.3%); with the smaller transverse dispersivity in the stream water in Scenario SW-GW-3b, the average detection time at the stream outlet is shorter (by 73 days, a decrease of 2.5%). The small difference in detection time at the compliance boundary explains why the plume geometries in the subsurface for the first realization of these scenarios are similar to that in the Base Case. The detection probabilities at the stream outlet and in the groundwater system for Scenario SW-GW-3a are identical to those shown in Figure 5.15 and 5.16 for the Base Case, 147  indicating that this difference in the longitudinal dispersivity in stream water does not change the detection probabilities at the stream outlet and in the subsurface units. The detection probabilities in the groundwater system for Scenario SW-GW-3b are also identical to those in Figure 5.16, indicating that the difference in the transverse dispersivity in stream water does not change the detection probabilities in the subsurface either. The results in Figure 5.37, however, show a smaller lateral dispersion leads to a lower detection probability at the sampling location 1m to the side of the centreline at the stream outlet, in comparison to the Base Case (Figure 5.15). The results show that the longitudinal dispersivity is less important than the transverse dispersivity in influencing the probabilities of detecting the contaminant at the stream outlet. This is likely because dilution of the plume (over one order of the magnitude) along the stream channel reduces the effect of longitudinal dispersion. However, with respect to monitoring the contamination along the stream channel before it reaches the stream outlet compliance boundary, the result in Figure 5.35 indicates that in general, smaller longitudinal dispersion tends to increase the probability of detection of the contaminant in the stream water, especially for a sampling location near the source zone in the upstream region of the channel and at the stream centreline in the downstream region of the channel. The results in Figure 5.36 indicates that smaller transverse dispersion tends to lower the probability of detection of the contaminant in the stream water, especially if a sampling location is at or near the stream banks in the lower reach.  5.7 Influence of detection threshold  In the earlier analysis of the groundwater monitoring network, the detection threshold was identified as a key factor influencing groundwater monitoring. To evaluate the influence 148  of the detection threshold on the integrated SW-GW monitoring network, Scenario SW-GW4 is designed with the detection threshold reduced to 1.0×10-4g/L (0.01% of the source zone concentration of 1.0g/L), one order of magnitude lower than that in the Base Case. The lower detection threshold applies at the compliance boundary in both the surface and subsurface domains, as well as at all the surface water observation points and the performance groundwater monitoring wells. Figure 5.38 shows the plume (concentrations above the detection threshold) in the surface water and groundwater system at the time the contaminant is detected at the compliance boundary in the first realization of this scenario. Plot A shows the plume in the stream water, Plot B shows the 3D envelop of the plume in the groundwater system, Plot C shows the plume in the unconfined aquifer and the bedrock in the cross-section at Y=70m (surface zone centre location), and Plot D shows the plume in stream water together with the plume beneath the streambed in the cross-section of Y=100m (the centreline of the domain). The contaminant is first detected at the stream outlet compliance boundary at the time of 4.8 years, about one half of the time (8.7 years) for the first realization of the Base Case. Recall that the detection at the compliance boundary is determined by comparing the solute concentration at any node directly to the detection threshold, and simulation of the SW-GW transport in each realization is terminated once the contaminant reaches the compliance boundary with a concentration above the detection threshold. With the lower detection threshold at the compliance boundary in this case, the contaminant is detected earlier at the stream outlet, and the transport simulation is terminated. In comparison to the Base Case (Figures 5.9 and 5.10), Figure 5.38 shows that at the time the contaminant is detected at the compliance boundary with the lower detection threshold, the concentrations of the plume in the stream are one order of magnitude lower 149  (Plot A), but the change in the geometry of the subsurface plume is relatively small in both the unconfined aquifer and the bedrock (Plot B, C, D). Also as shown in Figure 5.39, with the lower detection threshold, solute concentrations are about one order of magnitude lower in both the stream water and the streambed groundwater along the channel centreline. That’s because the shorter detection time at the compliance boundary results in less solute mass entering the stream, although the locations where the solute enters the stream from the subsurface are the same as in the Base Case (Figure 5.11). The pattern of solute concentration distribution across the stream in this case (Figure 5.40) is similar to the Base Case (Figure 5.12), but the concentrations in the stream are lower. Also as shown in Figure 5.40, the contaminant is detected at more locations at the stream outlet with concentrations above the lower detection threshold. The statistics of the maximum concentration in stream water and the detection time at the compliance boundary have converged in simulation of 100 realizations of coupled SW-GW flow and transport in Scenario SW-GW-4 (Figure F.1 in Appendix F). The results in Table 5.4 show that reducing the detection threshold by one order of magnitude results in a significant decrease in the average detection time at the compliance boundary (by 45.6%) with smaller variation. The results in Figure 5.41 show that in comparison to the Base Case, there are a greater number of observation points at the stream outlet with higher detection probabilities, which suggests that for stream water monitoring, a lower detection threshold will increase probabilities in detection of a contaminant, because a greater proportion of the plume is above the detection threshold, and target is therefore bigger. However, the results in Figure 5.42 show that the lower detection threshold leads to a relatively small change in detection probabilities in the unconfined aquifer (Plot A) and in the bedrock (Plot B), in comparison to 150  the Base Case (Figure 5.16). The results suggest that for integrated SW-GW monitoring, the probability of detection of a contaminant in stream water is more sensitive to the lower detection threshold than that in groundwater. Finally, because of the different flow system and boundary conditions, the results here are quite different from those with the lower detection threshold in the groundwater monitoring analysis (Scenario GW-5b in Section 3.6 of Chapter 3). Nevertheless, the results in both cases suggest that a lower detection threshold will increase the probabilities of detection of the contaminant in the bedrock.  5.8 Influence of source zone location  In the earlier analysis of groundwater monitoring, it was observed that the source zone location had a strong effect on the detection probabilities in the bedrock and the overlying unconfined aquifer. To evaluate the influence of the source zone location to integrated SWGW monitoring, Scenario SW-GW-5 is designed with the same parameters used in the Base Case, but with a different source zone location in the Z axis. In this scenario, the source zone is at Z=115m to Z=117m, located at the sediment – bedrock interface. Figure 5.43 shows the plume (concentrations above the detection threshold) in the surface water and groundwater system at the time the contaminant is detected at the compliance boundary in the first realization of this scenario. Plot A shows the plume in the stream water, Plot B shows the 3D envelop of the plume in the groundwater system, Plot C shows the plume in the unconfined aquifer and the bedrock in the cross-section at Y=70m (source zone centre location), and Plot D shows the plume in stream water together with the plume beneath the streambed in the cross-section of Y=100m (the centreline of the domain). 151  The results show that the contaminant is first detected at the stream outlet compliance boundary (Plot B). With the source zone at the sediment – bedrock interface, there is more solute transport in the bedrock (Plot C) at the time contaminant is detected at the compliance boundary, and the plume in the stream channel is slightly shorter with lower concentrations (Plot A). The results indicate that the distance from the source zone to the stream channel and to the bedrock is important in influencing the plume in the surface water and the groundwater. An enlarged figure of Plot D in Figure 5.43 shows that the solute mass enters the stream from the source zone side in the section between X=100m and X=240m along the channel, in comparison to the Base Case where the solute enters the stream in two locations (between X=80m and X=120m, between X=140m and X=220m). Also as shown in Figure 5.44, solute mass (with the highest concentration 1.4×10-2g/L in the streambed) enters the stream water (with the highest concentration 3.0×10-3g/L) at the location X=180m along the channel centreline, and this location is about 100m (downstream) from the location (X=80m) where the solute first entered the stream in the Base Case. This explains why the plume with concentration above the detection threshold is shorter and why the concentrations of the plume in the stream water are lower. The results demonstrate that the source zone location has a significant influence on the locations of the solute mass entering the stream from the subsurface, and the consequent concentrations in the stream water. The solute concentration profiles in cross-sections along the stream channel are shown in Figure 5.45. The results indicate that solute mass enters the stream water on the source zone side at the cross-sections at X=100m (Plot A) with lower concentrations, and at X=200m (Plot B) with higher concentrations. When the plume is carried by the stream flow to the downstream location at the cross-section X=300m (Plot C), the concentrations on the source zone side of the stream are lower than those on the other side of the stream, similar to what 152  was observed in Figure 5.12 for the Base Case. However, because of dilution, the concentrations at the downstream cross-sections at X=400m (Plot D) and X=500m (Plot E) are more evenly distributed due to the dilution by the flow augmentation along the stream. The statistics of the maximum concentration in stream water and the detection time at the compliance boundary have converged in simulation of 100 realizations of coupled SW-GW flow and transport in Scenario SW-GW-5 (Figure G.1 in Appendix G). The results in Table 5.4 show that in comparison to the Base Case, the deeper source zone location (at the sediment – bedrock interface) causes more variability in detection time at the compliance boundary and solute concentrations in the stream water. The results in Figure 5.46 show that in the case that the source zone is located at the sediment – bedrock interface, the detection probabilities in the bedrock are higher, and the area where the contaminant is likely to be detected in the bedrock is larger (Plot B). The detection probabilities in the overlying unconfined aquifer are also slightly higher (Plot A). The detection probabilities at the stream outlet are identical to those shown in Figure 5.15 for the Base Case, indicating that changing the source zone location to the sediment – bedrock interface did not make a difference in detection of the contaminant at the stream outlet, due to the effect of dilution by flow augmentation along the stream. However, as indicated in Figure 5.45, the influence of the distance between the source zone and the stream channel on the detection probabilities in the stream water depends upon where the samples are taken along the stream channel. For example, when the source is located at the sediment – bedrock interface, there would be no detection of the contaminant if the samples are taken from the locations at the cross-section of X=100m, but there would be 100% probability of detection of the contaminant if the stream water were sampled at the locations in the cross-section of X=200m. The source zone location has more influence on the stream water monitoring in the 153  upper reach than in the lower reach, and has more influence on the sampling locations on the source zone side of the stream than the other side.  5.9 Summary  Using the methodology presented in Chapter 4, detection probabilities at stream water observation points and performance monitoring wells of an integrated SW-GW monitoring network prior to detection at the compliance boundary were evaluated in a headwater stream catchment where an unconfined aquifer overlies fractured bedrock. Factors influencing plume detection in integrated SW-GW monitoring were discussed based on the results from simulation of coupled SW-GW flow and solute transport in 8 scenarios. For the assumed recharge rate on the land surface and given boundary conditions, a shallow stream with small flow rate is modeled. Due to the existence of fracture networks in the bedrock, the stream flow is non-uniform with Reynolds numbers close to the upper bound of laminar, and the distribution of stream flow velocities is non-symmetric and varies with location. The groundwater flow field is also irregular in both the bedrock and the overlying unconfined aquifer due to the fracture networks. The stream receives water from groundwater discharge along the entire reach. Solute mass from the continuous source enters the stream with groundwater discharge through the streambed. It is observed that dilution by flow augmentation along the stream channel has important effects on solute concentrations and detection probabilities in the stream water. Even if the contaminant enters the stream, it may still not be detected at the stream outlet compliance boundary. Among the factors examined in influencing flow and solute transport in the surface water and groundwater system and the integrated SW-GW monitoring, the permeability of the 154  unconfined aquifer has a significant impact. A higher permeability of the unconfined aquifer results in a deeper water table and a shorter stream length in the headwater basin. A higher permeability of the unconfined aquifer leads to more solute transport and higher detection probabilities in the subsurface (in both the unconfined aquifer and in the bedrock), but less solute transport and lower detection probabilities in the stream water. The influence of the fracture networks on flow and solute transport, and detection probabilities in the stream water is reduced when the unconfined aquifer has a higher permeability. Fracture apertures also have significant impacts, especially the horizontal fracture set. They play a major role in controlling the permeability of the bedrock, and in determining the detection probabilities in the groundwater system, especially in the bedrock. Fracture apertures also influence the locations of solute mass entering the stream water through the streambed, and hence the detection probabilities in the stream water. In the case with larger aperture for the horizontal fracture set, the contaminant enters the stream, but it is not detected at the stream outlet compliance boundary. The vertical fracture apertures have less influence on detection probabilities in the stream water than the horizontal fracture apertures in the case studied. In general, larger horizontal or vertical fracture apertures leads to an increase in the permeability of the bedrock and the detection probabilities in the subsurface (especially in the bedrock with larger horizontal fracture apertures), but a decrease in the stream flow and the detection probabilities in the stream water. The influence of the fractured bedrock on the flow, solute transport and detection probabilities in the stream increases when the fracture apertures are larger. The longitudinal and transverse dispersivities are important factors influencing solute transport and detection probabilities in the stream water, but they have an insignificant impact on detection probabilities in the groundwater system. A smaller longitudinal 155  dispersivity results in slower migration of the plume front, a longer detection time at the stream outlet, and higher detection probabilities along the stream, especially if sampling locations are on the source zone side bank in the upper reach. A smaller transverse dispersivity results in less lateral spreading of the plume (a narrow plume), a shorter detection time at the stream outlet, and lower detection probabilities when sampling locations are at or near the stream banks in the lower reach. However, the effects of dispersion in influencing the solute distribution and detection probabilities in the stream water are much less, in comparison to the effects of dilution due to flow augmentation. The dispersivities (or dispersion) have much less effects than the other factors examined. The detection threshold is a key parameter influencing the detection probabilities of the integrated SW-GW monitoring. Stream water monitoring is more sensitive to the lower detection threshold than the groundwater monitoring. A lower detection threshold will generally increase the detection probabilities in both the stream water and the groundwater (especially in the bedrock), and reduce the probabilities of failure in detection of the contaminant (especially in the stream water). The source zone location (or the distance from the source zone to the stream channel and to the bedrock) is important in influencing solute transport and plume detection in the stream and in the groundwater system. The source zone location has a strong effect on the locations where solute enters the stream from the subsurface, which will influence the probability of detection of the plume in the stream water. The influence of the source zone on the solute concentrations in the stream water decreases as the plume is carried by the stream flow farther away from the source zone and diluted by the flow augmentation along the channel. Therefore, sampling the stream water on the source zone side in the upper reach will generally yield a greater probability of detection. The probabilities of detection of the 156  contaminant in the bedrock are also affected by the source zone location. The closer the source zone is to the bedrock, the more likely the contaminant will be detected in the bedrock. Finally, the results in this study suggest that due to the uncertainties in association with the locations of solute mass entering the stream through the streambed in the subsurface because of the influence of fractured bedrock, and also the significant dilution of solute concentrations in the stream water, sampling both the stream water and the streambed groundwater will enhance the detectability of the contamination in the stream.  157  Table 5.1 Bulk hydraulic conductivity of the bedrock in each scenario  Scenario  Kx (×10-6m/s)  Ky (×10-6m/s)  Kz (×10-6m/s)  Average  Coefficient of Variation (%)  Average  Coefficient of Variation (%)  Average  Coefficient of Variation (%)  3.1  4.2  3.2  4.2  0.8  3.8  SW-GW-2a  7.7  5.2  8.1  5.3  1.3  13.6  SW-GW-2b  3.9  4.6  4.1  4.7  1.6  17.4  Base Case SW-GW-3a SW-GW-3b SW-GW-4 SW-GW-5 SW-GW-1  Table 5.2 Stream water flow in the first realization of each scenario Stream Length (m)  Water Depth at Stream Outlet (cm)  Total Discharge at Stream Outlet (m3/d)  Mean Velocity at Stream Outlet (m/d)  Reynolds Number  465  1.18  107.8  1141.9  555.0  SW-GW-1  210  1.24  119.1  1200.6  613.2  SW-GW-2a  408  1.07  96.5  1127.3  496.8  SW-GW-2b  451  1.17  106.2  1134.6  546.8  Scenario Base Case SW-GW-3a SW-GW-3b SW-GW-4 SW-GW-5  158  Table 5.3 Properties of stream water flow of all realizations in each scenario  Scenario  Water Depth at Stream Outlet (cm)  Total Discharge at Stream Outlet (m3/d)  Mean Velocity at Stream Outlet (m/d)  Average Reynolds Number  Average  Coefficient of Variation (%)  Average  Coefficient of Variation (%)  Average  Coefficient of Variation (%)  1.19  0.5  109.9  1.0  1154.6  0.6  565.9  SW-GW-1  1.24  0.2  119.1  0.3  1198.0  0.2  611.9  SW-GW-2a  1.10  1.2  99.7  2.3  1137.4  1.2  515.3  SW-GW-2b  1.17  0.5  107.0  1.1  1144.4  0.6  551.5  Base Case SW-GW-3a SW-GW-3b SW-GW-4 SW-GW-5  Table 5.4 Monte Carlo realizations, solute concentrations in stream water, detection time and location at compliance boundary in each scenario  Scenario  Max. Concentration in Stream Water (×10-3g/L)  Number of Realizations  Detection Time  Detection Location  (yrs)  Average  Coefficient of Variation (%)  Average  Coefficient of Variation (%)  Stream Outlet  Unconfined Aquifer  Fractured Bedrock  Yes  No  No  Base Case  200  16.2  22.0  7.9  9.1  SW-GW-3a  100  16.9  17.9  8.0  8.2  SW-GW-3b  100  17.4  25.9  7.7  9.6  SW-GW-4  100  2.4  20.7  4.3  8.3  SW-GW-5  100  5.3  40.0  4.3  39.5  SW-GW-1  150  2.1  19.7  14.5  5.8  Yes  Yes  Yes  SW-GW-2a  200  1.8  39.2  6.4  21.1  No  Yes  Yes  SW-GW-2b  200  9.6  52.0  9.1  17.4  Yes  No  No  -3  -4  * Detection threshold 1×10 g/L for all scenarios, except for SW-GW-4: 1×10 g/L; “Yes”, “No” represents detection of contamination occurs or not at the location.  159  Heads vs Simulation Time 129.036  Head (m)  129.032 129.028 GW  129.024  SW  129.020  Steady-State  129.016  Transient State  129.012 60  65  70  75  80  85  90  95  100  Time (Yrs)  Figure 5.1 Evolution of hydraulic heads in the stream water and in the streambed groundwater at observation point (at the node at X=100m, Y=100m) in the channel for the first realization of the Base Case  Stream Flux (Nodal Value) vs Simulation Time 1.0  3  Flux (m /d)  0.8 0.6 Steady-State  0.4 Transient State  0.2 0.0 60  65  70  75  80  85  90  95  100  Time (Yrs)  Figure 5.2 Evolution of stream flux at observation point (at the node at X=100m, Y=100m) in the channel for the first realization of the Base Case  160  Total Discharge along the Stream Channel  3  Discharge (m /d)  120 100 80 60 40 20 0 0  50  100  150  200  250  300  350  400  450  500  X (m)  Figure 5.3 Total discharge along the stream channel at the steady state for the first realization of the Base Case  Water Depth along the Stream Channel 0.014 0.012  Depth (m)  0.010 0.008 0.006 0.004 0.002 0.000 0  50  100  150  200  250  300  350  400  450  500  X (m)  Figure 5.4 Water depth along the stream channel at the steady state for the first realization of the Base Case  161  Mean Flow Velocity (V x ) along the Stream Channel 1400  Velocity (m/d)  1200 1000 800 600 400 200 0 0  50  100  150  200  250  300  350  400  450  500  X (m)  Figure 5.5 Flow velocity along the stream channel at the steady state for the first realization of the Base Case  162  A  Flow Velocity Profiles along Stream Channel with Bedrock 104  194.97  Y (m)  103  402.07 328.76  730.81  862.18  649.15  915.98  1080.61  1192.03  102  353.71  101  357.21  691.18  954.11  1135.97  1229.67  100  357.73  692.70  956.01  1137.16  1230.13  99  357.11  690.60  953.98  1134.57  1229.98  98  684.09  353.05  97  401.57  300  591.97  915.64 727.78  500  1223.97  1126.92 1076.90  1228.38 1194.43  885.99  700 V x (m/d)  X=200m  1129.61  946.48  648.06  X=100m  B  946.73  683.03  328.53  96 194.81 100  900  X=300m  1100  X=400m  1300  X=500m  Flow Velocity Profiles along Stream Channel with Equivalent Porous Medium 104  202.38  103  Y (m)  592.15  404.50  576.65  341.29  719.63  890.40  653.41  102  367.59  101  371.60  893.73  689.21  1088.14  923.95  696.85  1191.63  1113.48  931.12  1216.71  1120.15  1220.79  100  372.32  698.75  932.99  1121.93  1223.46  99  371.60  696.85  931.12  1120.15  1220.79  98 97 202.38 96 100  367.59  689.21  341.29  923.95  653.41 404.50  300 X=100m  576.65  500 X=200m  1113.48  893.73 719.63  1088.14  1216.71 1191.63  890.40  700 V x (m/d) X=300m  900 X=400m  1100  1300  X=500m  Figure 5.6 Flow velocity profiles at stream cross-sections at the steady state for the first realization of the Base Case. A Stream with the bedrock in subsurface, B stream with the equivalent porous medium.  163  A  B  C  D  Figure 5.7 Coupled steady-state surface water flow in the first realization of the Base Case. A Hydraulic head (m) and flow velocity vectors; B water depth (m); C surface flow velocity Vx (m/d); and D flux exchange with the subsurface (m/d). Dimension (m).  164  A  B  C  D  E  F  Figure 5.8 Coupled steady-state groundwater flow in the first realization of the Base Case. A Hydraulic head (m) in the 3D domain, B saturation profile, C head (m) in cross-section at Z=120m (in the unconfined aquifer), D head (m) in cross-section at Z=110m (in the bedrock), E head (m) and flow velocity vectors in cross-section at Y=100m (beneath the streambed), and F head (m) and flow velocity vectors in cross-sections at X=60m, X=250m and X=500m. Dimension (m).  165  A  B  Concentration  Concentration  C  D  Concentration Concentration  Figure 5.9 Evolution of plumes (concentrations above detection threshold) in surface water and groundwater before detection of the contaminant at the compliance boundary in the first realization of the Base Case. A 1500 days, B 2000 days, C 2500 days, and D 3184 days (detection at the stream outlet). Concentration (g/L), Dimension (m).  166  A  B Plume in stream  Source Zone Concentration  C  Concentration  D Source Zone Solute entering SW  Solute entering SW  Concentration  Concentration  Figure 5.10 Plume in the SW-GW system (concentrations above detection threshold) at time of detection at compliance boundary in the first realization of the Base Case. A plume in stream water (2D); B plume in groundwater (cross-section Y=70m, source zone); C solute entering stream water from groundwater through streambed (not through overland pathways), cross-sections at X=50m, 70m, 80m, 100m, 150m, 200m, 250m; and D solute entering stream water through streambed, cross-section Y=100m (domain centerline). Concentration (g/L), Dimension (m).  167  Solute Concentrations in Stream Water (SW) and Streambed (GW)  Concentration (g/l)  2.0E-02  1.5E-02 SW 1.0E-02  GW  5.0E-03  Solute in GW entering SW 0.0E+00 0  50  100  150  200  250  X (m)  300  350  400  450  500  Detection threshold: 1.0×10-3g/l  Figure 5.11 Concentrations along the stream centreline (in stream water and streambed) at the time of detection at the compliance boundary in the first realization of the Base Case  168  A Concentration (g/l)  Concentration Profiles at Cross-Sections in Stream Channel 1.60E-02 1.10E-02 Solute entering SW from source zone  6.00E-03  X=100m  1.00E-03 96  97  98  99  100  101  102  103  104  Y (m)  Concentration (g/l)  B  3.75E-03 3.65E-03  X=200m  3.55E-03 3.45E-03 96  97  98  99  100  101  102  103  104  C  Concentration (g/l)  Y (m)  2.00E-03 1.98E-03 X=300m  1.96E-03 1.94E-03 96  97  98  99  100  101  102  103  104  D  Concentration (g/l)  Y (m)  1.35E-03 1.33E-03 X=400m  1.31E-03 1.29E-03 96  97  98  99  100  101  102  103  104  E  Concentration (g/l)  Y (m)  1.02E-03 1.00E-03 X=500m  9.80E-04 9.60E-04 96  97  98  99  100 Y (m)  101  102  103  104  (Detection threshold: 1.0×10-3g/l)  Figure 5.12 Concentration profiles in stream water at the time of detection at the compliance boundary in the first realization of the Base Case. Cross-sections at A X=100m, B X=200m, C X=300m, D X=400m, and E X=500m.  169  A  Water Depth at Stream Outlet vs Number of Realizations  Depth (m)  0.0123 0.0121 0.0119 0.0117 0.0115 0  25  50  75  100  125  150  175  200  Realizations Minimum  B  Maximum  Average  3  Discharge (m /d)  Total Discharge at Stream Outlet vs Number of Realizations 114 112 110 108 106 0  25  50  75  100  125  150  175  200  Realizations Minimum  C  Maximum  Average  Mean Flow Velocity (V x ) at Stream Outlet vs Number of Realizations  Velocity (m/d)  1200 1175 1150 1125 1100 0  25  50  75  100  125  150  175  200  Realizations Minimum  Maximum  Average  Figure 5.13 Statistics of water depth, discharge and flow velocity at stream outlet in different number of realizations for the Base Case. A Water depth, B total discharge, and C mean velocity (Vx). “Minimum”, “Maximum” and “Average” represents the smallest, the largest and the arithmetic average values of each parameter, for each set of additional realizations.  170  A Concentration (g/l)  Maximum Solute Concentration in Stream vs Number of Realizations 2.5E-02 2.0E-02 1.5E-02 1.0E-02 5.0E-03 0.0E+00 0  25  50  75  100  125  150  175  200  Realizations Minimum  B  Maximum  Average  Detection Time (yrs)  Detection Time at Compliance Boundary vs Number of Realizations 12 10 8 6 4 2 0 0  25  50  75  100  125  150  175  200  Realizations Minimum  Maximum  Average  Figure 5.14 Statistics of the maximum solute concentration in stream at the time of detection at the compliance boundary, and detection time at compliance boundary in different number of realizations for the Base Case. A Maximum solute concentration, B detection time. “Minimum”, “Maximum” and “Average” represents the smallest, the largest and the arithmetic average values of each parameter, for each set of additional realizations.  171  Detection Probabilities at Stream Outlet 1.0  Probability  0.8 0.6 0.4 0.2 0.0 96  97  98  99  100  101  102  103  104  Y (m)  Figure 5.15 Detection probabilities at the stream outlet in the Base Case  A  B  Figure 5.16 Detection probabilities of groundwater monitoring wells in the Base Case. A Unconfined porous aquifer, B fractured bedrock. 172  A  B  C  D  E  F  Figure 5.17 Coupled steady-state SW-GW flow in the first realization of Scenario SW-GW-1 (higher Kuc). A Surface water depth (m), B surface flow velocity Vx (m/d), C flux exchange (m/d) between surface and subsurface, D saturation profile in the subsurface, E hydraulic head (m) in the subsurface domain, and F groundwater hydraulic head (m) and flow velocity vectors in cross-section at Y=100m (beneath the streambed). Dimension (m).  173  A  Depth (m)  Water Depth along the Stream Channel 0.014 0.012 0.010 0.008 0.006 0.004 0.002 0.000 0  50  100  150  200  250  300  350  400  450  500  400  450  500  X (m) Base Case  B  SW-GW-1  3  Discharge (m /d)  Total Discharge along the Stream Channel 140 120 100 80 60 40 20 0 0  50  100  150  200  250 300 X (m)  Base Case  C  350  SW-GW-1  Velocity (m/d)  Mean Flow Velocity (V x ) along the Stream Channel 1400 1200 1000 800 600 400 200 0 0  50  100  150  200  250 X (m)  Base Case  300  350  400  450  500  SW-GW-1  Figure 5.18 Comparison of stream water flow parameters for the first realizations of Scenario SW-GW-1 (higher Kuc) and the Base Case. A Water depth, B total discharge, and C mean velocity (Vx).  174  A Plume in stream  B  Concentration  C  Plume in SW & GW  Concentration  D  Source Zone Concentration  Solute entering SW  Concentration  Figure 5.19 Plume in the SW-GW system (concentrations above detection threshold) at time of detection at compliance boundary in the first realization of Scenario SW-GW-1 (higher Kuc). A Plume in stream water (2D), B plumes in stream water and groundwater (3D), C plume in groundwater (cross-section Y=70m, source zone), and D solute entering stream water through streambed (cross-section Y=100m, domain centerline). Concentration (g/L), Dimension (m).  175  Solute Concentrations in Stream Water (SW) and Streambed (GW)  Concentration (g/l)  8.0E-03 6.0E-03  Solute in GW entering SW  4.0E-03  SW  GW  2.0E-03 0.0E+00 0  50  100  150  200  250 X (m)  300  350  400  450  500  Detection threshold: 1.0×10-3g/l  Figure 5.20 Concentrations along the stream centreline (in stream water and streambed) at the time of detection at the compliance boundary in the first realization of Scenario SW-GW1 (higher Kuc)  176  A Concentration (g/l)  Concentration Profiles at Cross-Sections in Stream Channel 1.20E-03 1.00E-03 Solute in GW entering SW  8.00E-04  X=337m  6.00E-04 96  97  98  99  100  101  102  103  104  B  Concentration (g/l)  Y (m)  2.20E-03 1.80E-03 1.40E-03  X=350m  Solute in GW entering SW  1.00E-03 96  97  98  99  100  101  102  103  104  C  Concentration (g/l)  Y (m)  2.60E-03 2.40E-03 X=400m  2.20E-03  Solute in GW entering SW  2.00E-03 96  97  98  99  100  101  102  103  104  D  Concentration (g/l)  Y (m)  1.62E-03 1.58E-03 X=450m  1.54E-03 1.50E-03 96  97  98  99  100  101  102  103  104  E  Concentration (g/l)  Y (m)  1.03E-03 1.00E-03 X=500m  9.70E-04 9.40E-04 96  97  98  99  100 Y (m)  101  102  103  104  (Detection threshold: 1.0×10-3g/l)  Figure 5.21 Concentration profiles in stream water at the time of detection at the compliance boundary in the first realization of Scenario SW-GW-1 (higher Kuc). Cross-sections at A X=337m, B X=350m, C X=400m, D X=450m, and E X=500m.  177  Detection Probabilities at Stream Outlet 1.0  Probability  0.8 0.6 0.4 0.2 0.0 96  97  98  99  100  101  102  103  104  Y (m)  Figure 5.22 Detection probabilities at the stream outlet in Scenario SW-GW-1 (higher Kuc)  A  B  Figure 5.23 Detection probabilities of groundwater monitoring wells in Scenario SW-GW-1 (higher Kuc). A Unconfined porous aquifer, B fractured bedrock. 178  A  Depth (m)  Water Depth along the Stream Channel 0.014 0.012 0.010 0.008 0.006 0.004 0.002 0.000 0  50  100  150  200  250 X (m)  300  Base Case  B  350  400  450  500  400  450  500  SW-GW-2a  Total Discharge along the Stream Channel  3  Discharge (m /d)  120 100 80 60 40 20 0 0  50  100  150  200  250  300  350  X (m) Base Case  C  SW-GW-2a  Velocity (m/d)  Mean Flow Velocity (V x ) along the Stream Channel 1400 1200 1000 800 600 400 200 0 0  50  100  150  200  250 X (m)  Base Case  300  350  400  450  500  SW-GW-2a  Figure 5.24 Comparison of stream water flow parameters for the first realizations of Scenario SW-GW-2a (larger horizontal fracture apertures) and the Base Case. A Water depth, B total discharge, and C mean velocity (Vx).  179  A Plume in stream  B  Concentration  C  Source Zone  Plume in SW & GW  Concentration  D  Concentration  Solute entering SW  Concentration  Figure 5.25 Plume in the SW-GW system (concentrations above detection threshold) at time of detection at compliance boundary in the first realization of Scenario SW-GW-2a (larger horizontal fracture apertures). A Plume in stream water (2D), B plumes in stream water and groundwater (3D), C plume in groundwater (cross-section Y=70m, source zone), and D solute entering stream water through streambed (cross-section Y=100m, domain centerline). Concentration (g/L), Dimension (m).  180  Solute Concentrations in Stream Water (SW) and Streambed (GW)  Concentration (g/l)  4.0E-03 3.0E-03  Solute in GW entering SW  Solute in GW entering SW  SW  2.0E-03  GW  1.0E-03 0.0E+00 0  50  100  150  200  250 X (m)  300  350  400  450  500  Detection threshold: 1.0×10-3g/l  Figure 5.26 Concentrations along the stream centreline (in stream water and streambed) at the time of detection at the compliance boundary in the first realization of Scenario SW-GW2a (larger horizontal fracture apertures)  181  Concentration (g/l)  A  Concentration Profiles at Cross-Sections in Stream Channel 1.50E-04 1.30E-04 1.10E-04  X=160m  Solute in GW entering SW  9.00E-05 96  97  98  99  100  101  102  103  104  Y (m)  Concentration (g/l)  B  1.00E-03 8.00E-04  X=200m  6.00E-04  Solute in GW entering SW  4.00E-04 96  97  98  99  100  101  102  103  104  C  Concentration (g/l)  Y (m)  1.36E-03 1.32E-03 X=300m  1.28E-03 1.24E-03 96  97  98  99  100  101  102  103  104  D  Concentration (g/l)  Y (m)  Detection Threshold  1.03E-03 1.01E-03  X=400m  9.90E-04  Solute in GW entering SW  9.70E-04 96  97  98  99  100  101  102  103  104  E  Concentration (g/l)  Y (m)  7.85E-04 7.65E-04 X=500m  7.45E-04 7.25E-04 96  97  98  99  100 Y (m)  101  102  103  104  (Detection threshold: 1.0×10-3g/l)  Figure 5.27 Concentration profiles in stream water at the time of detection at the compliance boundary in the first realization of Scenario SW-GW-2a (larger horizontal fracture apertures). Cross-sections at A X=160m, B X=200m, C X=300m, D X=400m, and E X=500m.  182  Detection Probabilities at Stream Outlet 1.0  Probability  0.8 0.6 0.4 0.2 0.0 96  97  98  99  100  101  102  103  104  Y (m)  Figure 5.28 Detection probabilities at the stream outlet in Scenario SW-GW-2a (larger horizontal fracture apertures)  A  B  Figure 5.29 Detection probabilities in groundwater monitoring wells in Scenario SW-GW-2a (larger horizontal fracture apertures). A Unconfined porous aquifer, B fractured bedrock. 183  A  Depth (m)  Water Depth along the Stream Channel 0.014 0.012 0.010 0.008 0.006 0.004 0.002 0.000 0  50  100  150  200  250 X (m)  Base Case  B  300  350  400  450  500  400  450  500  SW-GW-2b  Total Discharge along the Stream Channel  3  Discharge (m /d)  120 100 80 60 40 20 0 0  50  100  150  200  250 300 X (m)  Base Case  C  350  SW-GW-2b  Velocity (m/d)  Mean Flow Velocity (V x ) along the Stream Channel 1400 1200 1000 800 600 400 200 0 0  50  100  150  200  250 X (m)  Base Case  300  350  400  450  500  SW-GW-2b  Figure 5.30 Comparison of stream water flow parameters for the first realizations of Scenario SW-GW-2b (larger vertical fracture apertures) and the Base Case. A Water depth, B total discharge, and C mean velocity (Vx).  184  A  B Plume in stream  Concentration  C  Plume in SW & GW  Concentration  D  Source Zone  Concentration  Solute entering SW  Concentration  Figure 5.31 Plume in the SW-GW system (concentrations above detection threshold) at time of detection at compliance boundary in the first realization of Scenario SW-GW-2b (larger vertical fracture apertures). A Plume in stream water (2D), B plumes in stream water and groundwater (3D), C plume in groundwater (cross-section Y=70m, source zone), and D solute entering stream water through streambed (cross-section Y=100m, domain centerline). Concentration (g/L), Dimension (m).  185  Solute Concentrations in Stream Water (SW) and Streambed (GW)  Concentration (g/l)  1.5E-02 Solute in GW entering SW  1.0E-02  SW GW  5.0E-03  0.0E+00 0  50  100  150  200  250  X (m)  300  350  400  450  500  Detection threshold: 1.0×10-3g/l  Figure 5.32 Concentrations along the stream centreline (in stream water and streambed) at the time of detection at the compliance boundary in the first realization of Scenario SW-GW2b (larger vertical fracture apertures)  186  Concentration (g/l)  A  Concentration Profiles at Cross-Sections in Stream Channel 5.75E-03 4.50E-03 3.25E-03  X=100m  Solute in GW entering SW  2.00E-03 96  97  98  99  100  101  102  103  104  B  Concentration (g/l)  Y (m)  3.25E-03 3.05E-03 2.85E-03  X=200m  Solute in GW entering SW  2.65E-03 96  97  98  99  100  101  102  103  104  C  Concentration (g/l)  Y (m)  2.12E-03 2.09E-03  X=300m  2.07E-03 2.04E-03 96  97  98  99  100  101  102  103  104  D  Concentration (g/l)  Y (m)  1.40E-03 1.38E-03 X=400m  1.36E-03 1.34E-03 96  97  98  99  100  101  102  103  104  Y (m)  Concentration (g/l)  E  1.02E-03 1.00E-03 X=500m  9.80E-04 9.60E-04 96  97  98  99  100 Y (m)  101  102  103  104  (Detection threshold: 1.0×10-3g/l)  Figure 5.33 Concentration profiles in stream water at the time of detection at the compliance boundary in the first realization of Scenario SW-GW-2b (larger vertical fracture apertures). Cross-sections at A X=100m, B X=200m, C X=300m, D X=400m, and E X=500m.  187  A  B  Figure 5.34 Detection probabilities in groundwater monitoring wells in Scenario SW-GW-2b (larger vertical fracture apertures). A Unconfined porous aquifer, B fractured bedrock.  188  Concentration (g/l)  A  Concentrations in Stream Water at Cross-Sections 1.55E-02 1.25E-02  Base Case  9.50E-03 6.50E-03 96  97  98  99  Concentration (g/l) Concentration (g/l) Concentration (g/l)  Concentration (g/l)  103  104 X=100m  Base Case  3.55E-03  SW-GW-3a  3.40E-03 96  97  98  99  100  101  102  103  104 X=200m  Y (m)  2.06E-03 2.01E-03  Base Case  1.96E-03  SW-GW-3a  1.91E-03 96  97  98  99  100  101  102  103  104  X=300m  Y (m)  1.38E-03 1.35E-03  Base Case  1.32E-03  SW-GW-3a  1.29E-03 96  97  98  99  100  101  102  103  104  Y (m)  1.5% of Increase in Average  E  102  3.70E-03  1.9% of Increase in Average  D  101  3.85E-03  3.1% of Increase in Average  C  100 Y (m)  6.2% of Increase in Average  B  SW-GW-3a  Solute in GW entering SW  X=400m  1.02E-03 1.00E-03  Base Case  9.80E-04  SW-GW-3a  9.60E-04 96  97  (Detection threshold: 1.0×10-3g/l)  98  99  100  101  102  103  104  X=500m  Y (m)  Figure 5.35 Comparison of the concentrations profiles in stream water at the time of detection at the compliance boundary in the first realization of Scenario SW-GW-3a (smaller αl in stream water) and the Base Case. Cross-sections at A X=100m, B X=200m, C X=300m, D X=400m, and E X=500m. 189  A Concentration (g/l)  Concentrations in Stream Water at Cross-Sections 1.60E-02 1.10E-02  Base Case  6.00E-03  Solute in GW entering SW  SW-GW-3b  1.00E-03 96  97  98  99  100  101  102  103  104 X=100m  Y (m)  Concentration (g/l)  B  3.75E-03 3.60E-03  Base Case  3.45E-03  SW-GW-3b  3.30E-03 96  97  98  99  100  101  102  103  104  C  Concentration (g/l)  Y (m)  X=200m  2.00E-03 1.96E-03  Base Case  1.92E-03  SW-GW-3b  1.88E-03 96  97  98  99  100  101  102  103  104  Y (m)  Concentration (g/l)  D  X=300m  1.35E-03 1.32E-03  Base Case  1.30E-03  SW-GW-3b  1.27E-03 96  97  98  99  100  101  102  103  104 X=400m  E  Concentration (g/l)  Y (m)  1.01E-03 9.90E-04  Base Case  9.70E-04  SW-GW-3b  9.50E-04 96  97  (Detection threshold: 1.0×10-3g/l)  98  99  100 Y (m)  101  102  103  104 X=500m  Figure 5.36 Comparison of the concentrations profiles in stream water at the time of detection at the compliance boundary in the first realization of Scenario SW-GW-3b (smaller αt in stream water) and the Base Case. Cross-sections at A X=100m, B X=200m, C X=300m, D X=400m, and E X=500m.  190  Detection Probabilities at Stream Outlet 1.0  Probability  0.8 0.6 0.4 0.2 0.0 96  97  98  99  100  101  102  103  104  Y (m)  Figure 5.37 Detection probabilities at the stream outlet in Scenario SW-GW-3b (smaller αt in stream water)  A  B Plume in stream  Concentration  C  Plume in SW & GW  Concentration  Solute entering SW  Concentration  D Source Zone  Concentration  Figure 5.38 Plume in the SW-GW system (concentrations above detection threshold) at time of detection at compliance boundary in the first realization of Scenario SW-GW-4 (lower detection threshold 1.0×10-4g/L). A Plume in stream water (2D), B plumes in stream water and groundwater (3D), C plume in groundwater (cross-section Y=70m, source zone), and D solute entering stream water through streambed (cross-section Y=100m, domain centerline). Concentration (g/L), Dimension (m). 191  Solute Concentrations in Stream Water (SW) and Streambed (GW) 1.4E-03  Concentration (g/l)  1.2E-03 Solute in GW entering SW  1.0E-03  SW GW  8.0E-04 6.0E-04 4.0E-04 2.0E-04 0.0E+00 0  100  200  300  X (m)  400  500 Detection threshold: 1.0×10-4g/l  Figure 5.39 Concentrations along the stream centreline (in stream water and streambed) at the time of detection at the compliance boundary in the first realization of Scenario SW-GW4 (lower detection threshold 1.0×10-4g/L)  192  Concentration (g/l)  A  Concentrations in Stream Water at Cross-Sections 1.40E-03 1.15E-03 X=100m  9.00E-04  Solute in GW entering SW  6.50E-04 96  97  98  99  100  101  102  103  104  B  Concentration (g/l)  Y (m)  3.75E-04 3.65E-04 X=200m  3.55E-04  Solute in GW entering SW  3.45E-04 96  97  98  99  100  101  102  103  104  C  Concentration (g/l)  Y (m)  2.03E-04 2.00E-04 X=300m  1.97E-04 1.94E-04 96  97  98  99  100  101  102  103  104  D  Concentration (g/l)  Y (m)  1.36E-04 1.34E-04 X=400m  1.32E-04 1.30E-04 96  97  98  99  100  101  102  103  104  E  Concentration (g/l)  Y (m)  1.01E-04 1.00E-04 X=500m  9.90E-05 9.80E-05 96  97  98  99  100 Y (m)  101  102  103  104  Figure 5.40 Concentration profiles in stream water at the time of detection at the compliance boundary in the first realization of Scenario SW-GW-4 (lower detection threshold 1.0×104 g/L). Cross-sections at A X=100m, B X=200m, C X=300m, D X=400m, and E X=500m.  193  Detection Probabilities at Stream Outlet 1.0  Probability  0.8 0.6 0.4 0.2 0.0 96  97  98  99  100  101  102  103  104  Y (m)  Figure 5.41 Detection probabilities at the stream outlet in Scenario SW-GW-4 (lower detection threshold 1.0×10-4g/L)  A  B  Figure 5.42 Detection probabilities in groundwater monitoring wells in Scenario SW-GW-4 (lower detection threshold 1.0×10-4g/L). A Unconfined porous aquifer, B fractured bedrock. 194  A  B  Plume in stream  Concentration  C  Plume in SW & GW  Concentration  D Solute entering SW  Source Zone Concentration  Concentration  Figure 5.43 Plume in the SW-GW system (concentrations above detection threshold) at time of detection at compliance boundary in the first realization of Scenario SW-GW-5 (source zone on bedrock). A Plume in stream water (2D), B plumes in stream water and groundwater (3D), C plume in groundwater (cross-section Y=70m, source zone), and D solute entering stream water through streambed (cross-section Y=100m, domain centerline). Concentration (g/L), Dimension (m).  195  Concentrations in Stream Water (SW) & Streambed (GW) with Difference Source Locations  Concentration (g/l)  1.8E-02 1.5E-02  Solute in GW entering SW in SW-GW-5  1.2E-02 9.0E-03 6.0E-03  Solute in GW entering SW in Base Case  3.0E-03 0.0E+00 0  SW-GW-5 (GW)  100  200  SW-GW-5 (SW)  X (m)  300  400  500  Detection threshold: 1.0×10-3g/l  Base Case (GW)  Base Case (SW)  Figure 5.44 Comparison of locations for solute entering stream water in Base Case and Scenario SW-GW-5 (source zone on bedrock) at the time of detection at the compliance boundary  196  B  Concentration (g/l)  Concentration (g/l)  A  Concentrations in Stream Water at Cross-Sections 1.90E-04 1.50E-04 X=100m  1.10E-04  Solute in GW entering SW  7.00E-05 96  97  98  99  100 Y (m)  101  102  103  104  3.50E-03 3.25E-03 3.00E-03  X=200m  Solute in GW entering SW  2.75E-03 96  97  98  99  100  101  102  103  104  C  Concentration (g/l)  Y (m)  2.01E-03 1.99E-03  X=300m  1.97E-03 1.95E-03 96  97  98  99  100  101  102  103  104  D  Concentration (g/l)  Y (m)  1.35E-03 1.33E-03  X=400m  1.31E-03 1.29E-03 96  97  98  99  100  101  102  103  104  E  Concentration (g/l)  Y (m)  1.02E-03 1.00E-03 X=500m  9.80E-04 9.60E-04 96  97  98  99  100 Y (m)  101  102  103  104  (Detection threshold: 1.0×10-3g/l)  Figure 5.45 Concentration profiles in stream water at the time of detection at the compliance boundary in the first realization of Scenario SW-GW-5 (source zone on bedrock). Crosssections at A X=100m, B X=200m, C X=300m, D X=400m, and E X=500m.  197  A  B  Figure 5.46 Detection probabilities in groundwater monitoring wells in Scenario SW-GW-5 (source zone on bedrock). A Unconfined porous aquifer, B fractured bedrock.  198  6 CONCLUSIONS  6.1 Conclusions: groundwater monitoring  The factors influencing the design of a groundwater monitoring network are examined for a catchment where a permeable unconfined aquifer overlies sparsely fractured bedrock. A contaminant source zone with a constant concentration is situated close to the ground surface. The unconfined aquifer is modeled as a homogeneous and isotropic porous medium with deterministic properties. The fractured bedrock is modeled as a lower permeability rock mass, but with a stochastic fracture network composed of three orthogonal fracture sets. Steady-state groundwater flow and non-reactive contaminant transport is simulated to represent a contaminant plume in the groundwater system, using the Monte Carlo approach to characterize the variability in flow and transport associated with the stochastic fracture networks. The probability of contaminant plume detection is defined in terms of the likelihood of detecting contamination in a performance monitoring network prior to its detection at a downstream compliance boundary. Detection is measured relative to a preset contaminant concentration threshold applied at both the compliance boundary and monitoring wells. The factors examined include the hydraulic conductivity of the unconfined aquifer, the fracture connectivity (fracture densities) and apertures in the bedrock, the rock matrix permeability, the detection threshold and the contaminant source zone location. The groundwater simulation results indicate: (1) The unconfined aquifer and the fractured bedrock are best viewed as an integrated hydrologic system for monitoring the contamination in the subsurface, and an understanding 199  of the large-scale hydraulic conductivities of both units is required in the design of a groundwater monitoring network. Monitoring the bedrock pathways may be advantageous even if the network of open fractures in the bedrock appears sparse. (2) The ratios of the hydraulic conductivity of the unconfined aquifer to the bulk hydraulic conductivity of the bedrock, coupled with contaminant detection threshold and source zone location play the major role in determining detection probabilities in the two units. With the source zone coinciding with the water table, and for the given detection threshold, it was observed that when the ratio of the hydraulic conductivity of the unconfined aquifer to that of the bedrock (Kuc/Kx) exceeds 10, detection at the compliance boundary occurs in the unconfined aquifer. When this ratio is less than 1, detection at the compliance boundary occurs in the bedrock. When the ratios are between 1 and 10, detection at the compliance boundary can be either in the unconfined aquifer or in the bedrock, or in both. (3) The hydraulic conductivity of the unconfined aquifer is significant in influencing plume detection in the groundwater system. The more permeable the unconfined aquifer, the more likely the detection of the contaminant occurs in the performance monitoring wells in the unconfined aquifer, and the less likely the detection occurs in the bedrock prior to detection at the compliance boundary. The influence of the fractured bedrock becomes less significant when the unconfined aquifer has a higher permeability. (4) For the case where transport in the bedrock is important in determining detection probabilities, fracture network connectivity (fracture densities) and fracture apertures have significant impacts on probability of detection of the contaminant in the groundwater system. Lower fracture network connectivity (smaller densities) will increase the tortuosity of the contaminant transport pathways in the bedrock, and cause more lateral spreading of the contaminant (detection of the contaminant in a wider area) in the bedrock prior to detection at 200  the compliance boundary. Larger fracture apertures lead to an increase in detection probabilities in the bedrock. The aperture of the horizontal fracture set is more important than that for the vertical fracture sets in determining if detection will occur in the bedrock at the compliance boundary. Matrix permeability has a minor influence on detection probabilities in the bedrock. Due to the fracture networks, groundwater flow field and plume geometries are irregular in both the bedrock and the unconfined aquifer. (5) The detection threshold is a key factor in influencing detection probabilities in both the unconfined aquifer and the fractured bedrock. A lower detection threshold will increase the probabilities in detecting the contaminant in the bedrock at the compliance boundary, and also increase detection probabilities in the performance monitoring wells in both subsurface units, especially in the bedrock, prior to arrival at the compliance boundary. (6) The source zone location is important for monitoring the contamination. Detection probabilities in the bedrock at the compliance boundary increase as the source zone is closer to the bedrock. For monitoring near the source zone, detection probabilities in the performance monitoring wells depend upon downgradient distance from source before the contaminant reaches the bedrock with concentrations above the detection threshold. Monitoring too close to the source may fail to detect the contaminant with the wells in the bedrock.  6.2 Conclusions: integrated surface water and groundwater monitoring  The factors influencing the design of an integrated surface water and groundwater (SWGW) monitoring network are examined in a headwater stream catchment where a permeable unconfined aquifer overlies sparsely fractured bedrock. The subsurface domain has the same 201  features and parameters as those in the groundwater monitoring analysis, except that the slope gradient of the land surface is larger, the unconfined aquifer is thicker, and the permeability of the rock matrix is two orders of magnitude higher. The surface domain consists of a V-shaped symmetric overland zone and a linear channel with a constant width. The contaminant source with a constant concentration is situated in the saturated zone of the unconfined aquifer. The integrated SW-GW monitoring network is composed of stream water observation points at the stream outlet (compliance boundary) and along the stream channel centerline, and performance monitoring wells in both subsurface units. The methods for calculation of detection probabilities and for evaluation of the monitoring network are the same as those in the groundwater analysis. Coupled steady-state SW-GW flow and nonreactive contaminant transport is simulated to represent a contaminant plume in the integrated SW-GW system, using the Monte Carlo approach to characterize the uncertainty of surface and subsurface flow and transport in association with the stochastic fracture networks. The factors examined include the dispersivities in the stream water, the hydraulic conductivity of the unconfined aquifer, the fracture apertures in the bedrock, the detection threshold and the contaminant source zone location. The coupled surface water and groundwater simulation results indicate: (1) Monitoring contamination in the stream water and in the groundwater system should be integrated, and understanding the interactions of water flow and solute contaminant between the surface and the subsurface at the catchment scale is essential for design of an integrated surface water and groundwater monitoring network. In particular, an understanding of the factors influencing the solute concentration distribution in the stream water, including the fracture networks in the bedrock, is important in the design of the monitoring networks. 202  (2) The fracture networks create an irregular flow field in the bedrock and also in the overlying unconfined aquifer. The shallow stream flow is non-uniform laminar with a nonsymmetric distribution of stream flow velocities due to the influence of the fracture networks. For the system geometry that is considered, the stream gains water from groundwater discharge along its entire reach, and solute mass enters the stream with groundwater discharge through the streambed. Dilution of solute concentrations by stream flow augmentation along the channel is observed to have a very significant influence on solute concentrations and detection probabilities in the stream water. Even if the contaminant enters the stream, it may still be first detected in unconfined aquifer or in the bedrock at the compliance boundary in some cases. (3) The hydraulic conductivity of the unconfined aquifer is a key factor in influencing flow, solute transport and detection probabilities in the integrated surface and subsurface water system. A higher permeability of the unconfined aquifer results in a shorter stream and a deeper water table in the headwater basin. A higher permeability of the unconfined aquifer leads to an increase in detection probabilities in the aquifer and also in the underlying bedrock, but a decrease in detection probabilities in the stream. A higher permeability of the unconfined aquifer will reduce the influence of fracture networks on flow, solute transport and detection probabilities in the stream. (4) For the case where the bedrock pathways are important, fracture apertures, especially for the horizontal fracture set, have a significant influence on both the surface water and groundwater flow and solute transport, and the SW-GW interactions. Larger fracture apertures lead to an increase in detection probabilities in the bedrock, but a decrease in detection probabilities in the stream (in case with larger aperture for horizontal fracture  203  set). The larger the fracture apertures, the more influence the fractured bedrock will have on the flow, solute transport and detection probabilities in the stream. (5) Dispersivities in the stream water are important parameters in calculation of solute concentrations and detection probabilities in the stream water, but they have an insignificant effect on detection probabilities in the groundwater system. A smaller longitudinal dispersivity results in slower migration of the plume front and higher detection probabilities along the stream. A smaller transverse dispersivity results in less lateral spreading of the plume and lower detection probabilities along the stream banks. In comparison to dilution due to the stream flow augmentation, dispersion has much less influence on the solute distribution and detection probabilities in the stream. Because of the effects of dilution, sampling both the stream water and the streambed groundwater will increase probabilities of detection of the contaminant in the stream. (6) The detection threshold is a key factor in influencing integrated SW-GW monitoring. A lower detection threshold generally leads to higher detection probabilities in both the stream water and the groundwater (especially in the bedrock), but detection probabilities in the stream water are more sensitive to the lower detection threshold. (7) The source zone location has a significant impact on solute transport and detection probabilities in both stream water and groundwater. The source zone location has a strong effect on the locations where solute enters the stream from the subsurface, which will influence detection probabilities and detection time at the stream outlet compliance boundary. Sampling the stream water on the source zone side in the upper reach will generally have a greater chance in detection of the contaminant in the stream. The probabilities of detection of the contaminant in the bedrock are also affected by the source  204  zone location. The closer the source zone is to the bedrock, the more likely the contaminant will be detected in the bedrock. In summary, the implications of the findings for hydrogeological practice include: (1) stream water monitoring and groundwater monitoring should be integrated for monitoring the contamination in a hydrological system; (2) An understanding of the SW-GW interactions in a catchment scale is essential for design of an integrated SW-GW monitoring network; (3) Fracture networks in the bedrock are important in influencing the solute transport and plume detection in both the surface and subsurface systems, even if they appear sparse; (4) The hydraulic conductivity of the unconfined aquifer is an important factor to be considered in design of an integrated SW-GW monitoring network (relative to the bulk hydraulic conductivity of the fractured bedrock), and it has a significant impact on both the stream water flow and groundwater flow; (5) Stream water monitoring is more sensitive to a lower detection threshold than groundwater monitoring; (6) Source zone location or the distance of the source zone to the stream is a key factor to be considered in influencing solute concentrations in the stream water, and sampling on the source zone side of the stream generally has more likelihood to detect the plume in the stream; (7) Due to significant dilution of solute concentrations by stream flow, sampling both stream water and streambed groundwater is advantageous in increasing probability of detection of the plume in the stream; (8) Dispersion has a much smaller effect relative to dilution by flow augmentation in influencing the solute concentrations in the stream water. Finally, it needs to be emphasized that in this study, both the advection and the effect of dispersion and diffusion are considered in evaluation of both the groundwater monitoring network and the integrated SW-GW monitoring network. The effects of dispersion and diffusion were typically ignored in the previous studies of groundwater monitoring (Jardine 205  et al., 1996, Storck et al., 1997). Also it needs to be noted that the conclusions are made with the assumptions of no vertical variation of flow velocities and solute concentrations in the stream, no change of solute concentrations in the stream due to rainfall events, and no borehole dilution in a performance monitoring well, and the performance monitoring wells provide perfect vertical resolution in the unconfined aquifer and the bedrock unit.  6.3 Suggestions for future research  This study evaluates the networks for groundwater monitoring and integrated SW-GW monitoring by examining detection probabilities in groundwater monitoring wells and at stream water observation points, and identifies the key factors influencing design of the monitoring networks in surface water and groundwater. However, there are certain limitations existing in the study due to the complexity of the issues investigated, the time and availability of the computer resources. The following factors are not considered and may have an impact on the results of this study, and should be investigated in the future: (1) Heterogeneity and stochasticity of the unconfined aquifer, (2) Heterogeneity of the bedrock (fractures and matrix) with a larger bedrock zone, (3) Fracture systems dominated by non-orthogonal and correlated fractures, (4) Undulating bedrock topography and a weathered zone at the bedrock-drift interface, (5) Complexity in stream channel geometry and streambed permeability, (6) Variability of vertical flow velocities and vertical solute mixing in stream, (7) Turbulent or transient (seasonal) flow assumed in the stream.  206  Therefore, it is suggested that the future study should examine whether and how these factors will influence solute concentration distribution and plume detection in groundwater and stream water, and design of groundwater monitoring or integrated SW-GW monitoring networks. In addition, the cost of sampling and well drilling was not considered in evaluation of the monitoring networks in this study. It is recommended that future research on evaluation of groundwater monitoring or integrated SW-GW monitoring networks in such a system be carried out with the aid of optimization analysis with the objectives such as minimizing cost of monitoring and maximizing detection probabilities to identify monitoring tradeoffs. Finally, with more powerful computation resources available in the future, it will be also interesting to investigate the factors influencing the integrated SW-GW monitoring network in a large and irregular watershed domain with a large catchment area and a variable surface (not a V-shape symmetric overland) for generation of the surface and subsurface flow in a porous and fractured system.  207  REFERENCES  Abdul A.S., 1985. Experimental and numerical studies of the effect of the capillary fringe on streamflow generation, Ph.D. Thesis, University of Waterloo, Waterloo, Ontario, Canada, pp. 1-210.  Adler P. M. and J. F. Thovert, 1999. Fractures and fracture networks, in Theory and Application of Transport in Porous Media, Volume 15, edited by J. Bear, Klumwer Academic Publishers, Dordrecht, the Netherlands, pp. 1-429.  Akan, A. O., 2006. Open Channel Hydraulics, Elsevier and Butterworth-Heinemann, Burlington, Massachusetts, ISBN: 0-7506-66857-1, pp.364.  Alexander R. B., J. R. Slack, A. S. Ludtke, K. K. Fitzgerald and T. L. Schertz, 1998. Data from selected U.S. Geological Survey national stream water quality monitoring networks, Water Resources Research, 34(9): 2401-2405.  American Society of Civil Engineers (ASCE), 2003. Long-Term Groundwater Monitoring – the State of the Art, Library of Congress Catalog Card No: 2003045245, ISBN 0-78440678-2, American Society of Civil Engineers, Reston, VA, USA, pp. 1-103.  Bauer, P., T. Gumbricht and W. Kinzelbach, 2006. A regional coupled surface water/groundwater model of the Okavango Delta, Botswana, Water Resources Research, 42, W04403: 1-15. 208  Berkowitz B., 2002. Characterizing flow and transport in fractured geological media: A review. Advances in Water Resources, 25:861-884.  Berkowitz B., O. Bour, P. Davy and N. Odling, 2000. Scaling of fracture connectivity in geological formations, Geophysical Research Letters, 27(14), 2061-2064.  Bodin J., F. Delay and G. de Marsily, 2003a. Solute transport in a single fracture with negligible matrix permeability: 1. fundamental mechanisms, Hydrogeology Journal, 11:418-433.  Bodin J., F. Delay and G. de Marsily, 2003b. Solute transport in a single fracture with negligible matrix permeability: 2. mathematical formalism, Hydrogeology Journal, 11:434-454.  Bonnet E., O. Bour, N. E. Odling, P. Davy, I. Main, P. Cowie and B. Berkowitz, 2001. Scaling of fracture systems in geological media, Review of Geophysics, 39(3), 347-383.  Chapman, S. W., B. L. Parker, J. A. Cherry, R. Aravena and D. Hunkeler, 2007. Groundwater–surface water interaction and its role on TCE groundwater plume attenuation, Journal of Contaminant Hydrology, 91: 203–232.  Chiles, J. P. and G. de Marsily, 1993. Stochastic models of fracture systems and their use in flow and transport modeling, in Flow and Contaminant Transport in Fractured Rock, 209  edited by J. Bear, C. F. Tsang, and G. de Marsily, Academic Press, San Diego, CA, USA, pp. 169-236.  Cieniawski, S. E., J. W. Eheart and S. Ranjithan, 1995. Using genetic algorithms to solve a multiobjective groundwater monitoring problem, Water Resources Research, 31(2), 399409.  Claessen, F. A. M., 1997. Comparing monitoring of surface and groundwater systems, European Water Pollution Control, 7(4), 27-35.  Conant J. B., J. Cherry and R. W. Gillham, 2004. A PCE groundwater plume discharging to a river: influence of the streambed and near-river zone on contaminant distributions, Journal of Contaminant Hydrology, 73: 249-279.  de Dreuzy, J. R., C. Darcel, P. Davy and O. Bour, 2004. Influence of spatial correlation of fracture centers on the permeability of two-dimensional fracture networks following a power law length distribution, Water Resources Research, Vol.40, w01502, doi:10.1029/2003wr002260.  di Giammarco, P., E. Todini, and P. Lamberti, 1996. A conservative finite element approach to overland flow: the control volume finite element formulation, Journal of Hydrology, 175, 267-291.  210  Dingman, S. L. 1994. Physical Hydrology, Macmillan College Publishing Company, New York, USA, pp575  Dixon, W. and B. Chiswell, 1996. Review of aquatic monitoring program design, Water Research, 30(9), 1935-1948.  Dougherty, D. E. and R. A. Marryott, 1991. Optimal groundwater management 1. simulated annealing, Water Resources Research, 27(10), 2493-2508.  Fischer H. G., E. J. List, R. C. Y. Koh, J. Imberger and N. H. Brooks, 1979. Mixing in inland and coastal waters, Academic Press, New York, pp. 40-147.  Hudak, P. F., 2002. Efficiency comparison of graphical approaches for designing contaminant detection networks in groundwater, Water Resources Research, 38(12), 18:1-5.  Hudak, P. F. and H. A. Loaiciga, 1993. An optimization method for monitoring network design in multilayered groundwater flow systems, Water Resources Research, 29(8), 2835-2845.  Hudak, P. F., H. A. Loaiciga and M. A. Marino, 1995. Regional-scale groundwater quality monitoring via integer programming, Journal of Hydrology, 164:153-170.  211  James, B. R. and S. M. Gorelick, 1994. When enough is enough: The worth of monitoring data in aquifer remediation design, Water Resources Research, 30(12), 3499-3513.  Jardine, K., L. Smith and T. Clemo, 1996. Monitoring networks in fractured rocks: a decision analysis approach, Ground Water, 34(3), 504-518.  Jones, J. P., E. A. Sudicky, A. E. Brookfield and Y.-J. Park, 2006. An assessment of the tracer-based approach to quantifying groundwater contributions to streamflow, Water Resources Research, 42, W02407: 1-15.  Jones, J. P., E. A. Sudicky and R. G. McLaren, 2008. Application of a fully-integrated surface-subsurface flow model at the watershed-scale: A case study, Water Resources Research, 44, W03407: 1-13.  Kollet, S. J. and R. M. Maxwell, 2006. Integrated surface–groundwater flow modeling: A free-surface overland flow boundary condition in a parallel groundwater flow model, Advances in Water Resources, 29, 945–958.  Kristensen, P. and J. Bøgestrand, 1996. European topic centre on inland waters: Surface water quality monitoring, National Environmental Research Institute, Denmark, ISBN 92-9167-001-4, pp. 82.  212  Li Q., A. J. A. Unger, E.A. Sudicky, D. Kassenaar, E. J. Wexler and S. Shikaze, 2008. Simulating the multi-seasonal response of a large-scale watershed with a 3D physicallybased hydrologic model, Journal of Hydrology, 357:317-336.  Loaiciga, H. A., R. J. Charbeneau, L. G. Everett, G. E. Fogg, B. F. Hobbs, and S. Roubani, 1992. Review of groundwater-water quality monitoring network design, Journal of Hydraulic Engineering, 118(1), 11-37.  MacQuarrie, K. T. B. and K. U. Mayer, 2005. Reactive transport modeling in fractured rock: A state-of-the science review, Earth-Science Review, 72:189-227.  Marryott, R. A., D. E. Dougherty and R. L. Stollar, 1993. Optimal groundwater management 2. application of simulated annealing to a field-scale contamination site, Water Resources Research, 29(4), 847-860.  Massmann, J. and R. A. Freeze, 1987a. Groundwater contamination from waste management sites: the interaction between risk-based engineering design and regulatory policy 1. methodology. Water Resource Research, 23(2), 351-367.  Massmann, J. and R. A. Freeze, 1987b. Groundwater contamination from waste management sites: the interaction between risk-based engineering design and regulatory policy 2. results. Water Resource Research, 23(2), 368-380.  213  Meyer, P. D. and E. D. Brill, 1988. A method for locating wells in a groundwater monitoring network under conditions of uncertainty, Water Resources Research, 24(8), 1277-1282.  Meyer, P. D., A. J. Valocchi and J. W. Eheart, 1994. Monitoring network design to provide initial detection of groundwater contamination, Water Resources Research, 30(9), 26472659.  Morisawa, S. and Y. Inoue, 1991. Optimum allocation of monitoring wells around a solidwaste landfill site using precursor indicators and fuzzy utility functions, Journal of Contaminant Hydrology, 7:337-370.  Neuman, S. P., 2005. Trends, prospects and challenges in quantifying flow and transport through fractured rocks, Hydrogeology Journal, 13:124-147.  Oxtobee, J. P. A. and K. Novakowski, 2002. A field investigation of groundwater/surface water interaction in a fractured bedrock environment, Journal of Hydrology, 269:169193.  Raven K. G., 1986. Hydraulic characterization of a small ground-water flow system in fractured monzonitic gneiss, Inland Waters Directorate, National Hydrology Research Institute, Environment Canada, NHRI Paper No. 30, IWD Scientific Series No. 149, Pages: 1-133.  214  Reed, P., B. Minsker and A. J. Valocchi, 2000. Cost-effective long-term groundwater monitoring design using a genetic algorithm and global mass interpolation, Water Resources Research, 36(12), 3731-3741.  Smith L. and Schwartz F. W., 1993. Solute transport through fracture networks, in Flow and Contaminant Transport in Fractured Rock, edited by J. Bear, C. F. Tsang, and G. de Marsily, Academic Press, San Diego, CA, USA, pp. 129-167.  Smith, R. E. and D. A. Woolhiser, 1971. Overland flow on an infiltrating surface, Water Resources Research, 7(4), 899-913.  Snow, D. T., 1969. Anisotropic permeability of fractured media, Water Resources Research, 5(6), 1273-1289.  Sophocleous, M., 2002. Interactions between groundwater and surface water: the state of the science, Hydrogeology Journal, 10:52-67.  Stephenson, K. M., K. Novakowski, E. Davis and G. Heron, 2006. Hydraulic characterization for steam enhanced remediation conducted in fractured rock, Journal of Contaminant Hydrology, 82:220-240.  Storck, P., J. W. Eheart and A. J. Valocchi, 1997. A method for the optimal location of monitoring wells for detection of groundwater contamination in three-dimensional heterogeneous aquifers, Water Resources Research, 33(9), 2081-2088. 215  Sudicky E. A., J. P. Jones, Y. J. Park, A. E. Brookfield and D. Colautti, 2008. Simulating complex flow and transport dynamics in an integrated surface-subsurface modeling framework, Geosciences Journal, Vol. 12, No. 2, pp. 107-122.  Therrien R. and E. A. Sudicky, 1996. Three-dimensional analysis of variably-saturated flow and solute transport in discretely-fractured porous media, Journal of Contaminant Hydrology, 23:1-44.  Therrien R., E. A. Sudicky and R. G. McLaren, 2003. FRAC3DVS: An efficient simulator for three-dimensional, saturated-unsaturated groundwater flow and density-dependent, chain-decay solute transport in porous, discretely-fractured porous or dual-porosity formations, Groundwater Simulations Group, University of Waterloo, Waterloo, Ontario, Canada, pp.146.  Therrien R., R. G. McLaren, E. A. Sudicky and S. M. Panday, 2005. HydroGeoSphere: A three-dimensional numerical model describing fully-integrated subsurface and surface flow and solute transport, Groundwater Simulations Group, University of Waterloo, Waterloo, Ontario, Canada, pp.322.  VanderKwaak, J. E., 1999. Numerical simulation of flow and chemical transport in integrated surface-subsurface hydrologic system, PhD thesis, Univ. of Waterloo, pp.217.  216  VanderKwaak, J. E., and K. Loague, 2001. Hydrologic-response simulations for the R-5 catchment with a comprehensive physics-based model, Water Resources Research, 37(4), 999-1013.  van Genuchten, M. Th., 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils, Soil Sci. Soc. Am. J., 44, 892-898.  Wealthall G. P., A. Steele, J. P. Bloomfield, R. H. Moss and D. N. Lerner, 2001. Sediment filled fractures in the Permo-Triassic sandstones of the Cheshire Basin: observations and implications for pollutant transport, Journal of Contaminant Hydrology, 50:41-51.  Weatherill, D., T. Graf, C. T. Simmons, P. G. Cook, R. Therrien and D. A. Reynolds, 2008. Discretizing the fracture-matrix interface to simulate solute transport, Ground Water, 46(4), 606-615.  Winter, T. C., 1999. Relation of streams, lakes, and wetlands to groundwater flow systems, Hydrogeology Journal, 7:28-45.  Yenigül N. B., A. M. M. Elfeki, J. C. Gehrels, C. van den Akker, A. T. Hensbergen and F. M. Dekking, 2005. Reliability assessment of groundwater monitoring networks at landfill sites, Journal of Hydrology, 308: 1–17.  217  APPENDICES Appendix A. Results of flow and solute transport in Scenario SW-GW-1 A  Water Depth at Stream Outlet vs Number of Realizations  Depth (m)  0.0126 0.0125 0.0124 0.0123 0  25  50  75 Realizations  Minimum  B  100  Maximum  125  150  Average  Total Discharge at Stream Outlet vs Number of Realizations Discharge (m3/d)  121 120 119 118 117 0  25  50  Minimum  C  75 Realizations  100  Maximum  125  150  Average  Mean Flow Velocity (V x ) at Stream Outlet vs Number of Realizations  Velocity (m/d)  1220 1200 1180 1160 0  25  50  75  100  125  150  Realizations Minimum  Maximum  Average  Figure A.1 Statistics of water depth, discharge and flow velocity at stream outlet in different number of realizations for Scenario SW-GW-1 (higher Kuc). A Water depth, B total discharge, and C mean velocity (Vx). “Minimum”, “Maximum” and “Average” represents the smallest, the largest and the arithmetic average values of each parameter, for each set of additional realizations. 218  A Concentration (g/l)  Maximum Solute Concentration in Stream vs Number of Realizations 4.0E-03 3.0E-03 2.0E-03 1.0E-03 0.0E+00 0  25  50 Minimum  B  75 Realizations  100  Maximum  125  150  Average  Detection Time (yrs)  Detection Time at Compliance Boundary vs Number of Realizations 20 15 10 5 0 0  25  50 Minimum  75 Realizations Maximum  100  125  150  Average  Figure A.2 Statistics of the maximum solute concentration in stream at the time of detection at the compliance boundary, and detection time at compliance boundary in different number of realizations for Scenario SW-GW-1 (higher Kuc). A Maximum solute concentration, B detection time. “Minimum”, “Maximum” and “Average” represents the smallest, the largest and the arithmetic average values of each parameter, for each set of additional realizations.  219  Appendix B. Results of flow and solute transport in Scenario SW-GW-2a  A  B  C  D  E  F  Figure B.1 Coupled steady-state SW-GW flow in the first realization of Scenario SW-GW-2a (larger horizontal fracture apertures). A Surface water depth (m), B surface flow velocity Vx (m/d), C flux exchange (m/d) between surface and subsurface, D saturation profile in the subsurface, E hydraulic head (m) in the subsurface domain, and F groundwater hydraulic head (m) and flow velocity vectors in cross-section at Y=100m (beneath the streambed). Dimension (m).  220  A  Water Depth at Stream Outlet vs Number of Realizations  Depth (m)  0.0117 0.0112 0.0107 0.0102 0  25  50  75  100  125  150  175  200  Realizations Minimum  B  Maximum  Average  3  Discharge (m /d)  Total Discharge at Stream Outlet vs Number of Realizations 108 104 100 96 92 0  25  50  75  100  125  150  175  200  Realizations Minimum  C  Maximum  Average  Mean Flow Velocity (V x ) at Stream Outlet vs Number of Realizations  Velocity (m/d)  1200 1160 1120 1080 0  25  50  75  100  125  150  175  200  Realizations Minimum  Maximum  Average  Figure B.2 Statistics of water depth, discharge and flow velocity at stream outlet in different number of realizations for Scenario SW-GW-2a (larger horizontal fracture apertures). A Water depth, B total discharge, and C mean velocity (Vx). “Minimum”, “Maximum” and “Average” represents the smallest, the largest and the arithmetic average values of each parameter, for each set of additional realizations.  221  A Concentration (g/l)  Maximum Solute Concentration in Stream vs Number of Realizations 5.0E-03 4.0E-03 3.0E-03 2.0E-03 1.0E-03 0.0E+00 0  25  50  75  Minimum  B  100 125 Realizations Maximum  150  175  200  Average  Detection Time (yrs)  Detection Time at Compliance Boundary vs Number of Realizations 16 12 8 4 0 0  25  50  75 Minimum  100 125 Realizations Maximum  150  175  200  Average  Figure B.3 Statistics of the maximum solute concentration in stream at the time of detection at the compliance boundary, and detection time at compliance boundary in different number of realizations for Scenario SW-GW-2a (larger horizontal fracture apertures). A Maximum solute concentration, B detection time. “Minimum”, “Maximum” and “Average” represents the smallest, the largest and the arithmetic average values of each parameter, for each set of additional realizations.  222  Appendix C. Results of flow and solute transport in Scenario SW-GW-2b  A  B  C  D  E  F  Figure C.1 Coupled steady-state SW-GW flow in the first realization of Scenario SW-GW-2b (larger vertical fracture apertures). A Surface water depth (m), B surface flow velocity Vx (m/d), C flux exchange (m/d) between surface and subsurface, D saturation profile in the subsurface, E hydraulic head (m) in the subsurface domain, and F groundwater hydraulic head (m) and flow velocity vectors in cross-section at Y=100m (beneath the streambed). Dimension (m).  223  A  Water Depth at Stream Outlet vs Number of Realizations  Depth (m)  0.0124 0.0120 0.0116 0.0112 0  25  50  75  100 125 Realizations  Minimum  B  Maximum  150  175  200  Average  3  Discharge (m /d)  Total Discharge at Stream Outlet vs Number of Realizations 115 112 109 106 103 0  25  50  75  100  125  150  175  200  Realizations Minimum  C  Maximum  Average  Mean Flow Velocity (V x ) at Stream Outlet vs Number of Realizations  Velocity (m/d)  1180 1155 1130 1105 1080 0  25  50  75  100  125  150  175  200  Realizations Minimum  Maximum  Average  Figure C.2 Statistics of water depth, discharge and flow velocity at stream outlet in different number of realizations for Scenario SW-GW-2b (larger vertical fracture apertures). A Water depth, B total discharge, and C mean velocity (Vx). “Minimum”, “Maximum” and “Average” represents the smallest, the largest and the arithmetic average values of each parameter, for each set of additional realizations.  224  A Concentration (g/l)  Maximum Solute Concentration in Stream vs Number of Realizations 2.5E-02 2.0E-02 1.5E-02 1.0E-02 5.0E-03 0.0E+00 0  25  50  75  100  125  150  175  200  Realizations Minimum  B  Maximum  Average  Detection Time (yrs)  Detection Time at Compliance Boundary vs Number of Realizations 12 8 4 0 0  25  50  75 Minimum  100 125 Realizations Maximum  150  175  200  Average  Figure C.3 Statistics of the maximum solute concentration in stream at the time of detection at the compliance boundary, and detection time at compliance boundary in different number of realizations for Scenario SW-GW-2b (larger vertical fracture apertures). A Maximum solute concentration, B detection time. “Minimum”, “Maximum” and “Average” represents the smallest, the largest and the arithmetic average values of each parameter, for each set of additional realizations.  225  Appendix D. Results of solute transport in Scenario SW-GW-3a  A Concentration (g/l)  Maximum Solute Concentration in Stream vs Number of Realizations 3.0E-02 2.0E-02 1.0E-02 0.0E+00 0  10  20  30  40  50 60 Realizations  Minimum  B  70  Maximum  80  90  100  Average  Detection Time (yrs)  Detection Time at Compliance Boundary vs Number of Realizations 12 8 4 0 0  10  20  30  40  50  60  70  80  90  100  Realizations Minimum  Maximum  Average  Figure D.1 Statistics of the maximum solute concentration in stream at the time of detection at the compliance boundary, and detection time at compliance boundary in different number of realizations for Scenario SW-GW-3a (smaller αl in stream water). A Maximum solute concentration, B detection time. “Minimum”, “Maximum” and “Average” represents the smallest, the largest and the arithmetic average values of each parameter, for each set of additional realizations.  226  Appendix E. Results of solute transport in Scenario SW-GW-3b  A Concentration (g/l)  Maximum Solute Concentration in Stream vs Number of Realizations 3.0E-02 2.0E-02 1.0E-02 0.0E+00 0  10  20  30  40  50  60  70  80  90  100  Realizations Minimum  B  Maximum  Average  Detection Time (yrs)  Detection Time at Compliance Boundary vs Number of Realizations 12 8 4 0 0  10  20  30  40  50  60  70  80  90  100  Realizations Minimum  Maximum  Average  Figure E.1 Statistics of the maximum solute concentration in stream at the time of detection at the compliance boundary, and detection time at compliance boundary in different number of realizations for Scenario SW-GW-3b (smaller αt in stream water). A Maximum solute concentration, B detection time. “Minimum”, “Maximum” and “Average” represents the smallest, the largest and the arithmetic average values of each parameter, for each set of additional realizations.  227  Appendix F. Results of solute transport in Scenario SW-GW-4  A Concentration (g/l)  Maximum Solute Concentration in Stream vs Number of Realizations 4.0E-03 3.0E-03 2.0E-03 1.0E-03 0.0E+00 0  10  20  30  40  Minimum  B  50 60 Realizations  70  Maximum  80  90  100  Average  Detection Time (yrs)  Detection Time at Compliance Boundary vs Number of Realizations 6 4 2 0 0  10  20  30  40 Minimum  50 60 Realizations Maximum  70  80  90  100  Average  Figure F.1 Statistics of the maximum solute concentration in stream at the time of detection at the compliance boundary, and detection time at compliance boundary in different number of realizations for Scenario SW-GW-4 (lower detection threshold 1.0×10-4g/L). A Maximum solute concentration, B detection time. “Minimum”, “Maximum” and “Average” represents the smallest, the largest and the arithmetic average values of each parameter, for each set of additional realizations.  228  Appendix G. Results of solute transport in Scenario SW-GW-5  A Concentration (g/l)  Maximum Solute Concentration in Stream vs Number of Realizations 1.5E-02 1.0E-02 5.0E-03 0.0E+00 0  10  20  30  40  50  60  70  80  90  100  Realizations Minimum  B  Maximum  Average  Detection Time (yrs)  Detection Time at Compliance Boundary vs Number of Realizations 12 8 4 0 0  10  20  30  40  50  60  70  80  90  100  Realizations Minimum  Maximum  Average  Figure G.1 Statistics of the maximum solute concentration in stream at the time of detection at the compliance boundary, and detection time at compliance boundary in different number of realizations for Scenario SW-GW-5 (source zone on bedrock). A Maximum solute concentration, B detection time. “Minimum”, “Maximum” and “Average” represents the smallest, the largest and the arithmetic average values of each parameter, for each set of additional realizations.  229  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

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

Comment

Related Items