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. 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 , E = -V<£ (4.3) and thus we have J = -o-V . (4.4) Using Ohm's law in the charge conservation equation V • J = 0 (4.5) we get V • (crV) = 0 (4.6) and in the presence of a current source we have V • (crV) = -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 , = -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 $n — {t = oo) is the secondary field ^. = ^ - ^ = :i[*JL_ (4.13) cr (1 - n) and the ratio of 4>a/t) is simply rj. Thus we can define apparent chargeabiUty—the chargeabiUty of a homogeneous earth corresponding to pair of , 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 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 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 t are integrated over a time window and normalized by r,. That gives data expressed in time units which are proportional to the chargeability. Other systems average , 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) = -I6(r-r.) ( 5 . 1 0 ) and the appropriate boundary conditions : d/dn = 0 at the surface and —• 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 {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 £ £ dln

• 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) = d + 9m = ||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.