Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Use of DC resistivity and induced polarization methods in acid mine drainage research at the Copper Cliff… Zudman, Yuval 1995

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

Item Metadata


831-ubc_1995-0177.pdf [ 10.24MB ]
JSON: 831-1.0052956.json
JSON-LD: 831-1.0052956-ld.json
RDF/XML (Pretty): 831-1.0052956-rdf.xml
RDF/JSON: 831-1.0052956-rdf.json
Turtle: 831-1.0052956-turtle.txt
N-Triples: 831-1.0052956-rdf-ntriples.txt
Original Record: 831-1.0052956-source.json
Full Text

Full Text

USE OF DG RESISTIVITY AND INDUCED POLARIZATION METHODS IN ACID MINE DRAINAGE RESEARCH AT THE COPPER CLIFF MINE TAILINGS IMPOUNDMENTS By Yuval Zudman B. Sc. (Geophysics) Tel-Aviv University, 1991 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF T H E REQUIREMENTS FOR T H E D E G R E E OF M A S T E R OF SCIENCE in T H E FACULTY OF GRADUATE STUDIES DEPARTMENT OF GEOPHYSICS AND ASTRONOMY We accept this thesis as conforming to the required standard T H E UNIVERSITY OF BRITISH COLUMBIA January 1995 © Yuval Zudman, 1995 In p resen t ing this thesis in partial f u l f i lmen t o f t h e r e q u i r e m e n t s f o r an advanced d e g r e e at the Univers i ty o f Brit ish C o l u m b i a , I agree that t h e Library shall make it f reely available f o r re ference and s tudy. I fu r ther agree that pe rmiss ion f o r ex tens ive c o p y i n g of this thesis f o r scholar ly pu rposes may be g ran ted by t h e head o f m y d e p a r t m e n t or by his o r her representat ives. It is u n d e r s t o o d that c o p y i n g o r pub l i ca t i on o f th is thesis f o r f inancial gain shall n o t be a l l o w e d w i t h o u t m y w r i t t e n permiss ion . D e p a r t m e n t o f ^r?0^f\C. <jjC^ OfY) &l $ fstflQ If) 0 M of Brit ish C o l u m b i a The Univers i ty Vancouver , Canada Date fiPRZL DE-6 (2/88) Abstract Oxidation of sulphide minerals contained in the mine tailings impoundments at Copper Cliff, Ontario generates acidic conditions and elevated concentrations of dissolved con-stituents in the pore water. The tailings groundwater migrates radially outwards and might pose an environmental hazard if reaches the nearby water systems. There is a need to estimate the dimension of the current problem and assess future prospects. That need prompted a combined DC resistivity and Induced Polarization (IP) survey along one of the major flowpaths in the tailings. The DC resistivity and IP data were inverted to produce the detailed electric conduc-tivity and chargeability structures of the cross section below the survey Une. Problems, generic of the inversion methodology and specific to the inversion of data from the tail-ings, are discussed and solutions are elaborated and demonstrated. The conductivity structure can be directly translated, through theoretical or empirical relations, to a map of the concentration of dissolved solids along the cross section and thereby provide in-sight about the current pore water quality. The sulphide minerals are the source of the IP response and thus the chargeability model can be used to estimate the amount and distribution of the sulphides. A new methodology to construct and invert IP decay curve information was devel-oped. Assuming Cole-Cole relaxation to model the IP phenomenon, data sets constructed through the inversion of IP data of different time windows were inverted using a damped least squares algorithm. Values of the resultant Cole-Cole parameters might contribute additional insight about the sulphides in the tailings. Lack of time and laboratory veri-fications prevented full development and comprehensive conclusions. i i The research in this thesis concentrated along one line of the survey data with the intention of showing the capability of the DC resistivity and IP methods and their po-tential contribution to estimation of location and concentrations of the sulphide minerals and their oxidation products. The conclusions should assist designers of a large scale sur-vey and prompt investigations of the relations between the physical properties mapped by these methods and the chemical parameters which are considered while assessing the damage and contemplating solutions. in Table of Contents Abstract ii List of Tables vii List of Figures viii 1 Introduction 1 2 Overview Of The Acid Mine Drainage Problem 4 2.1 Introduction 4 2.2 Acid Generation And Neutralization 5 2.3 Abatement And Control Measures 6 3 The Copper Cliff Tailings Disposal Area 8 3.1 Site Description 8 3.1.1 Tailings Composition And Mineralogy 9 3.1.2 Ground Water Hydrology 12 3.2 AMD Problem In The Tailings 12 3.3 Physical Properties 14 3.3.1 Physical Properties Of Rocks 14 3.3.2 Physical Properties Of The Disposal Area 16 3.4 Motivation For The Use Of The DC Resistivity And IP Methods In The AMD Research , , . . , 18 3.5 Other Geophysical Methods 19 iv 4 T h e D C Res i s t i v i ty A n d I n d u c e d P o l a r i z a t i o n M e t h o d s 22 4.1 Theoretical Background 22 4.1.1 DC Resistivity Method 22 4.1.2 Induced Polarization Method 24 4.2 DC/IP Measurements 26 4.3 Field Procedures 28 4.4 Traditional Data Plotting And Interpretation 29 4.5 The Geophysical Survey 32 4.5.1 Survey Lines 33 4.5.2 Field Procedures, Considerations And Available Data 34 4.6 Data Plots And Initial Interpretation 35 4.6.1 Data Plots 35 4.6.2 Initial Interpretation 36 5 Invers ion A n d Interpre ta t ion O f D C Res i s t iv i ty A n d I P D a t a 40 5.1 Inversion Methodology 40 5.1.1 Subspace Non-Linear Inversion 40 5.1.2 DC Resistivity Inversion 45 5.1.3 Induced Polarization Inversion 45 5.2 The DCIP2D Package 47 5.3 Typical Problems In The Inversion Of The Tailings Data Sets 48 5.3.1 High Conductivity Of The Upper Layer 48 5.3.2 Buried Metal Bodies 49 5.4 Assigning Errors To The Data 50 5.5 Inversion Results And Interpretation 56 5.5.1 Parameters For The Inversion 56 v 5.5.2 Inversions Results 59 5.5.3 Possible Improvements Of Results Through The Use Of Weights . 61 5.5.4 Interpretation 62 6 Inversion Of IP Decay Curve Data 75 6.1 Introduction 75 6.2 Modeling The IP Waveform With A Relaxation Model . 77 6.2.1 Equivalent Circuits . 77 6.2.2 Time Domain Behaviour Of The Cole-Cole Relaxation Model . . 79 6.2.3 Intrinsic Decay Time Constant And Frequency Dependence . . . . 80 6.3 Inversion 81 6.3.1 Inversion Methodology 81 6.3.2 IPINV2D Results As Input For The Inversion 81 6.3.3 Three Parameters Damped Least Squares Inversion 85 6.3.4 Inversion Results 86 6.3.5 Interpretation 87 7 Conclusions 90 7.1 DC Resistivity 90 7.2 Induced Polarization 91 Appendices 93 A Resistivity/Conductivity Conversion Table 93 References 94 vi List of Tables 5.1 Average concentrations of sulphides and maximum concentrations, (in mg/L), of oxidation products in sites along line A A . Sources: Coggans, (1994), Cherry et al, (1990) and King, (1994). The corresponding con-ductivities (S/m) and chargeabilities (Volt/Volt) of the upper metres of the inversion results 64 6.1 Parameters corresponding to figure 6.5 88 vii List of Figures 3.1 Location map of Sudbury, Northern Ontario 9 3.2 The Copper Cliff tailings disposal area. Survey lines AA and BB are designated by thick lines. The numbers along line AA are of sampling wells. Coordinates along the Une have their origin at the location of weU IN13 with the axis direction south to north. (Negative values towards south) 10 3.3 Cross-section of Une A A 11 3.4 ModeUed ftownet of cross-section AA. (After Coggans, 1992) 13 4.1 Responses of a non-polarizable and polarizable media to square-wave cur-rent 24 4.2 Phase lag 8 between the current and voltage cycles 25 4.3 Point current source and potential at the surface. The dipoles are grounded at infinity. 27 4.4 Common arrays, (a) Schlumberger, (b) Pole-Dipole, (c) Dipole-Dipole. . 29 4.5 Plotting points of apparent conductivity/chargeabiUty for dipoles n=l,3. 30 4.6 (a) Conductivity model — 0.07 S/m conductive layer overlying 0.005 S/m half-space and 10.0 S/m small anomaly, (b) ChargeabiUty model—0.01 Volt/Volt half-space, 0.05 Volt/Volt large anomaly and 0.5 Volt/Volt small anomaly, (c) Conductivity pseudosection, (d) ChargeabiUty pseudosection. 31 4.7 Conductivity pseudosection 38 viii 4.8 Chargeability pseudosection. Note the difference in the colour scaling of the two parts of the plot 39 5.1 Apparent chargeability data in V / V between locations 60 and -490 given separately, in a Une format, for each potential dipole. Plotting point is the midpoint of the potential electrodes at each position 53 5.2 (a) Assigned error as a function of datum value. Datum value of 1 is taken as the mean of the data set. The parameters used are A = 0.05, B = 0.1, C = 0.1, D = 0.1 and F equals 3.0 for positive data and 2.0 for negative ones, (b) The error relative to the datum as a function of the datum value. 55 5.3 Plots, (a) and (b) are maps of {df}' — d!jed^ jej plotted the same way dj" is plotted on a pseudosection. Plots (A) and (B) are the corresponding inversion results. Model A was produced using the error assigning scheme of equation 5.21. To produce model B the data were assigned errors which are 15 % of the average of the whole data set plus a contribution from the same exponential term of the more compUcated scheme 68 5.4 Conductivity inversion results of time domain data of Une A A. The colour coded values are logarithms to the base ten of the conductivity in S/m. . 69 5.5 ChargeabiUty inversion results of time domain data of Une A A . The colour coded values are chargeabiUty in Volt/Volt 70 5.6 Conductivity and ChargeabiUty inversion results of the phase data of Unes A A , AA-East and A A-West. Colour coded logarithm to the base ten of the conductivity in S/m. Phase change chargeabiUty values in miUradians. 71 ix 5.7 (a) Chargeability inversion results with misfits, (from top to bottom), 8N, ION and U N , where N is the number of data. Requesting lower misfit pushed the image of the pipe, along with other features, from the surface down, (b) Conductivity inversion results with misfits, (from top to bottom), IN, 2N and 4N. The chosen model has misfit 3N and is shown in figure 5;6 72 5.8 The conductivity borehole log at IN10 and the corresponding conductivity obtained from the inversion 73 5.9 Model and derivatives as function of depth. Model-solid line, first derivative-dashed line, second derivative-solid line and triangles, (a) Line A A model parameters at location -510. (b) Inversion of synthetic data forward modeled on a model of 20 metres conductive layer overlaying a resistive half-space. 74 5.10 The model corresponding to that of figure 5.4b but produced with tailored weights for the inversion 74 6.1 Curve of the vector sum of the in-phase and out-of-phase components of the voltage in the complex plane. The CR method provides many points along the curve—corresponding to many frequencies. The IP method uses only two separate frequencies. (After Zonge and Wynn, 1975) 76 6.2 Schematic section of mineralized rock, its equivalent circuit and the fre-quency and time domains IP responses, (After Pelton et al, 1978) 78 6.3 Schematic explanation of the the process of extracting intrinsic chargeabil-ity decay curves through IPINV2D inversions of IP time domain data. For clarity of the diagram only five time windows are considered in the figure. 83 6.4 Models (a), (b) and (c) are inversion results of the third, fifth and seventh time windows data correspondingly . 84 x Observed, (circles), and predicted, (solid line and crosses), chargeabili-ties as function of time. The corresponding model parameters and cell locations of each plot are given in table 6.1 xi C h a p t e r 1 I n t r o d u c t i o n Mine tailings have been continuously deposited since 1937 in the Copper Cliff tailings disposal area near Sudbury, Ontario. The tailings impoundments currently cover an area of 21 km2 to depth of up to 40 metres. Oxidation of sulphide minerals contained in the tailings occurs within the unsaturated zone and results with low pH conditions and high levels of dissolved solids in the pore water. Contaminated groundwater flow from the elevated tailings area recharges the local water systems and if not impeded and treated may pose an environmental hazard. The type and extent of the possible environmental impact of the contamination ori-ginating from the tailings has to be assessed. The first step in a comprehensive study is the estimation of the levels of the already existing oxidation products and of the total amount and distribution of the source of future oxidation—the sulphide minerals. Sparse in-situ measurements can provide only localized information and extensive drilling is not cost effective. Answers should be sought using surficial methods. The electric conductivity of pore water, (and of the porous medium it fills), is a phys-ical property very sensitive to the level of dissolved solids in the water. The oxidation products contribute significantly to that level and consequently to the conductivity of the medium. Chargeability is another electrical property of the earth. The physical mechanisms of the chargeability phenomenon are not fully understood but one of its known sources is sulphide minerals disseminated in the rock. DC resistivity is one of the geophysical methods capable of mapping conductivity of the earth. The Induced 1 Chapter 1. Introduction 2 Polarization is a method for detecting its chargeability. Data of both methods can be collected with the same instrumentation in a non-intrusive way. A combined D C / I P survey can provide information about earth conductivity and chargeability with relat-ively low costs and in a short time. Existing computerized inversion algorithms enable quantitative mapping of the conductivity and the chargeability using that information. The following work shows application of the DC resistivity and IP methods to the quantitative estimates of the conductivity and chargeability in the Copper Cliff tailings disposal site. The thesis begins with an overview of the general acid mine drainage prob-lem followed by a description of the Copper Cliff tailings disposal area and its specific acid mine drainage problem. It includes a discussion of the relevant physical properties to be considered and the geophysical methods for detecting them. The theoretical background to the methods used in the thesis—the DC resistivity and Induced Polarization—is given next along with the details of the survey and presentation of the data. Inversion and interpretation of the DC resistivity and IP data constitute the main part of the work. A package of computerized inversion algorithms was used to obtain the conductivity and chargeabiUty structures of the cross section below one survey Une. The inversion methodology and a short description of the computer programs package are given. Practical considerations in carrying out the inversions are discussed and solutions to the typical problems are proposed. This part concludes with demonstration of the inversion results and their interpretation. The pattern of the decay of the chargeabiUty with time contains additional inform-ation about the earth. In the last part of the thesis I propose a new methodology to study the IP decay as another tool of investigating the sulphide minerals in the taiUngs. Through the inversion of apparent chargeabiUty data of different time windows, IP decay curves in cells of the subsurface are constructed. The decay, at each cell, is assumed to be a Cole-Cole relaxation {Cole and Cole, 1941). Data points along the curves are Chapter 1. Introduction 3 inverted, using a damped least squares algorithm, to intrinsic Cole-Cole parameters of the relaxation. Time constraints and lack of laboratory data to verify the results limited that research and only initial results and interpretation are given. Chapter 2 Overview Of The Acid Mine Drainage Problem 2.1 Introduction Exposed to water and oxygen, sulphide minerals oxidize to soluble metals and sulphates. Products of the oxidation contribute to the acidity in the oxidizing site. In the absence of alkalinity to neutralize the acidity, pH level can drop quickly. Furthermore, low pH conditions encourage the activity of the acidophilic bacteria Thiobacillus ferrooxidans which catalyze the oxidation. The elevated solubility of metals in acidic water prevents their precipitation and results in high concentrations of dissolved metal and salts in the acid water. Water flushing the acid generation site causes an acid drainage into groundwater systems and downstream water bodies. The contaminated water has a negative impact on the environment. The damage is to the the quality of potable water reservoirs and the health of wildlife and plants as well as aesthetic damage. The leaching of sulphide oxidation products can happen naturally in the vicinity of sulphide-bearing rocks or in construction work sites where active sulphides have been exposed. The term Acid Rock Drainage (ARD) is generally used for all occurrences of acidic drainage due to sulphide oxidation. However, acid generation is usually associated with mining. Most of the metal deposits contain sulphide minerals in the ore body or the surrounding waste rock. The mining activity exposes large amounts of the sulphide minerals to water and oxygen. If no measures are taken to neutralize the acid generation, 4 Chapter 2. Overview Of The Acid Mine Drainage Problem 5 sooner or later, acid drainage occurs and is referred to as Acid Mine Drainage (AMD). A M D is currently experienced in many abandoned and operating mines around the world. With the exception of Alberta, Nova Scotia and Prince Eduard Island, all Ca-nadian provinces and territories have been identified as having acid generating mining sites. The environmental damage and financial liabilities from these sites and the con-cerns about future acid generation from mines yet to be developed are regarded as a major problem facing the Canadian mining industry today. 2.2 Acid Generation And Neutralization Sulphide minerals are compositions of sulphur with metals or semi-metals but no oxygen. Primarily found in rocks which lie beneath the upper soil and the water table the contact of these minerals with the atmospheric oxygen is minimal. When sulphides are exposed to air and humidity, chemical reactions take place which result with acid generation. The oxidation reactions of the mineral pyrrhotite (FeS) illustrate the acid generation process. The pyrrhotite is oxidized into dissolved iron and sulphate by FeS + 202 -> F e 2 + + SO\~ . (2.1) The iron ion F e 2 + may oxidize to F e 3 + and further react with pyrrhotite through 8 F e 2 + + 20 2 + SH+ -* 8Fe3+ + \H20 (2.2) FeS + 8 F e 3 + + AH20 -+ 9Fe2+ + SO\~ + 8H+ (2.3) FeS + F e 3 + + 7O2 + ^H20 -» 2Fe00H + S + 3H+ (2.4) 4 2 to produce one mole of H+ for every mole of consumed pyrrhotite. Reactions of other sulphide minerals with the same or other oxidants occurs in different rates but all result with the production of H+, S0\~ and dissolved metal ions. Chapter 2. Overview Of The Acid Mine Drainage Problem 6 The rate of the reactions is determined primarily by the oxygen concentration in the gas and aqueous phases. Other important factors are the pH level, temperature, surface area of the exposed sulphide grains and the activity of the bacteria Thiobacillus ferrooxidans which accelerates the ferrous-iron oxidation. Acid consuming minerals can be contained in the acid generating rock or encountered by the acid leachate while transported from the generating site. A common mineral like calcite (CaCOz) can neutralize the acidity through the reactions CaC03 4- H+ Ca2+ + ECOl (2.5) and CaC03 + 2H+ -* Ca2+ + H2C03 . (2.6) Other minerals with high acid consumption potential are siderite (FeCOs) and Magnesite (MgCOz)- The balance between the acid generation reactions and the neutralizing re-actions determines the the net acidity level of the water entering the environment. 2.3 Abatement And Control Measures The current approach to the abatement of the AMD problem defines three levels of control: 1. Prevention or minimization of acid generation. 2. Control of the AMD migration from the generation site. 3. Collection and treatment of contaminated drainage. Success of the primary control measure is the most effective and desireable. It prevents the need to use the other measures and assures negligible hazards to the environment. The prevention of acid generation can be achieved by exclusion of at least one of the essential ingredients: sulphide minerals, water or oxygen. Chapter 2. Overview Of The Acid Mine Drainage Problem 7 Water cover to reduce oxygen availability is considered the best solution where it can be implemented. An alternative way to minimize oxygen availability is with the use of soil covers. These must contain high moisture level to minimize oxygen diffusion. Till or clay covers can reduce access of both water and oxygen although not to the extent that acid generation would not occur at all. Other methods which were suggested involve excluding the sulphides or minimizing their effect. That can be done by conditioning of tailings or waste rock to remove sulphide minerals or controlling the pH level with base additives. The secondary control level intends to minimize the impact on the environment by controlling the migration of the A M D from the generation site. As the contaminants are entrained by water flow, prevention of water entry to the site is the only possible option. (Prevention of water outflow cannot be effective in the long term). Impermeable covers to prevent water infiltration, or diversion of the surface water are among the possible methods. These can be implemented only in relatively small or pre-planned disposal sites and success can never be total. Collection of the A M D and treatment to control water pH levels are the last resort if the first two control levels fail or do not achieve adequate results. An effective collection system must intercept the acid leachate which is then treated chemically to neutralize the acidity. The process is widely applied, sometimes in conjunction with other measures. The main problem is collecting all the A M D . Surface water from seeps at the base of a waste rock pile or tailings dam can be easily collected but there is always a ground flow which bypass the collection station and is very difficult to intercept. Another problem is the cost of the treatment process. The acidity is usually neutralized by chemical treatment—mainly by adding limestone or lime. With continuous acid generation, which might last for hundreds of years, chemical treatment means continuing costly operations in the mine for long time after ceasing production. Chapter 3 The Copper Cliff Tailings Disposal Area 3.1 Site Description The International Nickel Company (INCO Ltd.) operates a large complex of mining, milling, smelting and refining in the Sudbury district, Northern Ontario (figure 3.1). Ore extracted from the sublayer of the Sudbury Igneous Complex, containing up to 60 wt. % sulphides, is processed to produce nickel, copper, cobalt and precious metals. Since the 1930s the recovery is based on the metal-concentration process which results in the production of mill tailings. Ores from 10 mines within the area are upgraded at the Copper Cliff mills. The tailings—ground gangue minerals and ore fragments uneconomical to process further— are disposed off, since 1937, in near by shallow bedrock valleys. During the years, six impoundments were created by erecting dams between the bedrock ridges and dumping the tailings directly to the bottom of the valleys. Today the tailings cover an area of more then 21 km2 to thickness of up to 40 metres, (figure 3.2) The tailings disposal area is a huge artificial body elevated 20 to 30 metres above the surrounding terrain. The surface is fairly flat and barren. Few artificial ponds cover the lowest areas. In 1957 a re-vegetation program was initiated to stabilize the surface soil and prevent dusting. Today large parts of the CD, M and M l impoundments are covered by grass and small trees. 8 Chapter 3. The Copper Cliff Tailings Disposal Area 9 Figure 3.1: Location map of Sudbury, Northern Ontario 3.1.1 Tailings Composition And Mineralogy The ores are derived from noritic to gabroic rocks containing up to 60 wt. % sulphide minerals. The most common sulphides are pyrite (FcS^), pentlandite ((Ni, Fe)9Ss), chalcopyrite (CuFeSz) and mainly pyrrhotite (FeS). The residue of the floatation, magnetic separation and regrinding processes in the mill are the fine-grain to silty sands which are dumped with the liquid wastes of the processing in the form of blackish mud. The composition of the tailings varies due to different processes which took place during the years. Tailings solids are principally feldspars with 1-7 % sulphides and about 3 % carbonates content. Analysis of tailings cores reported by Coggans (1992) show pyrrhotite to be the principal sulphide mineral (about 3 wt. %) with smaller amounts of chalcopyrite, pentlandite and pyrite. About 1 wt. % of magnetite was found throughout the tailings dumped during the operation of an iron recovery plant. Carbonates, a main source of acid neutralizing form up to 3 wt. % of the total bulk tailings, mainly as calcite. Chapter 3. The Copper Cliff Tailings Disposal Area 10 Figure 3.2; The Copper Cliff tailings disposal area. Survey lines AA and BB are des-ignated by thick lines. The numbers along line AA are of sampling wells. Coordinates along the line have their origin at the location of well IN13 with the axis direction south to north. (Negative values towards south). Chapter 3. The Copper Cliff Tailings Disposal Area 11 Figure 3.3: Cross-section of line A A The bulk of the tailings is uniform fine-grained to silty sand, unoxidized and black in colour. The average porosity is 0.50 . From the surface down to 0.5-2.25 meters the tailings are highly oxidized and brownish-yellow in colour. Within the active oxidizing zone thin discontinuous hardpan lenses of tailings were found. These inhomogeneities are probably formed by the cementation effect of precipitated minerals like gypsum, jarosite and goethite. Chapter 3. The Copper Cliff Tailings Disposal Area 12 3.1.2 Ground Water Hydrology Precipitation rate at the tailings site is 0.78 m/a, 30 % of which recharge the groundwater system. The natural precipitation and artificial pumping of water into the ponds maintain the water table at an elevation of 20-30 metres above the surrounding terrain and 0-5 meters below the surface. Over most of the tailings area the water table follows the flat topography of the tailings with hydraulic gradients increasing sharply only near the tailings dams at the southern and southeastern parts of the impoundments (figure 3.3). Decreasing hydraulic head with depth, measured in all piezometer nests (Coggans, 1992), indicates the tailings impoundments as a regional groundwater recharge area. Cherry et al. (1990) identified seven major outwards radial flowpaths originating within the tail-ings. Modelled flow system of one of them— along line A A (figure 3.4)—suggests that the water flows downwards at an average velocity of 0.47-0.55 m/a and then laterally outwards along the tailings-bedrock boundary. Only minor flux penetrates the bedrock. Tailings water seeps occur on the surface at the base of the dams due to artesian con-ditions. De Vos (1992) found that a bedrock depression beneath the Pistol Dam at the southeast part of the tailings allow tailings pore water to flow outwards to the southeast through an overlying thin layer of glaciofluvial gravel and sand layer. It is possible in other locations around the impoundments, that tailings derived water extends outwards and finally will reach the many lakes and creeks in the area. 3.2 A M D Problem In The Tailings Sulphide oxidation reactions occur within the unsaturated zone of the tailings area. Ox-idation of the upper tailings is visibly evident by the yellowish colour of the surface. Mineralogical analysis of tailings core samples (Coggans, 1992) indicates almost com-plete depletion of sulphides in the upper 10 to 25 centimetres and various degrees of Chapter 3. The Copper CUff Tailings Disposal Area 13 Figure 3.4: Modelled flownet of cross-section A A . (After Coggans, 1992). oxidation down to 2.25 metres. As the water table level is 2 to 5 metres below the surface, significant sulphide oxidation potential still exists. Pore water pH level is as low as 3 near the surface and increases with increasing depth to neutral below 4 metres depth. High concentrations of dissolved metals and sulphate can be found down to a depth of 9 metres. The difference between the acidity and the metals contaminants fronts is explained by the acidity retardation effect of a sequence of neutralizing reactions which decrease the H+ concentration. The low pH water front, defined by less then pH=4.5 water, moves at velocity of 0.13 m/a—about a quarter of the velocity of the groundwater flow velocity. Of the dissolved metal cations, iron and nickel are the most abundant with concentrations of up to 3000 and 641 mg/L in the pore water. Other metals present in lower concentrations are cobalt, zinc, copper, lead, manganese and arsenic. The high concentrations of metals are toxic and highly undesirable. Modeling of the sulphide oxidation along line A A (Coggans, 1992) predicts that where the unsaturated zone is thin (less than 2 metres) oxidation will be completed within 5 to Chapter 3. The Copper Cliff Tailings Disposal Area 14 13 years, but where it is at its maximum—about 5 metres thick—the time required for complete oxidation is nearly 400 years. That suggests that were the water table allowed to fall and expose more tailings to oxygen, sulphide oxidation and the associated hazards will continue for hundreds of years. With the water table maintained at the current level, oxidation rate and contaminants production are decreasing. But a front of contaminated water is already propagating to-wards the Copper Cliff creek. If current hydrogeological conditions remain, contaminant water is expected to reach the creek in about 40 years and continue discharging for another 340 years. De Vos (1992) using coring, piezometers and DC resistivity sounding data identified a high conductivity, (0.25 S/m), plume extending from the Pistol Dam to the southeast. That plume is probably of tailings derived water and extends to yet unknown distance from the impoundments. 3.3 Physical Properties 3.3.1 Physical Properties Of Rocks The physical properties of interest for the purpose of the research are the conductivity, chargeability and magnetic permeability. Conductivity With the exception of clays, conductivity of rock-forming materials is very low. The rock matrix does not contribute to the electrical conductivity of the rock. Presence of conductive minerals may increase the bulk conductivity of a rock in dependence on the conductivity of the minerals and their geometrical distribution but the effect is very small. Relatively high concentrations of conductive minerals or high connectivity of conductive mineral veins are needed to elevate the bulk conductivity by significant amount. Chapter 3. The Copper CliiF Tailings Disposal Area 15 The dominant factors determining the rock conductivity are its water content and the pore water conductivity. Water appears in nature as an electrolyte and thus the pore water enables electrolytic conduction through the water filled pores. The bulk conductivity a of a rock is usually given by Archie's law (Archie, 1942) a = acrwti m (3.1) where <rw is the pore water conductivity, $ the volume fraction of water in the rock and a and m are parameters characteristic of the rock. In a homogeneous medium where a, m and t? are constants, any variations of the rock conductivity cr can be attributed to changes of the pore water conductivity aw. The water conductivity is a function of the number of dissolved ions and their mobilities—the terminal velocity of the ions under a constant electric field. Mobility is a complicated issue depending on viscous friction, temperature and the pore pres-sure. Most ions though, have, under equal conditions, the same mobility. Exceptions are the hydrogen ion H+ and the hydroxyl OH~ which have mobilities six times and three times higher then the average. The bulk rock conductivity depends on the pore water conductivity and thus, ultimately, on the total concentrations of dissolved solids in the water. Chargeability The main IP effect is caused by a build up of over-potential, opposing an external current flow, at the interface between electron-conducting minerals and ion-conducting electro-lyte. The presence of conducting-mineral grains which block or are next to the electro-lytic pore water paths is essential and it is expected that the IP response will increase with the concentration of mineral grains. The issue is more complicated though, because with increasing concentration, the more resistive blocked paths are short-circuited and thereby cause the response to decrease. Chapter 3. The Copper CHff Tailings Disposal Area 16 Most rocks materials show very low or nil polarization. Concentrations of conductive minerals like most of the sulphides, graphite and magnetite are usually the source of IP response. As the polarization occurs not inside the mineral grain but at it's surface, a response is expected when the minerals are disseminated throughout the medium. This is a great advantage over other geophysical methods which require connectivity of the conductors. The shape of the mineralization body does not play a role (Seigel, 1959). What do seem to be significant are the conductive mineral concentration and its grain size. Pelton et al. (1978), analyzing results of laboratory experiments carried out by Grisseman (1971), found that increasing concentration (up to 20 % level) increases the chargeability, whereas the effect of increasing the grain size is to decrease the chargeability. Permeability As magnetic susceptibility of most materials is very small the permeability of most rocks is very close to that of the free space. Sedimentary and metamorphic rocks have the lowest susceptibilities, igneous rocks usually show higher values but still quite insignificant. The source of magnetization of rocks is ferrimagnetic minerals like pyrrhotite and magnetite and the magnetization effect depends on the amount of these minerals in the rock. The electric conductivity, chargeability and magnetic permeability in the tailings area, their sources and expected values will be discussed in the next section. 3.3.2 Physical Properties Of The Disposal Area The tailings impoundments fill a few shallow valleys of hard and compact bedrock with possible thin layer of lacustrine sediments at their bottom. The hard, low porosity bedrock is expected to be very resistive. Conductivity measured at the surroundings of a new pyrrhotite disposal site is about 5 x 10 - 5 S/m. The tailings material consists of resistive waste rock and ore fragments with 50 % porosity. The highly acidic pore water makes the tailings very conductive. Borehole and Chapter 3. The Copper Chff Tailings Disposal Area 17 in-situ penetrating cone measurements on the tailings show conductivity values between 0.05-0.10 S/m under the water table. These values are in agreement with conductivities of acidic water measured by Pheme (1981) and Ebraheem et al. (1990) in acid contaminated sites. The conductivity of the dry layer above the water table is very low—below the minimum reading of the cone—and is increasing only in the vadose zone just above the water table. Polarization of the bedrock and of the main constituents of the tailings is negligible. IP anomalies should be attributed to the IP response of the sulphides, IP-like anomalies near artificial conductors or strong E M effects. Sulphide minerals are well known sources of IP. Induced Polarization surveys proved to be very useful in the search for sulphide ore bodies, and especially of disseminated sulphide mineralization. Laboratory measurements also indicate high IP responses pro-portional to the percentage of sulphide minerals in the rocks. According to Keller and Frischnecht (1965) mineralized rock with 2-8 % sulphide content, (the range expected in the tailings), have chargeabilities of 0.05-0.10 volt/volt. Spurious IP anomalies can occur due to non-linear resistance effects in the vicinity of good conductors like buried metal bodies. This "fence effect", (Wynn and Zonge, 1975), is very likely to occur in the man made tailings where buried pipe-lines and scrap metal are very common. Capacitance coupling in the electrical measuring system or E M coupling due to eddy currents in the earth are another troublesome sources of IP noise. Use of very low frequencies in a frequency domain system or delay between turn off of exciting current and the onset of secondary voltage measurements in a time domain system are the traditional ways to avoid these problems. In the conductive environment of the tailings, where displacement currents decay slowly and especially in the vicinity of good conductors these solutions might not work and the interpreter should bear the problem in mind. Chapter 3. The Copper Cliff TaUings Disposal Area 18 The magnetic permeability of most of the rocks is too small and has no significance in electric work. Most of the tailings material and the bedrock are not an exception. The pyrrhotite though, and especially magnetite which can also be found in the tailings, have exceptionally high magnetic permeabilities (Telford et al., 1990). In areas of high concentrations of these minerals any possible E M effects will be enhanced and might interfere with IP measurements. 3.4 Motivation For The Use Of The D C Resistivity And IP Methods In The A M D Research Knowledge of the total volume and distribution of sulphide minerals in the tailings is of primary importance for quantitative prediction of the A M D production potential of the tailings. Data from boreholes or in-situ measurements, although indispensable as control measures, cannot by themselves be used to build a model of sulphide distribution in the vast disposal area. An inexpensive, fast method should be employed to detect and map the sulphides. Off all the traditional geophysical surficial methods the IP method is the only one capable of detecting disseminated sulphides. The current ability to invert IP data can facilitate quantitative estimates of the amounts and distribution throughout the tailings and may be, as further step, prediction of A M D production. That makes obvious the need to test that method and its potential use. In order to invert IP data it is necessary to obtain the conductivity structure. Inver-sion of DC resistivity data is one of the ways to do it. As the DC potentials are collected anyway as part of the IP survey, it is an easy choice to use them for that purpose. Furthermore, in the final IP interpretation both chargeability and conductivity models should be at hand as doubts and uncertainties about one could be resolved by the other. Chapter 3. The Copper Cliff Tailings Disposal Area 19 Knowledge of the conductivity structure has also merits by itself. A conductivity model can roughly map the contaminant plume and suggest areas of high oxidation. (Ebraheem et al, 1990). Moreover, the tailings are relatively homogeneous in terms of soil structure, grain size and porosity. It may therefore be possible to establish a good correlation between the medium's conductivity and the conductivity of the pore water through Archie's law. With establishment of empirical or theoretical relationship between the pore water conductivity and the total dissolved solids concentration, quantitative mapping of the level of contamination and its spatial extent could be obtained. Such a result can be used as control measure on modeling of sulphide oxidation or contaminant plume flow. This work was not carried out here and it is an area for future research. A combined DC resistivity and IP survey has the potential to provide information about the amount and distribution of sulphide minerals in the tailings and contribute to the A M D research there in other ways. A typical drawback of this kind of survey— its cumbersome and slow field operation—is minimized in the mainly flat and treeless tailings area. All the above arguments support practical examination of the methods and their application to the A M D research on the tailings. 3.5 Other Geophysical Methods Before we turn to the description of the DC resistivity and IP surveys, it is important to mention other geophysical methods that could be employed. Electromagnetic methods are most suitable to detection and mapping of shallow con-ductors. Fast reconnaissance surveys by airborne E M or detailed quick and cheap pro-filing by ground methods are very common. Use of time domain or multi frequency equipment enables resolving conductivity variations with depth. These techniques have Chapter 3. The Copper Chff Tailings Disposal Area 20 been tried on the tailings but in locations not coinciding with the D C / I P survey and are not included in this thesis work. Pyrrhotite, the main constituent of the sulphide concentration in the tailings, is highly magnetized. That can suggest the use of magnetic methods. The relatively low concen-tration though, about 3 % wt. in average, might not be enough to be detected, not to mention mapped by magnetic methods. The presence of magnetite, which has higher susceptibility then pyrrhotite, in parts of the area, eliminates the possibility. On the other hand, a simple ground magnetic survey which can be carried out by one person is very efficient in locating buried metal bodies and might locate random magnetite con-centrations. Mapping these causes of IP noise might have helped in the interpretations of the IP results but the low quality of the available data did allow an effective use. Seismic wave propagation is not influenced by small variations of mineral concentra-tion in the medium but the tailings-bedrock boundary should be a good target for a seismic reflection/refraction survey. Mapping the depth to that boundary enables cal-culations of the total amount of tailings and help in the estimation of total amount of sulphides. Also, depth to the tailings-bedrock interface might be incorporated in the DC potentials and IP inversions since conductivity and chargeability are expected to change significantly at that boundary. To my best knowledge, seismic experiments were not carried out on the tailings but should be considered in future surveys. The IP method is sensitive to the presence of the possible sources of oxidation—the sulphides— but it can not provide any information about whether oxidation actually occurs, and if it does, about the oxidation rate. Natural electrical potentials in the earth are usually attributed to galvanic cells which are created in zones of oxidation. (Telford et ai, 1990). Self Potential measurements may contribute to the A M D research by locating areas of high oxidation rates and monitor remediation schemes to minimize them. SP data were acquired simultaneously with the D C / I P data. (The natural SP Chapter 3. The Copper CHff Tailings Disposal Area 21 must be measured and subtracted from the DC/IP measurements). The readings are noisy and do not repeat well and thus were not used. In summary, other geophysical surveys could have supplied additional information about the A M D problem and provide constraints upon the DC resistivity and IP work. These data do not exist or are inadequate for that purpose. In the next chapters only the DC resistivity and IP methods are considered. Chapter 4 The D C Resistivity And Induced Polarization Methods 4.1 Theoretical Background The DC resistivity and Induced Polarization methods have long been employed by geo-physicists for various purposes. With the DC resistivity method we look for the electrical conductivity of the earth using measurements of the potentials which are established by an input galvanic current. The method have been used in mineral exploration as well as in the search for geothermal reservoirs and in groundwater hydrology research. The Induced Polarization method detects the time dependent electrical polarization of earth materials and has been used mainly in conjunction with the DC resistivity in base metal exploration. Although concerned with totally different physical properties the field data of both methods are collected simultaneously with the same instrumentation and study of results usually involves cooperative interpretation. 4.1.1 D C Resistivity Method Despite the traditional name of "DC resistivity method", the following formulation will be in terms of electrical conductivity rather then resistivity with the relationship explained later. The basic governing equation is Ohm's law J = o-E (4.1) relating the current density vector J and the electric field E through the conductivity cr 22 Chapter 4. The DC Resistivity And Induced Polarization Methods 23 which is a second order tensor. In the isotropic medium we assume, cr is reduced to a scalar. In terms of resistivity we write E = p3 (4.2) where p is the resistivity. For an isotropic earth the conductivity cr and the resistivity p are simply reciprocals. A units conversion table is given in appendix A . The electrical field E is the gradient of a scalar potential 4>, E = -V<£ (4.3) and thus we have J = -o-V<j> . (4.4) Using Ohm's law in the charge conservation equation V • J = 0 (4.5) we get V • (crV<f>) = 0 (4.6) and in the presence of a current source we have V • (crV<f>) = -IS{r - r.) (4.7) where I is the current source at location rt. Equation 4.7 is a second order differential equation which relates the measurable potential to the electrical conductivity in the presence of a current source and it forms the basis for the DC resistivity method. With the formulation of the potential-conductivity relation and a given electrode array geometry we can define an apparent conductivity <ra. This is the conductivity of a homogeneous earth that would have given rise to a particular potential datum. Most often, for interpretation purposes, data are presented as apparent conductivities. A potentials data set or the associated apparent conductivity data set are interchangeable as input for inversion algorithms. Chapter 4. The DC Resistivity And Induced Polarization Methods 24 Injected current Non polarizable response Polarizable response o . Figure 4.1: Responses of a non-polarizable and polarizable media to square-wave current. 4.1.2 Induced Polarization Method DC resistivity method assumes the conductivity to be independent of time, but in fact many earth materials show time dependent conductivity. Looking at it from another point of view, with an alternating current source, the conductivity is dispersive as function of the frequency. As shown in figure 4.1, turning on an exciting direct current, we notice an instantaneous potential response followed by an increase with time. As the current is turned off the reverse phenomenon—an immediate drop followed by decay with time—is observed. Using an alternating current the phenomenon can be observed as a phase lag between the measured potential and the input current (figure 4.2), or by the change in maximum measured voltage as the frequency of the current is altered. As thorough understanding of the microscopic polarization is not known, a macro-scopic mathematical representation given by Siegel (1959) is commonly used. By his Chapter 4. The DC Resistivity And Induced Polarization Methods 25 Figure 4.2: Phase lag /3 between the current and voltage cycles. approach, polarization effects are represented by a volume distribution M of current di-poles anti-parallel to J a — which is the primary current due to external sources. At each point in the polarizable medium M = —rjJ^, where 77 is a property of the medium called the chargeability. The actual current density J , , is the sum of the primary current 3a and the polarization effects and we can write J , = 1.(1 - V ) . (4.8) In the presence of current dipoles, Ohm's law, in terms of the primary current, is written as = <r(l - V)E (4.9) and equation (4.7) takes the form V • (o-(l - 77)V^) = -16{r - r.) . (4.10) Consider an experiment in a homogeneous and isotropic polarizable medium. When Chapter 4. The DC Resistivity And Induced Polarization Methods 26 a direct current I is turned on, the instantaneous potential field at t = 0 is cf>, = -K . (4.11) cr This is the potential field due to the primary current density 3a. K is a geometrical factor, function of the electrode configuration. Build up of the potential field with time alters it to * . = • <4-12> The difference between $<r = $(* = 0) and <f>n — <j>{t = oo) is the secondary field ^. = ^  - ^ = :i[*JL_ (4.13) cr (1 - n) and the ratio of 4>a/<f>t) is simply rj. Thus we can define apparent chargeabiUty—the chargeabiUty of a homogeneous earth corresponding to pair of <j>, and <$>n measurements— to be the ratio Va = (4-14) The time domain apparent chargeabiUty na can be related to the chargeabiUty frequency domain measurements (Hallof, 1964) and all can serve as input for inversion algorithms looking for the intrinsic chargeabiUty 77. 4.2 D C / I P Measurements Basically the DC/IP experiment involves introducing a current source into the ground and measuring the resultant potentials at points on the surface, (see figure 4.3). Transmitter current of controUed form and magnitude is injected into the ground through two metal stakes. Voltages are measured by a receiver which is connected to a layout of cables and electrolyte filled porous-pot electrodes. In a time domain system the current is a square wave direct current with a fixed on-off time cycle. The potentials <f>v are measured during the "on" time and are used for Chapter 4. The DC Resistivity And Induced Polarization Methods 27 Figure 4.3: Point current source and potential at the surface. The dipoles are grounded at infinity. conductivity computations. After the current is shut-off, the residual secondary potential <f>t is measured during the "off" time. Ideally the measurement should be at the moment of current cut-off. Practically that cannot be done because of E M coupling effects that follow the cut-off. In most systems, values of <j>t are integrated over a time window and normalized by <f>r,. That gives data expressed in time units which are proportional to the chargeability. Other systems average <f>, values of a time window and the division by </>„ gives approximations of the apparent chargeabilities in units of volt/volt. A frequency domain system employs an alternating current as a source. Using two frequencies, (ideally DC and infinite frequency), a frequency effect can be defined by F E = P d c ~ p A C (4.15) PAC where pDC and pAC are the apparent resistivities at the low and high frequencies. The FE is proportional to the chargeability and in the limits of true DC current and an Chapter 4. The DC Resistivity And Induced Polarization Methods 28 infinite frequency, it can be shown (Hallof, 1994) that for relatively small polarization the relation FE ~ rj holds. A relative phase shift between the transmitter current and the receiver voltage can be measured using one frequency. The voltage lags the exciting current by an angle proportional to the magnitude of the polarization.(see figure 4.2). That angle expressed in radians is another measure of the polarization proportional to the apparent chargeability. Both time domain and frequency domain phase data have been acquired and used for the research in this thesis. 4.3 Field Procedures The current and potential electrodes can be set up in various array configurations. Factors like signal-to-noise ratio, E M coupling rejection, resolution, economical considerations and the nature of the target need to be incorporated into the array design. Most common in combined DC/IP surveys are the Schlumberger, Pole-Dipole and Dipole-Dipole arrays shown in figure 4.4. An array can be moved along a traverse line while looking for lateral variations, a procedure called profiling. Alternately, depth sounding can be carried out with expansion of the potential electrodes of an array about a fixed point. Because the fraction of current that flows at depth increases with expansion of the electrode spacing, the penetration depth of the array also increases. Both lateral and depth variation can be detected with soundings at several locations along a line or while profiling with an array which takes several potential measurements at each station. The electrode spread used for all the surface work of this research was the Colinear Dipole-Dipole array (figure 4.4 (c)). The electrodes are set along the line with a com-mon separation designated as the a-spacing distance. The number of potential dipoles Chapter 4. The DC Resistivity And Induced Polarization Methods 29 Figure 4.4: Common arrays, (a) Schlumberger, (b) Pole-Dipole, (c) Dipole-Dipole. is called the n-spacing number. The array is moved along the line in fixed intervals, (equal in length to the a-spacing), thus giving lateral coverage. The different distances of the potential dipoles from the current source result in sensitivity to conductivity and chargeability at different depths. The maximum penetration depth is a function of the length of the array. As the maximum n-spacing number is fixed by the receiver's capabil-ity, the a-spacing distance determines the maximum penetration depth—the larger the spacing, the deeper the penetration. Resolution is also dependent upon the a-spacing but it varies inversely to the a-spacing distance. This implies a trade off between resolution of the array and the maximum penetration depth in the choice of the a-spacing. 4.4 Traditional Data Plotting And Interpretation The traditional interpretation of DC resistivity and IP Dipole-Dipole data is based on data profiles or on pseudosection plots. The profiles are simply data points—apparent Chapter 4. The DC Resistivity And Induced Polarization Methods 30 Figure 4.5: Plotting points of apparent conductivity/chargeability for dipoles n=l,3. conductivity or apparent chargeability—plotted against the corresponding station loca-tions. Interpretation is done by matching the plotted data with derived theoretical curves assuming certain geological structure. Pseudosections display the data with the effect of variable electrode spacing illustrated on the plot. Apparent conductivity and apparent chargeability values taken at each station are plotted at the right angle vertex of an isosceles triangle, the other vertices of which are in the midpoints of current and potential dipole (figure 4.5). The values appear at points below the surface at increasing vertical depths corresponding to the n-spacing number of the potential dipoles. Interpretation of the pseudosections is done directly on the plots under the assumption that the image seen in the plot is a direct reflection of the geological structure. Pseudosections can also be interpreted with the aid of forward modeled synthetic data plots pre-derived analytically or numerically from simple shaped anomalies. Chapter 4. The DC Resistivity And Induced Polarization Methods 31 03 Figure 4.6: (a) Conductivity model — 0.07 S/m conductive layer overlying 0.005 S/m half-space and 10.0 S/m small anomaly, (b) Chargeability model—0.01 Volt/Volt half-space, 0.05 Volt/Volt large anomaly and 0.5 Volt/Volt small anomaly, (c) Con-ductivity pseudosection, (d) Chargeability pseudosection. Chapter 4. The DC Resistivity And Induced Polarization Methods 32 Examples of conductivity and chargeability pseudosections are given in figure 4.6. Figures 4.6c and 4.6d show apparent conductivity and chargeability generated by forward modeling of the synthetic models shown in figure 4.6a and 4.6b respectively. Values of the apparent chargeability and logarithms of the apparent conductivity to the base ten are colour coded and interpolated to yield smooth colour plots. The conductivity model consists of a conductive layer overlying a resistive half-space and a very conductive small body buried close to the surface. The chargeability model consist of a lairge% chargeable body at depth and a chargeable small body which coincides with the small conductive body close to the surface. The data are modeled as if they are collected by a Dipole-Dipole array, n-spacing equal six and a-spacing 5 metres. The pseudosections delineate quite well the lateral extent of the bodies. The conduc-tive layer, the big chargeable body and the small chargeable, conductive body all appear in the sections. Notice, though, that there is no quantitative relationship between the n-spacing number and the real depth and that the "apparent" values have no physical meaning. Exception is the conductive layer where the apparent conductivity equal to the real conductivity. 4.5 The Geophysical Survey As part of their efforts to identify appropriate geophysical methods for environmental applications, INCO Limited's Environmental Control Group and Alan King of INCO Exploration and Technical services Inc., initiated a project to assess the use of geophysical methods on A M D problems. During the summers of 1992 and 1993, airborne, ground and borehole E M measurement and combined DC resistivity and IP surveys were conducted on tailings disposal areas around Sudbury. The DC/IP survey was carried out in August, 1993 on lines A A and BB (figure 3.2), Chapter 4. The DC Resistivity And Induced Polarization Methods 33 the author taking part in the field work with the intention to invert and interpret the data as part of his M.Sc. work. Insufficiency of the data and their low quality which was revealed by subsequent analysis prompted INCO to redo and extend the survey. That was done during October-November, 1993. The second survey was conducted along guidelines suggested by us. Part of the experiment cost was financed by an NSERC grant. The rest of the section describes the survey and the available data of the second survey. 4.5.1 Survey Lines Line AA is approximately 3 km long and crosses the CD, M and M l tailings impound-ments where tailings were deposited from the 1940s through the 1960s. It follows one of the major flow-paths within the tailings and its physical and chemical hydrogeology are described in detail in Coggans (1992). It is expected that the cross section of the Une is typical in terms of physical and electrical properties. As the beginning of the Une goes through a lake, the DC/IP survey covers only the part of the Une which crosses the M and M l impoundments. That section is long, relatively straight and flat. The Waterloo Centre for Groundwater Research which conducts ongoing hydrogeological and geochemical research in the taiUngs area, installed a few sampUng weUs and piezometer nests along the section and their analyses results are available. All these facts make the section an exceUent site for testing the DC resistivity and IP methods. Line BB is short, about 500 metre long, which is perpendicular to the Pistol Dam. In the last years the area was the focus of the Waterloo research group which has already identified high sulphate concentrations in a plume extending below the dam. The Une is short and its topography is severe—elevation change is about 30 metres. These facts and the fact the dam was constructed with low sulphide tailings and that the water table is much deeper then throughout the impoundments make the area un-typical and not Chapter 4. The DC Resistivity And Induced Polarization Methods 34 very suitable for the purpose of the research. Data taken on that line are used only for reference and comparison. 4.5.2 Field Procedures, Considerations And Available Data The Colinear Dipole-Dipole array with n-spacing equal six was used for all the surface surveying. The depth of interest of the survey can exceed 40 metres. The a-spacing used in the first survey was 10 metres resulting with an eighty metres long array. It was felt that this might not be enough to sense information from the tailings-bedrock interface. It was also desireable to get an estimation of the bedrock properties to serve as reference values for the inversions. Thus it was suggested to survey line A A with both 10 and 20 metres spacings. That could also be beneficial in practical terms. Survey speed of a long line is almost doubled by doubling of the array length, ( and halving the number of stations). If a thorough survey of the whole tailings area is to be carried out, the economical consideration of total field work time would have to be taken into account. Conclusions about relative merits of different array spacing will be important while taking these decisions. The Colinear Dipole-Dipole configuration is sensitive to lateral variations along the line. The different n-spacing of the potential dipoles provide sensitivity to variations with depth. Naturally the exploration is mainly in two dimensions, but the conductivity and chargeability structure off line can influence the data. To assess the two-dimensionality of the problem, it was suggested to survey two lines parallel and close to line A A . Two different receivers were used, the Phoenix IPV-4 and the Scintrex IP-6. The first is a frequency domain instrument measuring the phase difference between the altern-ating current source and the resultant potentials. The latter is a time domain receiver. Although not intentional, the use of the two different systems enabled the comparison Chapter 4. The DC Resistivity And Induced Polarization Methods 35 between their outputs and added confidence in data reliability when agreement was ob-tained. The final decisions about the survey design were taken according to the considerations outlined above and the economical constraints. All the section of Une AA from the lake to the dam was surveyed with the time domain system at 10 metres spacings. The first 800 metres of that section and two parallel Unes 30 metres east and west of it, (AA-East and AA-West), were surveyed with the frequency domain system at 20 metre spacings. Line BB was surveyed with both systems at 20 metres spacings. In order to ground-proof, and if needed, caUbrate inversion results of the surface DC resistivity and IP data, in-situ measurements of conductivity and chargeabiUty were needed. It was suggested to take DC/IP borehole measurements in weUs along Une A A. Borehole measurements were taken only in one location, at weU IN-10. More borehole measurements and laboratory research of IP response of sulphide taiUngs were to be done by INCO later. The UBC In-Situ Testing group developed an electrical conductivity probing cone capable of providing a continuous record of conductivity with depth in addition to geotechnical parameters. Cone testing was carried out on the taiUngs in October 1993. Unfortunately none of the in-situ testing sites coincided with Une AA so the cone readings can be used only as general reference. 4.6 Data Plots And Initial Interpretation 4.6.1 Data Plots The time domain data set covers the whole length of the survey Une and its pseudosections do not significantly differ from the frequency domain data sets in the region where they overlap. Thus, for the purpose of the initial interpretation, only time domain data will be presented. Any discussions regarding the frequency domain data of Une AA and of Chapter 4. The DC Resistivity And Induced Polarization Methods 36 the parallel lines AA-West and AA-East will be carried out in the next chapter, with the inversion results at hand. Figure 4.7 shows plot of logarithms to the base ten of the apparent conductivity values. Figure 4.8 shows the apparent chargeabilities. Locations along the line are in metres from the zero point at the location of well IN13 (see figures 3.2 or 3.3). Positive locations are north of the well, negative ones, south of it. Apparent conductivities along the line, at least most these recorded with the n-spacing 1 to 3 are very high (0.05- 0.10 S/m). The main features in the conductivity pseudosection are more conductive region between locations 0 and -200 followed by a less conductive section extending to location -1150. Then come two very conductive regions followed by a very resistive one from location -1500 to the end of the line. There are also two small "pants-leg" conductive features at locations -290 and -1160, the latter is obscured somewhat by the high surrounding conductivity. Very low apparent chargeability values are most dominant along the chargeability pseudosection (figure 4.8). Slightly higher values appear at the northern (+60 to -100) and southern (-1680 to -1730) ends of the line. Strong "pants-leg" features coincide with the conductivity highs at locations -280 and -1160. The area between locations -1150 and -1500 contains varying very high values. Some of them, around location -1450 are so high that the second half of the section had to be plotted in a different colour scale in order not to obscure any other values. Still some small variations which occur between locations -350 and -1150 do not show up in the plot due to the colour scaling. 4.6.2 Initial Interpretation The generally high apparent conductivities are typical of an area with highly conductive acidic pore water. Some level of oxidation probably occured throughout the length of the line contributing to the high conductivities down to location -1500. The last 200 Chapter 4. The DC Resistivity And Induced Polarization Methods 37 metres of the line are actually On the down slope of the dam. As seen in figure 3.3, the water table is much lower there and that explains the low conductivity. It is difficult to determine from the pseudosection the depth of the higher conductivity layer. High conductivity values, recorded in all the n-spacing dipoles, at the northern part of the Une (60 to -200) and south of location -1140 where the Une enters the M l impound-ment, might be the result of higher dissolved soUds concentrations. Between locations -1400 and -1500 the conductivity values are exceptionally high. The length of that feature excludes the possibiUty of a metal body to be the source. A good possibiUty is that it is somehow related to a 30 centimeter thick layer of almost pure magnetite encountered at well IN12, (location -1409). If that layer extends to the south, it might be the cause of the high conductivity. The "pants-leg" features at locations -290 and -1160 are typical of small very con-ductive bodies. Their sources are metal pipes, reported by INCO, at locations -290 and between -1160 and -1170. These coincide with similar high chargeabiUty features at the same locations. The average sulphide concentration throughout the taiUngs is about 3 wt. %. The estimated concentration around weU IN10 (location -510), where taiUngs processed in the iron recovery plant were deposited, is less then 1 wt. % (King, 1994). Thus the generally low apparent chargeabiUties are not very surprising. The higher values at locations +60 to -100 and -1200 to 1400 can be attributed to higher sulphide concentrations. The higher chargeabiUty values at the end of the Une (-1700 to -1730) are surprising as this is the area of the dam itself which is known to be constructed from low sulphide taiUngs (King, 1994). It might be a small metal body buried at the base of the dam or possibly a topographic effect. Chapter 4. The DC Resistivity And Induced Polarization Methods 38 n-spacing n-spacing Figure 4.7: Conductivity pseudosection. Chapter 4. The DC Resistivity And Induced Polarization Methods 39 Figure 4.8: Chargeability pseudosection. Note the difference in the colour scaling of the two parts of the plot. Chapter 5 Inversion And Interpretation Of D C Resistivity And IP Data "The interpretation of geophysical data is a curious mixture of art, science, training, experience and intuition (or luck)." (Sumner, 1976). Probably training, experience and intuition will always be needed for geophysical data interpretation but the advances of inversion methods are intended to reduce the roll of luck in the geophysical work and make the art more scientific. The methodology for inversion of DC Resistivity and IP data given by Oldenburg and Li (1994) will be summarised at the beginning of this chapter. It will be followed by a short description of the program library DCIP2D which implements the methodology in a series of computerized algorithms. During the work, generic problems of the inver-sion methodology and ones typical to the taiUngs environment were encountered. The problems, ways to overcome them and results are given in the remainder of the chapter. 5.1 Inversion Methodology 5.1.1 Subspace Non-Linear Inversion The goal of an inversion is to find, for a given set of N observations dj ( j = 1,2,N), a model m(x) that reproduces the observations through the relation dj = Fj[m{x)} • (5.1) 40 Chapter 5. Inversion And Interpretation Of DC Resistivity And IP Data 41 The functional ,F/[m] is the forward mapping operator. It maps elements in H, the Hilbert model space, into elements in £, the N dimensional Euclidean data space. The technique used to achieve the inversion's goal defines an inverse mapping m(z) = Ff1[di] j = l,..,N. (5.2) Here DC resistivity and IP data inversions were carried out using the subspace inversion method of Oldenburg et al. (1993). The general methodology is the same whether the model is conductivity a or chargeability n and in the following both will be generically denoted by m. The method starts with partitioning the model domain into M rectangular cells as-suming constant model value mj, within each of them. The model is parametrized as M m = ]T) mklk (5.3) Jt=i where the fk a r e basis functions defined on the model space and are unity over the region occupied by the the cell and zero elsewhere. The problem thus, is reduced to finding a vector m = {m^mj, ...,mjvf} which reproduces the observations vector d°6* = {d^ 6*, dj6*, dw*} to certain accuracy. The partitioning is designed such that the number of model parameters M is much larger then the number of the data N. The large number of model parameters enables good representation of the earth by the discrete approximation and allows more thorough exploration of the model space than does a coarse parameterization. The degree to which the data are to be fitted is a problem which depends on the esti-mated errors of the data. Real data are never accurate, so exact fit ensures the recovered model will be wrong—fitting the additional noise demands an unwanted extra structure from the constructed model. On the other hand by under-fitting the observations we can miss features of the model contained in the data. Various measures of the misfit are possible. The one used in the following is the L2 norm. Chapter 5. Inversion And Interpretation Of DC Resistivity And IP Data 42 The misfit objective function is defined as Md,d°b')=\\Wd(d-d°»')\\2 (5.4) where d is data set predicted by forward modeling of m and Wd is an 7Y x N matrix. Assuming each datum dj to be contaminated with uncorrelated Gaussian noise with zero mean and standard deviation 6j, then Wd should be a diagonal matrix with the elements of the diagonal the reciprocals of the standard deviations ej. tpd is then distributed as a chi-squared variable with N degrees of freedom. Its expected value is N and its variance V2N. Fitting the data adequately does not ensure having the right model. The inverse problem is inherently non-unique, fitting of the data with one model means that there is an infinite number of models which can fit the data to the same degree. From all those models we wish to find one which has certain characteristics. Choosing the model which has the minimum norm is a reasonable option. Thus in addition to minimization of the misfit objective function, the minimization of a model objective function ipm(m,m0) = ||Wm(m - m 0)|| 2 (5.5) is sought. The matrix Wm is a model weighting matrix and m 0 is a reference model. Min-imization of ij)m will yield a model which is close to the base model mo in the Euclidean sense with the norm weighted by Wm. To achieve a model which fits the data adequately and has the desired characteristics, a general objective function V>(m) = Vm(m, m0) + p, (^(d, d"6') - j,*) (5.6) is to be minimized. The -0* is a desired target misfit and fi is a Lagrange multiplier. Both DC resistivity and IP problems are non-linear. A standard way to find a solution to a non-linear problem is by an iterative process. At each iteration, a model perturbation Chapter 5. Inversion And Interpretation Of DC Resistivity And IP Data 43 6m is sought such that when added to the last iteration model, (or to an initial guess in the first iteration), the forward modelled data fits the observed data better. Iterations are repeated until a model which reproduces the data with an acceptable fit is achieved. A proper perturbation can be found by first linearizing equation (5.6) about the current model. Differentiation of the linearized equation with respect to the model pa-rameters m,j and equating to zero yields an M x M system of equations to be solved for 6m. This must be done a few times at each iteration while searching for a value of p such that a linear target misfit is achieved. It must be emphasized that at each iteration a linear problem is solved to find a perturbation 6m which fits corresponding data perturbation to a target misfit. The choice of the iteration linear target misfit is quite arbitrary and determines only the form of the perturbation. The solution model of the non-linear problem is controlled by the misfit as defined by equation (5.4). As M becomes large the iterative procedure becomes computationally expensive. The subspace method avoids solving the M x M system by restricting the model perturbation to lie in a subspace of the model space spanned by q search vectors { V J } i = 1, q . The perturbation is expressed as the linear combination 6m = aw (5.7) »=i and the system of equations to be solved for the parameters a; is only q x q in dimension. In the subsequent inversions a typical value of M is about a thousand while the number of search vectors, (and parameters to solve for), is a few dozen. As solution of an M x M system involves about M 3 operations, the computation time saving is clear even though more iterations will be needed to reach the final model. The choice of the form of the model objective function of equation (5.5) is of prime importance. The objective function must be such that its minimization would yield a model with the desired characteristics required by our physical understanding and which Chapter 5. Inversion And Interpretation Of DC Resistivity And IP Data 44 reflects any a priori geological knowledge. Both conductivity and chargeability structures are expected to be smooth. In a two dimensional inversion, that requires minimization of the horizontal and vertical derivat-ives. Some measure of control on the model parameters' magnitude is also needed and may come with the requirement that the model will be close to a geologically reasonable reference model. These two requirements must be complimented by a flexibility in the objective function to produce different models which fit the data and the ability to in-corporate in the solution any a priori knowledge about the problem. The flexibility to produce different models is an obvious need due to the non-uniqueness of the solution. We want to be able to examine different models and verify that the main observed fea-tures are demanded by the data and not just artifacts of the mathematical process. The ability to incorporate a priori knowledge is important not only because it enables us to fix known model values but also because by doing so we drive the solution in the "right" direction and reduce the non-uniqueness. An objective function capable of accomplishing the above goals is i(}m(m, m0) = ft* f J ws(m — m0)2dx dz+ The constants /?,, 8m and 0Z weight the relative importance of the components of the model objective function. The functions wa, wx and wz are weighting functions defined on the model domain. By specifying certain forms for these functions the inversionist can get desired features in the model which fit his/her geological understanding. The discrete form of equation (5.8) is i M m , m0) = (ra - m 0 ) r {fi.WjW, + /WjW* + 0mW*W.} (m - m0) (5.9) Chapter 5. Inversion And Interpretation Of DC Resistivity And IP Data 4 5 where the matrix combination {B.WjW, + 8xWxrWx + 0MW?W,} is the matrix Wm of equation ( 5 . 5 ) . 5.1.2 D G Resistivity Inversion Definition of the forward mapping of the DC resistivity problem is by the equation V •{crV<f>) = -I6(r-r.) ( 5 . 1 0 ) and the appropriate boundary conditions : d<j>/dn = 0 at the surface and <j> —• 0 far away from the current source. The above boundary value problem is solved by the the finite difference method of Dey and Morrison ( 1 9 7 9 ) . This forward solution is incorporated into the subspace inversion algorithms. The model domain partitioning mesh is the same for the forward and the inverse problems and it is fine enough to create a strongly under-determined problem. In the inversion the model to be minimized is ln(cr) instead of the conductivity cr. That reduces the range of values the inversion algorithm needs to deal with and ensures positivity of the solution. 5.1.3 Induced Polarization Inversion Data for the IP inversion are the apparent chargeabilities defined by equation ( 4 . 1 4 ) . Using equations ( 4 . 1 4 ) and ( 5 . 1 ) one can define the IP forward mapping as F[cr{l - v)} -Va ~ - v)} ( 5 , 1 1 ) where . F [ C T ( 1 — 77)] and J-[cr} are the DC resistivity forward mapping operators in polar-izable and non-polarizable earth respectively. Oldenburg and Li ( 1 9 9 4 ) give three methods to invert the IP data. All the methods assume the conductivity structure to be known a priori or by inverting the DC resistivity data. Chapter 5. Inversion And Interpretation Of DC Resistivity And IP Data 46 One way is to solve the non-linear problem as described in the inversion method-ology subsection. Iteratively, the data equations are linearized about the last iteration chargeability model and a perturbation is found and added to the model. The second method involves manipulating the results of two inversions. The formal inverse mapping of equation (5.2) when applied to a <6», data set is a(l - n) = ^-1[4>v] . (5.12) Inversion of the (j>a data set is expressed as <r = • (5.13) Thus the chargeability n can be calculated by The third method assumes n <C 1. Neglecting high order terms, the potential measured over a medium composed of M different in polarization cells can be expressed as M Ql *, = 4>{o-(\ - r,)) = 4>{a) - £ j^rvm . (5.15) »=1 Using the above equation and equation (4.14), the apparent chargeabilities can be written as - E J£w<* ** = M \ % r*t. (5-16) and for small chargeabihties this can be approximated by £<nd4> £ dln<p ^ < 6 c V ' ' fri dinar In a matrix form the last equation can be written as M '7«i=]CJ^' j = 1 , J V (5.18) i=l Chapter 5. Inversion And Interpretation Of DC Resistivity And IP Data 47 (5.19) Equation (5.18) forms an approximate forward mapping which is incorporated in the 5.2 The DCIP2D Package The DCIP2D program library developed by the UBC-Geophysical Inversion Facility implements the above methodology in a series of five computer programs. The library can perform DC resistivity and IP forward modelings and inversions over two dimensional earth structures. Any electrode configuration along a Une can be handled. The forward modeUng programs DCF0R2D and IPF0R2D forward model conduc-tivity and chargeabiUty models suppUed by the user to DC potentials and apparent chargeabiUties data files. The data are generated using a finite difference code with the mesh and electrode positions defined by the user. The Inversion programs DCINV2D and IPINV2D invert DC potentials and apparent chargeabiUties data files to conductivity and chargeabiUty structures. The user suppUes the data files and parameters for the finite differences mesh with the same file format as the ones used for the forward modeUng programs. Special weighting schemes can be suppUed to modify the weighting matrices of equation (5.9) in another file. The IPINV2D and IPFOR2D programs use the method of Unearization of the data equations in the forward modelling of the IP data. The sensitivity matrix J of equation (5.19) is generated by the fifth program IPSENS2D. subspace formaUsm instead of equation (5.11). Chapter 5. Inversion And Interpretation Of DC Resistivity And IP Data 48 5.3 Typical Problems In The Inversion Of The Tailings Data Sets 5.3.1 High Conductivity Of The Upper Layer The water table throughout the tailings is 0-5 metres below the surface. Along Une AA it is about 2 metres deep. The dry layer above the water table is very resistive but the conductivity rises sharply below it to values of about 0.1 S/m. Because of the high conductivity below the water table the measured voltages are very small and to get a signal above the noise level, requires that the injected current be high. One of the problems of the first survey was that the transmitter was Umited to maximum output of 3 Amperes and long segments of the Une could not be surveyed. Even with maximum current of 10 Amperes, the DC voltages of the second survey data sets are relatively small. Small signals usuaUy tend to have low signal to noise ratio which is always undesirable. The problem is even worse with the IP time domain data as the secondary potentials used for the apparent chargeabiUty calculations are of the order of one hundredth of the primary potentials. Another issue related to the high conductive layer is that most of the current is channeled through it and thus the data are not sensitive to the bedrock conductivity. This is not a real problem as the focus of the research is on the taiUngs but the desire to get real bedrock conductivity in order to use it as a reference to subsequent inversions is probably destined to fail even with large electrode separations. Assuming smooth variation of earth conductivity, the inversion routines minimize the first derivatives of the conductivity model. The assumption breaks down in the inter-face between the highly conductive taiUngs and the bedrock. The very low hydrauUc conductivity of the bedrock causes the ground water to flow along that interface with minimal vertical penetration (Coggans, 1992). As the rock's conductivity is mainly due to the pore water conductivity, the result is a large and sharp conductivity change. The Chapter 5. Inversion And Interpretation Of DC Resistivity And IP Data 49 minimization of the vertical derivative causes "smearing" of the model in the vertical direction. The bottom interface of the taiUngs cannot be inferred accurately and con-ductivity values, especially near that boundary, are affected and do not reflect the true conductivity. 5.3.2 Buried Metal Bodies The taiUngs disposal area is a man-made structure. Metal bodies Uke pipes, tanks, fences and other metal debris are scattered throughout the area. At least two pipes are known to be in the vicinity of Une A A. These pipes have very characteristic responses in both the DC resistivity and IP data sets. In the pseudosections (figures 4.7 and 4.8) these responses dominate and mask the natural responses over long stretches of the Une. Inversion of the data can minimize that effect and locate the exact position of the metalUc body. The problem is that very fine model parameterization is needed to properly represent small bodies. Although the refinement can be done locally, it adds more ceUs and thus requires more computer time for the inversion. Even if that is done the results are not guaranteed to be good. With a fine parameterization a special weighting for the horizontal and vertical model derivatives minimization wiU be needed to prevent projection of the high conductivity values out of the conductive area. That might be very difficult as the pattern of the "smearing" of the conductivity is compUcated. A lot of trial and error will be needed to tailor a fine mesh and corresponding weighting scheme for each body and success is not guaranteed. If the fine parameterization solution is not used then the problem of big errors due to poor parameterization around very conductive small bodies arises. That problem and the way to deal with it are discussed in the next section. In addition to the poor parameterization, another problem concerns the inversion of chargeabiUty responses due to metal bodies. Wynn and Zonge (1975) discuss the Chapter 5. Inversion And Interpretation Of DC Resistivity And IP Data 50 so-called "pipeline effect" and suggest the effects of current focusing nonlinearities, elec-tromagnetic induction and complex electrode polarization as possible reasons. In the tailings very high IP responses are measured in the vicinity of metal pipes, and other unknown metal bodies, demanding high chargeabilities for the corresponding cells. The IPINV2D inversion program uses the method of linearization of the data equations. The basic assumption of that method, permitting the linearization, is that the chargeabilities of the model cells are small. That assumption does not hold in the vicinity of the pipes and the resultant chargeability structure can be distorted. The problem falls in the same category of the conductivity parameterization problem and is also dealt with in the next section. 5.4 Assigning Errors To The Data The goal of the inversion algorithm is to find a model which reproduces the observations to a certain degree of accuracy. The degree of fitness, or as it is usually called—the misfit, which we want to achieve should be determined by the data errors. Unfortunately usually we have no estimate of the amplitude of the errors nor of their distribution. The common assumption is that the data errors are Gaussian and independent. This is probably not true in most real world situations. Moreover, even using that assumption which enables us to decide about the degree to which the data are to be fit, we have yet to find the way to ascribe an error associated with each datum. Data errors should include both the measurement errors due to possible current leakage, telluric noise, E M effects, instruments inaccuracy etc . , and the errors arising due to the approximation of the 3D earth by a 2D, parametrized model. I chose to assume nothing about the physical or mathematical mechanisms contributing to the inaccuracies. Instead I take a very general approach and divide errors into three categories : Chapter 5. Inversion And Interpretation Of DC Resistivity And IP Data 51 a) Constant errors b) Proportional to the datum errors c) Random, erratic errors While assigning errors, each of these error types should be addressed in an appropriate way. The error assigning scheme, (and I have an automated one in my mind), can use only statistical information to accomplish that. Oldenburg and Li (1994) assigned a constant error—percentage of the average apparent chargeability—to IP phase data and percentage of each datum error to DC potentials data. Errors combined of both a percentage of the datum and a constant value were also tried. Inversion of DC potentials is usually straight forward. The basic conductivity struc-ture is reached with percentage error scheme that could be scaled later to fit borehole or in-situ measurements by adjusting the proportionality factor. We must realize though, that it has an inherent deficiency in it. Data values which are around the data average for each n-spacing and more likely the more accurate ones, are the most discriminated against. The small data values get assigned unrealistically small errors. Large data val-ues are assigned absolutely large errors but not large enough to be fitted to satisfactory degree. With such a scheme only the proportional to the datum errors are addressed. An error scheme which combines percentage of the datum value plus a constant addresses also possible constant errors. It is reasonable to choose a constant error which is a percentage of some average data value. As DC potentials decrease with distance from the current source, their relative accuracy is expected to worsen, hence requiring relatively larger assigned errors. Choosing a constant which is a percentage of the average of all the data achieves both targets. The assigned errors are proportional to the average and become relatively larger for the smaller data values. The error assigning scheme employed for all the DC potentials Chapter 5. Inversion And Interpretation Of DC Resistivity And IP Data 52 inversions was ERR = A • DATUM + B • MEAN (5.20) where ERR is the error assigned to DATUM, A and B are proportionality factors and MEAN is the average data value of the data set to be inverted. The proportional to the mean component is added to address possible constant errors. As such it should not reflect unusual data perturbations like the ones arising in the vicinity of metal debris. The MEAN value then, should be computed without considering outliers. The values of A and B as well as the definition of an outlier should be chosen considering the variance of the data values. The IP measurements are fractions of the potentials values—in the time domain—or small phase shifts—in the frequency domain. Thus, the IP data sets tend to be noisy. Moreover, apparent chargeabilities may be zero or negative. As an example see figure 5.1 which shows the time domain data between locations 60 and -490 in a line format. Zero or very small values are expected in a non-polarizable medium but may be erroneous otherwise and must be assigned appropriate errors in order to be fitted properly. Further complication arises with the negative apparent chargeabilities. These can occur in certain geoelectric layered earth situations (Nabighian, 1976) or near three-dimensional bodies and near-surface contacts (Siegel, 1959) but are quite unlikely. On the tailings the negative values are observed near small buried metal objects. All forward modeled data in the vicinity of these objects are prone to be fairly erroneous due to the poor parameterization and must be assigned relatively large errors. The negative values might be small in their absolute value but must be assigned large errors to let the inversion algorithm have freedom in exploring the model space in poorly parametrized regions. The way to address error assigning problem for small data values is to incorporate a constant value in the assigned error. Unlike the DC potential data, using the same Figure 5.1: Apparent chargeability data in V / V between locations 60 and -490 given separately, in a line format, for each potential dipole. Plotting point is the midpoint of the potential electrodes at each position. array, the apparent chargeabiUties values are a function of the chargeabiUty distribution only and are not expected to decrease if measured far from the current source. Along a few sections of Une AA, apparent chargeabiUties actually increase with the n-spacing. Thus, for the IP data it is better to compute averages of data recorded at the same dipole separately and use them as reference while assigning the errors. As most of the IP data are positive the average data value is positive. A simple way to systematically discriminate against the negative data is to add to the assigned error Chapter 5. Inversion And Interpretation Of DC Resistivity And IP Data 54 an amount that is proportional to the difference between the datum value and the mean data value. In a relatively homogeneous polarizable medium that component would be negligible but it would enlarge the errors assigned to any data which deviate from the average and negative data would be assigned larger errors then their equal in absolute value positive counterparts. There is still one more problem which needs to be addressed—that of huge chargeabil-ities, probably badly inaccurate which occur very close to the buried metal bodies. Al-though expected to be large and cannot be removed from the data set, these data, even modeled reasonably still can contribute a large and unwanted contribution to the Eu-clidean norm misfit. Outliers must be assigned very big errors but in a way which would not affect regular data. The natural choice is an exponential component to be added to the error assigning scheme. The final scheme used to invert the apparent chargeability data sets was ERR = A • \DATUM\ + B • MEAN + C • \DATUM - MEAN\ + D • D l • E (5.21) A,B,C and D are proportionality factors. Dl is the sum of the first three components and can be included in D. The last component is the exponential which may be formed in different ways. I chose where F is another factor determining from which data value the E component will become dominant. The assigned error of a datum which equals the mean value is pro-datum value gets further from the average its assigned error increases according to the other parameters. Figure 5.2a shows the assigned errors as function of data values where unity value is the mean. Figure 5.2b shows the errors relative to the corresponding data. (5.22) portional to the sum of A and B, (slightly bigger—because of the last term). As the Chapter 5. Inversion And Interpretation Of DC Resistivity And IP Data 55 Generally, the relative assigned error becomes larger as the datum value gets further from the mean. Around the zero datum value the relative error increases as the datum value approaches zero and negative data always get larger relative errors then their positive counterparts. a 1.2 [—1—1—r—[—T~~I—1—|—1—•—•—|—•—•—•—|—•—•—•—|—i i r 0.8 r N. / 0.6 : N. / 0.4 " \ . / 0.2 h Q E 1 \ I I I I L I I I I I I t _ _ l I I > • 1 • ' • " -6 -4 -2 0 2 4 6 DATUM VALUE Figure 5.2: (a) Assigned error as a function of datum value. Datum value of 1 is taken as the mean of the data set. The parameters used are A = 0.05, B = 0.1, C = 0.1, D — 0.1 and F equals 3.0 for positive data and 2.0 for negative ones, (b) The error relative to the datum as a function of the datum value. The above scheme is a bit complicated. Its great advantage is in facilitating relatively flat map of misfits between the observed and predicted data, normalized by the assigned errors. (See figures 5.3a and 5.3b). Achieving that, the effort can be directed to finding Chapter 5. Inversion And Interpretation Of DC Resistivity And IP Data 56 the right level of fit without worrying about improportionally big contributions to the misfit from certain data. A much simpler scheme which involves only a constant error value plus the exponent component works as well but the flexibility in exploring the model space is reduced and biased fit of certain features cannot be treated. Figure 5.3 shows inversion results and the corresponding misfit maps of the data between 60 and -460. Model 5.3A was produced using the error assigning scheme given in equation (5.21). To produce model 5.3B the parameters A and C where set to zero, B is 0.15 and the value of MEAN was set to be the average value of the whole data set. The corresponding misfit map in figure 5.3b is rougher. The model is quite similar to the one of 5.3A with slight differences around location 0 and addition of artifacts due to the influence the pipe at location -290. Only with a lot of experiencing and good ground proofing it will be possible to determine whether the improvements of the more complicated scheme are significant. 5.5 Inversion Results And Interpretation 5.5.1 Parameters For The Inversion Having the data plotted, inspected and providing an initial interpretation from the pseu-dosection is a preceding step to the inversion. That done, an initial inversion should be carried out. The resultant earth model serves as a basis for the following inversions which explore the model space and are guided by the physical and geological understanding of the problem. There are a few parameters which the DCIP2D inversion routines require from the user. These parameters are necessary for the programs to run and later on can be modified to produce more geologically reasonable model if required. First of these parameters is the mesh which is used by both the forward and inverse algorithms. The second parameter Chapter 5. Inversion And Interpretation Of DC Resistivity And IP Data 57 is the final desired misfit of the data. Third and fourth are the initial model from which the inversion begins and the reference model to which the model should be close. Last are the weighting matrices and the weights for the components of the model objective function. To these user controlled parameters we can add the data errors which were discussed in previous section and are supplied to the programs with the data files. Parameterization mesh The mesh structure is mainly determined by the spacing of the array used to collect the data, the depth of interest and computation time considerations. Ideally the mesh should have as many cells as possible to minimize discretization errors in the finite difference code, but cell dimensions much smaller then the array spacing are redundant as they go beneath the resolution of the data. I started the inversions with cells width, in the area of interest, equal to the electrode spacing. The cells thickness depends on the depth and features of the target and generally should increase with depth. The aspect ratio of the cells' dimensions is limited by the accuracy of the numerical calculation of the forward modeling. As the area of interest of the research is relatively shallow, I chose the cells thicknesses to be the minimum permitted by the aspect ratio limit down to 40 metres depth. The same mesh is used for both the resistivity and IP inversions. Target Misfit Because the final estimate of the error assigned to each datum is somewhat arbitrary, the choice of the target misfit is not absolutely determined. Our first choice is to set the target misfit equal to the number of data. Iteration in the inversion continues until the target misfit is reached and the model norm is stabilized. If the first condition cannot be achieved, a map of misfit between the predicted and observed data should be consulted. If there are unusual contributions from certain features, the parameters of the error assigning scheme should be adjusted. Chapter 5. Inversion And Interpretation Of DC Resistivity And IP Data 58 For all the survey data sets the same error assigning scheme yielded a relatively flat misfit map but in most cases the target misfit could not be achieved—estimation of the level of the data errors was too optimistic. The target misfit was then adjusted to a higher level and the inversion was repeated, iterating, until the target was achieved and model norm did hot decrease further from one iteration to the next. A few more inver-sions followed with increasingly larger values of target misfit and the most geologically reasonable model was chosen. Reference And Initial Models The model objective function is designed such that its minimization will result in a model which is close to a reference model. If a good idea of the geological true model is available it should be incorporated into the reference model of the inversion. Unfortunately none is available for the survey area and thus my reference models are half-spaces and serve as physically reasonable stabilizing values to prevent unreasonable values of the recovered models. Conductivities values which are intermediate between the conductivity of the taiUngs and the conductivity of the bedrock have been used for inverting the DC res-istivity data while zero chargeabiUty model has been used to invert the IP data. The purpose is to let the data constraints drive the solution to an accurate model stabiUzed around these reference values. The choice of an initial model is not critical. It is desireable to choose one which is close to the final model so a minimum number of iterations wiU be required. I usually choose the initial model to be the same as the reference model. Weights One of the advantages of the inversion methodology which was expounded in previous sections is the abiUty to drive the solution to one which agrees with the physical and geological understanding of the interpreter. A priori knowledge can be incorporated in the solution through the use of the weights of the components of the model objective function Chapter 5. Inversion And Interpretation Of DC Resistivity And IP Data 59 and the weighting matrices. With no known values of conductivity and chargeability, more weight should be given to the "flattest" model components. Thus the 8X and 8Z parameters were set to unity and 8S was given a small value of the order of 1/1000. The weighting matrices for the initial inversions were set such that equal weight was given to each cell. In following sections use of modification of these matrices to improve results will be explained and demonstrated. 5.5.2 Inversions Results Figures 5.4 and 5.5 show conductivity and chargeability inversion results of the time do-main data along line A A. The data set was divided into four sections in order to limit the size of data and model parameters for each inversion. Each of the four sections overlaps 150 metres of its adjacent sections so the final result is not affected. In figure 5.6 are the results of the inversions of the frequency domain data of lines A A , AA-East and AA-West. All the sections suffer from distortion of the results at the edges due to lack of data coverage at the ends of each section. The pattern of distortion depends upon the direction the array moved and the way it was set at the end of the line. Conductivities along the line are generally high, (about 0.07 S/m), with two exceptionally very high conductivity stretches at the northern and southern ends of the line. Higher chargeabil-ities, (up to 0.02 Volt/Volt), coincide with the higher conductivities. Everywhere else chargeabilities are negligible. Determination of the target misfit for each of the inversions proved to be an important aspect of the inversion. Fitting data better always results in additional structure to the resultant model and therefore possibly better resolution. However, overfitting the DC/IP surficial data can create artificial structures and tends to push the features of the model down (figure 5.7). That behaviour, noticed also by Oldenberg and Li (1994), might seem an undesirable trait but it can also be used to help choose an appropriate misfit level if Chapter 5. Inversion And Interpretation Of DC Resistivity And IP Data 60 the model contains an object of known depth. An example is given in figure 5.7b. The inversion of line AA-east was carried out four times with target misfits IN, 2N, 3N and 4N, where N is the number of data. Available to us were only the fact that the pipe at location -290 is very close to the surface and the borehole DC resistivity measurements of well IN10 at location -510 (figure 5.8). The chosen conductivity model is shown in figure 5.6 (target misfit 3N). The pipe appears on the surface, (with extension down due to minimization of the vertical derivative), and values around location -510—the location of well IN10—fit the ones of the borehole measurements. Depth to the image of an object like the pipe by itself cannot serve as reliable factor in determining the right level of fit. Theoretical relations between the level of fit and the amount by which an object is pushed down and develops its structure should be established in order to do so. However, as one can see in figure 5.7a and other examples, it can warn the interpreter about possibly overfitting the data. Combining depth to the pipes, IN10 borehole data, comparisons of time domain and phase data inversions, fit of overlapping sections and visual appearance of the models, all had some weight in the final decision about the choice of the conductivity models. Choosing chargeability models we could not benefit even from the meager ground proofing we had about the conductivity. Well IN10 is located at an area with very low sulphide concentration (King, 1994) and comparison of the negligible chargeabilities from the inversion and the borehole has no significance. With the reservation about using the pipe to scale the target misfit in mind, it is hard to confidently choose the right model from the ones of figure 5.7a. The choice in this case is based on the comparison to the phase data inversion results. It is clear, though, that a borehole measurement in well IN 13 (location 000) could solve the problem. Another potential difficulty is the use of a 2D algorithm to invert data which might be contaminated by 3D effects. To assess this, data were taken and inverted along lines Chapter 5. Inversion And Interpretation Of DC Resistivity And IP Data 61 AA-East and A A-West. The resultant conductivity models are similar to that of Une AA, (except in the vicinity of the pipe). There is also basic similarity in the chargeabiUty models although significant differences, which wiU be discussed later, exist. In general the two dimensionality assumption appears to be valid. The detailed reference to the inversion results is deferred to the interpretation sec-tion. A method which can possibly improve results at certain sections wiU precede that discussion. 5.5.3 Poss ible I m p r o v e m e n t s O f Resu l t s T h r o u g h T h e U s e O f W e i g h t s Two of the advantages of inversion results over the pseudosections are, the model pa-rameters being real and not apparent values and having these parameters defined on a real vertical coordinate. These advantages can faciUtate quantitative interpretations but they both suffer from the smearing of the model features in the vertical direction. The smearing, evident in all the models, is the result of the minimization of the vertical model derivative. With no data constraints in that direction we get a model which smoothly varies from the real values at the surface to the reference model values at depth. The depth to any possible horizontal interface cannot be weU defined and model parameters values are affected. For a layered earth with weU defined physical horizontal lower boundary an improved model can be produced using smaUer weights for the minimization at the mesh bound-aries corresponding to the physical interface. The boundary would be emphasized and model parameters would be less affected. The conductivity model of figure 5.4b is an example of a section which might better represent geology if the conductivity indicated a discontinuity rather than a gradual change of conductivity with depth. Whether or not a sharp conductivity changes exist in the taiUngs will be discussed in the next section. Assume, for now, that the earth model is a high conductive layer Chapter 5. Inversion And Interpretation Of DC Resistivity And IP Data 62 overlying a resistive half-space. Then, the depth to the interface can be inferred by analysis of the curvature of the smearing pattern in the inversion results. Consider the conductivity values of cells at the same horizontal location to be a function of depth. This function and its first and second derivatives are plotted in figure 5.9. Figure 5.9b shows the conductivity and the derivatives, at one horizontal location, of a synthetic data inversion result. The synthetic model from which the data was produced is a 20 metres conductive layer overlying a half-space. The maximum value of the second derivative is at 18 metres. That corresponds to the upper boundary of a model cell. Thickness of the cells, down to 38 metres, is 2 metres so the lower boundary of that cell is at 20 metres—the conductivity discontinuity of the synthetic model. Experiments on synthetic models using the relevant ranges of conductivity and depth show that the mesh point next to the one of the maximum of the second derivative coin-cides with the conductivity discontinuity of the synthetic model. The same process was carried out in all the horizontal locations of the inversion result shown in figure 5.4b. One example, at location -510, is given in figure 5.9a. It appears that the inversion model and the derivatives resemble the synthetic inversion result of figure 5.9b. Assuming that also the true earth model resembles the synthetic ones and thus the smearing pattern is the same, a smooth weighting scheme was tailored around the mesh boundary corresponding to the next to the maximums of the second derivatives. The result produced by an in-version using this weighting scheme is shown at the bottom of figure 5.10. The interface is better defined, sloping gradually from north to south with conductivity of the upper layer corresponds to the values of the borehole measurements at well IN10. 5.5.4 Interpretation The conductivity values in plots of figures 5.4 and 5.6 range from maximum of 0.15 S/m to 0.005 S/m which is the reference conductivity value. The inversion result of the frequency Chapter 5. Inversion And Interpretation Of DC Resistivity And IP Data 63 domain data of line AA agrees in general with that of the time domain data inversion. With the exception of the area around the pipe at location -290, the conductivity models of lines AA, AA-East and AA-West are similar in their main features. Small variations might be attributed to real small conductivity variations as well as to slightly different fit levels. The same typical conductivity values appear in the models of line BB. Conductivities of up to 0.12 S/m extend from the northern end of the line to location -250. Lower conductivities around 0.07 S/m are prevalent from that point to location -1140, where the line enters the M l impoundment. High and variable conductivities of up to 0.15 S/m characterize that area from -1140 to -1500. At that location the water table is dipping sharply below the dam and it is evident in the conductivity section as a thickening resistive layer. Except along these last 200 metres of the line, the dry and resistive layer above the water table cannot be seen. Nor do we see in the results the expected very low conductivities of the bedrock. The lower conductivity values are those of the chosen reference model and probably the sensitivity of the data to the bedrock parameters is very low. Chargeability values in the plots of figure 5.5, (time domain data inversions), range from zero to about 0.02. The phase inversion results, (in units of mrad), are proportional to the chargeability so comparison can be done but only qualitatively. The models of the three parallel lines, (figure 5.6), show a range values from zero to 3.5 mrad. Higher values of chargeability exist at the northern end,( between 100 to -200), of all three sections but the values are somewhat different, especially south of location -200. In the time domain inversion result (figure 5.5) the chargeabilities at the northern end are relatively low, in the order of 0.001 Volt/Volt, and appear to concentrate at depth of 12 metres. From location -150 to -1150 the values are even lower and with a small exception at -870 can be taken as zero. Only the fairly varying chargeabilities between locations -1150 and -1500 reach values of 0.01 Volt/Volt or higher. The highest values in that area appear close to Chapter 5. Inversion And Interpretation Of DC Resistivity And IP Data 64 the surface. Chargeabilities along Une BB are of the order of those near the north end of Une AA—about 0.003 Volt/Volt and are higher near the surface. Site Line Station Ave. Sulphides Fe Ni cr V N13 000 3% 2190 — 8700 0.1 0.002 IN10 -510 < 1% 10 38 3200 0.07 0.0004 IN11 -1020 1,5% 77 1.5 2200 0.05 0.0001 IN12 -1409 4% 1060 641 8700 0.3 0.02 Table 5.1: Average concentrations of sulphides and maximum concentrations, (in mg/L), of oxidation products in sites along Une AA. Sources: Coggans, (1994), Cherry et al., (1990) and King, (1994). The corresponding conductivities (S/m) and chargeabiUties (Volt/Volt) of the upper metres of the inversion results. There is a clear overlap of the main features in the conductivity and chargeabiUty models. This is not unexpected since the source of the high conductivity—higher total dissolved soUds level in the pore water-is Ukely the product of the sulphide oxidation. Table 5.1 shows the maximum ionic concentrations of iron, nickel and sulphate—the main products of oxidation in the taiUngs—at weU sites IN13, IN10, IN11 and IN12. Also in the table are the average sulphide concentrations in each site. The metal ions levels in IN12 and IN13 located at northern and southern sections of the Une, are one to three orders of magnitude higher then in IN 10 and IN 11 which are in the middle section of the Une. Sulphate concentrations are about three times higher. The association of higher concentrations of sulphides with more oxidation products is obvious from the table. That adds credibiUty to the inversion models which show that the higher conductivities are associated with the higher chargeabiUties and vice versa. It also should be noted that metals ions come only from the oxidation process. Sulphate is contained also in the mill process water and its maximum concentrations are controlled by the activity of gypsum (Coggans, 1994). That may explain why the concentrations of sulphate in IN10 and Chapter 5. Inversion And Interpretation Of DC Resistivity And IP Data 65 IN11 are lower, but in the same order of the ones of IN12 and IN13, and thus, why the conductivity values are high everywhere but span a narrow range. The fact that the higher chargeabilities are associated with higher oxidation products verify, at least qualitatively, our results and implies a possible use of them. What needs to be explained is why the chargeabilities between locations -1150 and -1500 are an order of magnitude higher then those at the northern section of the line. Another fact that should be noted is that higher chargeabilities at the north end are concentrated around the depth of ten metres. In the southern end the highest chargeabilities are on or very close to the surface. In the previous chapter I suggested that the layer of magnetite encountered in IN12 could be a source of both the higher conductivities and higher chargeabilities if it extends to the south. Unfortunately a definite answer cannot be given without more ground proofing. Another question, mentioned in the previous section, concerns the vertical extent of the conductive layer and chargeable bodies. Coggans, (1992) and Cherry et al, (1990) report higher concentrations of iron ions and sulphate down to 9 metres only. As these are the most abundant ions in the tailings pore water the high conductivities should not extend deeper than that. This depth is lower then the thickness of the conductive layer shown in figure 5.10 and in contrast to the assumption used to prepare the model shown in this figure. Other known facts contradict the one mentioned above. In support of the results here, borehole measurements in IN10, (see figure 5.8), show conductivity higher than 0.055 S/m down to the deepest measurement at 15 metres. Penetrating cone conductivity measurements at the C-D and M impoundments show high conductivities in the range of 0.04 to 0.2 S/m down to the bottom of the tailings. The peak of the tritium concentrations in the pore water which correlates with recharge water from 1963-64 is at 12 metres depth. Deposition in the M impoundment, through which most of line AA goes, ceased in the year 1960. If oxidation processes started close to the deposition time Chapter 5. Inversion And Interpretation Of DC Resistivity And IP Data 66 that implies that the advection front of the oxidation products should be found deeper then 12 metres and added diffusion could take these products even further. It is not in the scope of this thesis to solve the contradiction between the first fact and all the others but my belief, mainly based upon the borehole and in-situ measurements, is that highly conductive water may have reached the bottom of the taiUngs. In conjunction with the facts that only minor groundwater flux penetrates the bedrock and that the bedrock porosity is low, that impUes that the upper conductive layer probably has a sharp bottom boundary coinciding with the bottom of the taiUngs. The last conclusion corroborates the assumption made in the previous section. That assumption faciUtated possible solution to the question of the depth of the conductive layer in relatively homogeneous section. The method cannot be appUed to the chargeabil-ity sections as the more chargeable bodies are too narrow and their chargeabiUties vary in the horizontal direction. Thus the vertical extent of these bodies cannot be resolved without any additional information. Nor can the chargeabiUty values be trusted until verified. These values are smaller then the ones measured on artificial rocks with the same sulphide content as in the taiUngs. A laboratory investigation to find the expected chargeabiUty of taiUngs with different concentrations of sulphide was supposed to be done by INCO. If carried out, it wiU enable better understanding and hopefuUy verify our results. A few more remarks should be mentioned. There are two small chargeable bodies, one at location -870 and another towards the southern end of the Une. The first is in the midst of an area of very low chargeabiUties and is confined laterally. It is probably a small metal body. The chargeabiUties at the end of the Une do not seem to be due to small body but it is also not probable to get real IP response there because the dam was built with taiUngs that had small sulphide content. The body is relatively small and should not be of importance. Guesses about what it might be are left to the reader. Chapter 5. Inversion And Interpretation Of DC Resistivity And IP Data 67 The metal pipe at location -290 which appears in the time domain inversion results and in that of Une AA-East do not appear in the inversions of Une A A and A A-West. In the pseudosections of these Unes the pipe has a much weaker response but still big enough to be assigned large errors. The pipe appears in aU the sections when a very close fit is sought. It tends to fade with increasing target misfits. It might be interesting to investigate the mechanism by which, the model features are affected by the changes of the misfit, but as the pipes are not in the focus of the work the fact they do not show up in few of the final results should not be a problem. The DCIP2D package assumes two dimensionality of the problem at hand. In order to verify that this is the case with Une AA, the two parallel Unes AA-East and AA-West were estabUshed. As was mentioned before the conductivity sections are very similar to that of Une AA. The chargeabiUty sections are quite different one from the others especially south of location -150. My feeUng is that the low values of apparent phase data in that area might be very erroneous. Many negative values measured in the first and second n-spacing dipoles affected strongly the results and they cannot be considered very reUable. What is more certain is that whatever the source of the IP response is, there is not much of it along the kilometre south of location -150. Chapter 5. Inversion And Interpretation Ot DC Resistivity And IP Data 68 Figure 5.3: Plots (a) and (b) are maps of {df* - <%ed)/ej plotted the same way df1 is plotted on a pseudosection. Plots (A) and (B) are the corresponding inversion results. Model A was produced using the error assigning scheme of equation 5.21. To produce model B the data were assigned errors which are 15 % of the average of the whole data set plus a contribution from the same exponential term of the more complicated scheme. Chapter 5. Inversion And Interpretation Of DC Resistivity And IP Data 6 9 D S 5 g 2 S 2 s ~ ~ _i erf erf — a m — e s s s a W T T T I Q o z(m) z(m) Figure 5 . 4 : Conductivity inversion results of time domain data of line AA. The colour coded values are logarithms to the base ten of the conductivity in S/m. Figure 5.5: Chargeability inversion results of time domain data of line AA. The colour coded values are chargeability in Volt/Volt. Chapter 5. Inversion And Interpretation Of DC Resistivity And IP Data 71 Figure 5.6: Conductivity and Chargeability inversion results of the phase data of lines A A , AA-East and AA-West. Colour coded logarithm to the base ten of the conductivity in S/m. Phase change chargeability values in miliradians. Chapter 5. Inversion And Interpretation Of DC Resistivity And IP Data 72 Figure 5.7: (a) Chargeability inversion results with misfits, (from top to bottom), 8N, ION and U N , where N is the number of data. Requesting lower misfit pushed the image of the pipe, along with other features, from the surface down, (b) Conductivity inversion results with misfits, (from top to bottom), IN , 2N and 4N. The chosen model has misfit 3N and is shown in figure 5.6 Chapter 5. Inversion And Interpretation Of DC Resistivity And IP Data 73 Conductivity (mS/m) 30 40 50 60 70 80 90 * — - x - — x - - Inversion - o — ° — — o - Borehole Figure 5.8: The conductivity borehole log at IN10 and the corresponding conductivity obtained from the inversion. Chapter 5. Inversion And Interpretation Of DC Resistivity And IP Data 74 Figure 5.9: Model and derivatives as function of depth. Model-solid line, first derivat-ive-dashed fine, second derivative-solid line and triangles, (a) Line A A model parameters at location -510. (b) Inversion of synthetic data forward modeled on a model of 20 metres conductive layer overlaying a resistive half-space. I_0g a ^ Figure 5.10: The model corresponding to that of figure 5.4b but produced with tailored weights for the inversion. Chapter 6 Inversion Of IP Decay Curve Data 6.1 Introduction A generalization of the basic idea of the IP method is implemented by the Complex Resistivity method (CR). Earth responses to input alternating galvanic currents over a spectrum of frequencies are measured and analysed (figure 6.1). Zonge (1972) used the method for in-situ mineral discrimination. Wynn and Zonge (1975) and Zonge and Wynn (1975) showed its application for removal of E M coupling effects, host rock response identification and mineral discrimination. Pelton et al. (1978) showed another approach, by using relaxation models, to use the CR method for mineral discrimination and E M coupling removal. The time domain analogue of the CR method is the analysis of the IP curve, (or IP waveform as it is sometimes called), as a function of relaxation (or excitation) time. Modeling of the IP decay or excitation must include a parameter with time units—a "time constant"—to characterize the waveform. The nature of the waveform and value of its characteristic time constant does not necessarily correspond to the chargeabiUty, and their analysis may provide additional information about the subsurface. In the analysis of a time domain IP curve an apparent decay time constant can be defined. In analogy with the apparent conductivity and chargeabiUty, that time constant is the one corresponding to a homogeneous earth structure response. Pelton et al. (1978) used that parameter and found it to be, in laboratory measurements and in the field, 75 Chapter 6. Inversion Of IP Decay Curve Data 76 Figure 6.1: Curve of the vector sum of the in-phase and out-of-phase components of the voltage in the complex plane. The CR method provides many points along the curve—corresponding to many frequencies. The IP method uses only two separate fre-quencies. (After Zonge and Wynn, 1975). very sensitive to grain size and sulphide concentration. The chargeability according to their findings is much less sensitive to these parameters. As grain size distribution of the tailings is very homogeneous, the time constant might be used to estimate the concentration distribution of the sulphides in the disposal area. Instruments taking continuous measurements of the IP decay curve exist but were not used in the survey. This chapter will describe an attempt to learn about the IP decay and thus sulphide concentrations in the tailings using the available regular time domain chargeability data. Chapter 6. Inversion Of IP Decay Curve Data 77 6.2 Modeling The IP Waveform With A Relaxation Model 6.2.1 Equivalent Circuits A complex conductivity cr* can be represented by real and imaginary parts cr* = cr' + icr' ,n (6.1) where all components are a function of frequency. An IP experiment can be described as a linear, time invariant system, obeying Ohm's law, with current density input J , output electric field E and complex conductivity cr* as the transfer function. A convenient method to simulate such a system is by the use of equivalent electrical circuits. If a simple circuit can represent the IP phenomenon adequately, it can be used for its modeUng. Combinations of resistors and capacitors Uke the ones in figure 6.2b were used to describe the spectra of dielectric data (Cole and Cole, 1941) and Pelton et al. (1978) showed that they can be appUed to the analysis of CR data. The circuit of figure 6.2b is a simple network which exhibits an IP relaxation be-haviour Uke the one expected from the simple rock fraction shown in figure 6.2a. The complex impedance (iu}X)~c simulates a metalUc grain-electrolyte interface, the resistor Ri represents the resistance of the electrolyte in the blocked path and the resistor R0 represents the resistance of the un-blocked path. The total complex impedance of the circuit is where c is a parameter in the range 0 < c < 1 called the frequency dependence, u is the angular frequency of the current, r is a parameter given by (6.2) (6.3) Chapter 6. Inversion Of IP Decay Curve Data 78 Figure 6.2: Schematic section of mineralized rock, its equivalent circuit and the frequency and time domains IP responses, (After Pelton et al., 1978). Chapter 6. Inversion Of IP Decay Curve Data 79 and 7] is the chargeabiUty as per Siegel's definition, given by Figure 6.2c show the ampUtude and phase of the impedance as a function of frequency on a double logarithmic plot. The ampUtude is asymptotic to RQ in the lower frequencies where only conduction through un-blocked paths takes place. At higher frequencies the ampUtude value asymptotes to the value of Ro and Ri in parallel. Between those asymptotes there is a dispersive zone where the ampUtude decreases with the frequency. The phase curve is symmetric with slopes equal to — c and c around a maximum at the point of maximum ampUtude curvature. The above description of the IP phenomenon is very simpUstic but has been found to fit CR data adequately (Madden and Cantwell, 1967). The parameter r has time units and it characterizes the decay of chargeabiUty with time along with the parameter c. In the foUowing the behaviour of the Cole-Cole model in the time domain is adopted to model the IP decay. 6.2.2 Time Domain Behaviour Of The Cole-Cole Relaxation Model The transfer function of the Cole-Cole model is the Fourier transform of its impulse response. By transformation of equation (6.2) to the time domain and integration, Pelton et al. (1979) derived an expression for the step function response. The time domain decay is the subtraction of the step function response from unity and is given by oo / i \ n ( t \nc v(t)=vS°hEw^k <6-5> where Io is the current at t = 0, (time of the current cut-off). The series is rapidly convergent in the interval 0 < t/r < lit but convergence is very slow otherwise. For larger values of t/r Pelton et al derived, in a sUghtly different way, an alternative Chapter 6. Inversion Of IP Decay Curve Data 80 (6.6) Note that for frequency dependence c = 1, equation (6.5) degenerates to V(t) = nR0I0e-t/T . (6.7) With values of c less then unity the decay is slower then exponential—in agreement with laboratory (Madden and Cantwell, 1967) and in-situ studies (Pelton et al, 1978). The shown in figure 6.2d. Equations (6.5) and (6.6) can be used to model IP decay with time assuming that excitation time is infinite and that the earth response can be described by the Cole-Cole model. If the decay time constant is short compared with the duration of excitation then the first assumption is valid. The second assumption is surely too simplistic but, perhaps, it can be applied to small fractions of the earth. To apply this weaker assumption, the formulation of the problem must be modified. 6.2.3 Intrinsic Decay Time Constant And Frequency Dependence The small section of mineralized rock in figure 6.2b might not represent an earth model adequately. But if we go to the smaller scales, that illustration of the rock and the corresponding Cole-Cole model, becomes more realistic. In order to analyse the IP decay assuming Cole—Cole behaviour, we should consider the parameters 77, c and r to be properties characteristic of each point of the earth. The definition of the chargeability fits that approach. Defining also the decay time constant r and frequency dependence c for each point in the medium, we can use equations (6.5) and (6.6) to model the response of the earth at that point. time domain response corresponding to the frequency domain response of figure 6.2c is Chapter 6. Inversion Of IP Decay Curve Data 81 6.3 Inversion 6.3.1 Inversion Methodology The steady state voltage of the equivalent circuit is Vo = Rolo so using equation (6.5) we can write m=m=v{t=0)±tm^. (,8) Vo n = 0 1(1 +nc) In chapter 5 we derived the system of equations M Vai=2ZJJiVi, J. = 1 , J V (6.9) » = i where nai is the jth apparent chargeabiUty datum and J is the data sensitivity matrix. Incorporating equation (6.8) into that expression gives VaM = E = 0) E Vn 111 \ ' * = h~,N. (6.10) The last equation constitutes an inverse problem to be solved for 7?i(£ = 0), ci and r,-given data points rja(t) along apparent chargeabiUty decay curves. 6.3.2 IPINV2D Results As Input For The Inversion Good quaUty measurements of IP waveforms would be needed if implementation of the mentioned above inversion is to be attempted. These were not available to us. The only information available about the time variation of the chargeabiUty was from the IP6 time domain data. These include averages of apparent chargeabiUty measured in 10 separate time windows during the two seconds off time of the input square waved current. The data of the first time window are the largest in magnitude and seem to have the best quality. They were the ones which have been used in the previous inversions. One can consider the data in the different time windows, at each location, as discrete points along the IP waveforms. But each curve includes only 10 points and is usually Chapter 6. Inversion Of IP Decay Curve Data 82 very noisy. However, IPINV2D inversions of the data of the time windows other then the first are possible. Figure 6.4 shows IPINV2D inversion results of the data of the third, fifth and seventh time windows corresponding to the inversion of the data of the first time window of figure 5.5d. The same features appear in all the models but the amplitude decreases as expected. As shown in figure 6.3, overlapping the models produced from data from all the time windows we get for each cell a series of chargeability values at that cell as a function of time. These values can be considered as approximations of the intrinsic chargeabilities along the decay curve at the centre of the cell. Equation (6.8), at the discrete times of the time windows, can be used for a three parameter inversion where the consecutive chargeability values rjj, j = 1,.., 10 for each cell serve as a data set to be inverted to values of r)(t = 0), c and r of the cell. These data sets are far from ideal. They result from the two step IP inversion and, possibly, contain a lot of accumulated errors. Moreover, the field data for the IPINV2D are taken during the short time interval 0.080-1.84 sec. If the decay constant is short the missing data at the early times are the most important. If the decay constant is long, then we need data from time beyond two seconds. Despite these caveats, by carrying out the IPINV2D inversions for each time window we obtain enough information to recover, in principle, the intrinsic Cole-Cole parameters for every cell. Variations in r and c may provide additional information about the state of the sulphides and we therefore pursue the inversion described. Figure 6.3: Models (a), (b) and (c) are inversion results of the third, fifth and seventh time windows data correspondingly. Figure 6.4: Schematic explanation of the the process of extracting intrinsic chargeability decay curves through IPINV2D inversions of IP time domain data. For clarity of the diagram only five time windows are considered in the figure. Chapter 6. Inversion Of IP Decay Curve Data 85 6.3.3 T h r e e P a r a m e t e r s D a m p e d Least Squares Invers ion The three parameter inverse problem is solved by a standard damped least squares procedure. The discrete form of equation (6.8) can be written as ni(t) = Vs(t,ctTtVo) = ri(tjyp) j = l,2,..,10 (6.11) where T/0 = rj(t = 0) and p = {r?o,c, T}. AS the dependency of n on the parameters of p is non-linear we look for a model perturbation Sp and iterate until a satisfactory final solution is be found. At the nth iteration we can write the approximation Vi (P ( n ) + *P) = Vi ( P H ) + £ ~H L*Pi 3 = 1. 2' - 1 0 ( 6- 1 2) »=i °"» and set t]j ( p ^ + 5p) = nf', where nf" are the data points resulting from the IPINV2D inversions. Defining the data perturbation as SV = v°b'-V (P ( n )) (6-13) the last equations can be written in the already familiar matrix form JSp = Sn (6.14) where J is the sensitivity matrix Jji = dvj/dpi. The system should be weighted by a matrix Q = diag{j^,j^, ••,"•} where tj is the assigned error of the jth datum, and the parameters non-dimensionalized by a multipli-cation by a diagonal weighting matrix W. The elements of W are adjusted such that the numerical values of the model parameters will be of the same order. The linearized system is then written as Gm = d (6.15) Chapter 6. Inversion Of IP Decay Curve Data 86 where G = QJW-\ m = WSp and d = QSn. The model perturbation 6p is found by minimizing the objective function #m) = <j>d + 9<t>m = ||Gm - d||2 + 0||m||2 (6.16) where 0 is a ridge regression parameter regularizing the system of equations which results from the minimization. The minimization is done with respect to the perturbation para-meters and 6. The solution is constructed using SVD decomposition of the matrix which results from the differentiation with respect to the perturbation parameters. That is done few times using range of of values of 6 until a satisfactory perturbation is achieved. The resultant perturbation is added to the last iteration model and the procedure is repeated until a minimum misfit model is achieved. 6.3.4 Inversion Results Issues we were concerned about while inverting the DC/IP data, emerge also while invert-ing the cells chargeabiUty data sets. Assigned errors to the data, the initial model and the values of the elements of the weighting matrix W are all parameters which influence the final result. To my best knowledge, Cole-Cole parameters of mine taiUngs have never been investigated before so I had no idea of their expected values. Thus an exploration of the model space by adjustment of the inversion parameters mentioned above was carried out as the first step. The assigned errors were fixed, (as 5% of the datum plus 5% of the average of its data set), and the ranges of possible values of c and reasonable final misfits were investigated. With low values of c—less then 0.5—convergence could not be achieved and values of n0 tended to be very big. With high values of c good fit could be achieved but the values of rjQ were small. We expect T/O to be of the order of the datum of the early time window, (usually 0.001-0.02 V/V). Visually inspecting the data curves it seemed that it should Chapter 6. Inversion Of IP Decay Curve Data 87 be about 1.5-4 times higher. That limits the range of c values to the range 0.5-0.95. To get reasonable fit the corresponding r values must be 0.05-1.0 sec. Having an idea of the ranges of values of r and c the weighting matrix was fixed as W = diag{10.0,1.0,1.0} and each inspected data set was inverted to the lowest possible level of fit. Many data sets of cells in the areas of high chargeability were inverted. Four examples of observed data sets and the ones predicted by the inversion results are given in figure 6.5. Values of r, in most cases, were 0.2-0.4 sec. No distinction can be made by that parameter between the results in the different areas. On the other hand, the value of c seemed to be higher, (0.65-0.75), at locations between -1250 and -1450 than in the area between 0 and -150 where c value was in the range 0.90-0.95. 6.3.5 Interpretation Pelton et al. (1978) inverted results, reported by Grisseman (1971), of CR measurements on artificial rocks composed of cement, quartz and pyrite. They found the decay time constant of rocks with 2-5 % sulphide content and 0.2-0.3 mm grain size to be 2 x 10 - 4 sec to 5 x 10 - 4 sec correspondingly. The time constant of artificial rocks composed of smaller grains is expected to be even shorter. These results are totally different then the ones I got from the tailings which contain silty sand with w 0.01 mm grain size. Also the variation of the value of c and its range are not in agreement with the typical values Pelton et al. obtained in inversion of in-situ measurements. They found c to be ranging from 0.2 to 0.5 and the typical value is 0.25. These discrepancies are bothering but it must be remembered that the mineral deposits on which the in-situ measurements were taken are environments different than the tailings and might have different parameters. The estimated concentrations of sulphides along line AA are low and range from 3% to 4%. Assuming that r behaves like the time constant as per the definition of Pelton et al. we expect it to be sensitive to sulphide concentration, but probably it is not sensitive Chapter 6. Inversion Of IP Decay Curve Data 88 enough to detect the small concentration variations we encountered. The values of c above 0.9 at the area near location 000 are more typical to those of E M coupling. On the other hand the fact that the pseudosections and inversion results of the different time windows show similar features implies that the coupling is negligible. Also the variation of the c value are intriguing. I cannot say too much about the physical meaning of the values I obtain but it is clear that these values are different in the area at the north end of the line and at the south end. The unanswered question why the chargeability values at the southern end are an order of magnitude larger may be connected to and might find some clues in the differences of c values. Figure Cell location Cell depth Vo r c % Sulphides a -70 to -80 2-4 metre 0.0019 0.2855 0.9123 3 b -20 to -30 10-12 metre 0.0027 0.2880 0.9345 3 c -1350 to -1360 2-4 metre 0.0259 0.3396 0.6848 4 d -1270 to -1280 8-10 metre 0.0244 0.3124 0.6665 4 Table 6.1: Parameters corresponding to figure 6.5 Chapter 6. Inversion Of IP Decay Curve Data 89 IP DECAY CURVE IP DECAY CURVE < 0.0008 0.2 0.4 0. IP DECAY CURVE IP DECAY CURVE Figure 6.5: Observed, (circles), and predicted, (solid line and crosses), chargeabilities as function of time. The corresponding model parameters and cell locations of each plot are given in table 6.1 . C h a p t e r 7 Conc lus ions 7.1 D C Res i s t i v i ty The DC resistivity data inversion have produced a conductivity structure along line AA. The inversions of the different data sets along the line are in agreement. The two dimen-sionality, which the inversion routine assumes, appears to be valid except in the vicinity of buried metal bodies, (which are not in the focus of the research). The conductivity values of the upper conductive tailings are as expected and are in agreement with meas-urements in one borehole on the line and a few penetrating cone measurements in its vicinity. More borehole or in-situ measurements, reaching the bottom of the tailings, are needed to determine definitely the depth of the conductive layer. The additional ground proofing information can be incorporated in future inversions to minimize the inherent non-uniqueness of the results. The inversion results provide detailed conductivity structure along line AA. With a few more survey lines covered, a semi three-dimensional conductivity structure can be built to serve as a basis for a quantitative mapping of the TDS concentrations. To accom-plish that, correlation coefficients would have to be established between the conductivity and the TDS level. It should be feasible as the porosity and grain size throughout the tailings are relatively uniform and also interference from clays is expected to be minor. The DC resistivity data and the models resulting from their inversion are very reli-able. Time domain and frequency domain systems produced similar results. Also, data 90 Chapter 7. Conclusions 91 collected with the 20 metre spacing array seem to be adequate to construct accurate models. The corelation between them and the chargeability models can be exploited while interpreting the inversion results of the more problematic IP data. 7.2 Induced Polarization ChargeabiUty models were constructed by inverting both time domain and phase data sets. The correlation between chargeabiUty values and known concentrations of sulphides is manifested in all the results. Estimation of their accuracy is difficult as no ground proofing is available. Discrepencies between the models of the parallel Unes exist, mainly along the areas of known very low sulphides concentrations. It seems, though, that the two dimensional assumption is vaUd. The generally small chargeabiUties imply, in agreement with the current knowledge, that at least along Une AA, sulphide concentrations are relatively low. But as they are enough to produce significant amounts of unwanted oxidation products—evident in the conductivity models and in-situ measurements—there is a need to quantify their total amounts and distribution. Relationship between chargeabiUty values and sulphide concentrations must be estabUshed in order to do so. Estimation of the AMD production potential wiU require a further investigation of the oxidation activation processes. The accuracy of the chargeabiUty models and whether their accuracy will be adequate to serve as a basis to these calculations is yet to be determined. It seems, though, that in terms of prediction of AMD in the taiUngs they are superior to any other available source. The smearing of features—a major problem which might be difficult to solve for the IP inversions—affects mainly the chargeabiUty values at depth. The sulphide distribution in the few metres just below the water table is the more important in the AMD research and thus the smearing problem should not be overemphasized. Chapter 7. Conclusions 92 Unlike the conductivity models chargeabilities vary significantly along the line. Thus dense coverage provided by shorter spacing of the array is more essential to the construc-tion of chargeability models than for the conductivity models. The inversion of the 20 metre spacing phase data from line AA seems slightly inferior to those of the 10 metre spacing, time domain data inversions. The inversions of phase and time domain data taken on line BB, (both with 20 m spacing array), are very similar. If another survey is to be conducted on the tailings the array spacing would be an important factor to deter-mine. A definite conclusion as to whether less then 20 metres spacing is really necessary cannot be drawn from our results. It seems however that much larger a-spacing should be avoided if a detailed result is sought. The IP decay curve inversion was a very experimental endeavour and time constraints limited its development to a very crude stage. The values obtained for the model para-meters must be considered with care. One use of the procedure—mapping 7?0—which would have yielded a true chargeability sections cannot be accomplished without a robust method to invert the data sets from the different cells to the same level of fit. The values of the decay time constant, expected to be of use in mapping the sulphide concentrations, do not seem to reveal any additional contributions beyond what is obtained by the IPINV2D inversion. Directly relevant to the IP interpretation is only the fact that the frequency dependence values in one area significantly differ from the ones in another. If the reason for the difference can be determined then this parameter might be used to identify it in other cases. Fitting the IP decay curves using the simplest Cole-Cole relaxation to model the phenomenon is encouraging. Improvements can be achieved with the use of real IP waveform data and more comprehensive inversion routines. Appendix A Resistivity/Conductivity Conversion Table In this thesis electric conductivity a in units of Siemens per metre (S/m) was considered. As its values can span many orders of magnitude the inverse problem was solved and results were given in terms of logarithms to the base ten of the conductivity. Apparently, many people tend to think in terms of the resistivity p in Ohm-m which is the reciprocal of S/m. A units conversion table is given here as a quick reference for any confused readers. p (Ohm-m) a (S/m) 1.0 2.0 3.0 4.0 5.0 10.0 15.0 20.0 40.0 100.0 200.0 300.0 400.0 500.0 1000.0 2000.0 4000.0 1.000 0.500 0.333 0.250 0.200 0.100 6.66 ) 5.00 > 2.50 > 1.00 > 5.00 > 3.33 > 2.50 > 2.00 > 1.00 > 5.00 > 2.50 > x IO"2 x IO"2 x 10"2 x IO"2 x 10"3 x 10"3 x 10"3 x 10"3 x 10"3 x IO"4 x 10"4 0:000 -0.301 -0.477 -0.602 -0.699 -1.000 -1.176 -1.301 -1.602 -2.000 -2.301 -2.477 -2.602 -2.699 -3.000 -3.301 -3.602 93 References [1] Archie, G.E., 1942, The electrical resistivity log as an aid in determining some reservoir characteristics, Transactions of American Institute of Mineral Metallurgy Engineering, v. 146, 54-62. [2] Cherry, J.A., Robertson,W.D., Blowes, D.W., and Coggans, C.J., 1990, Hydro-geology and geochemistry of the INCO Ltd., Copper Cliff mine tailings impound-ment, Waterloo Center for Groundwater reaserch, University of Waterloo, Waterloo, Canada. [3] Coggans, C.J., 1992, Hydrology and geochemistry of the INCO Ltd., Copper Cliff, Ontario, mine tailings impoundments, M.Sc. Thesis, University of Waterloo, Water-loo, Canada. [4] Cole, K.S., and Cole, R.H., 1941, Dispersion and absorption in dielectrics, J. Chem. Phys, v. 9, 341. [5] Grisseman, C , 1971, Examination of the frequency-dependent conductivity of ore-containing rock on artificial models, Scientific rep. no. 2, Electronics Laboratory, University of Insbrook, Austria. [6] De Vos, K. J., 1992, Earth resistivity investigation of a plume of tailings-derived wa-ter, Copper Cliff, Ontario, B.Sc. thesis, University of Waterloo, Waterloo, Canada. [7] Dey, A., and Morrison H.F., 1979, Resistivity modeling for arbitrary shaped three-dimensional structures, Geophysics, v. 44, 753-780. 94 References 95 [8] Ebraheem, M.W., Bayless, E.R., and Krothe, N.C., 1990, A study of acid mine drainage using earth resistivity measurements, Ground Water, v. 28, 361-368. [9] Hallof, P.C., 1964, A comparison of the various parameters employed in the variable frequency induced polarization method, Geophysics, v. 29, 425-434. [10] Keller, G.V., and Frischnecht, F.C, 1966, Electrical methods in geophysical prospect-ing, Pergamon Press, London. [11] King, A, 1994, INCO Exploration And Technical Services Inc., Personal comunuca-tion. [12] Madden, T.R., and Cantwell, T., 1967, Induced polarization, a review, in Mining geophysics, v. 2, Theory, Tulsa, SEG, 373-400 [13] Nabighian, M.N., and Elliot, C.L., 1976, Negative induced-plarization effects from layered media, Geophysics, v. 41, 1236-1255. [14] Oldenburg, D.W. and Li, Yaoguo, 1994, Inversion of induced polarization data, Geophysics, v. 59, 1327-1341. [15] Oldenburg, D.W., McGillivray, P.R, and Ellis, R.G., 1993, Generalized subspace methods for large scale inverse problems, Geophys. J. Int., 114, 12-20. [16] Pelton, W.H., Smith, B.D., and Sill, W.R., 1979, Interpretation of complex resistiv-ity and dielectric data, part 1, Phoenix Geophysics Ltd., Toronto, Canada. [17] Pelton, W.H., Ward, S.H., Hallof, P.G., Sill, W.R, and Nelson, P.H., 1978, Min-eral discrimination and removal of EM coupling with multifrequency resistivity and dielectric data, Geophysics, v. 43, 588-609 References 96 [18] Pehem, P.E, 1981, Geophysical survey of the Nordic contamination plume, B.Sc. Thesis, University of Waterloo, Waterloo, Canada. [19] Siegel, H.O., 1959, Mathematical formulation and type curves for induced polariza-tion, Geophysics, v. 24, 547-565. [20] Sumner, J.S., 1976, Principles of induced polarization for geophysical exploration, Elsevier, Amsterdam. [21] Telford, W.M., Geldart, L.P, and Sheriff, R.E., Applied Geophysics, Cambridge Uni-versity Press, Cambridge. [22] Wynn, J.C, and Zonge, K.L., 1975, EM coupling, its intrinsic value, its removal and the cultural coupling problem, Geophysics, v. 40, 831-851. [23] Zonge, K.L., and Wynn, J . C , 1975, Recent advances and applications in complex resistivity measurements, Geophysics, v. 40, 851-864. 


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items