@prefix vivo: . @prefix edm: . @prefix ns0: . @prefix dcterms: . @prefix dc: . @prefix skos: . vivo:departmentOrSchool "Applied Science, Faculty of"@en, "Civil Engineering, Department of"@en ; edm:dataProvider "DSpace"@en ; ns0:degreeCampus "UBCV"@en ; dcterms:creator "Sheffer, Megan Rae"@en ; dcterms:issued "2009-08-14T19:45:21Z"@en, "2002"@en ; vivo:relatedDegree "Master of Applied Science - MASc"@en ; ns0:degreeGrantor "University of British Columbia"@en ; dcterms:description """Internal erosion is of great concern to owners and operators of earthfill dams worldwide. The migration of fine-grained soil particles from the dam core in the direction of seepage leads to the development of preferential zones of increased fluid flow, which can compromise the structural stability of the embankment. Standard dam safety monitoring techniques provide sparse sampling of subsurface hydraulic conditions and may not be sufficient to effectively detect the onset of internal erosion. Consequently, there is a clear need for a comprehensive monitoring tool that is sensitive to changing seepage conditions. The self-potential (SP) method is a non-invasive geophysical technique that responds directly to seepage through the phenomenon of streaming potential: the voltage gradient induced by the flow of water through a porous medium. However, current standard methods of SP data interpretation do not provide information about soil properties and seepage flow rates required by geotechnical engineers to assess dam performance. A three-dimensional numerical modelling tool is presented for predicting the SP response to fluid flow based on a comprehensive seepage analysis. Both the hydraulic regime and the resultant electrical potential distribution are calculated based on the distribution of hydraulic and electrical properties within the subsurface. Effective characterization of these parameters is fundamental in achieving an accurate numerical solution. An examination of the influence of internal erosion on hydraulic conductivity, cross-coupling conductivity and electrical conductivity is achieved through theoretical analyses based on published parametric data. The numerical procedure is validated against a theoretical closed-form solution, and is further verified through a comparison with measured values in a controlled physical model of an embankment dam. The capacity for threedimensional analysis and interpretation of the SP response to seepage is illustrated through preliminary models of a zoned embankment dam located in British Columbia."""@en ; edm:aggregatedCHO "https://circle.library.ubc.ca/rest/handle/2429/12220?expand=metadata"@en ; dcterms:extent "7408386 bytes"@en ; dc:format "application/pdf"@en ; skos:note "RESPONSE OF THE SELF-POTENTIAL METHOD TO CHANGING SEEPAGE CONDITIONS IN EMBANKMENT DAMS by MEGAN RAE SHEFFER B.Sc, Queen's University, 1995 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in THE FACULTY OF GRADUATE STUDIES Department of Civil Engineering We accept this thesis as conforming _ to the required standard THE UNIVERSITY OF BRITISH COLUMBIA April 2002 © Megan Rae Sheffer, 2002 In p resen t i ng this thesis in partial fu l f i lment of the requ i rements for an a d v a n c e d d e g r e e at the Univers i ty of Brit ish C o l u m b i a , I agree that the Library shall m a k e it f reely avai lable fo r re fe rence and s tudy. I fur ther agree that pe rm iss i on for ex tens ive c o p y i n g of th is thes is fo r scho lar ly p u r p o s e s may b e gran ted by the h e a d of m y d e p a r t m e n t o r by his o r her representat ives. It is u n d e r s t o o d that c o p y i n g o r pub l i ca t i on of this thesis for f inancial gain shall no t be a l l o w e d w i t hou t my wr i t ten p e r m i s s i o n . D e p a r t m e n t of T h e Un ivers i ty of Brit ish C o l u m b i a V a n c o u v e r , C a n a d a Da te D E - 6 (2/88) ABSTRACT Internal erosion is of great concern to owners and operators of earthfill dams worldwide. The migration of fine-grained soil particles from the dam core in the direction of seepage leads to the development of preferential zones of increased fluid flow, which can compromise the structural stability of the embankment. Standard dam safety monitoring techniques provide sparse sampling of subsurface hydraulic conditions and may not be sufficient to effectively detect the onset of internal erosion. Consequently, there is a clear need for a comprehensive monitoring tool that is sensitive to changing seepage conditions. The self-potential (SP) method is a non-invasive geophysical technique that responds directly to seepage through the phenomenon of streaming potential: the voltage gradient induced by the flow of water through a porous medium. However, current standard methods of SP data interpretation do not provide information about soil properties and seepage flow rates required by geotechnical engineers to assess dam performance. A three-dimensional numerical modelling tool is presented for predicting the SP response to fluid flow based on a comprehensive seepage analysis. Both the hydraulic regime and the resultant electrical potential distribution are calculated based on the distribution of hydraulic and electrical properties within the subsurface. Effective characterization of these parameters is fundamental in achieving an accurate numerical solution. An examination of the influence of internal erosion on hydraulic conductivity, cross-coupling conductivity and electrical conductivity is achieved through theoretical analyses based on published parametric data. The numerical procedure is validated against a theoretical closed-form solution, and is further verified through a comparison with measured values in a controlled physical model of an embankment dam. The capacity for three-dimensional analysis and interpretation of the SP response to seepage is illustrated through preliminary models of a zoned embankment dam located in British Columbia. ii TABLE OF CONTENTS Abstract ii Table of Contents iii List of Tables v List of Figures vi Acknowledgements viii 1 Introduction 1 1.1 Seepage Monitoring in Embankment Dams 1 / 1.2 Objectives and Thesis Outline 3 2 Theoretical and Conceptual Overview 6 2.1 Internal Erosion and Piping 6 2.2 The Streaming Potential Phenomenon 10 2.2.1 Introduction to Electrokinetic Phenomena 10 2.2.2 Properties of the Solid-Fluid Interface 11 2.2.3 Electrical Flow in Porous Media 13 2.2.4 Conceptual Model of Streaming Potential at the Pore Scale 16 2.3 Development of Flow Theory 18 2.3.1 Continuum Approach to Flow in Porous Media 18 2.3.2 Fundamental Equations 19 2.3.3 Coupled Flow Theory 21 2.3.4 Governing Equations: Assessing the Importance of Transient Effects 25 2.4 Summary and Conclusions 27 3 Forward Modelling 29 3.1 Previous Work 29 3.2 Finite Volume Formulation 3jl 3.2.1 Governing Equations and Generalized Approach 31 3.2.2 Finite Difference Equations 32 3.2.3 Numerical Solution Techniques 35 3.3 Seepage Flow Analysis 36 3.3.1 Mesh Geometry and Boundary Conditions 36 3.3.2 Free-Surface Approach 38 3.3.3 Solution Method 39 3.4 Electrical Flow Analysis 40 3.4.1 Mesh Geometry and Boundary Conditions 40 3.4.2 Solution Method 41 3.5 Validation of the Numerical Procedure 3DSP 42 3.6 Summary and Conclusions 46 4 Hydraulic and Electrical Properties 48 4.1 Physical Parameters 48 4.1.1 Hydraulic Conductivity 48 4.1.2 Electrical Conductivity 50 4.1.3 Cross-Coupling Conductivity 53 in 4.2 Sensitivity Analyses 58 4.2.1 Quantifying Internal Erosion 59 4.2.2 Effect of Internal Erosion in 1-D 63 4.2.3 Effect of Internal Erosion in 3-D 67 4.3 Summary and Conclusions 75 5 Numerical Analyses 77 5.1 Laboratory-Scale Embankment 77 5.1.1 Description of Laboratory Models and Physical Properties 77 5.1.2 Experimental and Numerical Results 81 5.2 Mica Dam, British Columbia 88 5.2.1 Project Description 88 5.2.2 Hydraulic Flow Analysis 90 5.2.3 Electrical Flow Analysis 93 5.2.4 Discussion of Results 96 5.3 Summary and Conclusions 98 6 Conclusions and Recommendations 100 6.1 Summary 100 6.2 Conclusions 101 6.3 Recommendations 102 References 104 Appendix A - Properties of the Solid-Fluid Interface 112 Appendix B - Development of the Continuity Equations 116 Appendix C - Compilation of Experimental Streaming Potential Data 132 Appendix D - Notation 139 iv L I S T O F T A B L E S Table 2-1: Typical ionic concentrations 12 Table 2-2: Ionic mobility in dilute solutions at 25°C 14 Table 2-3: Cation exchange capacity 16 Table 2-4: Primary and coupled flow phenomena 22 Table 4-1: Typical hydraulic conductivity values 49 Table 4-2: Typical fluid electrical conductivity values 51 Table 4-3: Typical electrical conductivity and resistivity values 52 Table 4-4: 1-D hydraulic and electrical properties 64 Table 4-5: 3-D isotropic hydraulic and electrical properties 68 Table 4-6: Unsaturated bulk conductivity of intact zoned embankment materials 74 Table 5-1: Physical properties derived from 1-D experimental data 79 Table 5-2: Measured bulk electrical conductivity and corresponding estimates of cross-coupling conductivity 80 Table 5-3: Numerical modelling input parameters: laboratory embankment 81 Table 5-4: Hydraulic head values reported at manometer port locations and volumetric seepage flow through embankment 83 Table 5-5: Numerical modelling input parameters: Mica dam 92 Table 5-6: Core piezometer data 93 Table A-l: Isoelectric points of geologic materials 112 v L I S T O F F I G U R E S Figure 2-1: Zoned embankment dam cross-section 6 Figure 2-2: Stern model of the solid-fluid interface 13 Figure 2-3: Convection and conduction currents in a saturated pore channel 17 Figure 3-1: Finite volume rectilinear grid configurations 33 Figure 3 -2: Hydraulic boundary conditions 37 Figure 3-3: Development of a seepage face on a free-outflow boundary 40 Figure 3-4: Injection well in a homogeneous subsurface 42 Figure 3-5: The method of images 44 Figure 3-6: Surface response from an injection well 45 Figure 4-1: Effect of fluid conductivity on measured values of the streaming potential coefficient for quartz sand 55 Figure 4-2: Variation of the effective pore fluid conductivity and streaming potential coefficient as a function of hydraulic radius in a quartz-water system 56 Figure 4-3: Influence of type and amount of fines content on hydraulic conductivity 61 Figure 4-4: Homogeneous and heterogeneous soil elements 63 Figure 4-5: Effect of internal erosion on the 1-D SP distribution over a homogeneous soil element 65 Figure 4-6: Effect of internal erosion on the 1-D SP distribution over a heterogeneous soil element 66 Figure 4-7: Transverse cross-section of the field-scale dam model with erosion channel .67 Figure 4-8: Transverse cross-section of hydraulic head and flux 68 Figure 4-9: Longitudinal cross-section of flow within saturated zone beneath crest 69 Figure 4-10: Distribution of sources of streaming current flow sources in longitudinal cross-section 71 Figure 4-11: Conduction current flow in longitudinal cross-section 72 Figure 4-12: SP profiles over dam surface at a longitudinal distance of 50m 73 Figure 4-13: SP difference contours over dam surface 73 Figure 4-14: Effect of unsaturated bulk conductivity on the surface SP distribution 75 Figure 5-1: Laboratory embankment dimensions with manometer and electrode measurement locations 78 Figure 5-2: Grain size distribution of embankment materials 78 Figure 5-3: Predicted hydraulic head distribution 82 Figure 5-4: Measured and modelled surface SP distribution over the homogeneous dam. 85 Figure 5-5: Measured and modelled surface SP distribution over the dam with upstream defect 86 Figure 5-6: Modelled SP distribution through a cross-section of the embankment at the 22.5cm reservoir level 87 Figure 5-7: Plan view of Mica dam site 89 Figure 5-8: Mica dam reservoir elevation 90 Figure 5-9: Mica dam seepage conditions at high pool 91 Figure 5-10: Modelled surface SP distribution at high pool 95 Figure 5-11: Modelled surface SP transverse profiles 95 vi Figure 5-12: Modelled surface SP difference data 96 Figure 5-13: Continuous SP response recorded from installed electrodes 98 Figure A-l: Stern model of the electrical double layer at the solid-fluid interface 114 Figure A-2: Change in potential away from a solid surface with potential y 0 115 vii A C K N O W L E D G E M E N T S Funding for this research was awarded through the BC Hydro - UBC Professional Partnership Program, with partial funding from a Natural Sciences and Engineering Research Council of Canada grant, and is greatly appreciated. The cooperation of Powertech Labs Incorporated and BC Hydro in providing experimental and in-situ data is gratefully acknowledged. Special thanks are extended to Dick Campanella for his assistance, John Howie for his guidance and encouragement in pursuing this area of research, Doug Oldenburg and Eldad Haber in the Earth and Ocean Sciences Department for their help and collaboration in the development of the numerical code, and Bob Corwin for sharing his enthusiasm, knowledge and experience with the SP technique. viii 1 INTRODUCTION 1.1 Seepage Monitoring in Embankment Dams Seepage flow monitoring is of great importance to owners and operators of embankment dams worldwide. Permeable construction materials coupled with hydraulic gradients imposed by the upstream reservoir result in a limited amount of seepage flow through the body of an earth dam. Well-designed structures typically consist of an impervious central core, which is flanked by an outer shell of higher permeability materials, including transition filters and drains, that serve to stabilize the core and resist the thrust of retained water to the dam foundation. Deviations from expected or predicted seepage flow rates may be indicative of preferential seepage zones or piping through the core or foundation, which can severely compromise the structural stability of the embankment. Internal erosion, or the loss of fine-grained material due to seepage forces, is one of the most prevalent failure modes in embankment dams (ICOLD, 1995; Sherard, 1963). Complete failure can result in catastrophic land damage and loss of life due to the collapse of the dam structure and mass release of water from the reservoir. The majority of failed dams studied by the International Congress on Large Dams did not have a functional monitoring system to assess dam performance (ICOLD, 1995). Regular inspection and monitoring programs are therefore paramount in ensuring dam safety. Standard instrumentation is used to monitor seepage flow from the dam drainage system, displacements and pore water pressure within the dam body at discrete measurement points. These systems may be installed at the time of construction, but their existence and efficiency is dependent on the design and age of the dam. Standard methods provide sparse sampling of subsurface hydraulic conditions and are often not sufficient to effectively detect the early onset of internal erosion. A striking example of this is the 1995 sinkhole incident at Whiteman's dam, located outside of Canmore, Alberta. At the time of discovery and investigation of the two damaged zones beneath the crest, which were attributed to piping of the fine-grained core material, absolutely no other signs of dam distress were evident and no anomalous changes in seepage or piezometer levels were recorded (Courage et al, 1997). Small increases in seepage flow due to piping will tend not to exceed the maximum flow capacity of downstream drainage systems, which are typically designed with a factor of 1 safety of 5 to 10 (Cedergren, 1993), and consequently will not generate a significant response in piezometers located in the downstream shell. Moreover, the proximity of these piezometers or observation wells to piped or eroded zones will dictate their sensitivity to changing seepage conditions. Increased seepage may be detectable at drainage collection sites or weirs, but these give no indication of the source or location of leakage. Once seepage flows exceed the drainage capacity of the structure, there is an impending risk of dam failure. Consequently, there is a clear need for a comprehensive monitoring system capable of detecting internal erosion in the very early stages of development. The application of geophysical methods developed primarily for resource exploration to geotechnical and hydrologic problems has grown as different and more comprehensive techniques are sought to characterize the near subsurface. Geophysics provides a means of characterizing the properties of the ground through the application of physical theories and measurements. A number of techniques are sensitive to variations in soil and fluid properties, and enable imaging of specific physical property contrasts that result from changing seepage conditions. These non-invasive methods may be employed at surface or down existing wells, which is preferable to conventional investigation methods that involve drilling since they avoid the risk of structural damage to the core. The information garnered from geophysical survey techniques can be used in conjunction with standard methods of investigation to generate a more comprehensive picture of subsurface flow conditions. The self-potential (SP) method is one such technique currently being used to assess seepage conditions in earth dams. Originally developed for application in the mining industry, this geophysical method detects natural or 'spontaneous' voltages in the ground generated by chemical, thermal or hydraulic processes. The SP method has found broad application ranging from geothermal exploration to engineering and environmental problems (Corwin, 1990). In seepage investigations, the SP method responds directly to seepage flow through the mechanism of streaming currents or streaming potential. Electrical currents, which are characterized by potential or voltage differences, are generated by the flow of water in the subsurface and across zones of differing soil properties. This electrokinetic coupled-flow process is the complementary phenomenon to electro-osmosis, which has been studied extensively in geotechnical engineering (Casagrande, 1983). 2 At present, the SP technique is predominantly a qualitative tool. The distribution of electrical potential, or voltage, can be characterized on the surface of the embankment and in the reservoir through a series of electrode measurements. Data are collected at different reservoir levels, typically at high and low pool, and compared to isolate the SP response to changes in the pattern of seepage. Areas of increased seepage appear as voltage anomalies relative to background values, and estimates of the location and depth of these zones can be determined using geometric interpretation techniques (Corwin, 1990). The SP method has been used to successfully identify anomalous areas in embankment dams, which were subsequently confirmed as leakage zones during remediation (Corwin, 1991). However, the current state of practice does not enable a direct quantitative assessment of seepage conditions from the measured SP response. In order for the SP method to develop into an accepted standard monitoring tool for seepage investigations, geotechnical engineers must be able to extract information on physical properties, such as hydraulic conductivity and flow rate, that influence the measured SP response. With an analytical means of predicting the SP response to seepage, the interpreter can attribute anomalies to internal erosion with greater confidence, and relations between changes in electrical potential response and changes in physical properties can be developed. With this capability, we can begin to investigate the sensitivity of the SP method to changing seepage conditions due to internal erosion, and assess its viability as a long-term monitoring tool for evaluating dam performance. 1.2 Objectives and Thesis Outline The purpose of the present research is to study the effect of changing seepage conditions on the self-potential (SP) response. This thesis is focused on the following three primary objectives. 1. Develop a 3-D numerical forward modelling procedure to predict the self-potential response to seepage flow based on a comprehensive and realistic seepage analysis. 2. Examine the interdependence of hydraulic conductivity, cross-coupling conductivity and electrical conductivity, and assess the influence of internal erosion on these three model input parameters for soils typically used in earth dam construction. 3 3. Study the sensitivity of self-potential measurements to changing seepage conditions through numerical analyses and comparison with measured results. In order to understand the effect of changing seepage conditions on the SP response, an overview of the physical phenomena occurring within the embankment is required. The development of internal erosion and piping in embankment dams, and the mechanism responsible for the SP response to seepage are discussed in Chapter 2. An understanding of the electrokinetic phenomenon of streaming potential is garnered through a study of flow at both microscopic and macroscopic levels. A conceptual model of hydraulic and electrical flow at the pore scale is introduced, along with a discussion of the electrical properties of the solid-fluid interface and the mechanisms of electrical current conduction in a porous medium. The governing macroscopic equations for hydraulic and electrical flow, namely Darcy's law and Ohm's law, are linked through the theory of coupled flow. These equations provide the basis for a numerical analysis of the streaming potential phenomenon. For the quantitative study of the SP response to seepage in embankment dams, numerical methods must be applied to determine both the hydraulic and electrical flow regimes. Previous studies have not always focused on a comprehensive seepage analysis, and consequently the SP response is based on a solution that may not be representative of actual seepage conditions through an embankment dam. Chapter 3 outlines the development of the three-dimensional forward modelling scheme, which involves discretization of the domain and coupled flow equations, and application of the boundary conditions using the finite volume method. A discussion of numerical solution techniques includes the application of Visual MODFLOW for a rigorous seepage analysis, and coded algorithms in MATLAB for a comprehensive solution to the electrical flow problem. Verification of the numerical procedure 3DSP is achieved through a comparison with the analytical solution to a simple 3-D flow problem. To predict the effect of changing seepage conditions on the SP response, the physical properties that characterize the porous medium must be estimated, and changes in these properties due to internal erosion must be quantified. Chapter 4 discusses the required model input parameters, which include hydraulic conductivity, electrical conductivity, and cross-coupling conductivity, and reveals their strong interdependence. The effect of internal erosion on soil structure, fluid and electrical properties are discussed, and an attempt is made 4 to quantify these changes from theoretical and laboratory studies in the literature. Finally, the sensitivity of the self-potential response to physical property variations is assessed through simple numerical analyses. Examples of the SP response to changing seepage conditions are presented in Chapter 5 with the application of the numerical procedure to actual seepage problems. Numerical models of a laboratory-scale embankment, over which SP measurements were made, enable verification of the numerical procedure using input parameters derived from measured data. Preliminary models of Mica Dam illustrate the three-dimensional character of the SP response over the surface of the embankment, and provide a basis for further modelling studies and data interpretation. 5 2 THEORETICAL AND CONCEPTUAL OVERVIEW How does the self-potential method respond to changing seepage conditions caused by internal erosion? The first step in answering this question requires understanding the process of internal erosion and the mechanism responsible for the SP response to seepage, the electrokinetic phenomenon of streaming potential. Section 2.1 discusses the causes and consequence of internal erosion in earth dams, and the streaming potential phenomenon is described in Sections 2.2 and 2.3 through a conceptual and theoretical study of flow at both microscopic and macroscopic scales. 2.1 Internal Erosion and Piping A typical zoned embankment constructed using modern methods of design consists of a central impervious core, which is flanked by transition zones that grade towards a coarse outer shell, as shown in Figure 2-1. The core is often constructed of a broadly graded material, such as glacial till, which contains a significant amount of silt- and/or clay-sized particles. Transition zones are designed to act as filters to prevent any seepage-induced loss of fine particles from the core or retained fill. Drainage layers intercept seepage through the dam body and foundation, and provide a downstream outlet, keeping the bulk of the downstream shell in a dry, and consequently more stable, state. Maximum controlled pool elevation 830 Height = 620 ft Acceptable foundation Core' Grout curtain (projected) Figure 2-1: Zoned embankment dam cross-section (from Burke et al, 1972) 6 Seepage through an embankment is not strictly uniform, with concentrated seeps or leaks through the impervious core developing as a result of poor construction control, through the uneven compaction of layers or substandard compaction around structures buried within the embankment. Weak bonds between the embankment and foundation or abutments, and cracks in the core caused by differential settlement, hydraulic fracture or earthquake loading can also lead to concentrated leakage. Due to a lack of confinement, erosive forces are largest where water emerges from these outflows in the downstream section of the dam (Sherard et al, 1963). Any removal of the surrounding soil by erosion only increases the seepage intensity, which can lead to further particle movement. When the erosive seepage forces exceed those resisting erosion, fine-grained soil particles within the soil matrix will be transported in the direction of flow. This process of internal erosion may lead to piping, or the development of preferential seepage pathways that migrate in an upstream direction, which can compromise the stability of the embankment and cause its ultimate failure. The basic role of filters and drains in embankment dams is to retain core material and finer-grained fill while effectively discharging seepage flow through the embankment or foundation without the accumulation of excess pore pressure (ICOLD, 1994). The role of the filter adjacent to the core is to prevent major loss of fines should leaks develop within the dam core. Coarser transition or drainage zones located downstream serve to protect this finer-grained filter. If piping develops in the core, the transition zone must assume the role of primary seepage barrier. This is accomplished through the process of self-filtration. Internal erosion in broadly graded soils results when seepage forces induce a migration of fine particles from the core out through the filter soil. Particles exceeding the maximum opening size of the filter, which is dictated in part by grain size distribution, are stopped by the filter and form an arched matrix of coarse particles. This retained layer filters finer particles migrating through the core, which in turn filters even finer ones, until the self-filtration action is complete and no further particles escape the core. Since the self-filtration process is initiated with seepage flow, finite migration of fine-grained particles from the core to the downstream shell is expected upon first filling of the reservoir, and minor settlements may occur. If internal erosion progresses and piping develops, increased seepage can result from an increase in porosity and permeability due to 7 the loss of fines and deposition of coarser-grained particles in the eroded zone. The seepage increase may be detectable at the toe of the dam, in the form of clear or muddy water flowing at a constant rate or at a slightly increasing rate for a number of years prior to failure (Sherard et al, 1963). Collapse, or self-repair of the eroded zone by settlement of the overlying material may interrupt this period of increased seepage, and can eventually translate to visible sinkholes or depressions on the dam crest in the later stages of development. Although modern design methods include the use of transition zones between the impervious core and pervious outer shell materials, filter criteria have developed considerably in more recent time, particularly with regard to the retention or piping criterion. As the core acts as a main barrier to seepage, the erosion resistance of the core material, and the compatibility of this material with those comprising the upstream and downstream filters are key factors in ensuring a low risk of piping failure. The erosion resistance of a given soil is linked to its internal stability (Sherard, 1979). A broadly graded material is considered internally unstable if the fine fraction is able to migrate freely within the primary structure of coarser particles. Kenney and Lau (1985) have shown that the internal stability of a broadly graded soil is a function of the coefficient of uniformity and shape of the gradation curve. In designing filters using standard Terzaghi criteria, which rely on representative grain sizes, the internal stability of broadly graded core or filter materials is not taken into account, as these criteria were initially developed for uniformly graded soils (Cedergren, 1973). Consequently, the development of more comprehensive methods of filter design has been the focus of many researchers (Kenney et al, 1985; Honjo and Veneziano, 1989; Jansen, 1988; Lafleur et al, 1989; Sherard and Dunnigan, 1985, 1989; Vaughan and Soares, 1982). A good review of the evolution of filter design and the current state of practice has been published by the International Commission on Large Dams (ICOLD, 1994). Since many dams were constructed in the mid-1900s, the behaviour of older dams that do not contain filters, and those in which filters and drains do not meet modern standards must be monitored carefully for internal erosion. The most dangerous, and likely the primary source of, failures due to piping are those caused by water flowing freely through open cracks in the core (Cedergren, 1993). Unless cracks develop, typical seepage velocities through intact core material are too small to 8 initiate the migration of fine particles, even in broadly graded soils classified as internally unstable (Sherard, 1979). In cohesionless soils, cracks tend to collapse upon saturation as capillary bonds are destroyed (Vaughan and Soares, 1982). However, clay cores, or cores constructed of broadly graded soils containing clay, are able to sustain small cracks due to inter-particle bonds. Clay cores of low plasticity that are compacted at low water contents are most susceptible to cracking due to differential settlement (Sherard et al, 1963). Cracks pose a serious threat to dam safety since they limit the capacity for self-filtering behaviour. The seepage velocity and pressure within the crack, along with the size and shape of the crack, influence the amount of segregation and deposition of soil particles, and consequently the degree to which self-filtration occurs at the core-filter interface (Vaughan and Soares, 1982). As a crack enlarges due to erosion, the increased seepage velocity will promote self-filtering as larger particles are transported downstream. Subsidence will occur once the size of the crack exceeds that which can be sustained by inter-particle forces. A filter design suitable for an intact clay core may not be sufficient to prevent internal erosion along a crack through the core, as was shown in the study of the Balderhead Dam failure (Vaughan and Soares, 1982). The dam was constructed in 1965 but suffered internal erosion damage in 1967 due to poor performance of the downstream crushed limestone filter drain. Low seepage velocities caused segregation of eroded particles within the crack, and did not appear to initiate self-filtering behaviour at the core-filter interface. The filter was too coarse to prevent the erosion of fines due to cracking of the core, which was constructed of an erosion-resistant broadly-graded glacial till with clay matrix, compacted dry of optimum. If piping does develop within an embankment, the loss of fine-grained material from the crack or eroded zone results in an increase in porosity and water content, and increased seepage flow arises from an increase in hydraulic conductivity. As a result of these structural changes, the electrical properties of the embankment will vary, and these in turn will influence the self-potential response. The validity of any numerical modelling of the SP response to seepage relies on realistic estimates of the relevant physical properties. An examination of these properties and their variation with the onset of internal erosion is the focus of Chapter 4. 9 2.2 The Streaming Potential Phenomenon In general terms, electrokinetic phenomena involve the flow of electricity and the lateral movement of one material phase with respect to another (Overbeek, 1952). In a geotechnical context, these phenomena describe the coupled behaviour of hydraulic and electrical flow within a saturated soil or porous rock formation. Streaming potential is just one electrokinetic process that can develop in a soil saturated or partially saturated with an electrolytic fluid such as water. An investigation into this phenomenon requires an understanding of the relation between hydraulic and electrical flow systems in a porous medium. This investigation begins at the pore scale: a brief discussion of electrochemistry at the solid-fluid interface leads to a conceptual model of streaming potential and methods of electrical current conduction in porous media. 2.2.1 Introduction to Electrokinetic Phenomena Relative movement between the electrically charged surface of a solid particle and free ions in solution define different electrokinetic processes, such as electrophoresis, sedimentation potential, electro-osmosis, and streaming potential. The solid surface is ordinarily either a solid particle suspended in solution or the medium through which the solution flows. Common applications of electrokinetic techniques in geotechnical engineering include the use of electrophoresis for dewatering dredged sediments, such as clay suspensions and slurry deposits, and the use of electro-osmosis for the consolidation and strengthening of soft clayey soils (Mitchell, 1991). Both of these processes involve the application of an external electrical current, or direct current (dc) electric field, to induce particle movement or fluid flow. Electrophoresis and sedimentation potential are phenomena that involve the forced movement of charged solid particles in solution. An applied electric field induces the movement of colloidal particles, such as suspended clay particles, through the process of electrophoresis, whereas a voltage difference, termed sedimentation potential, is generated by the settlement of particles in a gravitational field. Electro-osmosis, described as the movement of a fluid through a soil matrix in an applied dc electric field, has been studied extensively as a means of increasing the strength of 10 clays and fine-grained soils (Casagrande, 1983). Streaming potential is the complementary phenomenon: the generation of electrical current flow or an electric field, characterized by a voltage or potential gradient, by the movement of fluid through a porous medium from an imposed hydraulic gradient. 2.2.2 Properties of the Solid-Fluid Interface Electrokinetic phenomena are restricted to the capillaries of a porous medium and develop at the interface between the pore fluid and the solid particle. Consequently, investigations into the streaming potential phenomenon are aided by an understanding of the electrical character of this interface. The solid surfaces within a saturated soil display a net electrical charge, either positive or negative, which is controlled by surface chemistry. Factors influencing the polarity of this charge include mineralogy or chemical composition of the solid particles; pH of the pore fluid; concentration and type of dissolved ions in the pore fluid; imperfections in the molecular or lattice structure of the solid, specifically in clay minerals; and adsorption of organic material (Stumm, 1992). Most geologic materials exhibit a negative surface charge when saturated with natural waters displaying typical pH and ion concentrations. It is important to note that in order for electrokinetic effects to occur, the saturating solution must be electrolytic. That is, the water must contain free, or dissociated ions available to balance the charged surface of the soil grain. As most natural surface and soil waters contain free ions in solution due to the process of geologic weathering and other dissolved chemical species, this condition is easily satisfied. Table 2-1 illustrates average ionic concentrations as a function of fluid volume for fresh surface water sources, expressed in terms of both milligrams of solute and milliequivalents [meq], or millimoles of electrical charge. Although fluid chemistry is highly site-specific, the data given in Table 2-1 provide an order of magnitude approximation of reservoir waters in contact with glacial deposits existing in geologic and climatic zones similar to the examples provided. The solid-fluid interface in a saturated soil is characterised by an electrical double layer due to the attraction of ions in the pore fluid to the negative grain surface. The thickness of this double layer, termed the Debye length A, is controlled by the ionic concentration, dielectric properties and temperature of the pore fluid, such that a decrease in 11 dissolved ions results in an increase in thickness of the electrical double layer (Morgan, 1989). Table 2-1: Typical ionic concentrations (from 'Drever, 1988; 2 R.L.&L. Environmental Services Ltd, 1999). South Cascade Lake, Brilliant Dam Headpond, North Cascade Mountains, Kootenay River system, Ion Washington, USA 1 British Columbia, Canada2 Average River Water' mg/L meq/L mg/L meq/L mg/L meq/L Na + 0.41 0.018 2.1 0.091 5.15 0.224 K + 0.82 0.021 0.6 0.015 1.3 0.033 Ca 2 + 2.92 0.146 20.9 1.043 13.4 0.669 Mg 2 + 0.36 0.030 4.8 0.395 3.35 0.276 Cl\" 0.47 0.013 0.9 0.025 5.75 0.162 S0 4 2\" 1.83 0.038 11.0 0.229 8.25 0.172 HCO3\" 10.43 0.171 78 1.278 52 0.852 pH (mean) 6.62 7.52 The most accepted model of the double layer charge distribution is the Stern model, which distinguishes between an adsorbed layer of positive ions and a more loosely bound diffuse outer layer (Adamson, 1990). Figure 2-2 shows a simple schematic of the double layer, with a Helmholtz plane defining the boundary between the adsorbed Stern layer and the diffuse outer layer, and a slipping plane delineating the shear zone between the fixed and free pore fluid. The electrical potential at this shear plane is termed the zeta potential, ^ , which can be derived from experimental measurements and is used to approximate the true surface potential or surface charge of the solid. Under static conditions, the system is electrically neutral with ions in the double layer completely balancing the negative surface charge of the soil grain. Upon initiation of fluid flow, electrokinetic phenomena result from relative movement between the solid and fluid along the shear plane separating the Stern layer from the diffuse zone within the electrical double layer. A more detailed discussion of surface chemistry and structure of the electrical double layer can be found in Appendix A. 12 Figure 2-2: Stern model of the solid-fluid interface (adapted from Hibbert, 1993). 2.2.3 Electrical Flow in Porous Media Electrical current flow in porous media is controlled by the presence of water. Electrolytic conduction, or the flow of current by means of free ions in solution, predominates in most geologic systems since the pore water typically provides a path of lower resistance, as compared to the minerals comprising the soil matrix. Most solid geologic materials, with the exception of metals and semi-conducting minerals, can be considered solid electrolytes, such that electrical current is transferred primarily by means of ionic conduction within the minerals that comprise the rock or soil matrix. This type of conduction is facilitated by inherent and thermally induced imperfections in the molecular structure of the particles. The electrical conductivity, or resistance to current flow, is controlled by the mobility of different ions through this structure as well as the number of ions available for conduction (Keller and Frischknecht, 1966). However, electrolytic conduction in solid materials is very small compared to 13 electrolytic conduction in a fluid. Only the concentration and mobility of dissolved ions controls the electrical conductivity of a fluid, which is dependent on water temperature. Natural surface and ground water may contain a wide range of dissolved species at varying concentrations, as was illustrated in Table 2-1. Table 2-2 shows the mobility of typical dissolved ions in dilute solution at 25 degrees Celsius. Ionic mobility may be defined as the velocity of an ion in a potential gradient of 1 V/cm (Hem, 1985) and is directly influenced by fluid temperature. An increase in temperature decreases the viscosity of the fluid, thus increasing the mobility of dissolved ions. Equation 2-1 describes the electrical conductivity of a fluid, oy [S/m] as the product of ionic concentration c [eq/m3, or meq/L] and mobility m [m2/s.V] for each ion in solution, where F is Faraday's constant [C/eq] (Keller and Frischknecht, 1966). Table 2-2: Ionic mobility in dilute solutions at 25°C (from Keller and Frischknecht, 1966). Ion Mobility (m2/s.V) H + 36.2 x 10\"8 OH\" 20.5 x 10\"8 S0 4 2 \" 8.3 x 10\"8 N a + 5.2 x 10\"8 c r 7.9 x 10\"8 K + 7.6 x 10\"8 N 0 3 \" 7.4 x 10\"8 L i + 4.0 x 10\"8 HCO3\" 4.6 x 10\"8 a l = F^c,m, Equation 2-2 describes the increase in fluid conductivity with temperature T [°C] for a given solution, where p\\ is the temperature coefficient of the dissolved compound. The temperature coefficient of a typical natural electrolytic solution may range from 0.022, which 14 c o r r e s p o n d s t o a d i l u t e N a C l s o l u t i o n , t o a m a x i m u m o f 0 . 0 2 5 ( K e l l e r a n d F r i s c h k n e c h t , 1 9 6 6 ) . T h e s e v a l u e s t r a n s l a t e t o a 2 . 2 % to 2 . 5 % i n c r e a s e i n e l e c t r i c a l c o n d u c t i v i t y p e r ° C . ° / c r ) = ° / (25«c) t1 + AT ( T - 25)] ( 2 - 2 ) I n a d d i t i o n t o t h e i n f l u e n c e o f f l u i d c o n d u c t i v i t y a n d t h e c o n d u c t i v i t y o f t h e s o i l m a t r i x o n t h e b u l k c o n d u c t i v i t y o f a s o i l , s u r f a c e c o n d u c t i v i t y e f f e c t s at t h e s o l i d - f l u i d i n t e r f a c e i n finer-grained s o i l s c a n b e s i g n i f i c a n t . P a r t i a l d i s s o c i a t i o n o f a d s o r b e d c a t i o n s w i t h i n t h e e l e c t r i c a l d o u b l e l a y e r p r o m o t e s i o n e x c h a n g e a n d c o n t r i b u t e s t o t h e i o n i c c o n d u c t i v i t y o f the f l u i d i n t h e f o r m o f a c o n d u c t i o n p a t h i m m e d i a t e l y a d j a c e n t t o t h e p a r t i c l e s u r f a c e . T h i s s u r f a c e c o n d u c t i v i t y e f fect i s d o m i n a n t i n c l a y p a r t i c l e s d u e t o t h e i r l a y e r e d s t r u c t u r e a n d h i g h a v e r a g e s u r f a c e a r e a t o v o l u m e r a t i o , b u t c a n a l s o a r i s e i n n o n - c l a y e y s o i l s i f t h e a v e r a g e p o r e d i a m e t e r o f a g i v e n s o i l i s s u f f i c i e n t l y s m a l l . S u r f a c e c o n d u c t i v i t y i s s i g n i f i c a n t w h e n t h e s o i l p r o p e r t i e s s a t i s f y E q u a t i o n 2 - 3 , w h e r e A i s t h e D e b y e l e n g t h , a i s t h e a v e r a g e p o r e d i a m e t e r , a n d d i s s o m e r e p r e s e n t a t i v e g r a i n s i z e ( M o r g a n et a l , 1 9 8 9 ) . - < 1 0 , w h e r e a«0.2d to 025d ( 2 - 3 ) A A l t h o u g h c o h e s i o n l e s s s o i l s e x h i b i t f i n i t e c a p a c i t y f o r c a t i o n e x c h a n g e , t h e s e s o i l p a r t i c l e s a r e e s s e n t i a l l y n o n - c o n d u c t i n g a n d c a n b e c o n s i d e r e d i n s u l a t o r s i n t h e a b s e n c e o f a n y m o i s t u r e . C o n s e q u e n t l y , t h e c o n d u c t i v i t y o f t h e s o i l m a t r i x i n a s a n d o r g r a v e l i s n e g l i g i b l e s u c h t h a t p o r o s i t y a n d c o n n e c t i v i t y o f t h e p o r e s p a c e , d e g r e e o f s a t u r a t i o n , a n d e l e c t r i c a l c o n d u c t i v i t y o f t h e p o r e f l u i d g o v e r n e l e c t r i c a l f l o w . S i n c e a f i n i t e l a y e r o f f l u i d s u r r o u n d s i n d i v i d u a l s o i l g r a i n s e v e n i n a n u n s a t u r a t e d state , t h e t h i c k n e s s a n d c o n n e c t i v i t y o f t h i s w a t e r l a y e r w i l l u n d o u b t e d l y i n f l u e n c e t h e b u l k e l e c t r i c a l c o n d u c t i v i t y , a n d c o n s e q u e n t l y t h e f l o w o f e l e c t r i c a l c u r r e n t t h r o u g h t h e p o r o u s m e d i u m . I n c l a y s , a h e i g h t e n e d c a p a c i t y f o r c a t i o n e x c h a n g e at t h e p a r t i c l e s u r f a c e a n d w i t h i n t h e l a t t i c e s t r u c t u r e c o n t r i b u t e s to t h e b u l k c o n d u c t i v i t y . T a b l e 2-3 l i s t s t h e c a t i o n e x c h a n g e c a p a c i t y ( C E C ) o f a c l a y e y s a n d a n d t y p i c a l r a n g e s f o r d i f f e r e n t c l a y s , e x p r e s s e d i n t e r m s o f t h e a m o u n t o f a v a i l a b l e c h a r g e p r e s e n t p e r v o l u m e o f s o i l . 15 Table 2-3: Cation exchange capacity (from Drever, 1988; Keller and Frischknecht, 1966; and McNeill, 1980). Soil type CEC (meq/lOOg) Sand (1.7% organic, 7% clay) 6.3 Kaolinite 1 - 10 Montmorillonite 80-150 Illite 10 -40 Vermiculite 120-200 Chlorite <10 In fine-grained soils saturated with very fresh water, or fluid of low conductivity, surface conductivity effects govern the bulk electrical conductivity. Since the thickness of the electrical double layer is inversely proportional to ionic concentration, the net effect of an increase in pore fluid conductivity is to decrease the contribution of surface conductivity effects to the bulk electrical conductivity. Further discussion of the electrical properties of soils is contained in Chapter 4 along with typical ranges in conductivity values for various soil types. 2.2.4 Conceptual Model of Streaming Potential at the Pore Scale Under static no-flow conditions, the saturated system is in equilibrium with a balance of electrical charge across the solid-fluid interface. Upon initiation of fluid flow, this charge balance is disrupted giving rise to electrical currents through the phenomenon of streaming potential. As water passes through the pore channels of a soil, positive ions are pulled in the direction of fluid flow within the diffuse outer layer of ions. This physical movement of charged ions defines an electrical current within the pore channels, which is known as the streaming current or convection current, JConvection> displayed schematically in Figure 2-3. As positive ions are mobilized with the flow of water, a net negative charge develops in the upstream end of the pore due to a charge imbalance across the soil-fluid interface. Downstream, the accumulation of positive ions results in a net positive charge within the pore channel. Streaming potential is this separation of charge, or voltage difference, parallel to the direction of flow that defines the streaming or convection current. The voltage 16 gradient induces a secondary electrical current, called the conduction current, which flows in an opposite direction to the convection current. Figure 2-3: Convection and conduction currents in a saturated pore channel. The convection and conduction currents balance each other under steady-state fluid flow conditions, and create a closed current loop. Although the convection current is strictly limited to within pore channels experiencing fluid flow, the conduction current permeates the entire medium, as shown by the dashed arrows in Figure 2-3. More specifically, the conduction current will follow the path of least resistance, dictated by the electrical properties of the medium. In saturated, predominantly cohesionless soils, in which the individual soil grains do not display significant electrical conductivity, the conduction current will tend to flow through the pore fluid and along the solid-fluid interface by means of electrolytic conduction, as discussed in Section 2.2.3. Ultimately, only the conduction current can be measured. Therefore, the magnitude of the electrical potential defining this current is used to characterize the streaming current. 17 2.3 Development of Flow Theory In order to study the coupled behaviour of fluid and electrical current flow, it is essential to examine the common fundamental laws that govern each individual flow system at a macroscopic level. Mathematical relations essential to a numerical analysis of the streaming potential phenomenon are outlined using the theory of coupled flow. 2.3.1 Continuum Approach to Flow in Porous Media In both hydraulic and electrical flow systems, the physical quantity under study is regarded as a continuum, and is usually quantified volumetrically in terms of mass or density. However, the fundamental physical process takes place at the molecular level. In a hydraulic system, the basic quantity is the ion; in an electrical system, it is a single point charge. Ions can be grouped into discrete molecules, but for practical purposes, water composed of these molecules is treated as continuous and is quantified in terms of mass density. In the same fashion, individual protons or electrons can be grouped into discrete integer multiples of these basic charge units, but because of their minuscule size, charge is more effectively treated as a microscopic continuum, described in terms of charge density. The fluid continuum is bounded by the solid surfaces in the porous medium. Fluid flow within the void space o f the medium is dependent on pore structure, which includes the connectivity of the pore space and tortuosity of the flow path. However, the analysis of flow systems at this microscopic scale is impractical, since it is almost impossible to define the geometry and boundary conditions of the solid-fluid interface throughout the entire soil matrix under study. Consequently, a macroscopic approach is required. The effects of pore structure on flow can be explained using coefficients that represent the properties of an averaging volume, termed a control volume, or representative elemental volume (REV) . These macroscopic properties; which include volumetric porosity, and all volumetric flux terms and conductivity coefficients used in the flow equations defined in the following sections; are independent of the size of the R E V and define average values that can be applied throughout the domain for which this control volume is representative. The R E V is a unit that can be considered homogeneous and so must be large enough to average microscopic heterogeneities, but small enough to avoid incorporating macroscopic heterogeneities into the homogeneous approximation. The characteristic dimension of this 18 control volume must be much larger than that typifying microscopic structure, e.g. grain size, and much smaller than that of the flow domain under study. 2.3.2 Fundamental Equations The primary theory governing flow is the principle of continuity, which ensures that any physical quantity, such as mass or energy, is conserved in a given system. In a global sense, the physical quantity itself can neither be created nor destroyed. On a local scale, the principle of continuity ensures the balanced transfer of this physical quantity. In the same manner that the law of mass conservation is applied in the development of the continuity equation for fluid flow in a porous medium, the law of charge conservation is applied in the study of electrical current flow. Equations 2-4 and 2-5 describe the net flow rate of fluid or electrical charge across the boundaries of a unit control volume by means of an accumulation term and an external volumetric flow term. V .q = S v Q Transient hydraulic continuity (2-4) dt V. J = - / Transient electrical continuity (2-5) dt Equation 2-4 describes transient hydraulic flow conditions in a given control volume, where q is volumetric fluid flux [m3/s.m2], Ss is specific storage [m3/m3.m], h is hydraulic head [m of water], and Q is volumetric fluid flow [m /s.m ]. The accumulation term, Ss.dhldt, represents the time rate of change of the volume of fluid stored within the control volume, which is controlled by the compressibility of the soil matrix and pore fluid. The volumetric flow term, Q, describes any external source of fluid injected or withdrawn from the control volume, and is expressed as a function of this volume. Equation 2-5 describes transient electrical flow, where J is volume current density [A/m ], pe is volumetric charge density [C/m3], and / is current flow per unit volume [A/m3]. The accumulation term, dpjdt, describes the time rate of change of charge density or, more appropriately, the creation or dissipation of electrical current within the control volume. This charge accumulation or generation of current is caused by a variation in ionic mobility within the control volume, 19 which can also be described by a frequency-dependent (ac) bulk electrical conductivity. Any external current sources imposed on the control volume are represented by the volumetric flow term, /. A soil-fluid system subjected to an imposed gradient is forced out of static equilibrium, and macroscopic flow laws can be used to describe the irreversible flow phenomena that arise as a consequence of this gradient. Incorporation of these linear flow laws, namely Darcy's law and Ohm's law further develops the transient forms of the continuity equations. Equations 2-6 and 2-7 define hydraulic and electrical flow, where K is hydraulic conductivity [m/s], cr is electrical conductivity [S/m = 1/Ohm.m] and V is electrical potential [V = J/C]. V.(K.Vh) = Sx Q Transient hydraulic flow. (2-6) dt V.(o.VV) = —— - / Transient electrical flow (2-7) dt Under steady state flow conditions, no storage of fluid or charge can occur within the control volume, and the hydraulic and electrical accumulation terms must reduce to zero in order to uphold the principle of continuity. These conditions apply in the hydraulic case if it is assumed that the soil matrix is rigid and the pore fluid incompressible. In the study of electrical flow, steady-state conditions require no net accumulation of charge in the volume, which translates to an assumption of dc, or frequency-independent electrical conductivity. External sources of fluid or current may still be injected or withdrawn from the volume, provided that the net influx completely balances the net outflow. Equations 2-8 and 2-9 define steady-state hydraulic and electrical flow. V.(K.Vh) = -Q Steady-state hydraulic flow (2-8) V.(o.VV) = - / Steady-state electrical flow (2-9) If no external sources or sinks of fluid or electrical flow exist, the steady-state relations reduce to Equations 2-10 and 2-11. 20 V . (K .Vh) = 0 V.(o.VV) = 0 Steady-state hydraulic flow (no external sources) (2-10) Steady-state electrical flow (no external sources) (2-11) In order to solve the transient equation for a given flow system, boundary conditions that dictate flow behaviour at the periphery of the volume, and initial conditions that specify the potential distribution within the volume, must be defined. In a steady-state approach, only the conditions at the boundary need to be specified. A discussion of the fundamental relations and complete derivations of the continuity equations describing hydraulic and electrical flow can be found in Appendix B . 2.3.3 Coupled F low Theory In the same manner that Darcy's law may be used to describe fluid flow in a porous medium due to a hydraulic gradient, similar macroscopic laws can define flow induced by other gradients existing within the system. The flow of heat, and ionic compounds are described using linear flux equations commonly known as Fourier's law and Fick 's law, respectively. These flow systems, along with that described by Ohm's law, do not operate independently, and their interaction can be explained through the theory of coupled flow. Temperature, chemical or electrical gradients can induce fluid flow, and a hydraulic gradient may in turn induce the flow of heat, chemical compounds or electrical current. Consequently, the generalized macroscopic flow law can be broadened to include the influences of other potential gradients, VOy, where the cross-coupling coefficients, L,y, link these secondary gradients to the primary flux term, T„ as shown in Equation 2-12 (de Groot, 1951). Table 2-4 provides a summary of primary and secondary flows that may develop in a given system. O f the coupled flows listed, chemical and electro-osmosis, isothermal heat transfer, electrophoresis, and streaming current are of particular interest in the study of saturated soil systems (Mitchell, 1991). (2-12) 21 Table 2-4: Primary and coupled flow phenomena (adapted from Mitchell, 1991). FLOW, r GRADIENT, V Hydraulic Potential Temperature Electrical Potential Chemical Concentration ^ Fluid Darcy's Law (hydraulic conduction) Thermo-osmosis Electro-osmosis Chemical osmosis Heat Isothermal heat transfer Fourier's Law (thermal conduction) Peltier Effect Dufour Effect Current Streaming current (streaming potential) Seebeck Effect Ohm's Law (electrical conduction) Diffusion and membrane potentials Ionic (chemical) Streaming current Soret Effect Electrophoresis Fick's Law (diffusion) The importance of each component on the resultant flux is dependent on the relative magnitude of the gradient and the corresponding conductivity coefficient, such that for sufficiently slow processes, flow may result from not only the primary gradient, but also from the secondary gradients. In the study of coupled flow, the gradient terms are considered the thermodynamic forces and must be formulated properly in terms o f potential to uphold the principles of nonequilibrium thermodynamics (de Groot, 1951). Consequently, the hydraulic gradient must be formulated in terms of hydraulic potential O/, [N.m/kg], which describes the total mechanical energy per unit mass of fluid. Neglecting the effect of fluid velocity and assuming the fluid incompressible, hydraulic potential can be expressed in terms of hydraulic head as shown in Equation 2-13, where g is acceleration due to gravity [m/s ], z is elevation [m], u is pore fluid pressure with respect to atmospheric [Pa], and pw is fluid density [kg/m ] (Freeze and Cherry, 1979). 22 V < D / l = V ( g z + — ) = V ( g h ) ( 2 - 1 3 ) I n t h e a b s e n c e o f s i g n i f i c a n t t h e r m a l o r c h e m i c a l c o n c e n t r a t i o n g r a d i e n t s , t w o c o u p l e d f l o w e q u a t i o n s m a y b e u s e d t o d e s c r i b e e l e c t r o k i n e t i c p h e n o m e n a , as s t a t e d i n E q u a t i o n s 2 - 1 4 a n d 2 - 1 5 . T h = - L h . V O A - L f c . V O e Coupled hydraulic flow ( 2 - 1 4 ) r e = - L eh . V O h - L e . V O e Coupled electrical flow ( 2 - 1 5 ) T h e c r o s s - c o u p l i n g c o n d u c t i v i t y c o e f f i c i e n t s , L/,e a n d Le/,, l i n k t h e s e c o n d a r y g r a d i e n t s t o t h e p r i m a r y f l o w t e r m s , s u c h that t h e t e r m - L / , e V ® e d e s c r i b e s e l e c t r o - o s m o t i c f l u i d f l o w a n d -Le/,V 0 or — ->0 ( 2 - 1 6 ) ( 2 - 1 7 ) 2 3 The coupled flow equations may be re-stated using more familiar engineering parameters, by recognizing that hydraulic flux E/, represents the mass flux of fluid, and electrical flux Ye represents the flux of charge. Equations 2-18 and 2-19 describe volumetric fluid flux q and volume current density J, where kg [m/s per V/m] is the coefficient of electro-osmotic permeability and L [A/m ] is the streaming current cross-coupling conductivity coefficient. These coefficients are typically derived from one-dimensional laboratory tests specific to the study of electro-osmosis or streaming potential, and are discussed in further detail in Chapter 4. Appendix B contains a complete derivation of the coupled flow relations expressed in Equations 2-18 and 2-19, including a discussion of other forms of these equations commonly reported in the literature. q = ^ = - K . V h - k e . V V (2-18) J = T e = - L . V h - c .VV (2-19) In the study of the streaming potential phenomenon, the hydraulic flow equation may be decoupled and seepage analysed directly using Darcy's law. In the absence of external current sources, the hydraulic gradient drives fluid flow and induces electrokinetic effects, such that the electro-osmotic flow term represents a counter reaction to the primary flow of fluid. Mitchell (1991) reported negligible contribution of the electro-osmotic flow term in materials with K > 10\"9 m/s. Hence, fluid flow in cohesionless soils and most broadly graded glacial tills is governed by the primary hydraulic gradient, and Equation 2-18 is reduced to Darcy's law. For most geologic materials, the quantities described by Equation 2-17 are order 103 to 107 (Fitterman, 1978). Therefore, the contribution of the primary flow term to the overall fluid flow will be much greater than that of the coupled flow term, regardless of the hydraulic conductivity of the material. In the study of embankment dams, where large hydraulic gradients are imposed on the low permeability materials that comprise the core, an uncoupled approach is further justified. According to the convection current model proposed by Sill (1983), coupled electrical current flow can be described by a convection current or streaming current, induced by the hydraulic gradient, and a conduction current, induced by the resultant separation of 24 charge within the porous medium. Equation 2-19 can be expressed in terms of these two current densities, as shown in Equations 2-20 to 2-22. J = J . + J « * ( 2\" 2°) where, J 6 Y ; m , = - L . V h Convection (streaming) current (2-21) J c o n d = - o . V V Conduction current (2-22) For electrical continuity, the divergence of the total volume current density in the domain must equal the sum of the divergence of convection and conduction current densities. Consequently, Equations 2-23 and 2-24 describe the transient continuity relations for hydraulic flow and coupled electrical flow. V . (K .Vh)- S, ^--Q (2-23) ot V.(L.Vh) + V.(o\\VV) = ^ - / (2-24) dt 2.3.4 Governing Equations: Assessing the Importance of Transient Effects A steady-state analysis produces the maximum response obtained in a transient analysis once equilibrium is reached. The importance of transient terms can be assessed for a given problem by formulating a dimensionless form of the transient equation, and resolving the time scale or relaxation time (Beckie, 2000). The relaxation time, teq, for one-dimensional hydraulic flow is defined in Equation 2-25, where Ss,avg and Kavg describe specific storage and hydraulic conductivity averaged over a given dimension in the domain, and x is a representative length of that dimension. = Simg x civg 25 For times greater than the relaxation time, transient effects may be ignored and the problem solved using a steady-state analysis. Although certain hydraulic flow problems are more accurately represented by the transient continuity equation, in practice, most seepage analyses of embankment dams beyond the first filling of the reservoir are performed using a steady state approach, provided the reservoir elevation is fluctuating within a normal range of operation for given time period. This approach does not consider the time lag of the seepage front through the dam caused by a change in reservoir elevation. However, it provides a first approximation or \"snapshot\" of the hydraulic regime at a given reservoir elevation for the purposes of the coupled electrical flow analysis. Regardless of the nature of the hydraulic flow regime, steady state behaviour may be assumed in the analysis of electrical flow. Wurmstich (1995) concluded that electrical conduction current flow can be studied as a steady-state process, independent of the conditions governing the hydraulic system since hydraulic relaxation times are orders of magnitude larger than electrical relaxation times. In a physical sense, this means that the time required to reach equilibrium in a time-varying hydraulic system, considering the effects of consolidation for example, is much larger than the time required to reach equilibrium in an electrical system. The error in assuming steady-state conditions in the electrical problem can be assessed by means of the electrolytic polarization or induced-polarization (IP) effect. When an electrical current is imposed across a saturated soil resulting in a charge separation, the time required for the ions to return to their original positions can be quantified through time-dependent voltage measurements, or electrical conductivity measurements at different frequencies. Marshall and Madden (1959) reported a maximum frequency effect of less than 2% between ac and dc conductivity in studying the impact of electro-osmotic and streaming potential effects for a range of geologic materials. This result further substantiates the conclusion that neglecting the time-dependent accumulation of charge and assuming steady-state behaviour renders an acceptable solution to the electrical flow problem. The relations governing any steady-state flow analysis of the self-potential response to seepage are given in Equations 2-26 to 2-28. 26 V.(K.Vh) -Q (2-26) or, V.(K.Vh) 0 (2-27) V.(o.VV) -V.(L.Vh) (2-28) 2.4 Summary and Conclusions Chapter 2 introduced two seepage-induced physical phenomena occurring within embankment dams. The migration of soil particles through the process of internal erosion was explored to provide a basis for evaluating changes in soil and fluid properties that result from the onset of piping. The initiation of electrical current flow through the phenomenon of streaming potential was introduced as an electrokinetic coupled flow process. This phenomenon was explained in a physical sense at the pore scale, and mathematically in terms of linear macroscopic flow equations to reveal the link between hydraulic and electrical flow within a porous medium, so to provide a basis for a quantitative analysis of the phenomenon. Key points and conclusions include the following: Internal erosion is a prevalent mode of failure in embankment dams. Internal erosion and piping result in a loss of fine-grained material, which changes the structure of the soil matrix and can lead to an increase in porosity, water content, and seepage. These physical changes can influence the electrical properties of the eroded zone. A hydraulic gradient can induce electrical currents in a saturated porous medium through the electrokinetic phenomenon of streaming potential. Electrokinetic phenomena rely on the presence of an electrical double layer at the interface between soil particles and pore fluid, in which positive ions from solution balance the negative surface charge of the solid. In cohesionless soils, electrical current flow is concentrated within the pore fluid through the process of electrolytic conduction. Surface conductivity effects facilitate current flow at the solid-fluid interface, and predominate in clays, finer-grained soils, and soils saturated with fresh water. 27 Coupled flow theory is used to describe the link between hydraulic and electrical flow, and defines the streaming potential phenomenon in terms of a linear relation between streaming current and hydraulic gradient. In the study of the streaming potential phenomenon, the hydraulic flow equation may be decoupled, such that fluid flow may be completely described using Darcy's law. Neglecting the time lag induced by transient hydraulic boundary conditions, steady-state conditions may be assumed for the electrical flow problem, irrespective of whether the hydraulic flow problem is resolved in steady state or transient form. 28 3 FORWARD MODELLING Prediction of the self-potential response to seepage flow in a domain with known hydraulic and electrical parameters is a forward modelling problem. In the study of embankment dams, numerical methods are required to solve the coupled flow equations developed in Chapter 2 to tackle the complexities inherent in these structures. A mathematical model, comprised of governing flow equations and boundary conditions, and discretization of the flow domain are required for each numerical modelling application, given system-specific input parameters. Chapter 3 describes the numerical procedure 3 DSP developed to examine the SP response to seepage. Section 3.1 outlines previous modelling studies, which have taken different approaches to the coupled flow problem. In the present study, a three-dimensional finite difference approach is applied. The commercially available software Visual MODFLOW is used in conjunction with a series of coded algorithms to solve the coupled flow problem. Section 3.2 describes the finite volume method of solution, and Sections 3.3 and 3.4 detail model considerations and boundary conditions specific to the analysis of the hydraulic and electrical flow systems. 3DSP is verified against an analytical solution to a three-dimensional flow problem in Section 3.5. 3.1 Previous Work The self-potential method can be used to characterize naturally occurring electrical currents in the ground generated by chemical, thermal or hydraulic flow. These applications have led to the development of interpretive modelling techniques that rely on coupled flow theory to describe the resulting electrical potential distribution. Two coupled flow relations describing fluid and electrical current flow were stated as Equations 2-18 and 2-19 in Section 2.3.3, pertaining to the present investigation. Previous studies of the streaming potential phenomenon have focused on electrical flow coupled with fluid or thermal flow, where the primary hydraulic or thermal flux was solved using an uncoupled approach, and electrical flow analysed using a general form of Equation 2-19, which is stated here as Equation 3-1 with secondary potential gradient representing the hydraulic or temperature gradient. 29 J = - L 2 1 . V ^ - o.VV Coupled electrical flow (3-1) Two main approaches to the coupled electrical flow problem described by Equation 3-1 have been explored. Nourbehecht (1963) and later Fitterman (1978) introduced the total potential approach, which combines the two potentials into a single pseudo potential, \\> to describe electrical flow, as shown in Equation 3-2. J = -c.Vvj/ where, \\|/ = V + (3-2) a The total potential method solves for \\\\i at the surface of the earth, and relates it to electrical potential by assuming the driving potential, (f> is constant at the surface. This specified potential boundary condition is acceptable if represents temperature, but is not realistic in the study of hydraulic flow where a zero fluid flux boundary is more appropriate for a saturated subsurface. A second approach to the coupled electrical flow problem proposed by Sill (1983) describes electrical flow in terms of a convection current induced by the primary gradient, and a resultant conduction current, corresponding to the two terms on the right-hand side of Equation 3-1. This method, which was introduced in Equations 2-21 and 2-22, approximates our understanding of the physical phenomenon more effectively and results in the direct solution of electrical potential, V. The convection current approach is implemented in the present study as it correctly models the electrical potential anomaly at surface in response to a single point primary fluid flow source, such as a pumping well, in a uniform subsurface (Sill, 1983). This numerical result is equivalent to the physical response measured in practice (Sato et al, 1998), which is not predicted correctly using the total potential approach (Fitterman, 1978). The original algorithm developed by Sill (Sill and Killpack, 1982) solves the 3-D steady-state problem described by Equations 2-26 and 2-28, formulated in terms of hydraulic pressure, using a nodal-based 2-D finite difference approach in the program SPPC (Wilt and Butler, 1990). Discrete volumetric fluid flow sources are required as input for solution to the 30 hydraulic problem, and zero flux boundaries assumed at the edge of the domain, such that positive and negative 3-D flow sources must completely balance one another within the flow domain. Consequently, a seepage regime that satisfies Equation 2-27 in a given earth structure may not be realistically approximated given that specified potential boundaries cannot be imposed. SPPC is further limited by restrictions to mesh size and geometry, which can contribute significant error in the approximation of complex heterogeneous systems or problems of irregular geometry. The convection current approach has since been applied in most numerical investigations of the streaming potential phenomenon, which have focused primarily on geothermal applications (Yasukawa, 1993; Ishido and Tosha, 1998). Wurmstich (1995) developed a 3-D finite difference algorithm for analysis of streaming currents generated by oil well pumping (Wurmstich and Morgan, 1994) and performed a 2-D finite element analysis of a simple case of controlled seepage through a model embankment to resolve the cross-coupling coefficient, L of the fill material (Wurmstich et al, 1991). Berube (2000) introduced a 2-D finite element code for analysis of the streaming potential phenomenon due to seepage. The present investigation differs from these previous studies in applying a three-dimensional approach to the coupled hydraulic and electrical flow problem in embankment dams, with emphasis on an accurate seepage analysis and the characterization of physical properties from a geotechnical perspective. 3.2 Finite Volume Formulation 3.2.1 Governing Equations and Generalized Approach The continuity equations relevant to the study of streaming potential were derived in Chapter 2 using the theory of coupled flow and the convection current approach proposed by Sill (1983). The rationale for using a steady-state approach to both the seepage and electrical flow problems was discussed in Section 2.3.4, and the governing relations repeated here in Equations 3-3 and 3-4. V.(K.Vh) = -_> or V.(K.Vh) = 0 (3-3) 31 V.(c.VV) = - / _ , , , „ where, IamvMlon = V.(L.Vh) (3-4) In order to solve these equations numerically, we must apply discrete approximations to the equations and the flow domain, assign system-specific boundary conditions and specify the distribution of physical properties. The flow domain is characterized by a series of three-dimensional discrete volumes, each of which is assigned physical properties required for solution, namely hydraulic, cross-coupling and electrical conductivities. Hydraulic head and electrical potential are resolved from the discrete equations, which are satisfied in each discrete volume comprising the flow domain. Discrete volumes situated at physical boundaries internal or external to the domain are subject to conditions governing flow at the boundary, which are incorporated into the corresponding discrete flow equation. With specified hydraulic boundary conditions and hydraulic conductivity values, the hydraulic head distribution can be resolved using Equation 3-1. If external sources of fluid flow are imposed on the domain, in-situ examples of which include a pumping well or surface recharge, the steady-state equation containing the volumetric flow term Q is applied to the appropriate discrete volumes encompassing these sources. Otherwise, Q is equal to zero throughout the domain. Based on the computed distribution of head and specified cross-coupling conductivity values, convection or streaming current sources are calculated within each discrete volume throughout the domain using the relation V.(L.Vh). In a steady-state analysis, these streaming current source terms serve as volumetric conduction current sources for the solution to the conduction current, or dc electrical resistivity problem, described by Equation 3-4. Resolution of the electrical potential in each discrete volume is achieved with a specified distribution of electrical conductivity and restrictions to conduction current flow imposed at the external boundary of the flow domain. 3.2.2 Finite Difference Equations The formulation of a numerical solution to the coupled flow problem requires a discrete approximation to the continuous flow equations shown in Equations 3-3 and 3-4. The method of finite differences was chosen based on the availability and compatibility of software and computer code, and the capacity for modelling three-dimensional flow problems through a cell-centred finite volume approach. 32 The finite volume method relies on the principle of continuity applied to a series of control volumes that comprise the flow domain. The flow domain is discretized to form a grid or mesh, as outlined in Figure 3-1, which shows the finite volume cell and node configurations for two- and three-dimensional systems. Nodes are located at the cell centres and are representative of the entire grid cell volume. The conductivity, potential, and volumetric flow source term for a given cell are defined at the node, whereas gradient and flux terms are represented at the cell boundaries. Figure 3-1: Finite volume rectilinear grid configurations: a) finite volume grid cell; b) 2-D nodal arrangement; c) 3-D nodal arrangement. Discrete equations are derived from the continuous flow relations by integrating over a discrete volume, using the central-difference approach for a second order approximation to the integrals and first derivatives. The discrete approximation of the one-dimensional form of Equation 3-3 is illustrated in Equation 3-5, and is shown for the three-dimensional case in Equation 3-6. f | < * , £ ) & = - J\"\" fi& =* K J ^ ± + K J ^ = -QAc, (3-5) J*wccfc dxc J j c w C Axwc Axec 33 A*wc ^ e c AV„C ^ . ^ A x c A , c + ^ c . ^ A x c A , c + ( 3 . 6 ) A j s c A z f c ^ b c - ^ T ^ ^ c = - 2 -Ax; c Ay c Az c The discrete equations must be applied to each cell in the flow domain, and are therefore more efficiently solved in matrix form. Equation 3-5 describing 1-D fluid flow can be stated as a matrix equation, as illustrated in Equation 3-7. This relation can be stated in more general terms to describe both hydraulic and electrical flow in 3-D as shown in Equation 3-8, where D and G are difference matrices that impose the central difference approximations, representing the divergence and gradient operators. Hydraulic or electrical conductivity terms are contained in the conductivity matrix M, and the potential vector [u] and source matrix S represent hydraulic or electrical potential and volumetric source terms, respectively, for each cell within the domain. -1/Ax 1/Ac K wc 0 0 1/Axw c 1/Axw c 0 0 -1/Ax 1/Ax, K K h. DMG[u] = S (3-7) (3-8) One benefit of the cell-centred finite volume approach is the ease with which it handles physical property discontinuities. The conductivity at the interface between two adjacent cells is described by the harmonic average of the independent cell coefficients, which is computed by taking the integral of the inverse conductivity, or resistivity term over the volume spanning the interface. Equation 3-9 defines the interface hydraulic conductivity between cells 'w' and 'c', as displayed in Figure 3-1. 34 2KwKcAx wc Interface conductivity (3-9) A number of references deal with the subject of numerical analysis and the finite volume method. The reader is directed to publications by Remson (1971) and de Marsily (1986) on the application to hydrogeology, and by Gerald and Wheatley (1989) and Golub and Ortega (1992) for a generalised approach. 3.2.3 Numerical Solution Techniques The discrete flow equations and imposed boundary conditions together form a system of equations for all cells in the domain, which are organized in matrix form, as shown in Equation 3-7 and simplified here as Equation 3-10. The sparse coefficient matrix A, which depends on the conductivities and dimensions of the grid cells, links the potential vector [u], to the source term vector S. This linear system of equations is solved for both hydraulic and electrical flow systems. Direct methods of solution, such as Gaussian elimination or LU decomposition, provide an \"exact\" solution through back-substitution and are useful for small or medium scale problems, but computational memory constraints render these methods inefficient for larger scale multi-dimensional systems. Iterative solution techniques use very little memory, but raise the issue of convergence of the calculated and exact solutions with consecutive iterations. Stationary iterative techniques, such as the Jacobi and Gauss-Seidel methods, successive over relaxation (SOR) and symmetric successive over relaxation (SSOR) methods, require many iterations for convergence and therefore may not be more efficient than direct methods. However, these methods may be used for pre- or pre-post conditioning of the matrices prior to solution using faster non-stationary or Krylov space iterative techniques, such as the conjugate gradient and bi-conjugate gradient methods. A [u] = S (3-10) 35 Regardless of the numerical technique used to solve Equation 3-10, an inexact solution [u] results due to truncation and round-off errors, which arise as a consequence of the finite approximations to infinite equations and numbers, respectively. Decreasing the grid cell spacing can minimize truncation error resulting from the finite difference approximations. Round-off error can be reduced using double precision arithmetic, and through the use of iterative methods over direct methods of solution, which prevent the accumulation of round-off error with successive approximations. However, the condition number of the matrix A dictates the sensitivity of the solution to round-off error, such that for large condition numbers, A is ill-conditioned and the problem is more difficult to solve. Consequently, matrix preconditioning can be performed to lower the condition number, thereby reducing the influence of round-off error and the number of iterations required to reach convergence. In absolute terms, convergence is achieved once the residual or difference between solutions from two successive iterations is within a specified tolerance. If the magnitude of the potential distribution is unknown for a given problem, a relative residual term is used to assess convergence. Further information on direct and iterative solution methods can be found in the references listed previously in Section 3.2.2. 3.3 Seepage Flow Analysis 3.3.1 Mesh Geometry and Boundary Conditions The flow domain dimensions and the grid cell mesh are designed in Visual MODFLOW version 2.7.2, a software package developed by Waterloo Hydrogeologic Inc. (Guiguer and Franz, 1997). Symmetric and irregular grids can be generated, and although a rectilinear mesh is required, zones of inactive cells can be delineated to approximate a domain with an irregular geometry. Typically, the finite element method is preferred over finite differences for problems with complex geometries, since it provides more flexibility in grid design. However, the advantages of finite element models can be at the expense of fundamental issues, such as satisfying the law of mass conservation (Beckie, 2001). In dam seepage studies, sloped physical property boundaries must be approximated by a stepped surface. If a mesh spacing is chosen such that cell dimensions are smaller than 36 the R E V o f a given soil unit, the discrete approximation should not affect the final solution. The use of non-uniform mesh spacing enables more detailed resolution in zones where large changes in flow conditions are expected. The maximum representative cell size for a given zone can be determined by comparing solutions obtained on a series of coarser grids. Substantial error between two successive solutions implies that the mesh is too coarse to adequately model flow conditions. Although hydraulic boundary conditions for a given domain are problem-specific, several common boundaries apply to the study of earth dams. If the foundation of the embankment is considered impermeable, a type of specified flux, or Neumann boundary, may be imposed to restrict flow across this interface. The segment ae in Figure 3-2 represents the base o f an embankment, where the no-flow boundary dh/dn = 0 is applied. Figure 3-2: Hydraul i c boundary conditions. Seepage analyses require consideration of the hydraulic gradient imposed by the reservoir, and the downstream outflow boundary, or seepage face. A Dirichlet or specified-potential boundary condition is used to define hydraulic head along the upstream face corresponding to a given reservoir level h i , as shown by segment ab in Figure 3-2. In M O D F L O W , this type of boundary condition is imposed at the cell node. Fluid flow resulting from this imposed gradient must exit the dam by means of a free-outflow boundary, coinciding with the downstream face of the dam or the interface with a downstream toe drain. This seepage face, illustrated by segment cd in Figure 3-2, is defined by means o f a double or mixed boundary condition, where specified head and flux conditions, h = z and dh/dn < 0 , dictate fluid outflow in a direction perpendicular to the interface. If standing 37 water exists downstream of the dam, a Dirichlet boundary may be used to specify the tailwater level, as shown by segment de where h = I12. The use of drain cells in MODFLOW is the most appropriate boundary condition for simulation of a seepage face (Anderson and Woessner, 1992). Cells lining the outflow boundary are assigned as drains, and are activated if they are found to lie below the computed phreatic surface. The volumetric flow rate is governed by a conductance term, described in Equation 3-11, where the conductance c [m /s] is defined in terms of vertical hydraulic conductivity Kz. The dimensions Av and Az define the cross-sectional area of the seepage path entering the drain, which are equivalent to the cell dimensions for a drain elevation corresponding to the base of the grid cell, and Ax describes the seepage length, typically defined as half the cell length. The validity of the seepage face approximation is verified by ensuring that the calculated hydraulic head is equivalent to elevation in each cell along the seepage face. In this way, the law of mass conservation is ensured across the domain. Fluid entering the system from the reservoir is allowed to exit at a rate controlled by the hydraulic conductivity of the soil unit at the downstream toe. K.AyAz c= - / (3-11) Ax 3.3.2 Free-Surface Approach In an unconfined or saturated-unsaturated flow system, the position of the phreatic surface with respect to outflow boundaries, such as the downstream seepage face in a dam, must be determined. The phreatic surface, indicated by segment be in Figure 3-2, delineates the boundary along which pore fluid pressure is atmospheric, such that hydraulic head is equal to elevation. A zone exhibiting negative pore pressures, or capillary fringe, extends above the phreatic surface, across which flow may occur. However, in many engineering applications, fluid flow is assumed to dominate below the phreatic surface and a saturated flow or free-surface approach is applied. With this assumption, the phreatic surface represents the top flow line, and a mixed boundary condition of h = z and dh/dh = 0 may be applied. 38 In seepage analyses performed with transient boundary conditions, such as a fluctuating reservoir, the rate of response of the phreatic surface is highly dependent on the unsaturated soil properties and a free-surface approach to solution can lead to erroneous results (Freeze, 1971). However, in a steady-state analysis, the effect of unsaturated fluid flow can be neglected such that a fully saturated free-surface model may realistically simulate flow conditions (Romano et al, 1999). MODFLOW calculates the position of the phreatic surface using Dupuit-Forchheimer theory in a free-surface approach. Central to the theory are the assumptions that the hydraulic gradient is equivalent to the slope of the phreatic surface and does not vary with depth, and that the vertical flow component can be ignored, such that the 3-D Darcy flow equation can be resolved using a 2-D form of the equation, as shown in Equation 3 - 1 2 . These assumptions are valid if the slope of the phreatic surface is small, which is satisfied in a finite difference analysis if many layers are incorporated into the mesh such that only the uppermost grid cells require analysis using Dupuit-Forchheimer theory. However, near steeply dipping or vertical outflow boundaries where the vertical flow component figures prominently, the Dupuit assumptions may lead to an inexact approximation of the location of the phreatic surface. An illustration of the analysis of unconfined flow from both a saturated-unsaturated and a saturated approach is shown in Figure 3 - 3 . KMx, y)^ + Kyh(x, y)^ = ~Q ( 3 - 1 2 ) 3.3.3 Solution Method Within Visual MODFLOW, the user may specify several solution options, including solver type and convergence criteria. The iterative bi-conjugate gradient stabilized method using Stone incomplete decomposition for preconditioning was chosen for all seepage analyses in MODFLOW. Upon convergence of the solution, zero mass balance was verified for the steady state system, and the hydraulic head distribution exported for use in the electrical analysis. 3 9 •>-x Figure 3-3: Development of a seepage face on a free-outflow boundary: a) saturated-unsaturated flow net; b) free-surface flow net; c) Dupuit-Forchheimer flow net (from Freeze and Cherry, 1979). 3.4 Electrical Flow Analysis 3.4.1 Mesh Geometry and Boundary Conditions The study of steady-state hydraulic flow using the free-surface approach was discussed in Section 3.3.2. Since streaming current sources are induced by fluid flow through porous media, convection current flow is limited to the saturated zone beneath the phreatic surface, neglecting the effects of flow in the unsaturated zone. Consequently, convection current flow is constrained by means of the hydraulic no-flow boundary at the free surface, and volumetric convection current source terms are calculated for each grid cell 40 within the saturated zone of the fluid flow domain, which is illustrated as the area bound by abce in Figure 3-2. Conduction current flow is the direct result of streaming or convection current source terms within the saturated subsurface, as described by Equation 3-4. However, unlike the flow of streaming current, conduction current flow is unlimited within the fluid flow domain. The electric field, or voltage gradient that defines the conduction current, decays with distance from the distribution of positive and negative volumetric convection current sources within the hydraulic domain. Consequently, a zero-flux Neumann boundary can be applied to restrict electrical flow at the external edge of the domain, such that dWIdh = 0 on the edges of grid cells delineating the boundary. However, this boundary condition must not influence the electrical flow solution within the hydraulic domain. In order to avoid boundary effects, the finite volume mesh must be expanded to define an electrical flow domain. Since the electrical potential matrix, represented by [u] in Equation 3-10, is under-determined, a known potential, or V=0 at the external boundary, may also be specified in one grid cell. The electrical domain is constructed using MATLAB code by padding the original hydraulic mesh with cells at an incremental expansion rate of 1.3. The appropriate size of the electrical domain should be established by trial and error; a satisfactory solution exists when the electrical potential distribution within the bounds of the hydraulic domain remains constant with increasing domain size. A general guideline limiting the size of the padded grid is that the cell aspect ratio, or the ratio of cell dimensions, should not exceed 10 (Haber, 2001). 3.4.2 Solution Method Although similar iterative solution techniques are used in the fluid and electrical flow problems, the hydraulic system is analysed using Visual MODFLOW, and the electrical system is solved using a series of algorithms programmed using the MATLAB interface. Conductivity matrices of cross-coupling conductivity L and electrical conductivity cr are generated for the fluid domain mesh designed in MODFLOW. Streaming current volumetric source terms are calculated directly for each saturated cell using Equation 3-8, where the potential vector [u] and conductivity matrix M describe hydraulic head and cross-41 coupling conductivity values, respectively. The electrical flow analysis is then performed on the padded mesh to resolve electrical potential, where the coefficient matrix A in Equation 3-10 defines electrical conductivity terms, and the source term vector S defines conduction current sources, which are equivalent to convection current source terms of inverse sign. The iterative bi-conjugate gradient method in conjunction with an S S O R preconditioner was used to solve the electrical problem. Coded algorithms specific to the generation of the difference matrices, harmonic average calculation and iterative solver were provided courtesy of Eldad Haber of the Department of Earth and Ocean Sciences at the University of British Columbia. 3.5 Validation of the Numerical Procedure 3DSP Agreement between the numerical and analytical solutions to a three-dimensional coupled flow problem provides verification that the numerical procedure 3DSP can adequately model the steady-state SP response to seepage. The problem of a point fluid flow source within a homogeneous medium, which represents a single injection well screened over a small finite portion of its length, was chosen for its 3-D character and simplicity. In the hypothetical field-scale problem, observation points are defined by their radial distance from the screened portion of the well, as shown in Figure 3-4. obse rva t i on point it: z X Fluid flow It Figure 3-4: Injection well in a homogeneous subsurface. 42 Equation 3-13 gives the potential distribution at a radial distance, r, from a point flow source in a homogeneous medium (Telford et al, 1990). The same equation can be used to evaluate both hydraulic and electrical potential distributions, where u represents hydraulic head h or electrical potential V; S represents volumetric fluid flow rate Q, or streaming current or conduction current source strength /; and C represents hydraulic conductivity K or electrical conductivity a. Since the medium is homogeneous, the streaming current source term can be formulated in terms of the hydraulic flow source through Poisson's equations describing fluid and convection current flow, as shown in Equation 3-14, where L is the cross-coupling conductivity coefficient. u = — (3-13) 47trC I , =-I =— (3-14) conduction convection ^ v ' Due to the presence of a boundary between the subsurface and atmosphere, the hydraulic and electrical flow fields will be distorted near surface. To account for this effect in the analytical solution, the method of images must be employed. The potential distribution is now a function of both the buried source and an image source of equal magnitude positioned an equal distance above the ground surface, as shown in Figure 3-5. The effect of the image source can be added, as shown in Equation 3-15, where the reflection coefficient kr describes the influence of the boundary on the potential distribution. When observations are made at the ground surface, kr is unity and Equation 3-15 reduces to Equation 3-16. u = S ( l k ^ where k=^-^- (3-15) 4TIC v r r' j C , + C 2 u z = 0 2nr C (3-16) 43 Figure 3-5: The method of images. A fully saturated sand layer was modelled under steady-state flow conditions with a hydraulic conductivity of 5x l0\" 4 m/s, an electrical conductivity of 5x l0\" 4 S/m (5 pS/cm), and 5 2 a streaming current cross-coupling coefficient of 7.5x10\" A / m . The fluid injection rate was chosen as 0.005 m 3/s (432 m 3/day) at a depth of 15 m. Hydraulic head and self-potential distributions were resolved in the fluid domain using a non-uniform 37x37x37 cell mesh. Axisymmetric profiles o f the analytical and numerical solutions for hydraulic head and self-potential along the ground surface in the x-z plane of the well are shown in Figure 3-6. The maximum error between the two solutions is less than 1%. 44 5 analytical Distance from well (m) a) hydraulic head distribution -7 -I , , , , , , , ! -20 -15 -10 -5 0 5 10 15 20 Distance from well (m) b) electrical potential distribution Figure 3-6: Surface response from an injection well. 45 3.6 Summary and Conclusions Chapter 3 satisfied the first objective of this thesis by introducing the 3-D finite difference procedure 3DSP as a numerical tool to predict the self-potential response to fluid flow based on a comprehensive and realistic seepage analysis. A forward modelling approach is employed to calculate the hydraulic head distribution resulting from specified boundary conditions and the distribution of hydraulic properties across the flow domain. The subsequent steady-state electrical flow analysis relies on specified electrical flow properties, and iteratively evaluates the distribution of electrical potential throughout the flow domain based on the computed hydraulic head values. The mathematical electrical flow model describes the generation of convection current flow due to laminar fluid flow in a porous medium, and the resulting feedback flow of conduction current. Key conclusions and contributions include the following: • A comprehensive analysis of steady-state hydraulic flow conditions offers a solution to Equation 3-3, where flow may originate from zones of specified hydraulic head and/or unique volumetric fluid source terms within the domain. • The cell-centred finite volume method enables a 3-D representation of flow conditions and physical properties, and effectively represents physical boundaries through the process of harmonic averaging. • Careful mesh design and the correct use of boundary conditions are important in achieving an accurate and representative numerical solution. • The incorporation of the flow software Visual MODFLOW allows flexibility in grid design, and enables transient or steady-state hydraulic flow analyses to be performed. • Possible limitations of MODFLOW include the saturated flow approach to solution, which neglects the influence of flow in the unsaturated zone. The use of Dupuit-Forchheimer theory to calculate the position of the phreatic surface may lead to inaccurate results near steeply sloped or vertical seepage boundaries, since the vertical component of flow is neglected. • The flow domain of interest must be expanded for the purposes of the electrical flow analysis to impose the required Neumann and/or Dirichlet boundary conditions. 46 3DSP is a self-contained numerical tool that generates the SP response to seepage, requiring the user to specify only grid geometry and a mesh spacing that define the hydraulic flow domain under study, hydraulic boundary conditions, and the distribution of hydraulic and electrical properties within the domain. Verification of the numerical procedure 3DSP was achieved through excellent agreement with a closed-form solution. 47 4 HYDRAULIC AND ELECTRICAL PROPERTIES Valid characterization of the physical properties that control hydraulic and electrical flow in the subsurface is the main limitation to any forward modelling procedure. Three physical parameters are required to solve the steady-state or transient continuity equations stated in Chapter 2. Hydraulic conductivity K, bulk electrical conductivity cr, and streaming current cross-coupling conductivity L, must be specified for each discrete cell in the flow domain as input to the numerical modelling procedure described in Chapter 3. Hydraulic conductivity and electrical conductivity have been studied extensively through in-situ and laboratory methods of measurement. Investigation of the cross-coupling conductivity coefficient is limited to experiments focused on electrokinetic behaviour, and consequently a small data set exists. Section 4.1 summarizes typical values of each physical property for a range of soil types, and examines several common empirical and theoretical relations that may be used to estimate these parameters. In order to study the effect of changing seepage conditions on the self-potential response, the influence of internal erosion on the physical properties of the embankment material must be assessed. Stemming from the Chapter 2 discussion of structural changes that occur through the process of internal erosion, Section 4.2 seeks to quantify changes in the model input parameters resulting from the onset of internal erosion, and assess the sensitivity of the SP response to these approximations. 4.1 Physical Parameters 4.1.1 Hydraulic Conductivity Hydraulic conductivity is a property of both the soil particles and pore fluid comprising the porous medium, and is dependent on particle size and shape, void ratio, amount and mineralogy of the fine fraction, degree of saturation, and temperature (USACE, 1986). Hydraulic conductivity can be described for a saturated uniform soil of grain diameter d [m] by Equation 4-1, where p. is dynamic fluid viscosity [Pa.s], k is the intrinsic permeability of the soil [m2], and c is a coefficient of proportionality that describes the grain size distribution, roundness and packing. Hydraulic conductivity values for soils range 48 approximately from 10\"2 m/s for coarse gravels to 10\"11 m/s for clays, with cohesionless soils alone spanning seven orders of magnitude. Table 4-1 lists hydraulic conductivity values for a range of earth materials, along with unit conversion factors. K = in [m/s], where k = cd2 (4-1) Table 4-1: Typical hydraulic conductivity values (from Domenico and Schwartz, 1998). Hydraulic Conductivity, K Material (m/s) Gravel 3 x l 0 - 4 - 3 x l 0 - 2 Coarse sand 9 x l 0 \" 7 - 6 x l 0 \" 3 Medium sand 9 x l 0 \" 7 - 5 x l 0 \" 4 Fine sand 2 x l 0 _ 7 - 2 x l 0 - 4 Silt, loess I x l 0 \" 9 - 2 x l 0 \" 5 Glacial till I x l 0 \" 1 2 - 2 x l 0 ' 6 Clay l x l O ' n - 5 x l ( V 9 Unweathered marine clay 8x lO\" ' 3 -2x lO\" 9 Limestone, dolomite I x l 0 \" 9 - 6 x l 0 ' 6 Sandstone 3 x l 0 \" ' ° - 6 x l 0 \" 6 Unfractured igneous and metamorphic rocks 3 x l 0 \" ' 4 - 2 x l 0 ' 1 0 To convert m/s to: Multiply bv: cm/s 102 ft/s 3.28 gal/day.ft2 2.12xl0 6 cm 2 IO\"3 Characterization of the hydraulic conductivity distribution in a particular formation should be performed using direct measurement techniques. However, it is of interest to examine the dependency of this parameter on soil properties, such as grain size distribution and soil fabric, in order to calculate estimates based on other measured parameters. Efforts 49 to relate hydraulic conductivity to grain size and/or porosity have resulted in several empirical relations, such as Hazen's formula, the Masch and Denny method, and the Kozeny-Carman equation. These formulae are suitable only for saturated cohesionless soils and provide general estimates for a given grain size distribution (USACE, 1986). Hazen's formula, shown here as Equation 4-2, was derived from experiments conducted on relatively uniform sands exhibiting a coefficient of uniformity, C u < 5, and an effective grain size 0.01cm « dio « 0.03cm (Taylor, 1948). This formula makes no allowance for changes in void ratio or the shape of soil particles, and the coefficient represents the average of many values ranging from 81 to 117. Consequently, Hazen's formula is best used to provide estimates for uniform sands, but may be applied to give rough estimates of K for soils in the fine sand to gravel range defined by d\\o [cm] (Taylor, 1948). A: [cm/s] =100d,20 (4-2) To account for variations in grain size distribution, Masch and Denny (1966) defined an inclusive standard deviation, 07 to measure the spread of the gradation curve. This parameter is applied in conjunction with the mean grain size, t/50 of the soil to estimate hydraulic conductivity using graphical curves derived from tests performed on dense sands. As porosity is a function of the grain packing arrangement, many researchers have attempted to define the coefficient c in Equation 4-1 in terms of porosity. The Kozeny-Carman equation (Bear, 1972), shown in Equation 4-3, relates hydraulic conductivity K [m/s] to a representative mean grain size d„, [m], which is not strictly defined. This relation was derived from experiments conducted on uniform sands at various densities. K = PWS P n3 180 (4-3) 4.1.2 Electrical Conductivity A variety of surface and downhole electrical methods, such as resistivity and electromagnetic techniques, may be used to characterize the distribution of electrical 50 conductivity in the subsurface, with resistivity methods most commonly employed in the analysis of laboratory samples. Bulk electrical conductivity is a function of porosity and the connectivity of the pore space; degree of saturation and fluid conductivity; and the electrical conductivity, mineralogy and grain-size distribution of the soil matrix. The influence of saturation and pore connectivity has guided many researchers to attempt correlation between measurements of electrical conductivity or resistivity and hydraulic conductivity in a given formation (Jones and Buford, 1951; Kelly, 1977; Mazac et al, 1985). Pore fluid electrical conductivity is a primary factor influencing the bulk conductivity of porous media. The principle of electrolytic conduction was introduced in Section 2.2.3, where the number of ions present in solution and the mobility of these ions were stated as the two factors controlling fluid conductivity. Increasing temperature results in a 2.2% to 2.5% increase in fluid conductivity per °C for natural water, as shown by Equation 2-1 in Section 2.2.2. Natural surface and ground waters exhibit a wide range of electrical conductivities, from <10 pS/cm for fresh water to 50,000 u.S/cm for sea water (Telford et al, 1990). Table 4-2 lists a range of fluid conductivities for water in contact with sediments, along with typical values from reservoirs in British Columbia. Table 4-2: Typical fluid electrical conductivity values (from Telford et al, 1990; 2R.L.&L. Environmental Services Ltd., 1999; 3BC Hydro, 2001). Electrical Conductivity, ay Fluid (uS/cm) Surface water (sediments) 100- 1000 Ground water (sediments) 100- 10,000 Distilled water <1 Brilliant Dam headpond (all seasons)2 127-172 Mica Dam reservoir (October)3 139 Terzaghi Dam reservoir (August)3 85 To convert uS/cm to: Multiply by: S/m, or 1/Ohm.m IO\"4 In addition to the influence of fluid conductivity, surface conductivity effects at the solid-fluid interface and, to a lesser extent, the conductivity of the soil matrix contribute to 51 t h e o v e r a l l e l e c t r i c a l c o n d u c t i v i t y . S u r f a c e c o n d u c t i v i t y e f f e c t s p r e d o m i n a t e i n f i n e r - g r a i n e d m a t e r i a l s a n d i n s o i l s s a t u r a t e d w i t h v e r y f r e s h w a t e r , w h i c h e x h i b i t s a l o w c o n d u c t i v i t y r a n g e d u e t o s m a l l i o n i c c o n c e n t r a t i o n s . T h e p r e s e n c e o f c l a y m i n e r a l s f u r t h e r c o n t r i b u t e s to t h e b u l k c o n d u c t i v i t y o f a s o i l b y i n c r e a s i n g t h e c o n d u c t i v i t y t h e s o i l m a t r i x . T y p i c a l r a n g e s o f e l e c t r i c a l c o n d u c t i v i t y f o r d i f f e r e n t s o i l t y p e s w i t h v a r y i n g d e g r e e s o f s a t u r a t i o n a n d w a t e r c o n d u c t i v i t y are s h o w n i n T a b l e 4 - 3 , b o t h i n u n i t s o f e l e c t r i c a l c o n d u c t i v i t y a n d r e s i s t i v i t y . Table 4-3: Typica l electrical conductivity and resistivity values (adapted from Culley et al, 1975). Bulk Electrical Bulk Electrical Conductivity, a Resistivity, p Material ((iS/cm) (Ohm.m) Clay 100- 10,000 1 --100 Loam 140 --1400 7 -70 Organic silt 33 - 1500 7--300 Clayey soils 15- -100 100 -670 Sandy soils 1 -14 710- 10,000 Loose sands 0.1 - 10 1000- 100,000 Alluvial sands and gravel 1 - 100 100- 10,000 Glacial till 1 - 1000 10- 10,000 Limestone 2 - 140 70- -5000 Sandstone 1 - 300 33- 10,000 Crystalline rocks 0.01 - 10 1000- 1,000,000 To convert uS/cm to: Multiply by: S/m, or 1/Ohm.m IO-4 A n u m b e r o f m o d e l s h a v e b e e n d e v e l o p e d f o r e s t i m a t i n g t h e b u l k e l e c t r i c a l c o n d u c t i v i t y o f a p o r o u s m e d i u m b a s e d o n the e l e c t r i c a l c o n d u c t i v i t y o f t h e s a t u r a t i n g f l u i d a n d g e o m e t r i c a l f a c t o r s . T h e s e m o d e l s are a l l b a s e d o n t h e e m p i r i c a l r e l a t i o n k n o w n as A r c h i e ' s l a w , w h i c h d e s c r i b e s t h e b u l k c o n d u c t i v i t y o f a f u l l y s a t u r a t e d c o h e s i o n l e s s d e p o s i t , i n w h i c h t h e c o n d u c t i v i t y o f t h e s o i l m a t r i x i s a s s u m e d n e g l i g i b l e ( A r c h i e , 1 9 4 2 ) . C o n s e q u e n t l y , t h i s r e l a t i o n i s m o s t s u i t a b l e f o r m e d i u m t o c o a r s e r - g r a i n e d s o i l s . E q u a t i o n 4 -4 i s a r e p r e s e n t a t i o n o f A r c h i e ' s l a w , w h e r e cr i s b u l k e l e c t r i c a l c o n d u c t i v i t y , n i s p o r o s i t y , oy 5 2 is the fluid conductivity, and m is the cementation exponent, which describes the connectivity of the pore space and is equivalent to 1.3 for soils. For partially saturated sands, Archie's law can be expressed as Equation 4-5, where S is the degree of saturation (Telford et al, 1990). o~ = nmaf (4-4) CJ = S2nmCTf (4-5) The Waxman-Smits model is the most widely accepted approach used to estimate the electrical conductivity of shaley sands, which is equivalent to a classification of clayey sands in unconsolidated deposits. This model depends on both pore fluid conductivity and the conductivity of clay particles within the matrix, attributed to the high cation exchange capacity of clay (Waxman and Smits, 1968). Equation 4-6 shows the Waxman-Smits relation for bulk electrical conductivity a [S/cm], where crs is the solid matrix conductivity, represented by the equivalent cation conductance, B [S.cm /meq] and the concentration of exchange cations associated with the clay, Qv [meq/mL]. Representative values of Qv range from 0 to 0.5 for soils containing kaolinite and illite, and 0.3 to 1.5 meq/mL for soils containing montmorillonite. The equivalent cation conductance is described as a function of oy [S/cm], as shown in Equation 4-7. cr = (as + Gf)nm where, cr, = BQV (4-6) B = 0.046 [l-0.6e\"77ff/] (4-7) 4.1.3 Cross-Coupling Conductivity The streaming current cross-coupling conductivity coefficient describes the physical character of the electrical double layer at the solid-fluid interface in a soil. Although many laboratory and in-situ estimates of hydraulic and electrical conductivity exist for a wide range of geologic materials, a similar database of streaming current cross-coupling conductivity values has yet to be established. 53 A s w i t h a l l f l o w - r e l a t e d p h y s i c a l p r o p e r t i e s , the f u n d a m e n t a l c h a r a c t e r o f t h e c r o s s -c o u p l i n g c o e f f i c i e n t i s m o s t e f f e c t i v e l y e x a m i n e d t h r o u g h o n e - d i m e n s i o n a l f l o w s t u d i e s . T h e o n e - d i m e n s i o n a l f o r m o f the s t e a d y - s t a t e c o u p l e d e l e c t r i c a l f l o w r e l a t i o n i n E q u a t i o n 2 -2 8 r e d u c e s t o E q u a t i o n 4 - 8 , s u c h that t h e c r o s s - c o u p l i n g t e r m L, c a n b e e x p r e s s e d as a f u n c t i o n o f b o t h a s t r e a m i n g p o t e n t i a l c o u p l i n g c o e f f i c i e n t , C a n d t h e b u l k e l e c t r i c a l c o n d u c t i v i t y o f t h e m e d i u m , as s h o w n i n E q u a t i o n 4 - 9 . I f t h e h y d r a u l i c g r a d i e n t i n t h e c o u p l e d f l o w e q u a t i o n s i s e x p r e s s e d i n t e r m s o f h y d r a u l i c h e a d , L a n d C h a v e u n i t s [ A / m 2 ] a n d [ V / m ] , r e s p e c t i v e l y . I f h y d r a u l i c g r a d i e n t i s s t a t e d i n t e r m s o f p r e s s u r e , as i s o f t e n f o u n d i n t h e l i t e r a t u r e , L h a s u n i t s o f [ A . m / N , o r V . S / P a . m ] i f C i s s t a t e d i n t e r m s o f [ V / P a ] . T h e f o r m u l a t i o n o f t h e c o u p l e d f l o w e q u a t i o n s i n t e r m s o f p r e s s u r e i s o u t l i n e d i n A p p e n d i x B . A V L A h a L = -Ca ( 4 - 9 ) T h e s t r e a m i n g p o t e n t i a l c o u p l i n g c o e f f i c i e n t C i s d e t e r m i n e d e x p e r i m e n t a l l y , a l t h o u g h i n - s i t u m e t h o d s h a v e b e e n s u g g e s t e d ( F r i b o r g , 1 9 9 6 ) . T h e e x p e r i m e n t a l p r o c e d u r e i n v o l v e s m e a s u r i n g the v o l t a g e d i f f e r e n c e a l o n g t h e l e n g t h o f a s a t u r a t e d s o i l e l e m e n t s u b j e c t e d t o a g i v e n h y d r a u l i c g r a d i e n t . U n d e r l a m i n a r f l o w c o n d i t i o n s , t h e l i n e a r i n c r e a s e i n v o l t a g e d i f f e r e n c e w i t h i n c r e a s i n g p r e s s u r e h e a d d i f f e r e n c e d e f i n e s C , w h i c h i s o f t e n r e p o r t e d i n u n i t s o f [ m V / c m ] o r [ m V / k P a ] . D e c e p t i v e l y s i m p l e , s t r e a m i n g p o t e n t i a l m e a s u r e m e n t s r e q u i r e c a r e f u l l y d e s i g n e d e q u i p m e n t a n d p r o c e d u r e s t o e n s u r e d a t a q u a l i t y ( M o r g a n , 1 9 8 9 ) . T h e r e l a t i o n b e t w e e n C a n d f l u i d c o n d u c t i v i t y c a n b e s h o w n t h r o u g h a n e x a m i n a t i o n o f m e a s u r e d v a l u e s o f t h e s t r e a m i n g p o t e n t i a l c o u p l i n g c o e f f i c i e n t f o r u n i f o r m s a n d . A c o l l e c t i o n o f p u b l i s h e d v a l u e s o f C m e a s u r e d o n s a m p l e s o f p u r e q u a r t z s a n d s a t u r a t e d w i t h f l u i d s o f v a r y i n g c o n d u c t i v i t y i s s h o w n i n F i g u r e 4-1 ( A h m a d , 1 9 6 4 ; O g i l v y et a l . , 1 9 6 9 ; I s h i d o a n d M i z u t a n i , 1 9 8 1 ; F r i b o r g , 1 9 9 6 ) . A l t h o u g h the d a t a are l i m i t e d , a d e c r e a s e i n t h e a b s o l u t e v a l u e o f C w i t h i n c r e a s i n g f l u i d c o n d u c t i v i t y i s e v i d e n t . A s b u l k e l e c t r i c a l c o n d u c t i v i t y i s c o n t r o l l e d b y p o r e f l u i d c o n d u c t i v i t y i n c o h e s i o n l e s s s o i l s , t h i s t r e n d 5 4 I reinforces the importance of monitoring and reporting values of bulk conductivity during measurements of C for a given soil sample, such that L may be calculated. -0 .25 -0.2 ~ -0 .15 E > O -0.1 -0.05 0 X \\ \\ \\ \\ \\ \\ \\ s \\ X 1 < X \\ < 1 10 100 1000 a f ( uS /cm) Figure 4-1: Effect of fluid conductivity on measured values of the streaming potential coefficient for quartz sand. The concentration of different ions in solution control pore fluid conductivity, and have been shown to influence C. Morgan et al (1989) reported the ratio C/cy for Westerly granite to exhibit average values of 4 x 10\"4 A/m2 for 1:1 electrolytes, and 2 x 10\"4 A/m2 for 2:1 electrolytes. These values are based on experiments conducted on crushed samples saturated with solutions of NaCl and CaCb, respectively. The results suggest that L may be influenced by the distribution of ions in solution. Since natural waters contain primarily mono- and divalent ions in varying concentrations, it is important to perform measurements of the streaming current cross-coupling coefficient using in-situ fluids, or water displaying similar chemical, and consequently electrical, properties. No unique relation between Clcrj and grain size was evident in the data reported in Figure 4-1. Reports on the variation of C or L with grain size vary, such that no firm conclusions can be drawn (Mitchell, 1991; Ahmad, 1964; Ogilvy et al., 1969). However, C can be considered independent of porosity in cohesionless soils if surface conduction effects are neglected (Wurmstich and Morgan, 1994). In fine-grained soils, surface conductivity 55 effects increase the effective conductivity of the pore fluid, such that C appears to decrease with a smaller pore radius. Figure 4-2 illustrates the theoretical effect of changing pore radius on the streaming potential coefficient and effective fluid conductivity of crushed quartz saturated with a 250 pS/cm solution (pH = 7, 10\"3 N KN03, T = 45°C). In this analysis, the streaming potential coefficient is calculated using a form of the H-S equation described later in this section. Effective fluid conductivity cxey/ is calculated in terms of fluid conductivity and surface conductance, such that creff = oy +o~Jr, where r represents the hydraulic radius of the pore channels, and the surface conductance as is assumed equal to lxlO\"8 S (Ishido and Mizutani, 1981). From Figure 4-2, the effect of pore size appears to influence both C and aejf below a hydraulic radius of 0.01 mm, which translates to a representative grain size of 0.08 to 0.1 mm, using the approximation given in Equation 2-3. This range in grain size corresponds to a soil classification of fine sand. Figure 4-2: Variation of the effective pore fluid conductivity and streaming potential coefficient as a function of hydraulic radius in a quartz-water system (adapted from Ishido and Mizutani, 1981). Published data from one-dimensional streaming potential and electro-osmosis experiments are summarized in Appendix C. As many researchers report only the conductivity or composition of the saturating fluid, a value for L could not be derived directly from the measured data, but instead was calculated using an estimated value of bulk 56 c o n d u c t i v i t y . F o r a l l c o h e s i o n l e s s s a m p l e s , t h e b u l k c o n d u c t i v i t y w a s e s t i m a t e d u s i n g A r c h i e ' s l a w , a s s u m i n g a n a v e r a g e p o r o s i t y o f 0 . 3 5 u n l e s s o t h e r w i s e n o t e d . I n s a n d a n d s i l t s a m p l e s s a t u r a t e d w i t h l o w c o n d u c t i v i t y f l u i d s , s u r f a c e c o n d u c t i v i t y e f f e c t s i n f l u e n c e t h e b u l k c o n d u c t i v i t y s u c h that A r c h i e ' s l a w a l o n e p r o v i d e s t o o l o w a n e s t i m a t e . A s a r e s u l t , a r a n g e o f L v a l u e s are r e p o r t e d f o r t h e s e s o i l s , c o r r e s p o n d i n g t o a b u l k c o n d u c t i v i t y v a l u e b o u n d e d b y t h e p r e d i c t e d v a l u e u s i n g A r c h i e ' s l a w a n d t h e f l u i d c o n d u c t i v i t y . T h r o u g h t h e O n s a g e r r e c i p r o c a l r e l a t i o n s , v a l u e s f o r L w e r e c a l c u l a t e d f r o m t h e e l e c t r o - o s m o t i c p e r m e a b i l i t y c o e f f i c i e n t , ke u s i n g t h e r e l a t i o n s h o w n i n E q u a t i o n 4 - 1 0 , w h e r e yw i s t h e u n i t w e i g h t o f w a t e r [ N / m ] . L = y v K ( 4 - 1 0 ) F r o m t h e c o m p i l a t i o n o f e x p e r i m e n t a l d a t a c o n t a i n e d i n A p p e n d i x C , t h e s t r e a m i n g c u r r e n t c r o s s - c o u p l i n g c o e f f i c i e n t L a p p e a r s t o v a r y b e t w e e n 8 x l 0 \" 6 a n d 8 x l 0 \" 5 A / m 2 ( 8 x l 0 \" 4 a n d 8 x l 0 \" 3 m V . S / k P a . m ) f o r m o s t s o i l s . T h e m e a n o f t h i s r a n g e i s r e m a r k a b l y c l o s e t o t h e a v e r a g e v a l u e o f 5 x l 0 ' 5 A / m 2 r e p o r t e d b y C a s a g r a n d e ( 1 9 4 9 ) b a s e d o n h i s e x p e r i e n c e i n e l e c t r o - o s m o s i s t e s t i n g . I n u n i f o r m s a n d s , L c a n b e s e e n t o r a n g e u p t o a n a v e r a g e v a l u e o f 5 x l 0 \" 4 A / m 2 . T h i s u p p e r l i m i t m a y b e e x c e s s i v e l y h i g h d u e t o t h e e s t i m a t e d r a n g e o f b u l k c o n d u c t i v i t y v a l u e s u s e d i n t h e c a l c u l a t i o n o f L d e s c r i b e d a b o v e . S i n c e u n i f o r m s a n d s a r e n o t t y p i c a l l y u s e d i n e m b a n k m e n t d a m c o n s t r u c t i o n , a m o r e r e a s o n a b l e u p p e r l i m i t f o r t h e s e m a t e r i a l s m a y b e e s t i m a t e d as l x l O \" 4 A / m 2 , a c c o r d i n g to t h e v a l u e s r e p o r t e d i n A p p e n d i x C . C a l c u l a t i n g a t h e o r e t i c a l u p p e r b o u n d to L u s i n g E q u a t i o n 2 - 1 7 s u p p o r t s t h e v a l i d i t y o f a c r o s s - c o u p l i n g c o n d u c t i v i t y v a l u e a t t r i b u t e d t o a g i v e n s o i l u n i t , p a r t i c u l a r l y w i t h r e s p e c t t o m a t e r i a l s o f l o w h y d r a u l i c p e r m e a b i l i t y . F r o m t h e d e f i n i t i o n s o f Lh, Le a n d Lhe r e p o r t e d i n A p p e n d i x B , t h e m a x i m u m v a l u e o f L f o r a g i v e n s a m p l e m a y b e d e s c r i b e d i n t e r m s o f h y d r a u l i c c o n d u c t i v i t y a n d b u l k c o n d u c t i v i t y , as s h o w n i n E q u a t i o n 4 - 1 1 . A l l s o i l s i n A p p e n d i x C f o r w h i c h K i s l i s t e d s a t i s f y t h i s c o n d i t i o n . £max = JKyw a ( 4 - 1 1 ) 5 7 A theoretical description of the streaming potential coupling coefficient C [V/Pa] is given by the Helmholtz-Smoluchowski (H-S) equation, which was derived for flow in a porous plug from a capillary flow model (Overbeek, 1952). The H-S equation defines the ratio of electrical potential and pressure gradients as shown in Equation 4-12, where 8 is the electric permittivity of the fluid [F/m], ^ is zeta potential [V], ju is fluid viscosity [Pa.s], and oy is the electrical conductivity of the fluid [S/m]. The zeta potential is the electrical potential at the shear plane within the electrical double layer, which approximates the true surface potential of the solid material. Consequently, the H-S equation is of particular use in electrochemical studies of the solid-fluid interface, which aim to characterize surface potential. ^ AF s£ C = 7 ^ = ( 4- 1 2) AP jUCT, The H-S equation is valid only under laminar flow conditions and for porous media in which the curvature of the pore radius is much larger than the thickness of the double layer, such that surface conductivity effects can be neglected and bulk conductivity is a function only of fluid conductivity. In materials such as fine-grained or clayey soils where surface conductivity effects cannot be neglected, the H-S equation can be modified to account for surface conductivity as shown in Equation 4-13, where as is the surface conductivity and r is the average pore radius. A comprehensive derivation of Equations 4-12 and 4-13 is contained in the listed reference by Overbeek (1952). r AV eC C = 7IT = - 2cO ( 4 \" 1 3 ) AP af + • 4.2 Sensitivity Analyses To effectively model the SP response to seepage, accurate characterization of the hydraulic and electrical input parameters is required. The effect of internal erosion on these 58 physical properties must also be quantified to enable an assessment of the sensitivity of the method to anomalous seepage patterns. Several assumptions central to the numerical analysis contribute a margin of error to the modelling results, and are summarized below in point form. • Properties of a given soil unit are assumed homogeneous and, although the effect of anisotropy in K may be included in the hydraulic flow analysis using M O D F L O W , isotropic L and 2500 g/m2 of fines from the sample, which is equivalent to a minimum 1% change in mass due to the loss of fines. Characterization of disturbed and undisturbed sections of the core during the WAC Bennett Dam sinkhole incident enabled an assessment of the effects of internal erosion on soil properties in-situ. The core is described as a broadly graded, non-plastic silty sand with some gravel (Stewart and Watts, 2000). Investigations at Sinkhole 1 showed approximately a 2% loss in silt content and an increase in void ratio of between 0.09 and 0.12, which corresponds to an increase in porosity of 0.05 to 0.07. Estimates of hydraulic conductivity from in-situ techniques showed little to no variation between disturbed and undisturbed zones (Watts et al., 2000). 0 3 6 9 12 15 FINES, PERCENT a) concrete sand 1 61 3 6 9 12 15 FINES, PERCENT b) sand-gravel mixture1 0 2 4 6 8 10 FINES, PERCENT c) uniform fine sand2 Figure 4-3: Influence of type and amount of fines content on hydraulic conductivity (from 'Barber and Sawyer, 1952; 2Fenn, 1966). 62 4.2.2 Effect of Internal Erosion in 1 -D In a homogeneous soil element subjected to 1-D steady state flow conditions, the coupled electrical flow equation simplifies to Equation 4-8, which is the same relation used to define the streaming potential coefficient in lab experiments. As a result, numerical modelling in 1-D is of limited use, as it will simply predict an SP distribution that was originally used to define the streaming potential coefficient C. However, it is useful to x model a single soil element to illustrate the effect of a change in porosity on L and cr. A heterogeneous soil element is modelled to examine the effect of a change in hydraulic conductivity, and consequently a change in the hydraulic head distribution, on the SP response. Figure 4-4 displays the dimensions of the modelled homogeneous and heterogeneous saturated sand elements, each subjected to steady-state flow under a hydraulic gradient, Vh = 0.2. From Section 4.2, we will assume the loss of 4% fines fraction due to internal erosion results in an order of magnitude increase in K and a 0.07 increase in porosity between soil units 1 and 2. ffpipfi ^ ^ ^ ^ ^ ^ Q <=^> m i l p i l l \" V - S O I L 1/2 Q c = > S O I L 1 S O I L 2 0 5 m 1 m Figure 4-4: Homogeneous and heterogeneous soil elements. The SP response to changes in K, L or a caused by internal erosion may be evaluated by altering each parameter independently. Assumed base values of porosity, hydraulic conductivity, bulk electrical conductivity, and cross-coupling conductivity are represented by soil 1 and listed in Table 4-4. Bulk electrical conductivity is calculated as the sum of the value predicted using Archie's law, assuming a fluid conductivity of 100 pS/cm, and an 63 estimate of surface conductivity. Values of surface conductivity for soils 1 and 2 were derived from laboratory measurements performed on glacio-fluvial sands, which exhibited a decreasing trend from 10 pS/cm to 1.6 pS/cm corresponding to hydraulic conductivities of 0.01 cm/s and 0.1 cm/s, respectively (Huntley, 1986). The streaming potential coefficient was derived using Equation 4-9 and is listed for each soil type in Table 4-4. Table 4-4: 1-D hydraulic and electrical properties. K a Uo L -C Soil n (cm/s) (pS/cm) (Ohm.m) (A/m 2) (V/m) 1 0.3 0.01 31 322 8.0 x IO\"5 0.0258 2a 0.37 0.1 37 270 8.0 x IO\"5 0.0216 2b 0.37 0.1 29 345 8.0 x IO\"5 0.0276 2c 0.37 0.1 31 322 8.9 x 10\"5 0.0287 2d> 0.37 0.1 31 322 8.0 x 10\"5 0.0258 Soils 2a and 2b represent a change in bulk conductivity with porosity assuming L remains constant. An increase of 6 pS/cm in bulk conductivity is predicted using Archie's law assuming no change in surface conductivity, as represented by soil 2a. Soil 2b represents a net decrease in bulk conductivity of 4 pS/cm, which accounts for a change in as due to the loss of fines. Consequently, the streaming potential coefficient must vary in order to maintain a constant value of L, as illustrated by the change in magnitude of the SP profiles shown in Figure 4-5, where all SP values are reported with respect to the upstream end of the sample. A 19% increase or 7% decrease in bulk conductivity due to internal erosion, as described by soils 2a and 2b, correspond to a -16%) or +7% change in the amplitude of the SP response, and similarly the magnitude of the streaming potential coefficient. If a change in porosity is assumed to affect the contribution of pore fluid conductivity and surface conductivity equally such that a change in bulk conductivity is not expected, L is solely influenced by the streaming potential coefficient. If the absolute value of C is assumed to increase with increasing pore radius as per Figure 4-2, any increase in the magnitude of C would result in a greater SP response for a given bulk conductivity. This translates to a SP profile plotting above the soil 1 response, such as that shown for soil 2c in 64 F i g u r e 4 - 5 , w h i c h r e p r e s e n t s a n a r b i t r a r y 1 0 % c h a n g e i n C . I f s u r f a c e c o n d u c t i v i t y e f f e c t s are a s s u m e d n e g l i g i b l e s u c h that C d o e s n o t v a r y c o n s i d e r a b l y w i t h p o r o s i t y , i n t a c t a n d e r o d e d s o i l s w o u l d e x h i b i t a s i m i l a r S P r e s p o n s e w h e n s u b j e c t e d to a g i v e n h y d r a u l i c g r a d i e n t , as e x h i b i t e d b y S o i l s 1 a n d 2 d i n F i g u r e 4 - 5 . T h e e f f e c t o f a s s u m i n g C i n d e p e n d e n t o f p o r o s i t y r e s u l t s i n a S P p r o f i l e e q u i v a l e n t to t h e b a s e c a s e e x h i b i t e d b y s o i l 1, r e g a r d l e s s o f t h e b u l k c o n d u c t i v i t y v a l u e . T h e c h a n g e i n s t r e a m i n g c u r r e n t s o u r c e s t r e n g t h i s b a l a n c e d b y t h e c h a n g e i n b u l k c o n d u c t i v i t y , s u c h that t h e S P r e s p o n s e a l o n g t h e l e n g t h o f t h e s a m p l e r e m a i n s t h e s a m e . Figure 4-5: Effect of internal erosion on the 1-D SP distribution over a homogeneous soil element. A l t h o u g h s i n g l e e l e m e n t tests c a n p r o v i d e i n f o r m a t i o n a b o u t t h e h y d r a u l i c a n d e l e c t r i c a l p r o p e r t i e s o f a s o i l , t h e i n f l u e n c e o f i n t e r n a l e r o s i o n o n h y d r a u l i c f l o w p a t t e r n s are n o t a c c o u n t e d f o r w h e n t h e h e a d d i s t r i b u t i o n t h r o u g h t h e s a m p l e i s l i n e a r . T h e c a p a c i t y o f t h e S P m e t h o d t o d e t e c t t h e o n s e t o f i n t e r n a l e r o s i o n m a y b e d e m o n s t r a t e d i n a h e t e r o g e n e o u s e l e m e n t , w h e r e the e f f e c t o f i n t e r n a l e r o s i o n o n h y d r a u l i c c o n d u c t i v i t y , a n d c o n s e q u e n t l y t h e h y d r a u l i c r e g i m e i s i n c l u d e d i n t h e m o d e l . F i g u r e 4 - 6 s h o w s t h e S P r e s p o n s e t o i n t e r n a l e r o s i o n i n a n e l e m e n t c o m p o s e d o f s o i l 1 a n d s o i l 2 , w h e r e s o i l 2 d r e p r e s e n t s t h e e r o d e d z o n e w i t h a n i n c r e a s e i n K, a s s u m i n g L a n d 6 5 cr remain constant with the change in porosity. The non-uniform hydraulic head distribution resulting from the discontinuity is mirrored by the SP response, which exhibits a maximum anomalous increase of 82% from the homogeneous soil 1 SP profile. The heterogeneity corresponds to an equivalent increase in seepage flow rate of 82%> from the base case. 0 . 0 5 0 . 1 5 0 . 2 5 0 . 3 5 0 . 4 5 0 . 5 5 0 . 6 5 0 . 7 5 0 . 8 5 0 . 9 5 S a m p l e l e n g t h ( m ) Figure 4-6: Effect of internal erosion on the 1-D SP distribution over a heterogeneous soil element. The analysis of a saturated soil element illustrates that the SP response in 1-D is governed by C and K. In a homogenous element subjected to a given hydraulic gradient, an increase or decrease in C dictates the electrical potential across the element, since C = AV/Ah. If C remains constant, the effect of a variation in bulk electrical conductivity due to changes in porosity will not affect the 1-D SP response since convection and conduction current flow are constrained by the same electrical flow boundaries. Any change in streaming current source strength is balanced by a change in bulk conductivity, such that the resultant SP response across the element is unchanged. However, a significant SP response results from a non-uniform hydraulic gradient caused by a change in K, as shown by the heterogeneous soil element portrayed in Figure 4-6. Even if C remains constant, the effect of a porosity change on K, and consequently the hydraulic gradient, appears to have a significant influence on the SP response. 66 4.2.3 Effect of Internal Erosion in 3-D A multi-dimensional analysis enables an approximation of the SP response to a zone of increased seepage at depth. Deviations in hydraulic flow, streaming current flow and conduction current flow occur in response to changes in physical properties. These variations may arise at a boundary between different geologic or soil units, or within a given zone due to the onset of internal erosion. Consequently, it is of interest to study the effect of internal erosion on the hydraulic and electrical flow regimes within a zoned embankment, and to characterize the SP response across the surface of a dam. A schematic of the field-scale 3-D numerical model embankment is shown as Figure 4-7 in transverse cross-section at the mid-point of the lOOm-long crest. The erosion channel exhibits a cross-sectional area of 200 m2, penetrating the core at a depth of 65 m. This zone was modelled with an increased hydraulic conductivity to represent the loss of approximately 7% fines from the core to the downstream shell. Properties of the core, eroded zone, and outer shell materials are listed in Table 4-5, corresponding to a fluid conductivity of 100 pS/cm chosen for the reservoir. The saturated bulk electrical conductivity of the core is representative of glacial till or well-graded silty sand, and is assumed to decrease with the loss of fines through the process of internal erosion. The streaming potential coefficient of the core material is assumed constant, such that the streaming current cross-coupling coefficient of the piped zone is smaller than that of the intact core. 30 < — — J > 20 350 (dimensions in m) Figure 4-7: Transverse cross-section of the field-scale dam model with erosion channel. 67 Table 4-5: 3-D isotropic hydraulic and electrical properties. K L Csat ^uns at Material (cm/s) (A/m2) (S/m) (S/m) Shell 1 x 10\"3 3.0 x 10\"5 3 x 10\"J 2.0 x 10\"4 Core 1 x 10\"6 2.0 x 10\"5 8 x 10-3 2.5 x 10\"3 Piped zone 1 x lO\"4 1.0 x 10\"5 4 x 10\"3 2.0 x IO- 4 The effect of the presence of an erosion channel on the hydraulic regime is shown in Figure 4-8, where fluid flow is channelled through the more conductive piped zone. This channelling effect in both the hydraulic and electrical flow regimes is further illustrated in Figure 4-9, which depicts a longitudinal cross-section parallel to and directly beneath the crest. Fluid flow is concentrated within the more conductive eroded zone, while streaming current flow predominates in the intact core beneath the phreatic surface, due to the higher cross-coupling conductivity of this material. b) core with piped zone Figure 4-8: Transverse cross-section of hydraulic head (m) and flux. 68 0 20 40 60 80 10£ a) hydraul ic head (m) and flux 100 | 1 1 1 1 1 1 1 90 -80 -70 -i / / 4 \\ \\> \\ \\ i / / / \\ \\ \\ i / / / s. \\ \\ i / / \\ \\ / / S \\ \\ / / / if \\ \\ / / / / 1 \\ \\ \\ \\ / / / / \\ \\ \\ \\ / / \\ t \\ \\ I / / / \\ 1 \\ \\ \\ 1 / / / X \\ \\ \\ 1 1 / \\ \\ \\ ! \\ 10 20 30 40 50 60 70 80 90 Longitudinal distance (m) b) streaming current flow F igure 4-9: Long i tud ina l cross-section of flow wi th in saturated zone beneath crest. 6 9 In the absence of single point sources of fluid flow in the subsurface, Equation 3-4 indicates that streaming or convection current source terms are concentrated where there exist both gradients in L , such as at physical property boundaries, and significant hydraulic gradients. The distribution of streaming current source amplitude, ICOnvection, is shown in Figure 4-10 for a cross-section through the centre of the core, both intact and with the erosion channel. Streaming current flows uniformly from the sources lining the upstream shell-core boundary to the negative sink at the downstream toe of the core in Figure 4-10a. The large hydraulic gradient between the erosion channel and downstream toe of the core results in an increase in streaming current flow between the lower boundary of the erosion channel and the downstream core shell boundary, as shown by the distribution of current sources in Figure 4-10b. As electrical current is most likely to follow the path of least resistance, conduction current flow is concentrated in the saturated core material linking conduction current sources, which are of equal magnitude but opposite polarity to the convection current sources shown in Figure 4-10. The conduction current flow pattern is shown for cross-sections of intact core and core containing the erosion channel in Figure 4-11. The SP response at the surface of the embankment is presented both in absolute and difference form. SP profiles corresponding to the intact dam and the dam containing the erosion channel are shown in Figure 4-12 along an upstream-downstream line crossing the midpoint of the crest. All potentials are referenced to a base point at the surface of the reservoir above the upstream toe. Evident in the figure is the net increase in SP with distance downstream, which reflects the distribution of conduction current source terms and consequently conduction current flow in the subsurface. A net decrease in the magnitude of the SP response results from the presence of the erosion channel. The difference between the two data sets over the entire surface of the dam is presented in Figure 4-13. These data represent the change in SP response due to the presence of the erosion channel. A maximum difference anomaly of -11 mV, or -21% of the response to the intact core, occurs on the downstream face of the dam, and corresponds to an 11% increase in seepage. 70 x10\"3 b) core with piped zone Figure 4-10: Distribution of streaming current flow sources (I) in longitudinal cross-section. 71 a) intact core Transverse distance (m) b) core with piped zone Figure 4-11: Conduction current flow in longitudinal cross-section. 72 > E, CL if) 160 140 120 100 80 60 40 20 0 -20 • intact core •core with piped zone 50 100 150 200 Transverse distance (m) 250 300 350 Figure 4-12: SP profiles over dam surface at a longitudinal distance of 50m. 50 100 150 200 Transverse distance (m) 250 300 Figure 4-13: SP difference contours (mV) over dam surface (crest delineated by dotted lines). The sensitivity of the SP response to the presence of the erosion channel may be compared to the sensitivity of the SP response to the physical properties of the eroded zone. The influence of the hydraulic and electrical parameters can be assessed as in the 1 -D case through an independent analysis of each parameter. Although an increase in hydraulic conductivity with the loss of fines is expected from the published laboratory data shown in Section 4.2.1, a negligible difference in hydraulic conductivity was reported for zones of weak and intact core during the Bennett dam sinkhole 73 i n v e s t i g a t i o n . M a i n t a i n i n g a c o n s t a n t v a l u e o f K f o r b o t h i n t a c t a n d e r o d e d z o n e s r e s u l t s i n a h y d r a u l i c d i s t r i b u t i o n i d e n t i c a l t o that s h o w n i n F i g u r e 4 - 8 a . T h e d e c r e a s e i n c r o s s - c o u p l i n g c o n d u c t i v i t y a n d s a t u r a t e d b u l k c o n d u c t i v i t y v a l u e s , as d e t a i l e d i n T a b l e 4 - 5 f o r t h e p i p e d z o n e , i n f l u e n c e t h e f l o w o f s t r e a m i n g a n d c o n d u c t i o n c u r r e n t b e t w e e n t h e u p s t r e a m a n d d o w n s t r e a m b o u n d a r i e s o f the c o r e . H o w e v e r , t h e s u r f a c e S P d i s t r i b u t i o n r e s u l t i n g f r o m t h e s e c h a n g e s i n e l e c t r i c a l p r o p e r t i e s i s e q u i v a l e n t t o t h e o r i g i n a l c a s e w i t h i n t a c t c o r e . S i m i l a r l y , e l e c t r i c a l f l o w a n a l y s e s o f t h e e r o d e d c o r e h y d r a u l i c m o d e l , p e r f o r m e d to a s s e s s t h e e f f e c t o f a r e a l i s t i c c h a n g e i n C, asat o r L o f t h e p i p e d z o n e , p r o d u c e < 1 % d e v i a t i o n f r o m t h e s u r f a c e S P d i s t r i b u t i o n o f t h e p i p e d z o n e c a s e s h o w n i n F i g u r e 4 - 1 2 . A l t h o u g h c h a n g e s i n t h e e l e c t r i c a l p r o p e r t i e s o f t h e e r o d e d z o n e a p p e a r t o h a v e l i t t l e e f f e c t o n t h e S P d i s t r i b u t i o n at s u r f a c e , e f f e c t i v e c h a r a c t e r i z a t i o n o f asat a n d L t h r o u g h o u t t h e e m b a n k m e n t i s i m p o r t a n t to t h e p r e d i c t i o n o f a n a c c u r a t e S P d i s t r i b u t i o n . A s w a s s h o w n b y t h e 1 -D s a t u r a t e d f l o w a n a l y s e s p r e s e n t e d i n S e c t i o n 4 . 2 . 2 , t h e s e p a r a m e t e r s c o n t r o l t h e m a g n i t u d e o f t h e S P r e s p o n s e . T h e i m p o r t a n c e o f c h a r a c t e r i z i n g t h e e l e c t r i c a l c o n d u c t i v i t y d i s t r i b u t i o n o f t h e s u b s u r f a c e m a y b e f u r t h e r i l l u s t r a t e d t h r o u g h a s t u d y o f t h e e l e c t r i c a l c o n d u c t i v i t y o f t h e u n s a t u r a t e d z o n e . T h e i n f l u e n c e o f t h e b u l k u n s a t u r a t e d e l e c t r i c a l c o n d u c t i v i t y o n c o n d u c t i o n c u r r e n t f l o w a n d c o n s e q u e n t l y t h e S P d i s t r i b u t i o n i s e a s i l y a s s e s s e d i n a m u l t i - d i m e n s i o n a l n u m e r i c a l a n a l y s i s , w h e r e t h e m a t e r i a l a b o v e a n d b e l o w t h e p h r e a t i c s u r f a c e m u s t b e c h a r a c t e r i z e d . T h e o r i g i n a l S P d i s t r i b u t i o n o v e r t h e s u r f a c e o f t h e d a m w i t h i n t a c t c o r e i s c o m p a r e d t o s i m i l a r m o d e l s i d e n t i f i e d as c a s e s A a n d B i n F i g u r e 4 -14, w i t h u n s a t u r a t e d b u l k c o n d u c t i v i t i e s as o u t l i n e d i n T a b l e 4 - 6 . T h e e f f e c t o f t h e u n s a t u r a t e d z o n e e x h i b i t i n g t h e s a m e e l e c t r i c a l c o n d u c t i v i t y as t h e s a t u r a t e d z o n e i s i l l u s t r a t e d b y c a s e A , w h i c h s h o w s a m a x i m u m 6 5 % d e c r e a s e i n S P a m p l i t u d e . A n o r d e r o f m a g n i t u d e d e c r e a s e i n the u n s a t u r a t e d c o n d u c t i v i t y v a l u e s f r o m t h o s e r e p o r t e d i n T a b l e 4 - 5 f o r m s c a s e B , w h i c h r e s u l t s i n a m a x i m u m a m p l i t u d e i n c r e a s e o f 2 4 % . Table 4-6: Unsaturated bulk conductivity of intact zoned embankment materials. Original intact core Case A Case B Material Olmsat (S/m) CTu„Sat (S/m) CTunsat (S/m) Shell 2.0 x 10\"4 3.0 x 10° 2.0 x 10\"5 Core 2.5 x 10\"3 8.0 x 10\"3 1.0 x 10'4 7 4 200 350 Transverse distance (m) Figure 4-14: Effect of unsaturated bulk conductivity on the surface SP distribution. T h e 3 - D a n a l y s i s o f a f i e l d - s c a l e e m b a n k m e n t i l l u s t r a t e s h o w t h e d i s t r i b u t i o n o f p h y s i c a l p r o p e r t i e s i n f l u e n c e s h y d r a u l i c a n d e l e c t r i c a l f l o w i n t h e s u b s u r f a c e , a n d d i c t a t e s t h e o v e r a l l m a g n i t u d e o f t h e S P r e s p o n s e to s e e p a g e . A c o m p a r i s o n o f t h e h y p o t h e t i c a l e m b a n k m e n t s w i t h i n t a c t c o r e a n d a c o r e e x h i b i t i n g a n e r o d e d z o n e d e m o n s t r a t e t h e i n f l u e n c e o f t h e h y d r a u l i c r e g i m e o n t h e c h a r a c t e r o f t h e S P d i s t r i b u t i o n . T h e a n t i c i p a t e d c h a n g e i n e l e c t r i c a l p a r a m e t e r s d u e t o t h e o n s e t o f i n t e r n a l e r o s i o n h a s l i t t l e e f f e c t o n t h e r e s u l t a n t S P r e s p o n s e c o m p a r e d to t h e e f f e c t o f a c h a n g e i n h y d r a u l i c c o n d u c t i v i t y , as s h o w n b y t h e s e n s i t i v i t y a n a l y s i s o f e r o d e d z o n e p a r a m e t e r s . T h e v a r i a t i o n i n h y d r a u l i c h e a d g r a d i e n t s t h a t r e s u l t f r o m a c h a n g e i n h y d r a u l i c c o n d u c t i v i t y a p p e a r t o h a v e t h e greates t i n f l u e n c e o n t h e S P r e s p o n s e , a n d d i c t a t e w h e t h e r a z o n e o f i n c r e a s e d s e e p a g e m a y b e d e t e c t e d as a n S P a n o m a l y at s u r f a c e . 4.3 Summary and Conclusions C h a p t e r 4 f u l f i l l e d t h e s e c o n d t h e s i s o b j e c t i v e w i t h a n e x a m i n a t i o n o f t h e t h r e e m a c r o s c o p i c p a r a m e t e r s r e q u i r e d i n a n u m e r i c a l a n a l y s i s o f h y d r a u l i c a n d e l e c t r i c a l f l o w : h y d r a u l i c c o n d u c t i v i t y , e l e c t r i c a l c o n d u c t i v i t y a n d s t r e a m i n g c u r r e n t c r o s s - c o u p l i n g c o n d u c t i v i t y . T h e c o m p l e x i t y o f t h e s e p a r a m e t e r s w a s i n t r o d u c e d w i t h e s t a b l i s h e d t h e o r e t i c a l m o d e l s , a n d a c o m p i l a t i o n o f p u b l i s h e d d a t a i l l u s t r a t e d t h e r a n g e o f v a l u e s 75 possible for a given material and the importance of proper characterization. A series of numerical analyses were performed in an effort to quantitatively assess the influence of internal erosion on subsurface properties and the resultant SP response. While the electrical parameters were shown to influence the total magnitude of the SP response, the distribution of hydraulic head in the subsurface was shown to dictate the character of the SP distribution. Key conclusions and contributions include the following: • Predictions of bulk electrical conductivity using Archie's law quantify differences in saturating fluid properties, porosity and saturation, but do not account for surface conductivity effects that may be significant in cohesionless soils in fresh water environments. Consequently, these approximations of the actual bulk conductivity distribution should not replace representative in-situ measurements. • The hydraulic conductivity, streaming current cross-coupling coefficient, and saturated and unsaturated bulk conductivity must be characterized for all soils at given site in order to accurately interpret the magnitude of the SP response due to seepage. • The direct relation between the distribution of electrical potential and the distribution of hydraulic potential, or hydraulic head, defined in the mathematical model derived from coupled flow theory suggests the pattern of the SP response is dictated by the hydraulic regime. • The effect of internal erosion on hydraulic conductivity appears to be significantly larger in magnitude than the estimated variation in electrical properties. Consequently, for a given set of parameters, the degree of variation in the hydraulic head distribution resulting from a change in hydraulic conductivity dictates the amplitude of the anomalous SP response. • The quantitative effect of internal erosion on hydraulic and electrical soil properties under in-situ conditions, in particular the electrical conductivity, streaming potential coefficient and cross-coupling conductivity, must be further investigated to confirm the anticipated variations and perhaps re-evaluate their influence on the SP response. 76 5 NUMERICAL ANALYSES Evaluation of the numerical procedure 3DSP through application to physical seepage problems and a comparison with measured data is the focus of Chapter 5. Section 5.1 presents a comparison of laboratory and numerical data sets that describe the SP response over a controlled physical model of an embankment, with a discussion on physical property characterization. Preliminary modelling of the SP response to uniform seepage through Mica dam in British Columbia is described in Section 5.2, and illustrates the 3-D character of hydraulic and electrical flow through an embankment. 5.1 Laboratory-Scale Embankment A concurrent laboratory research program designed to investigate the SP response to seepage through embankment materials enabled a comparison between measured and numerical results. All laboratory data were collected and provided courtesy of Powertech Labs Incorporated. 5.1.1 Description of Laboratory Models and Physical Properties Two laboratory-scale embankments were designed and constructed by Powertech Labs Inc. to study the influence of a pervious zone on the measured SP response to steady-state seepage flow. The embankments were constructed of pure quartz silica sand (Ottawa sand) graded according to ASTM C33-93 and compacted in layers at approximately 5% moisture content within a Plexiglas tank, with dimensions as shown in Figure 5-1. Filter mesh and gravel were placed above the surface of embankment at the downstream toe to prevent sloughing and migration of the sand down a Plexiglas collector trough designed to guide seepage water downstream of the embankment. This extended drainage system was an improvement on an original design, which consisted of a vertical outflow pipe located directly beneath the downstream toe. Seepage was induced by means of an upstream reservoir maintained at constant levels of 22.5 cm and 18.0 cm. Two embankment configurations were tested at these reservoir levels: a homogeneous dam and a dam containing an upstream defect. A pervious zone extending the full 15 cm crest length of the embankment was introduced using a 3.75 cm layer of concrete sand placed 4 cm above the base. The defect extended from the upstream face of the dam to a distance of 60 cm from 77 the upstream toe, as indicated in Figure 5-1. Grain size distributions of both the Ottawa sand and concrete sand mixtures are plotted in Figure 5-2. 6 cm Delineates pervious defect zone Figure 5-1: Laboratory embankment dimensions with manometer and electrode measurement locations. i i r i i I III / / / / — O t t a w a Sc C o n c r e t e and s a n d li / / / 0.01 0.1 1 1 0 Particle size (mm) Figure 5-2: G r a i n size distribution of embankment materials. The hydraulic head distribution within each embankment was monitored by means of six manometers connected to ports located at the side of the tank, as indicated in Figure 5-1. Electrical SP measurements were recorded automatically using a series of gel-filled Ag-AgCl electrodes placed in contact with the soil along the upstream and downstream face of the 78 dam. Hydraulic and electrical readings were monitored until the system stabilized, and were then recorded for a period of five days under steady-state flow conditions. One-dimensional testing was carried out to characterize the SP coupling coefficient of the Ottawa sand and concrete sand (Sasitharan et al, 2001). Although the electrical conductivity, temperature and pH of the fluid were monitored, the bulk electrical conductivity of each sample was not measured directly, and the unidirectional flow source was uncontrolled such that fluid properties varied during testing. Consequently, accurate determination of the streaming current cross-coupling coefficient is difficult given that the electrical conductivity of the 1 -D test samples must be approximated, and the fluid used in the 1 -D testing was not always chemically or electrically equivalent to that used in the model embankment experiments. Table 5-1 summarizes the physical properties of Ottawa sand and concrete sand measured and calculated from the results of the 1-D testing. Estimates of hydraulic conductivity were derived from measured volumetric flow rates, and porosity was calculated based on bulk density measurements reported in Sasitharan et al (2001). The SP coupling coefficient and fluid conductivity were measured for each test. To account for the large variation in fluid conductivity, Archie's law was used to approximate the bulk conductivity of both 1-D test samples. Since the tap water used in testing exhibited a low electrical conductivity, the bulk conductivity values reported in Table 5-1 include estimates of surface conductivity, obtained from tests performed on soils of similar grain-size distribution reported by Huntley (1986). The range of values reported for the streaming current cross-coupling coefficient was calculated from the measured SP coupling coefficient and calculated bulk conductivity values using Equation 4-9. Table 5-1: Physical properties derived from 1-D experimental data. K n -C C\\v crsat {estimated} L Material (m/s) (V/m) (uS/cm) (pS/cm) (A/m 2) Ottawa sand 4.5 x 10\"4 0.40 0.172 14-16 14-15 2 .4-2.6 x IO\"4 Concrete sand 5.5 x IO\"3 0.45 0.074 2 0 - 2 8 12-15 8.9 - 11 x IO\"5 79 Bulk conductivity measurements of saturated and partially saturated Ottawa sand, and saturated concrete sand were performed on unique soil samples using a Miller soil box and Nilsson resistivity meter (ASTM G57-78). The samples were saturated with fluid exhibiting a conductivity of 28 pS/cm, which is in the range reported for the 1-D testing of concrete sand, but higher than that used in conjunction with the 1-D tests of Ottawa sand. Table 5-2 lists the measured bulk conductivity values and shows the corresponding cross-coupling coefficient for each material, calculated using the measured SP coefficient reported in Table 5-1 and the average measured saturated conductivity, assuming the measured data are representative of the electrical character of the 1-D test samples. Table 5-2: Measured bulk electrical conductivity and corresponding estimates of cross-coupling conductivity. o- s a t{S= 100%} a;,„sa. {S= 10%} Cunsat {S - 5%} L Material (US/cm) (LiS/cm) (LiS/cm) (A/m 2) Ottawa sand 29 (27.7-30.8) 15.4 10.4 5.0 x 10\"4 Concrete sand 10 - - 7.4 x 10 - 5 A comparison of the cross-coupling coefficients presented in Tables 5-1 and 5-2 with the laboratory data compiled in Appendix C shows that the reported values reach, and tend to exceed, the upper range of measured or derived data for most soils. In particular, the value of 5.0 x 10\"4 A/m2 listed in Table 5-2 for Ottawa sand appears high relative to the range reported for sand. The magnitude of this coefficient is due in part to the measured average saturated conductivity of 29 mS/cm, which appears incongruous considering the reported saturating fluid conductivity value of 28 itS/cm. However, a theoretical analysis of bulk conductivity may not accurately represent the surface conductivity effects associated with the low fluid conductivity coupled with the fine-grained nature of the sand. Since the actual bulk conductivity distribution of the embankments is unknown, the data presented in Tables 5-1 and 5-2 provide upper and lower bounds to the electrical input parameters for the purposes of numerical modelling. Two sets of input parameters were applied in the numerical analysis of the embankments and are listed in Table 5-3. The measured hydraulic conductivity and 80 streaming current cross-coupling coefficient derived from the 1-D testing of each material were used directly in the numerical models. The Ottawa sand used to construct the embankments was compacted to an average density of 1700 kg/m , which corresponds to an average porosity of 0.36. This approximation was used along with the measured reservoir water conductivity of 28 pS/cm and visual observations of saturation conditions to derive a theoretical estimate of the electrical conductivity distribution within each embankment. These calculated estimates form case A and are listed in Table 5-3 a. For comparative purposes, measured bulk conductivity values were also used to characterize the electrical conductivity distribution, and the corresponding parameters in Table 5-3b define the case B analysis. Table 5-3: Numerical modelling input parameters: laboratory embankment. a) Case A - estimated electrical conductivity K L Csat \"iinsat Material (cm/s) (A/m 2) (uS/cm) (uS/cm) Ottawa sand 4.5 x IO\"4 2.5 x IO - 4 17 15 Concrete sand 5.5 x 10'3 1.0 x IO\"4 14 -b) Case B - measured electrical conductivity K L Csat \"unsat Material (cm/s) (A/m 2) (uS/cm) (uS/cm) Ottawa sand 4.5 x lO\"4 5.0 x IO\"4 29 15 Concrete sand 5.5 x IO\"3 7.4 x IO - 4 10 -5.1.2 Experimental and Numerical Results A 141><5x31 cell mesh was used to model the laboratory dams under steady-state flow conditions using 3DSP. Measured and modelled hydraulic head data were compared at each port location, as shown in Table 5-4 for the two embankment configurations at both reservoir levels. Figure 5-3 displays the results of the 3-D numerical seepage analyses taken through a central cross-section of each embankment. As the embankments and pervious defect were designed to test the SP response to seepage in 2-D, these cross-sections are representative of the numerical solution in all cells parallel to the crest. 81 a) homogeneous dam; 22.5 cm reservoir level b) homogeneous dam; 18.0 cm reservoir level c) dam with upstream defect; 22.5 cm reservoir level d) dam with upstream defect; 18.0 cm reservoir level Figure 5-3: Predicted hydraulic head distribution. 82 Table 5-4: Hydraulic head values reported at manometer port locations (cm of water) and volumetric seepage flow through embankment. 22.5 cm reservoir level 18.0 cm reservoir level Port Homogeneous Upstream defect Homogeneous Upstream defect measured modelled measured modelled measured modelled measured modelled 1 22.5 22.1 22.2 22.2 - - - -2 18.5 19.0 20.2 19.7 - - - -J 4 22.8 22.5 22.1 22.4 18.2 18.0 17.5 18.0 5 18.6 18.7 20.2 19.4 15.1 14.4 15.0 15.6 6 10.1 8.8 9.8 9.2 8.6 7.7 7.8 7.9 Q (cm3/s) 0.917 1.39 0.667 1.74 0.508 0.926 0.435 0.995 Table 5-4 reveals strong agreement between the measured and modelled hydraulic head distributions, but poor correlation between measured and modelled volumetric seepage flow rates through the embankment. Although the modelled seepage quantities do not correspond to the measured values, it is encouraging to note that the decrease in seepage through the homogeneous dam corresponding to a decrease in reservoir level is equal to 45% in both measured and modelled data. Discrepancies between measured and modelled hydraulic data may be caused in part by the characterization of each soil unit with a homogeneous and isotropic hydraulic conductivity value. The modelled SP distribution through the body of each embankment is presented in profile form to show the SP response at the surface of the dam. Figures 5-4 and 5-5 display profiles of both the measured and modelled SP data along the upstream and downstream face of each embankment, with all electrical potential values reported with respect to a base electrode position at the centre of the crest. In each plot, the numerical solution resulting from the Case A input parameters listed in Table 5-3 is represented by a solid line, and the solution based on the Case B parameters is represented by a dotted line. Figure 5-4 displays the surface SP profiles over the homogeneous dam at two reservoir levels. Evident in the plots is the net increase in SP with distance downstream, with an absolute response of almost 30 mV at the 22.5 cm reservoir level, and 20 mV at the 18.0 cm reservoir height. This decrease in magnitude of the SP response with reservoir elevation 83 is revealed in both measured and modelled data. Excellent agreement exists between the measured data and the Case A modelling results over the entire surface of the homogeneous embankment at both reservoir levels, with the predicted response falling within the 2 mV absolute error range of most measured data points. The discrepancy between the measured data and the Case B modelling results emphasizes the importance of proper parameter characterization. Although the shape of the SP response remains constant, an increase in amplitude results from the application of the Case B physical properties in the numerical model. Excellent agreement between measured and Case A modelled data along the downstream face of the dam with upstream defect is illustrated in Figure 5-5a. A positive trend in the measured data down the upstream face of the dam corresponds to up to a 5 mV departure from modelled results at this higher pool level. Although a positive trend along the upstream face is also suggested in Figure 5-5b, the measured data corresponding to the 18.0 cm reservoir level show significantly more scatter than the other data sets. The positive deviation of the measured data with respect to the predicted results along the upstream face in Figure 5-5a, and to a lesser extent in Figure 5-5b, suggest that additional flow behaviour may exist in the upstream portion of the laboratory dam, which is not captured in the numerical model. The discrepancy may be caused by external sources of error influencing the measured results, such as the proximity of the reservoir re-supply conduit, or inhomogeneities in the soil surrounding the defect, which may cause further variations in the hydraulic distribution and consequently change the shape of the SP distribution upstream. The numerical model predicts a disturbance in the hydraulic regime due to the upstream pervious defect, which is substantiated by the measured hydraulic head data. However, due to the more moderate hydraulic gradients observed in the upstream portion of the dam, the impact of the defect on the hydraulic, and consequently the electrical potential distributions is slight. The perturbation in electrical flow translates to a change on the order of 1 mV within the dam, as shown in Figure 5-6 for the high reservoir level, but does not correspond to significant variation along the downstream face. Consequently, the SP signatures on the surface of the homogeneous and upstream defect dams at a given reservoir level are very similar. 84 a) 22.5 cm reservoir level b) 18.0 cm reservoir level Figure 5-4: Measured and modelled surface SP distribution over the homogeneous dam. 85 30 -10 1 — , — , — , — , — , — , — : — , — > — > — i — i — i — i — i — i — i — • — . — i — i — i — • — i — ' — • — ' — l 0 0.2 0.4 0.6 0.8 1 1.2 1.4 Distance downstream (m) a) 22.5 cm reservoir level b) 18.0 cm reservoir level Figure 5-5: Measured and modelled surface SP distribution over the dam with upstream defect. 86 b) dam with upstream defect Figure 5-6: Modelled SP distribution (mV) through a cross-section of the embankment at the 22.5 cm reservoir level. 87 5.2 Mica Dam, British Columbia A 3-D seepage and electrical flow model of Mica dam illustrates the application of the numerical procedure to an actual dam, and provides a first approximation of the change in the SP signature at the surface of the embankment between representative high and low reservoir elevations. Construction details, hydraulic readings and self-potential survey data were acquired from internal reports and reports delivered by external consultants to BC Hydro. 5.2.1 Project Description Mica dam is a zoned earthfill dam located on the Columbia River in British Columbia, Canada. The embankment, which has been in operation since 1973, is comprised of a near-vertical glacial till core flanked by inner sand and gravel shells and an outer rockfill shell. Mica dam has a maximum height of 244 metres and a crest length of approximately 813 metres, and is situated in a narrow section of the river valley, with abutments exhibiting an average slope of 40°. The sound bedrock underlying the core is overlain by fluvial overburden of varying thickness in other sections of the dam. Figure 5-7 shows a plan view of the dam site, which includes an underground power generation station. 88 89 5.2.2 Hydraulic Flow Analysis A 95x41x45 non-uniform grid was generated to represent Mica dam in a numerical analysis of the SP response to seepage induced by the upstream reservoir. Two steady-state seepage analyses were performed at representative high and low pool elevations to set upper and lower bounds on the hydraulic head distribution for the purposes of the electrical flow modelling. As shown in Figure 5-8, the Mica reservoir level fluctuates and achieves a peak high pool and minimum low pool elevation annually. Although these fluctuations can give rise to transient flow conditions within the low-permeability core, transient effects are ignored and the steady-state behaviour is deemed representative of the maximum possible hydraulic response at a given pool level, for the purposes of preliminary modelling. A high pool level of 754.4 metres was chosen to represent the maximum normal level of operation, and a low pool level of 731.4 metres was chosen to correspond with the reservoir elevation reported during an on-site field SP survey in May 1998. 760 -, \\\\\\\\'\\ I I I I _ r T _ _ . r . 1111 1 1 ] 1 1 j [ ~r~t rr\" 1 1111 1111 /DO • 750 7 A P. I ! I I I 1 1 1 1 I M 1 1 1 1 x i ; ; 1 1 1 1 1 ; : : i : : : i I I I I I 1111 \\ 1 1 1 1 1 ! ii i i \\ i [' ' 1 \\ j\\ i i iV 111 1 f i \\ : i i i i i ; : : : IHO -; ; ; ; ; / 1111 1 1 1 1 \\ 1 1 1 1 N 1 I ] i i i i I I I I ion (n n c V i l l i 11111 v 1 ' 1 f \\ i i i I I I ! \\ 1 1 i: i: \\ 1 1 1 1 levat 0 c l \\ I 11 11111 1 \\ l l/l 11111 \\ 1 1 1 1 1 I\\I i j 11111 \\! 1 /; i :'; LU ' o u I M i I I I i i iV i 1 1 1 1 1 III ; I 1 I I ! > 11111 / yon I I I I I 11111 —1—t—t—1—1— 1 I 1 1 I 111: i ; 1 1 1 1 11111 1 £AJ -7 1 £ I I I I I 11111 1 1 1 1 1 11111 1 1 1 1 1 1 1 1 1111 1 1 1 | 11111 1 I A'J / 1 0 -710 - ! ! ! 1 1 1 11111 1 1 1 1 1 11 11 01-Jan- 02-Jul- 31-Dec- 01-Jul- 30-Dec- 30-Jun- 29-Dec- 28-Jun- 27-Dec- 27-Jun- 26-Dec-97 97 97 98 98 99 99 00 00 01 01 Figure 5-8: Mica dam reservoir elevation. The geometry of the embankment and its influence on the flow regime is illustrated in Figure 5-9, which shows cross-sectional plots of the hydraulic head distribution through the 90 dam at high pool. A transverse cross-section through the centre of the embankment is displayed in Figure 5-9a, which intersects the crest approximately at station 25+00 ft in Figure 5-7. Longitudinal cross-sections extend from the western abutment at 350 metres (station 11+66 ft) to the eastern edge of the dam at 1170 metres (station 38+34 ft). These sections are taken beneath the crest, paralleling the position of the phreatic surface within the core in Figure 5-9b, and intersecting the tailwater at approximately elevation 573 metres directly beneath the downstream edge of the crest in Figure 5-9c. The curvature of the crest is accounted for in the numerical model, as shown by the irregular zone boundaries in the longitudinal sections. Hydraulic conductivity values for all zones defined in the numerical model were derived from measured data or reported estimates, and are listed in Table 5-5. A B i i i i A\" B' a) transverse cross-section b) longitudinal cross-section A-A' 9 1 c) longitudinal cross-section B-B' Figure 5-9: Mica dam seepage conditions at high pool (754.4 m). Table 5-5: Numerical modelling input parameters: Mica dam. K L °sat u^nsat Zone (cm/s) (A/m 2) (uS/cm) (uS/cm) Shell 1 x 10° 1.8 x 10\"5 18 6 Core 1 x icr 6 2.0 x 10\"5 100 27 Overburden 1 x 10\"3 1.8 x 10 -5 18 6 Spoil fill (berm) 1 x 10\"3 1.8 x 10\"5 18 6 Transition 5 x If/3 1.8 x 10\"5 18 6 Bedrock 1 x lO\" 1 0 0 3 3 Reservoir water - - arf= 150 -A comparison of hydraulic piezometer readings with modelled hydraulic head data provides a simple means o f confirming that the numerical model approximates the hydraulic head distribution through the core of the dam. This comparison by no means constitutes calibration of the model, nor does it take into account transient effects, flow from the abutments, or the potential influence of piezometer lead trenches on the measured data. Given the sparse sampling offered by piezometers within the dam, an accurate and comprehensive depiction o f hydraulic conditions is difficult to deduce from the observed data alone. 92 T a b l e 5 - 6 d i s p l a y s m e a s u r e d a n d m o d e l l e d h y d r a u l i c h e a d d a t a at h i g h a n d l o w p o o l f o r a s e r i e s o f p i e z o m e t e r s l o c a t e d at d i f f e r e n t e l e v a t i o n s w i t h i n t h e c o r e . T h e m e a s u r e d d a t a are a p p r o x i m a t e a n d r e p r e s e n t the a v e r a g e o f d i s c r e t e r e a d i n g s c o r r e s p o n d i n g to t h e s p e c i f i e d p o o l l e v e l . P r e f e r e n c e w a s g i v e n to d a t a c o m p i l e d w i t h i n t h e p a s t d e c a d e to m o s t c l o s e l y a p p r o x i m a t e \" s t e a d y - s t a t e \" c o n d i t i o n s a r i s i n g f r o m t h e c u r r e n t l i m i t e d f l u c t u a t i o n i n p o o l l e v e l as c o m p a r e d t o t h e f i r s t f i l l i n g o f t h e r e s e r v o i r . A l t h o u g h p l a n c o o r d i n a t e s are n o t s t a t e d i n t h e t a b l e , t h e e q u i v a l e n t d i s t a n c e p a r a l l e l t o t h e c r e s t i s s t a t e d i n t e r m s o f l o n g i t u d i n a l s t a t i o n i n b o t h m e t r e s a n d feet. T h e s i m i l a r i t y b e t w e e n m e a s u r e d a n d m o d e l l e d p i e z o m e t r i c d a t a s u g g e s t s that t h e n u m e r i c a l m o d e l a p p r o x i m a t e l y r e p r e s e n t s g e n e r a l s e e p a g e c o n d i t i o n s . Table 5-6: Core piezometer data (all values in metres unless otherwise noted). Longitudinal Hi£ ,h pool Low pool Piezometer Elevation station on crest {754 to 754.4 m} {730-731.4 m} (m/ft) Measured modelled measured modelled PHI 590.7 686/22+50 656 654 653 645 PH2 590.7 686/22+50 648 633 646 627 PH6 650.7 686/22+50 695 700 694 689 PH7 650.7 686/22+50 693 686 692 678 PH9 677 930/30+50 697 695 696 689 PH10 677 503 / 16+50 697 697 694 691 PH12 718 930/30+50 732 720 724 720 5 . 2 . 3 E l e c t r i c a l F l o w A n a l y s i s T h e S P d i s t r i b u t i o n w a s c a l c u l a t e d t h r o u g h t h e d a m f o r b o t h h i g h p o o l a n d l o w p o o l l e v e l s u s i n g e s t i m a t e d i n p u t p a r a m e t e r s d e f i n e d i n T a b l e 5 - 5 . T h e s e p a r a m e t e r s r e p r e s e n t a s i m p l e a p p r o x i m a t i o n t o t h e t r u e d i s t r i b u t i o n , a n d p r o v i d e a n i n i t i a l a t t e m p t at c h a r a c t e r i z i n g t h e s u b s u r f a c e . A l t h o u g h t h e r e s e r v o i r f l u i d c o n d u c t i v i t y h a s b e e n o b s e r v e d to v a r y s e a s o n a l l y , a n a v e r a g e v a l u e o f 1 5 0 p S / c m w a s c h o s e n f o r t h e a n a l y s i s . C r o s s - c o u p l i n g c o n d u c t i v i t y v a l u e s w e r e a p p r o x i m a t e d f r o m p u b l i s h e d d a t a c o n t a i n e d i n A p p e n d i x C . S a t u r a t e d a n d u n s a t u r a t e d b u l k c o n d u c t i v i t y v a l u e s f o r t h e c o r e w e r e d e r i v e d f r o m d c r e s i s t i v i t y s u r v e y d a t a o b t a i n e d a l o n g t h e c r e s t o f t h e d a m . 9 3 Figure 5-10 displays the SP distribution over the upstream and downstream face of the dam at high pool level, with data reported with respect to a base electrode position at a longitudinal distance of 705 metres (23+12 feet) on the upstream edge of the crest. The significant variation in SP response perpendicular and parallel to the crest direction is apparent, especially on the downstream face of the dam. Figure 5-11 displays the surface SP data at both high and low reservoir elevations in transverse cross section, intersecting the crest at a distance of 700 metres (22+96 feet), as indicated by section C-C in Figure 5-10. The purpose of modelling the SP response to seepage through the embankment at representative high and low reservoir elevations is to isolate variations in SP caused by the change in seepage conditions. The low pool data is subtracted from the high pool data to generate an SP difference data set, which can then be analysed for seepage-related anomalies with external sources of noise minimized. In the numerical analysis of Mica dam, complexities inherent to the dam site, such as instrumentation and buried structures, are not included in the model. Consequently, SP difference data is most effective for comparing numerical and measured results. A plot of the modelled SP difference may be compared to that calculated from field survey data acquired at high and low pool levels, measured with respect to a base electrode stationed at the same position assumed in the numerical modelling. Measured and predicted SP difference profiles are presented in Figure 5-11 in transverse cross-section. Measured data are the product of field surveys performed in May 1998 (731.4 m) and August 2000 (749.4 m), and represent an 18-metre change in reservoir elevation. The modelled difference profile predicts a maximum SP difference of 40 mV corresponding to a 23-metre change in reservoir elevation. Although the magnitude of the measured difference profile is less than one half the magnitude of the modelled difference, the general trend with distance downstream is very similar. The modelled SP difference data over the entire surface of the dam is displayed in Figure 5-12 to show the significant variation in response in both longitudinal and transverse directions. 94 Transverse distance (m) Figure 5-10: Modelled surface SP distribution at high pool (754.4 m). 140 -I -60 -I 1 1 1 1 . . 1 1 . 1 1 1 1 1 1 1 1 1-0 100 200 300 400 500 600 700 800 900 Transverse distance (m) Figure 5-11: Modelled surface SP transverse profiles (section C-C in Figure 5-10). 9 5 50 Figure 5-12: Modelled surface SP difference data (754.4 m - 731.4 m). 5.2.4 Discussion of Results The comparison of measured and modelled hydraulic head and SP data at M i c a dam suggests that a 3-D numerical model may be used to successfully approximate observed behaviour in an embankment. However, inconsistencies in these results highlight the implications of certain assumptions in the numerical model, which may require further consideration to enable a comprehensive assessment of field conditions. The magnitude of the predicted SP response is largely influenced by the input parameters chosen to represent the physical properties of each zone. The physical parameters used in the numerical modelling were chosen to reflect average measured values, i f available, assuming homogeneous and isotropic conditions within each soil unit. A more comprehensive characterization of physical property distribution could lead to dramatically 9 6 different results, by defining the overall magnitude of the predicted response and perhaps the range in response expected due to the inherent non-uniformity of given soil unit. The hydraulic head distributions used in the electrical flow analysis were not the product of a fully calibrated seepage assessment, and may not represent actual flow conditions within the embankment, abutments and foundation. Although the modelled hydraulic regime within the core appears to compare with piezometric data, hydraulic conditions may not be correctly represented by a free-surface approach. The influence of an accurate characterization of hydraulic conditions on the character of the SP response was illustrated in Chapter 4, and may help to explain some of the.discrepancy between measured and modelled results. The elevated magnitude of the modelled SP difference profile relative to the observed difference profile shown in Figure 5-11 may be attributed in part to the greater difference in pool elevation used in the numerical analysis, but is likely also a function of the steady-state approximation of hydraulic flow conditions. Neglecting transient effects assumes a maximum hydraulic response corresponding to a given reservoir elevation. Consequently, it is of interest to estimate the lag time of a change in reservoir level at a point within the dam core. Applying Equation 2-25 in Section 2.3.4 with a length scale of 25 to 50 metres representing the maximum width variation of the core, translates to a time to equilibrium on the order of several months to over one year. Therefore, at the time of a field survey, the hydraulic conditions within the core of the embankment are likely consistent with an earlier pool level. This is substantiated by continuously monitored SP data collected from two distinct pairs of electrodes oriented parallel to and orthogonal to the crest. The crest electrodes are placed 200 metres apart, as indicated by points Tl and T2 in Figure 5-7, and the orthogonal set are 130 metres apart on the downstream face of the dam, as indicated by points T3 and T4. Figure 5-13 displays the measured SP data and the reservoir elevation plotted versus time. Although a distinct time lag is evident between the profile of reservoir elevation and the SP profiles, it is difficult to quantify the delay without further examination. This may include an analysis of piezometric data, a longer period of measured SP data, and information on the variability of the electrical conductivity of the reservoir water and embankment materials. 97 Preliminary SP distributions calculated over the surface of M i c a dam at high pool and low pool levels illustrate the general behaviour of the observed response, but further analysis is required to accurately characterize the complexities of a field site. 140 120 100 > 80 E CL w 60 40 20 0 01-May 98 Reservoir level 31- 31- 30- 02- 01- 01- 31- 02- 01- 01- 31- 03- 02-Jul- Oct- Jan- May- Aug- Nov- Jan- May- Aug- Nov- Jan- May- Aug-98 98 99 99 99 99 00 00 00 00 01 01 01 Figure 5-13: Continuous SP response recorded from installed electrodes. 5.3 Summary and Conclusions Chapter 5 presented further verification of the numerical procedure 3 D S P through a comparison with measured values in a controlled physical model of an embankment dam. Two model configurations, consisting of a homogeneous dam and a dam with upstream defect, were studied at two reservoir levels. Significant changes in the SP response resulted from a decrease in reservoir level, which was mirrored in both the measured and predicted surface profiles. However, inclusion of the defect did not translate to a significant change in the SP distribution along the downstream face of the embankment. A coupled flow analysis of M i c a dam at representative high and low reservoir elevations illustrated the application of the numerical procedure to an actual embankment, and highlighted the advantage of 3 -D modelling to capture the variability in the SP response at surface in directions both parallel and perpendicular to the crest. Calculation of the SP difference between high and low pool enabled a comparison with measured data, which 9 8 showed encouraging results that suggest the procedure is capable of modelling the general flow behaviour within a field embankment. Key conclusions and contributions include the following: • Effective characterization of physical properties and the hydraulic flow regime are key to a successful prediction of the SP response to seepage using a forward modelling procedure. Even within a controlled laboratory experiment, accurate representation of input parameters is challenging. • Physical property measurements performed on samples not representative of actual conditions due to sample size, saturation and preparation techniques may lead to unrealistic estimates of field behaviour. • Variations apparent in the numerical modelling results, caused by a change in input parameters, highlight the importance of measuring the electrical conductivity along with 1-D SP coupling coefficient measurements, and replicating in-situ soil and fluid conditions. • Measured and modelled results over the laboratory embankment exhibit remarkably similar trends and order of magnitude response. The predicted SP difference signature over the surface of Mica dam shares similar trends with measured data, despite the limited complexity of the numerical model. • The interpretation of the SP response to seepage over an actual embankment requires consideration of transient effects, and a comprehensive assessment of seepage conditions. Further evaluation of the influence of model assumptions on the predicted response may be required to accurately represent field behaviour. 99 6 CONCLUSIONS AND RECOMMENDATIONS 6.1 Summary The main objective of this thesis was to explore the effect of changing seepage conditions on the self-potential response to flow in porous media. This was accomplished with the development of a 3-D numerical modelling procedure, a theoretical analysis of the effect of internal erosion on the physical parameters controlling the SP response, and application of the numerical procedure to actual seepage problems. Chapter 2 provided an overview of the physical processes central to the study of the streaming potential phenomenon in soils. Included were a description of the physical process of internal erosion, an introduction to electrokinetic processes and the phenomenon of streaming potential, and a theoretical discussion of hydraulic and electrical flow in porous media. These provided the framework for an analysis of the effect of changing seepage conditions on streaming current flow and the response of the self-potential method to the phenomenon. The link between streaming current and the hydraulic gradient and electrical parameters was established in Chapter 2 with the development of the coupled flow equations. These governing equations are central to the mathematical model on which the numerical procedure 3 D S P is based. The development of the three-dimensional finite difference forward modelling procedure described in Chapter 3 enables characterization of the SP response to complex hydraulic systems. This capability was illustrated in a series of numerical analyses presented in Chapters 4 and 5. These included: • a sensitivity analysis of the SP response to internal erosion in a field-scale hypothetical zoned embankment with piped zone through the core; • verification of the numerical procedure through a comparison with measured SP data obtained in a laboratory-scale physical model of a dam; and • application of 3 DSP to characterize the surface SP distribution in preliminary models of Mica dam at high and low reservoir elevations. 1 0 0 6.2 Conclusions Forward modelling of the SP response to seepage relies on an accurate representation of in-situ hydraulic flow conditions. Since the mathematical model describing streaming current flow is defined by the hydraulic gradient, confidence must be placed in the calculated hydraulic head distribution used as input to the electrical flow analysis. Multidimensional analyses in Chapters 4 and 5 demonstrate that the finite volume method can be used to effectively model seepage in embankment dams. Only software or code that can be shown to accurately model the hydraulic head distribution should be used in a coupled electrokinetic flow analysis. The sensitivity analyses presented in Chapter 4 , and the results of the laboratory dam study illustrate the impact of the hydraulic potential distribution on the character of the SP distribution, and the influence of electrical properties on the magnitude of the SP response. Changes in hydraulic conductivity caused by the onset of internal erosion that influence the hydraulic flow pattern appear to have the greatest impact on the SP response. This was shown in the 3-D sensitivity analysis, where the inclusion of a hypothetical erosion channel influenced the hydraulic and electrical flow regimes to result in a SP anomaly at surface. The inclusion of an upstream pervious defect in the laboratory embankment resulted in only a slight perturbation in the hydraulic flow pattern, and consequently did not translate to a significant change in the SP distribution at surface. However, effective characterization of the resistivity distribution was shown to be essential in predicting an accurate SP response, and is consequently key to any meaningful interpretation of SP data. The preliminary modelling of Mica dam illustrates the value of 3-D modelling to effectively characterize the SP distribution in an embankment. More comprehensive analyses and further assessment of the impact of model assumptions on the predicted SP response will only improve our capacity to interpret field results. Due to the influence of the hydraulic regime on the SP distribution, a calibrated 3-D analysis may benefit from an attempt to capture not only the seepage patterns caused by the upstream reservoir, but other flow sources that may influence the SP response, such as flow from the abutments and the foundation, and the effect of precipitation and flow through the unsaturated zone. Although the contribution of these flow sources may be minor in the assessment of volumetric seepage, the SP response to these flows may prove significant. 101 6.3 Recommendations Self-potential (SP) measurements are being used to identify possible zones of anomalous seepage in earth dams. The main benefit of this method over other techniques is its simplicity and suitability for long term, continuous monitoring. However, in order for the SP method to be truly practical for geotechnical surveillance, quantitative interpretation of soil properties and seepage flow is required. Numerical modelling will only improve our ability to interpret SP data, but its application is limited by our knowledge of the hydraulic and electrical properties of the subsurface. The three main physical parameters controlling the SP response; namely hydraulic conductivity, electrical conductivity, and the cross-coupling conductivity; are each a function of fundamental soil and fluid properties. While significant research has been conducted in an effort to characterize and study hydraulic and electrical conductivity for a wide-range of geologic materials, a similar collection of information on the streaming current cross-coupling coefficient has yet to be compiled. Furthermore, the effect of internal erosion on all three parameters is not well understood. This lack of data limits our ability to apply numerical modelling techniques to field-scale problems. In order to determine the viability of SP as a standard technique in seepage flow monitoring, further research into the nature and generation of self-potential anomalies within the framework of earth dam configurations must be conducted. A physical understanding of streaming potential and the influence of material characteristics on the SP response could be assessed through a laboratory analysis involving electrical measurements on unconsolidated sediments used in earth dam construction. The ultimate objective would be to establish a quantitative link between measured SP readings and actual seepage conditions, and to assess the viability of SP as a monitoring tool to detect internal erosion in earth dams. Experimental research focused on streaming potential measurements should proceed in both of the following areas: 1. Dependence of the cross-coupling coefficient, SP coupling coefficient, and bulk electrical conductivity on soil and fluid properties. Properties of the fluid, properties of the solid, and activity at the solid-fluid interface have an effect on the streaming current and measured streaming potential. Variability of these parameters and their effect on the measured cross-coupling coefficient is of prime 102 concern in the accuracy and reproducibility of measurement for a given geologic material and saturating fluid. The sensitivity of the electrical parameters to changes in soil structure should be studied using soils and saturating fluids representative of in-situ conditions, and wi l l help to ascertain the validity of theoretical estimates. 2. Characterization of the SP cross-coupling coefficient for a variety of soils used in earth dam construction. Soil samples representative of those used in dam construction could be created in the laboratory based on design criteria, or acquired directly from drilling programs at existing projects. The development of a database of cross-coupling coefficients for different soils would provide further insight into the degree of variability and range of values expected, and determine the importance of accurate characterization of this parameter for use in field SP data interpretation. If accurate characterization of this parameter is considered significant for field applications, the development of in-situ methods of measurement should ensue. Experimental research focused on the characterization of electrical input parameters wi l l lead to a better understanding of the streaming potential phenomenon in soils, and may provide insight to means of improving the mathematical coupled flow model. Further analyses of embankments with well-characterized physical properties w i l l help to ascertain whether the current model is capable of accurately characterizing flow behaviour. Forward modelling is useful for predictive analyses, and can provide insight into the magnitude and behaviour of the SP response in different geologic and hydrologic settings. However, for numerical modelling to develop into an efficient and effective tool to aid in the interpretation of measured data, the application of inverse modelling to the SP coupled flow problem in embankment dams should be also be pursued. 103 REFERENCES Adamson, A.W. 1990. Physical chemistry of surfaces, 5m ed. John Wiley and Sons, New York. Ahmad, M.U. 1964. A laboratory study of streaming potentials. Geophysical Prospecting, 12: 49-64. Anderson, M.P., and Woessner, W.W. 1992. Applied groundwater modeling: Simulation of flow and advective transport. Academic Press, San Diego. Archie, G.E. 1942. The electrical resistivity log as an aid in determining some reservoir characteristics. Transactions of the Society of Petroleum Engineers of the American Institute of Mining, Metallurgical and Petroleum Engineers (AIME), Vol. 146, p.p. 54-67. Arnold, M. 1973. Laboratory determination of the coefficient of electro-osmotic permeability of a soil. Geotechnique 23: 581-588. ASTM. 1994. Standard method for field measurement of soil resistivity using the Wenner four-electrode method (G57-78). In 1994 Annual book of ASTM standards, Sect. 3, Vol. 03.02. American Society for Testing and Materials, Philadelphia, p.p. 218-222. ASTM. 1994. Standard specification for concrete aggregates (C33-93). In 1994 Annual book of ASTM standards, Sect. 4, Vol. 04.02. American Society for Testing and Materials, Philadelphia, p.p. 10-16. ASTM. 1994. Standard test method for measuring the soil-geotextile system clogging potential by the gradient ratio (D5101-90). In 1994 Annual book of ASTM standards, Sect. 4, Vol. 04.09. American Society for Testing and Materials, Philadelphia, p.p. 180-186. Barber, E.S., and Sawyer, CL. 1952. Highway subdrainage. Proceedings, Highway Research Board, p.p. 643-666. BC Hydro. 2001. Internal reports, British Columbia Hydro and Power Authority. Bear, J. 1972. Dynamics of fluids in porous media. Dover Publications, New York. Beckie, R. 2001. Personal communication. Beckie, R. 2000. Course notes, GEOL565: Advanced groundwater hydrology. Department of Earth and Ocean Sciences, University of British Columbia. 104 Berube, A.P. 2000. Finite element modelling of coupled flows for the interpretation of the streaming potential phenomenon in embankment dams. Report, Lulea University of Technology, Sweden. Burke, H.H., Content, C.S., and Kulesza, R.L. 1972. Current practice in abutment and foundation treatment. Journal of the Soil Mechanics and Foundations Division, American Society of Civil Engineers, 98(SM10): 1033-1052. Casagrande Consultants. 1982. Laboratory investigation on effectiveness of electro-osmosis to improve weak foundation soils at proposed Greens Creek tailings dam, Admiralty Island, Alaska. Report for Steffen, Robertson and Kirsten (B.C.) Inc., Vancouver, B.C. Casagrande, L. 1983. Stabilization of soils by means of electro-osmosis - state of the art. Journal of the Boston Society of Civil Engineers, American Society of Civil Engineers, 69(2): 255-302. Casagrande, L. 1949. Electro-osmosis in soils. Geotechnique 1: 159-177. Cedergren, H.R. 1993. Improving seepage safety of earth dams. Proceedings, 1993 Canadian Dam Safety Association Conference, St. John's, Nfld., September 19-23, p.p. 151-164. Cedergren, H.R. 1973. Seepage control in earth dams. In Embankment dam engineering, Casagrande volume. John Wiley and Sons, New York. Corwin, R.F. 1991. Evaluation of effects of cutoff wall construction on seepage flow using self-potential data, east embankment Wells dam. Report for Douglas County Public Utility District No. 1, Washington. Corwin, R.F. 1990. The self-potential method for environmental and engineering applications. Geotechnical and Environmental Geophysics, Society of Exploration Geophysics, 1: 127-145. Courage, L.J.R., Peterson, T.W.P., Rothbauer, A.L., and Pelz, B.A. 1997. Whiteman's earthfill dam internal erosion. Transactions, 19th Congress on Large Dams, Florence, Vol. 2. International Commission on Large Dams, p.p. 9-30. Culley, R.W., Jagodits, F.L., and Middleton, R.S. 1975. E-Phase system for detecting buried granular deposits. Innovations in Subsurface Exploration of Soils, Transportation Research Record 581. Transportation Research Board, National Research Council, Washington. De Groot, S.R. 1951. Thermodynamics of irreversible processes. Selected topics in modern physics, Volume 3. Edited by J. de Boer, H. Brinkman, and H.B.G. Casimir. North Holland Publishing Company, Amsterdam. 105 Domenico, P.A., and Schwartz, F.W. 1998. Physical and chemical hydrogeology, 2na ed. John Wiley and Sons, New York. Drever, J.I. 1988. The geochemistry of natural waters, 2nd ed. Prentice Hall, New Jersey. Fenn, CD. 1966. The effects of montmorillonite and kaolinite dispersions on the permeability of porous media. M.S. thesis, Mississippi State University, State College, Mississippi. Fitterman, D.V. 1978. Electrokinetic and magnetic anomalies associated with dilatant regions in a layered earth. Journal of Geophysical Research, 83(B12): 5923-5928. Freeze, R.A., and Cherry, J.A. 1979. Groundwater. Prentice-Hall, New Jersey. Freeze, R.A. 1971. Influence of the unsaturated flow domain on seepage through earth dams. Water Resources Research, 7(4): 929-941. Friborg, J. 1996. Experimental and theoretical investigations into the streaming potential phenomenon with special reference to applications in glaciated terrain. Ph.D. thesis, Lulea University of Technology, Sweden. Gerald, C.F., and Wheatley, P.O. 1989. Applied numerical analysis, 4th ed. Addison-Wesley Publishing Company, Reading. Golub, G.H., and Ortega, J.M. 1992. Scientific computing and differential equations: an introduction to numerical methods. Academic Press, Boston. Gray, D.H., and Mitchell, J.K. 1967. Fundamental aspects of electro-osmosis in soils. Soil Mechanics and Bituminous Materials Research Laboratory, University of California, Berkeley. Guiger, N., and Franz, T. 1997. User's manual for Visual MODFLOW. Waterloo Hydrogeologic Inc., Waterloo. Haber, E. 2001. Practical aspects of running EH3D. Technical Note, Department of Earth and Ocean Sciences, University of British Columbia. Hameiri, A. 2000. Soil/geotextile filtration behaviour under dynamic conditions of vibration and cyclic flow. Ph.D. thesis, University of British Columbia, Vancouver. Hem, J.D. 1985. Study and interpretation of the chemical characteristics of natural water, 3rd ed. U.S. Geological Survey Water-Supply Paper 2254, Department of the Interior, United States Government Printing Office. Hibbert, D.B. 1993. Introduction to electrochemistry. Macmillan Physical Science Series. Macmillan Press, London. 106 Honjo, Y., and Veneziano, D. 1989. Improved filter criterion for cohesionless soils. Journal of Geotechnical Engineering, 115(1): 75-94. Huntley, D. 1986. Relations between permeability and electrical resistivity in granular aquifers. Ground Water, 24(4): 466-474. ICOLD. 1995. Dam failures: statistical analysis. International Commission on Large Dams, Bulletin no. 99. ICOLD. 1994. Embankment dams: granular filters and drains. International Commission on Large Dams, Bulletin no. 95. Ishido, T., and Mizutani, H. 1981. Experimental and theoretical basis of electrokinetic phenomena in rock-water systems and its application to geophysics. Journal of Geophysical Research, 86(B3): 1763-1775. Ishido, T., and Tosha, T. 1998. Feasibility study of reservoir monitoring using repeat self-potential measurements. Transactions, Geothermal Research Council Meeting, San Diego, September 20-23, Vol. 22, p.p. 171-178. Jansen, R.B., ed. 1988. Advanced dam engineering for design, construction and rehabilitation. Van Nostrand Reinhold, New York. Jones, P.H. and Buford, T.B. 1951. Electric logging applied to ground-water exploration. Geophysics 16(1): 11-139. Jouniaux, L., and Pozzi, J.P. 1995. Streaming potential and permeability of saturated sandstones under triaxial stress: consequences for electro telluric anomalies prior to earthquakes. Journal of Geophysical Research 100(B6): 10197-10209. Keller, G.V., and Frischknecht, F.C. 1966. Electrical methods in geophysical prospecting. International series of monographs in electromagnetic waves, Vol. 10. Edited by A.L. Cullen, V.A. Fock, and J.R. Wait. Pergamon Press, New York. Kelly, W.E. 1977. Geoelectric sounding for estimating aquifer hydraulic conductivity. Ground Water, 15(6): 420-425. Kenney, T.C., Chalal, R, Chiu, E., Ofoegbu, G.L, Omange, G.N., and Ume, CA. 1985. Controlling constriction sizes of granular filters. Canadian Geotechnical Journal, 22(1): 32-43. Kenney, T.C. and Lau, D. 1985. Internal stability of granular filters. Canadian Geotechnical Journal, 22(2): 215-225. Kovacs, G. 1981. Seepage hydraulics. Developments in water science, Vol. 10. Elsevier, Amsterdam. 107 Lafleur, J., Mlynarek, J., and Rol l in , A . L . 1989. Filtration of broadly graded cohesionless soils. Journal of Geotechnical Engineering, 115(12): 1747-1768. Marshall, D.J . and Madden, T.R. 1959. Induced polarization, a study of its causes. Geophysics, 24(4): 790-816. Marsily, G . de 1986. Quantitative hydrogeology: groundwater hydrology for engineers. Academic Press, Orlando. Masch, F .D . and Denny, K . J . 1966. Grain size distribution and its effect on the permeability of unconsolidated sands. Water Resources Research, 2(4): 665-677. Mazac, O., Kel ly , W . E . , and Landa, I. 1985. A hydrogeophysical model for relations between electrical and hydraulic properties of aquifers. Journal of Hydrology, 79: 1-19. M c N e i l l , J .D. 1980. Electrical conductivity of soils and rocks. Technical Note TN-5 , Geonics Limited. Mil ler , D . G . 1960. Thermodynamics of irreversible processes - the experimental verification of the Onsager reciprocal relations. Chemical Reviews, 60(1): 15-37. Mitchell , J. K . 1991. Conduction phenomena: from theory to geotechnical practice. Geotechnique, 43(3): 299-340. Morgan, F .D . 1989. Fundamentals of streaming potentials in geophysics: laboratory methods. In Detection of subsurface flow phenomena. Edited by G.-P. Merkler, H . Mili tzer, H . Hotzl , H . Armbruster, J. Brauns. Lecture notes in earth sciences, V o l . 27, Springer-Verlag, Berlin, p.p. 133-144. Morgan, F .D. , Will iams, E.R. , and Madden, T.R. 1989. Streaming potential properties of Westerly granite with applications. Journal of Geophysical Research, 94(B9): 12449-12461. Nourbehecht, B . 1963. Irreversible thermodynamic effects in inhomogeneous media and their applications in certain geoelectric problems. Ph.D. thesis, Massachusetts Institute of Technology, Cambridge. Ogilvy, A . A . , Ayed, M . A . , and Bogoslovsky, V . A . 1969. Geophysical studies of water leakages from reservoirs. Geophysical Prospecting, 17: 36-62. Olsen, H . W . 1969. Simultaneous fluxes of liquid and charge in saturated kaolinite. Soil Science Society of America Proceedings 33(3): 338-344. Onsager, L . 1931. Reciprocal relations in irreversible processes, I. Physical review, 37: 405-426. 108 Overbeek, J.Th.G. 1952. Electrokinetic phenomena. In Colloid science: irreversible systems, Volume 1. Edited by H.R. Kruyt. Elsevier Publishing Company, Amsterdam. Parks, G.A. 1967. Aqueous surface chemistry of oxides and complex oxide minerals. In Equilibrium Concepts in Natural Water Systems, Advances in Chemistry Series, Vol. 67. Edited by R.F.Gould, p.p. 121-160. Parks, G.A. 1965. The isoelectric points of solid oxides, solid hydroxides, and aqueous hydroxo complex systems. Chemical Reviews, 65: 177-198. Pengra, D.B., Li, S.X., and Wong, P. 1999. Determination of rock properties by low-frequency AC electrokinetics. Journal of Geophysical Research, 104(B12): 29485-29508. Remson, I., Hornberger, G.M., and Molz, F.J. 1971. Numerical methods in subsurface hydrology. John Wiley and Sons, New York. R.L.&L. Environment Services Ltd. 1999. Brilliant expansion project - summary of aquatic inventory data of the Kootenay system near Brilliant dam. Report prepared for Columbia Power Corporation, Castlegar, B.C. Romano, C.G., Frind, E.O., and Rudolph, D.L. 1999. Significance of unsaturated flow and seepage faces in the simulation of steady-state subsurface flow. Ground Water 37(4): 625-632. Sasitharan, S., Sheffer, M.R., and Gaffran, P.C. 2001. An investigation of streaming potential and its application to earthfill dams. Proceedings, Canadian Dam Association Annual Conference, Fredericton, N.B., September 30 - October 4, p.p. 135-142. Sato, H., Shima, H., and Toshioka, T. 1998. A trial application of streaming potential to detecting underground water movement while pumping up. Proceedings, European Association of Exploration Geophysicists 60th Conference and Technical Exhibition, Leipzig, Germany, June 8-12, Expanded Abstract 10-08. Schuch, M., and Wanke, R. 1967. Stromungsspannungen in einigen torf- und sandproben. Zeitschrift fur Geophysik 33: 94-109. Sherard, J.L. 1979. Sinkholes in dams of coarse, broadly graded soils. Transactions, 13th Congress on Large Dams, New Delhi, India, Vol. 2. International Commission on Large Dams, p.p. 25-35. Sherard, J.L., and Dunnigan, L.P. 1989. Critical filters for impervious soils. Journal of Geotechnical Engineering, 115(7): 927-947. Sherard, J.L., and Dunnigan, L.P. 1985. Filters and leakage control in embankment dams. Proceedings, Symposium on Seepage and Leakage from Dams and Impoundments, American Society of Civil Engineers, May 1985, p.p. 1-30. 109 Sherard, J.L., Woodward, R.J., Gizienski, S.F., and Clevenger, W.A. 1963. Earth and earth-rock dams. John Wiley and Sons, New York. Sill, W.R. 1983. Self-potential modeling from primary flows. Geophysics, 48(1): 76-86. Sill, W.R., and Killpack, T.J. 1982. SPXCPL: Two-dimensional modeling program of self-potential effects from cross-coupled fluid and heat flow. User's guide and documentation for version 1.0. Earth Sciences Laboratory, University of Utah Research Institute, Salt Lake City. Report for U.S. Department of Energy, Division of Geothermal Energy. Sprunt, E.S., Mercer, T.B., and Djabbarah, F. 1994. Streaming potential from multiphase flow. Geophysics 59(5): 707-711. Stewart, R.A., and Watts, B.D. 2000. The WAC Bennett dam sinkhole incident. Proceedings, 53rd Canadian Geotechnical Conference, Montreal, P.Q., October 15-18, p.p. 39-45. Stumm, W. 1992. Chemistry of the solid-fluid interface: processes at the mineral-water and particle-water interface in natural systems. John Wiley and Sons, New York. Taylor, D.W. 1948. Fundamentals of soil mechanics. Wiley, New York. Telford, W.M., Geldart, L.P., and Sherriff, R.E. 1990. Applied geophysics, 2nd ed. Cambridge University Press, New York. Tuman, V.S. 1963. Thermo-telluric currents generated by an underground explosion and other geologic phenomena. Geophysics 28(1): 91-98. USACE. 1986. Engineering and design - seepage analysis and control for dams. United States Department of the Army, Corps of Engineers, Engineer Manual no. 1110-2-1901. Vaughan, P.R., and Soares, H.F. 1982. Design of filters for clay cores of dams. Journal of the Geotechnical Engineering Division, American Society of Civil Engineers, 108(GT1): 17-31. Watts, B.D., Gaffran, P.C, Stewart, R.A., Sobkowicz, J.C, and Kupper, A.G. 2000. WAC Bennett dam - characterization of sinkhole no. 1. Proceedings, 53r Canadian Geotechnical Conference, Montreal, P.Q., October 15-18, p.p. 67-75. Waxman, M.H., and Smits, L.J.M. 1968. Electrical conductivities in oil-bearing shaly sands. Transactions of the Society of Petroleum Engineers of the American Institute of Mining, Metallurgical and Petroleum Engineers (AIME), Vol. 243, Part 2, p.p. 107-122. 110 Wilt , M . J . , and Butler, D . K . 1990. Numerical modeling of SP anomalies: documentation of program SPPC and applications. In Geotechnical applications of the self potential (SP) method, Report 4. Department of the Army, U S Army Corps of Engineers, Technical Report R E M R - G T - 6 . Wurmstich, B . 1995. 3-D self-consistent modeling of streaming potential responses: theory and feasibility of applications in earth sciences. Ph.D. thesis, Texas A & M University. Wurmstich, B . , and Morgan, F .D . 1994. Modeling of streaming potential responses caused by oi l well pumping. Geophysics 59(1): 46-56. Wurmstich, B . Morgan, F .D. , Merkler, G.-P., and Lytton, R . L . 1991. Finite element modeling of streaming potentials due to seepage: study of a dam. Proceedings, Society of Exploration Geophysicists Annual Meeting, Expanded Technical Abstracts. Yasukawa, K . 1993. A coupled self potential (SP), fluid, and heat flow model for subsurface flow systems. M . S . thesis, University of California, Berkeley. Yeung, A . T . , and Mitchell , J .K. 1993. Coupled fluid, electrical and chemical flows in soil. Geotechnique, 43(1): 121-134. 111 APPENDIX A - Properties of the Solid-Fluid Interface The solid surfaces within a porous medium saturated with an electrolytic solution display a net electrical charge, either positive or negative, which is controlled by surface chemistry. Factors influencing the polarity of this charge include mineralogy or chemical composition of the solid particles; pH of the pore fluid; imperfections in the molecular or lattice structure of the solid, specifically in clay minerals; concentration and type of dissolved ions in the pore fluid; and adsorption of organic material (Stumm, 1992). Simple oxide minerals such as quartz exhibit a limited capacity for ion exchange, which arises from the presence of a pH-dependent surface charge. The chemical functional groups existing at the surface of silicate minerals, such as quartz and feldspar, are hydroxyl groups (e.g. Si-OH, Al-OH), which tend to exhibit a positive charge at low pore fluid pH, and a negative charge at high pH. Other types of minerals behave similarly, but employ different surface functional groups (e.g. Ca-COsH in carbonate minerals). The fluid pH defining zero surface charge is termed the point of zero charge (PZC), and differs for each mineral. The isoelectric point (IEP) is an equivalent parameter, which describes the fluid pH defining zero charge at the shear plane within the electrical double layer, an interface that approximates the actual surface of the solid. A comprehensive listing of PZC and IEP values is contained in research conducted by Parks (1965), and IEP values obtained for several rocks and minerals are shown in Table A-l. Most oxides exhibit a negative surface charge when in contact with natural waters, which typically have a pH less than 7 (Parks, 1967). Table A- 1: Isoelectric points of geologic materials (from Ishido and Mizutani, 1981). IEP Material Room temperature 45°C Quartz 2.6 2.0 Orthoclase feldspar 2.5 1.0 Granite 1.5 Kaolinite =1 Montmorillonite =2 112 Under typical subsurface conditions, clay minerals also tend to exhibit a negative surface charge, which is amplified due to their layered structure. Although clay particles exhibit siliceous surface functional groups, the surface charge distribution is more complex and is enhanced by the weak ionic bonds between layers. A s a result, these complex oxide minerals exhibit a much higher capacity for ion exchange. Imperfections in the internal clay structure result from the substitution of positive ions in the lattice with other dissolved cations of lower valency, which imparts a net negative charge to the particle structure. This negative charge is then balanced by weakly bound positive cations at the particle surface, which are also prone to substitution and contribute to the high surface electrical conductivity typical in clays. Different types of dissolved ions and organic substances in the pore fluid may also influence surface charge. The presence of non-coordinating ions increases the ionic strength, or electrical conductivity of the fluid and, although this does not affect the P Z C , may serve to increase the magnitude of the surface charge. Significant amounts of other ions, such as metals and ligands, can affect the polarity or magnitude of surface charge. Surface charge w i l l be decreased i f there is a significant amount of dissolved organic carbon, which wi l l adsorb to the grain surfaces. A n electrical double layer forms when a charged solid surface comes into contact with an electrolytic solution. The most accepted model of the double layer charge distribution is the Stern model, which distinguishes between an adsorbed layer of positive ions and a more loosely bound diffuse outer layer (Adamson, 1990). Under static flow conditions, the system is electrically neutral with ions in the double layer completely balancing the surface charge of the soil grain. In a porous medium, the velocity of fluid flow tapers from the centre of the pore toward the edges of the capillary such that a thin layer of fluid with a thickness of one to several molecules can be considered stationary at the solid-fluid interface. This adsorbed layer is termed the Stern, or Helmholtz layer and is separated from the diffuse outer layer, or Gouy layer of ions by the outer Helmholtz plane (OHP), as shown in Figure A - l . The Stern layer is characterized by adsorbed water molecules at the solid surface, delineated by the Inner Helmholtz Plane (IHP), and attracted positive cations, which are typically surrounded by a sheath of water molecules, which sit at the Outer Helmholtz Plane (OHP). Irrespective of the polarity of the solid surface charge, negative 113 anions and neutral organic molecules may also specifically adsorb to the solid surface, the anions losing their solvation sheath of water molecules in the process. in diffuse layer Figure A - l : Stern model of the electrical double layer at the solid-fluid interface (from Hibbert, 1993). The Debye length A describes the thickness of the electrical double layer, which is a function of the electric permittivity s [F/m], temperature T [°K] and ionic strength I [mol/m ] of the fluid, as shown in Equation A-l, where k& is Boltzman's constant [1.380622 x 10\" J/K] and e is electronic charge [C] (Morgan et al, 1989). The ionic strength of the fluid is defined by the valency n and molecular concentration c [mol/m3] of each ionic species in solution. This relation shows that the thickness of the double layer is inversely proportional to the ionic concentration of the fluid, which controls electrical conductivity. ^ ek^T where, I = - V n 2 c (A-l) V 2e2I 2 ^ 114 The electrical charge on the OHP can be positive or negative, depending on the distribution of ions balancing the negative surface charge of the solid. The slipping plane between the Helmholtz layer and the diffuse outer layer of ions is the boundary between the fixed and flowing fluid. This shear plane may or may not coincide with the OHP, and will tend to sit at a slightly greater distance from the grain surface. The electrical potential at this interface is termed the zeta potential, C,, which can be derived from experimental measurements and is used to approximate the true surface potential of the solid. Ishido and Mizutani (1981) discuss the dependence of zeta potential on fluid chemistry and pH in porous media. A graph of electrical potential from the grain surface into the pore fluid is shown in Figure A-2, where the surface potential is considered to be positive. The upper curve in the diagram shows the drop in potential between the surface (i|/0), the OHP (vj/e), and finally the shear plane (Q. If strong contact adsorption of negative anions exists in the Stern layer, the OHP and zeta potential could change sign with respect to the surface potential, as shown in the lower curve of Figure A-2. Solvation .layer Figure A-2: Change in potential away from a solid surface with potential \\ j / 0 (from Hibbert, 1993). 115 APPENDIX B - Development of the Continuity Equations B-l Fundamental Equations A general form of the continuity equation describes the conservation of a physical quantity in terms of the divergence of a flux vector, T, and a volumetric accumulation or storage term, S, as follows: — V X = Continuity Equation (B-l) where, r = (rx \\, ry j , rz k), flux vector with units quantity 1 time volume _ Flux is defined as the rate of flow of a physical quantity across a surface. When a unit volume, or control volume, of material is considered, the flux at the surface bounding this volume is of interest. The divergence of flux is simply the net flux of a physical quantity within a unit volume where inflow is positive and outflow negative. Divergence is a measure of how much a vector spreads out, or diverges, from a given point. Consequently, the accumulation term represents the time rate of change of the total physical quantity contained within the control volume. If the influx of the physical quantity to the control volume is equal to the flux out of the control volume, the storage term S is zero and the system is steady-state. Alternately, a transient state results if there is a net storage or net loss of the physical quantity. If the inflow to the control volume is greater than the outflow, the divergence is negative, S is positive, and the volume is considered a flow sink. If the divergence is positive, S is negative and the volume is a flow source. It is important to note that although flux is defined at the boundary of the control volume, the divergence of flux is defined its centre. The storage term is therefore representative of the entire control volume and must be expressed as a function of this volume. Phenomenological laws, such as Darcy's law and Ohm's law have been developed through experimentation to describe irreversible flow as a linear relation between flux and a quantity 1 time area drx ery drz , . u . + + — - , scalar source with units dx dy dz 116 driving force. Flux can be expressed mathematically as the product of a conductivity term L , and the gradient of a scalar potential as seen in the following equation: r = -L.(VO) Phenomenological Law (B-2) where, L = (Lx i , Ly j , Lz k), conductivity vector VO = ( — i, — j,—k), gradient vector of the scalar quantity, j (B-26) A soil-fluid system subjected to an imposed gradient is forced out of static equilibrium, and the irreversible flow phenomena that arise as a consequence of this gradient can be studied quantitatively using the theory of non-equilibrium thermodynamics. In the study of coupled flow, the gradient terms are considered the thermodynamic forces and must be formulated properly in terms of potential to uphold thermodynamic principles (de Groot, 1951). Miller (1960) and later Yeung and Mitchell (1993) outlined the basic steps required for the development of coupled flow equations that satisfy the laws of non-equilibrium thermodynamics, assuming the macroscopic flow laws apply to a given system, and perturbations from equilibrium are small such that a state of local equilibrium can be assumed for infinitesimal volumes that comprise the larger system. The formalism of non-equilibrium thermodynamics is briefly summarized here as 4 distinct steps: 1) Define the flux and gradient terms through a dissipation function, which describes the internal entropy production of an irreversible system at a given temperature; 2) Formulate the coupled flow equations using Equation B-26; 3) Apply the Onsager reciprocal relations, L/,e = heh to define the cross-coupling conductivity terms (Onsager, 1931); 4) Relate all phenomenological coefficients L, and L,y to common geotechnical engineering parameters. The basic entropy equation is stated in Equation B-27, where S is entropy [J/°K], Tj is heat flux, T is temperature, and oSint is the entropy produced internally during the irreversible process. The flux and gradient terms can then be defined through equality with a dissipation function, /cr[N/s.m ], as shown in Equation B-28. The entropy production A : is always positive according to classic thermodynamic principles. 127 where, St 5V (B-27) (B-28) In the study of electrokinetic phenomena, two coupled flow equations emerge from Equation B-26 in the absence of significant thermal or chemical concentration gradients: r h = - L h . V O h - L he .V O e Cowp/erf hydraulic flow (B-29) r e = - L e / l .V O - L e . V O e Coupled electrical flow. (B-30) In order to satisfy Equation B-28, the primary hydraulic and electrical flux and gradient terms in Equations B-29 and B-30 take the form: kg m Tf, - pw q, mass flux with units . . \\_m s.m cXD, - 3, - SO, - ,. VO / j=( i ,—-j ,—-k) , hydraulic potential gradient with units dx dy dz \"1 d~ ' N~ m kg Ml T e = pe v = J, charge flux or volume charge density with units C m3 m 3 s.m2 dV - 5V - dY -and VO =VV = (— i, — j, — k), electrical potential gradient with units dx dy dz \"1 J~ \"V\" _m C The hydraulic potential O/, [J/kg] describes the total mechanical energy per unit mass of fluid through the sum of potential, kinetic and elastic energy terms, as shown in Equation B-31, where g is acceleration due to gravity [m/s2], z is elevation [m], v is fluid velocity 128 [m/s], p is fluid pressure [Pa] and p0 represents atmospheric pressure. Since pore fluid velocities are small in porous media, the energy lost during fluid flow through the kinetic energy term is negligible, and Equation B-31 is reduced to Equation B-32, which describes the hydraulic potential of a compressible fluid, assuming pw is only a function of fluid pressure. If the fluid is assumed incompressible, hydraulic potential can be described using Equation B-3 3. Bernoulli's equation (B-31) Compressible pore fluid (B-32) Incompressible pore fluid (B-3 3) The primary conductivity terms, L/, [kg.s/m ] and L e [S/m] are derived from substitution of the primary flux and gradient terms into the macroscopic flow equations, as shown in Equations B-34 and B-35, where k is permeability [m ] and p: is dynamic fluid viscosity [Pa.s, or Poise]. k pi T h = ^ . V O ^ Primary hydraulic flow (B-34) r e = - t T . V O e Primary electrical flow (B-35) In order to derive the coupled flow conductivity terms using the same procedure applied to the primary flow terms, electrical charge must be quantified in terms of ionic mass. This is accomplished using Faraday's constant, F = 96,485 C/eq, and the atomic weight fraction of charge-carrying ions, M [kg/eq], such that: r L he = - ^ ^ , electro-osmotic cross-coupling conductivity with units . v a> p , p> p P-Po P w As2 129 and L eh = - , streaming current cross-coupling conductivity with units which satisfies Onsager's reciprocal relation, indicating that the flux and gradient terms are formulated properly according to the principles of non-equilibrium thermodynamics. Equations B-29 and B-30 can be simplified by defining hydraulic potential in terms of hydraulic head. With the assumption of an incompressible pore fluid, which is inherent in an assumption of steady-state flow conditions, hydraulic potential can be defined using Equation B-3 3, which simplifies to O a = gh . Therefore, the coupled flow equations can be re-stated as: pwq = - .V (gh) - L he .V V Coupled hydraulic flow (B-36) M p e v = - L eh. V (gh) - cr. V V Coupled electrical flow, (B-37) which reduce to Equations B-38 and B-39 when stated in terms of volumetric fluid flux and volume current density or volumetric charge flux, where ke = [m/s per V/m, or m2/s.V] is the coefficient of electro-osmotic permeability and L = Le/,.g [A/m ] is the streaming current cross-coupling conductivity coefficient. q = - K . V h - k e . V V Coupled hydraulic flow (B-38) J = pgy = - L . V h - cr.VV Coupled electrical flow (B-39) By defining the flow equations in terms of volumetric flux, the conductivity coefficients are macroscopic properties of the total average volume of a given porous medium. This is in contrast to previous studies of coupled contaminant transport, which have considered the migration of dissolved chemical constituents in terms of a diffusive flux, and expressed fluid flow in terms of pore fluid velocity (Mitchell 1991; Yeung and Mitchell 1993). Statement of the coupled flow equations in terms of volumetric flux per unit pore 130 volume through the relation v = — or —, where v is pore fluid velocity, assumes that all flow n n is strictly limited to within the pore space, and requires that the conductivity coefficients be divided by porosity. The most common form of the coupled flow equations reported in the literature express the hydraulic gradient in terms of pressure P [Pa], such that VP = p„,gVh . Using this convention, the coupled flow relations take the form of Equations B - 4 0 and B - 4 1 , where the streaming current cross-coupling conductivity in terms of a pressure gradient becomes L = ^ [A.s2/kg]. Pw k q = . V P - k e . V V Coupled hydraulic flow ( B - 4 0 ) P J = pe\\ = - L . V P - cr.VV Coupled electrical flow ( B - 4 1 ) The reciprocal nature of the coupled flow equations enables the calculation of cross-coupling conductivity from electro-osmotic permeability. If the hydraulic gradient is expressed in terms of hydraulic head, as shown by Equations B - 3 8 and B - 3 9 , L can be derived from reported values of ke using the relation shown in Equation B - 4 2 . If the coupled flow equations are formulated in terms of pressure, as in Equations B - 4 0 and B - 4 1 , L is equivalent to k g , as shown in Equation B - 4 3 . L v h = y w k . ( B - 4 2 ) L V P = k e ( B - 4 3 ) 131 APPENDIX C - Compilation Of Experimental Streaming Potential Data Published data pertaining to the study of streaming potential is presented in tabular form for a variety of geologic media. In cases where the bulk electrical conductivity of the laboratory samples was not reported, a range of values for the cross-coupling coefficient is calculated based on theoretical estimates, with an upper bound imposed by the electrical conductivity of the fluid. Footnotes to Tabulated Data * Data translated from Friborg (1996) 1. All values recorded at room temperature unless otherwise noted. 2. Bulk conductivity estimated using Archie's law in cohesionless soils, and the Waxman-Smits relation in silty clay samples. Unit Conversions Streaming Potential Coefficient, C To convert V/m to: Multiply by: raV/m 103 V/Pa Y,,,\"1 « (9810)\"' « 1.0194xl0 - 4 mV/kPa (Y, ; ' )x l0 6 ~ 101.94 mV/atm 1.0326xl04 Electrical Conductivity, a To convert S/m (l/Ohm.ml to: Multiply by: pS/cm 104 Cross-Coupling Coefficient, L To convert A/m to: Multiply by: A.s2/kg or A.m/N y,/ 1 « (9810)\"' * 1.0194x10 132 a, (uS/cm) 8 8 £ *— *— T T 8 § ! 1 8 8 t T - TJ- 8 8 8 8 i l l y § 8 8 8 Fluid conductivity a, (S/m) measured calculated 1.26E-02 1.26 E-02 1.26E-02 1.26E-02 1.26E-02 1.26E-02 1.3E-01 1.3E-01 1.3E-01 1.3E-01 1.3E-01 1 3E-02 1.3E-02 1.3E-02 1.3E-02 Fluid conductivity a, (S/m) measured calculated 3.70E-04 8.00E-04 1.47E-03 2.35E-03 4.17E-02 4.17E-02 1.5E-03 4.17E-02 6.25E-03 9.09E-02 9.09E-02 9.09E-02 4.17E-02 1 .OE-02 1 OE-02 1 .OE-02 1 .OE-02 --C (V/m) 2.29E-01 2.00E-01 1.45E-01 7.4E-02 5.11 E-02 4.59E-02 4.17E-03 5.12E-03 8.1 OE-02 6.67E-02 1.72E-01 5.71 E-03 1.23E-02 9.80E-02 2.10E-01 9.81E-01 5.40E-02 6.41 E-02 6.30E-03 6.5E-03 4.4E-03 7.0E-03 9.0E-03 5.5E-04 5.7E-04 5.9E-04 6.2E-04 6.5E-04 3.49E-03 4.69E-03 5.44E-03 6.00E-03 In 5.50E-03 4.50E-04 2.80E-04 9.70E-04 7.80E-06 5.30E-06 2.31 E-09 1.16E-08 4.40E-10 1.81 E-09 3.71 E-09 6.55E-09 3.70E-10 Ahmad (1964) Ahmad (1964) Ahmad (1964) Sasitharan et al (2001) Ogilvyetal (1969) Ogilvyetal (1969) Ahmad (1964) Ahmad (1964) Ogilvyetal (1969) Ogilvyetal (1969) Sasitharan et al (2001) Ahmad (1964) Friborg (1996) Schuch and Wanke (1967)* Schuch and Wanke (1967)* Ogilvyetal (1969) Ogilvyetal (1969) Ahmad (1964) Friborg (1996) Friborg (1996) Friborg (1996) Gray and Mitchell (1967) Gray and Mitchell (1967) Gray and Mitchell (1967) Gray and Mitchell (1967) Gray and Mitchell (1967) Gray and Mitchell (1967) Gray and Mitchell (1967) Gray and Mitchell (1967) Gray and Mitchell (1967) Casagrande (1982) Material 1. Coarse sand 2. Medium sand 3. Fine-medium sand clayey 4. Fine sand 5. Glacial till poorly-graded sandy till well-graded till well-graded silty till 6. Silty clay soft, sensitive silty clay 133 rn \"O \"O \"O . ^ 1 3 \"° trj D) ro ro O) O) a l l CD II X C L co co co o) CD a> E E E i— LO LO J , i n i n 0 co co ^ d ° S II II I- £= c E O O o o CO CO II X a . p\" f-co CO d d T f CO d CO CO d d O O c ro ca eg 5 5 | r> CTJ o o O II T f T f CO CO CO CO d d II II c c O O b d O o ° o ° ° o 2 2 3 £ t II II II E E E E CM 3 CO * II ° o II - c n <\"> csi <*? M CM N o o o .co \"I 2 I « I »_ 'to & > (0 JO . y o o o o o o o o o ^ : ronirorororororococ: Lfj >i :N .N .N M n n n Ip - d b d d d d d d d „ I sS # -S sP _s cS sf?\" # se 9 o o o ^ o o c D o a > T - < u C M n r f i n S C M O ' f L O ' O JS? II n II II II n II II II II Lo g g g g g g g g g _? IS T f T f o o T f T f o o LU LU LO CO LU LU LU LU LU LU T - i ^ T f co r~ N (D LO T - (N O O O O O CO T f \"3- T f cp p Cp o LU LU LU LU O ( D t T— CO CM CM T f T f T f o cp o 1 LU 00 LU LU CD CD 00 CM LO LO LO O O O L O T f T f o o o LO LO o o LU LU LU LU LU < - * • * I D in i n W T - T -L U L U T f LO Tf' LO CM CM CO CO LO CO CO 9 9 9 [ii ii) LU o co cn CM CM T f T f T f LD 9 9 9 LU iii LU O) CO f -^ CM CD LU LU LU cn o I S -CM CO CD L O L O L O L O L O L O L O L O p c p c p c p c p c p c p c p L U L U L U L U L U L U L U L U ^ r c n L O C M C M c n c D C N T— T—• CN CO* • T—• CM* CO E o CO CM T f CD c o c o ° 2 CM CM CO CO CM CM N t 03 N o Lo co o cn CO CM CO T f T f c o j - c o T f o CO T f T f LO co 12. CO CO CM CM o o o o CO CO T f o o o LU LU LU LU LU O (N N S S O CN CM O O C n c N C O C D C O C O T - r -LU LU LU LU CM CN co CM CM CO O CO CO CO LU LU LU O CM CM CO CO CO CO CO CM 9 9 9 LU LU LU CO CO N -LO LD O CO CO ^ C N C N C N C M C M C O C O C O C O e p p o o o o o o o LUUJLULULULULJJLULLI N- r~ r- r-0 ) 0 1 0 ) 0 ) T f T f ^ T f p p o o o CO CO CO CO CO CO CO CO CO o o o o LU LU LU LU 0) T f LO LO T - O) CO 00 ' CM '— *— CM CM CM CM p o o o 111 L U L ' r - T f c LO CO C CM CO ' LU LU LU LU CO r - CO T f CO T f Tj- L O 1 3 4 | § § | § 8 8 8 8 8 8 8 8 8 8 8 8 8 CN CN T - s $ 18900 18900 18900 18900 18900 18900 Fluid conductivity o, (S/m) I 9 9 9 9 9 9 9 9 9 9 9 1.3E-01 1.3E-01 1.3E-01 1.3E-01 1.3E-02 1.3E-02 1.3E-02 1.5E-02 1.5E-02 Fluid conductivity o, (S/m) measured 1.79E-04 1.79E-04 4.17E-02 00+368' I 00+368' I 00+368' I 00+368' l 00+368' I 00+368' I -C (V/m) 9.78E-04 1.16E-03 1.22E-03 1.29E-03 1.43E-03 2.11 E-03 2.42E-03 3.05E-03 3.58E-03 3.78E-03 1.65E-04 2.52E-04 2.73E-04 2.96E-04 1.8E-04 4.9E-04 6.0E-04 3.94E-01 1.10E+00 3.94E-03 2.35E-03 3.82E-02 5.38E-05 7.28E-05 7.22E-05 9.29E-05 8.28E-05 7.40E-05 K (m/s) 3.90E-10 8.30E-10 1.52E-09 4.85E-09 9.70E-09 2.60E-10 3.40E-10 7.70E-10 2.61 E-09 6.57E-09 3.40E-10 1.49E-09 1.27E-09 6.12E-09 2.17E-05 9.6E-06 5.7E-08 6.6E-06 3.8E-07 1.4E-08 Source Gray and Mitchell (1967) Gray and Mitchell (1967) Gray and Mitchell (1967) Gray and Mitchell (1967) Gray and Mitchell (1967) Gray and Mitchell (1967) Gray and Mitchell (1967) Gray and Mitchell (1967) Gray and Mitchell (1967) Gray and Mitchell (1967) Olsen (1969) Gray and Mitchell (1967) Gray and Mitchell (1967) Gray and Mitchell (1967) Gray and Mitchell (1967) Gray and Mitchell (1967) Gray and Mitchell (1967) Gray and Mitchell (1967) I I < Ahmad (1964) Ahmad (1964) Ishido and Mizutani (1981) Pengra and Wong (1999) Pengra and Wong (1999) Pengra and Wong (1999) Pengra and Wong (1999) Pengra and Wong (1999) Pengra and Wong (1999) Material I 8. Illitjcclay 10. Glass beads 6mm (coarse gravel) 3mm (medium gravel) 11. Quartz 12. Sandstone Fontainebleau Berea Bandera 135 TJ __ ro 0) B (0 T J o E E E E E to fo 0) _) ai O O O o o o o ra ro ro ro Z 2 2 2 _>•_•_> o o o o o o o o O o ro ra 2 Z _» _ i b b ro cP fN ro oi ro T i -ll II 5 g 5P 5P Sp T f T f T l -d t j co ro T i -ll II ll S g g _2 SP cn TJ-T f LO LO CO II II g g 2 E =P LO £f cr Sp O CM s _> s s , rN :N ro ro o b b b b SP sP ;=p £ CD O CO o oo oo r-- T f f- cn co co II n II II g g g g LO -g ll co\" -2 ° 3 CM CO II X SP £ CO 00 SP SP cP O) LO CD •\"\"I CM O CO CM CM CM II II II II o o o o o o ro ro ro ro ro ro Z Z Z Z Z — 2 -5 S 2 \"S CM CM CM CM CM CM b b b b b b E l < T f O O LU LU LO r-CO LO LO o LO O so- CD O T f o 8E-.OE-LU CM OE-5E-T\" LO T f T-CD CD CD LO CD CD O O O O O p LU LU UJ LU UJ O) O) O LO T f CO IO T - N LO LO LO LO p O O p LU LU LU LU CD LO CD *-O T - CM CD i/i id s CO LO LO LO LO LO LO LO O p O O O O O LU LU LU LU LU LU LU 1^ CO LO CD CO CO O LO LO T f CO LO T f CD CM CO T f LO LO CD E co CD T f CO CM CO O CD O T f LO LO CD O O O O r-, O O 1 0 9 M T - S O O LO co co T - ^ co co co ro g T f CO 2 CM T f (Q LO co S cn LO ~ cn jo co cn E b 1 c o LO LO CM p o p LU LU LU CO CD CD LO LO O T f T f ^ LU LU CO CO CO CO CO CO CN CM CM CN O O O O c o CO T f o o I I LU LU CO r -CO CO CM O T— CN CM T— CM CM O p O O O O LU UJ LU LU LU 111 CO T f CO O CM T f CD LO CO CO cn LO 136 o p p g o o o o o o o o cn cn co co oo co I r— '1 o 8 U J LLI O O CM' csi CO CO CO CO CO CO N I— O O CO CO LLI LU CO CO > O LU LU LU LU T - CO CM CO CO LD r - CN T - I O ) i d ^ ^ CN CO CO 00 CO N- CO CN LU UJ CO O) T -CM o CO U J LU CM O co* i r i LO cn cn co co co o> co O) O ) , O l O ) ^ T - *r •-' — cn D) CO O) C C T -O O § 5 « T3 T3 0) C C „ 00 CO C CU QL 0. o VI '1 cu 'E ro 5 TT VI cu 137 Material 12. Sandstone 13. Limestone Whitestone Indiana 14. Granite All soils (general estimates) Additional information 1 o ro 0.2 M NaCI; n = 29% 0.2 M NaCI; n = 15% NaCI crushed Westerly granite (grain size ~ uniform fine-med sand); pH=5.5, NaCI k e = 5e-5 cm^/s.v average range of C values L (A/m2 ) estimated/calculated 2.5E-07 to 9.8E-07 1.7E-04 to 6.5E-04 9.0E-06 to 3.5E-05 2.2E-05 to 8.4E-05 5.4E-05 2.9E-06 2.4E-06 4.0E-05 5.8E-05 5.8E-05 6.4E-05 4.9E-05 measured o- (uS/cm) ro co in in ^ 1100 508 4547 4547 8 8 itivity a (S/m) estimated 2 2.55E-04 5.11E-04 5.11E-04 4.55E-01 1.10E-01 5.08E-02 4.55E-01 4.55E-01 J.01 b-02 3.01 E-02 Bulk conduc measured 138 APPENDIX D - Notation A Debye length c zeta potential [V or mV] °f fluid electrical conductivity [S/m or pS/cm] c ionic concentration [eq/m3 or meq/L; where #equivalents = #mol xvalence] m ionic mobility [m2/s.V or m/s per V/m] F Faraday's constant (96,489 C/eq) fh temperature coefficient T temperature [°C or °K] a average pore diameter d representative or effective grain size q 3 2 volumetric fluid flux vector [m /s.m ] q volumetric fluid flux scalar of vector component j volume current density vector [A/m ] j volume current density scalar or vector component ss specific storage [m /m .m] h hydraulic head [m] t time [s] Q 3 3 volumetric fluid flow [m /s.m ] Pe volumetric charge density [C/m ] I volumetric current flow [A/m ] K hydraulic conductivity vector [m/s or cm/s] K hydraulic conductivity scalar or vector component cr bulk electrical conductivity vector [S/m or u.S/cm] cr bulk electrical conductivity scalar or vector component V electrical potential [V or mV] r, flux vector potential cross-coupling coefficient vector 139 hydraulic potential [N.m/kg] Pw fluid density [kg/m ] U pore pressure with respect to atmosphere (gauge pressure) [Pa] 8 acceleration due to gravity [m/s2] xyz positional coordinates [m] S entropy [J/°K]; volumetric source; degree of saturation L streaming current cross-coupling coefficient vector [A/m ] L streaming current cross-coupling coefficient scalar or vector component ke electro-osmotic permeability coefficient vector [m /s.V] ke electro-osmotic permeability coefficient scalar or vector component Jconv convection (streaming) current vector Jcond conduction current vector teq relaxation time potential M> total or pseudo potential D divergence matrix M conductivity matrix G gradient matrix [u] potential vector S source matrix A coefficient matrix c conductance [m2/s]; coefficient of proportionality K reflection coefficient P dynamic viscosity [Pa.s = N.s/m2 = kg/m.s = 1 Poise, P] k intrinsic permeability [m ] n porosity P bulk electrical resistivity [Ohm.m] os surface conductivity [S/m or S/cm] or surface conductance [S] B equivalent cation conductance [S.cm2/meq] Qv concentration of exchange cations [meq/mL] c streaming potential coupling coefficient 140 unit weight of water, defined as: yw = pw g [N/m3] s electric permittivity, defined as: s= ks0 [F/m] k dielectric constant [no units] permittivity of free space (8.85 x 10\" F/m) 0~sat saturated bulk electrical conductivity [S/m or pS/cm] Q~unsat unsaturated bulk electrical conductivity [S/m or pS/cm] e electronic charge [C] kb Boltzman's constant (1.380622 x 10\"23 J/°K) I ionic strength [mol/m ] n valency 141 "@en ; edm:hasType "Thesis/Dissertation"@en ; vivo:dateIssued "2002-05"@en ; edm:isShownAt "10.14288/1.0063692"@en ; dcterms:language "eng"@en ; ns0:degreeDiscipline "Civil Engineering"@en ; edm:provider "Vancouver : University of British Columbia Library"@en ; dcterms:publisher "University of British Columbia"@en ; dcterms:rights "For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use."@en ; ns0:scholarLevel "Graduate"@en ; dcterms:title "Response of the self-potential method to changing seepage conditions in embankment dams"@en ; dcterms:type "Text"@en ; ns0:identifierURI "http://hdl.handle.net/2429/12220"@en .