Newer results of Monte Carlo inversion of IP data in water base protection and ore exploration

The paper presents the latest results of Monte Carlo inversion of IP data in the areas of water base protection and ore exploration. The method of determining the time constant spectrum using Monte Carlo inversion and the parameters characterizing the degree of environmental contaminations and ore deposition are presented. Among the field applications in Hungary, the investigation of the ionic pollution of the Ráckeve water base and the characterization of the Felsőtelekes waste dump of the Rudabánya iron ore mine are presented. In addition to these, the paper presents the investigation of ore deposition in the area of a gold mine in Mongolia using the Monte Carlo inversion of multi-electrode pole-dipole IP method.


Introduction
The induced polarization (IP) method is a geophysical tool discovered by Conrad Schlumberger in a mining region of France. It has primarily been used for mineral and ore exploration (Bleil 1953;Pelton et al. 1978;Wait 1959;Keller and Frischknecht 1966;Sumner 1976;Fink et al. 1990). Later detailed studies of measurable IP effects associated with nonmetallic minerals (e.g., Olhoeft 1985;Vanhala and Soininen 1995) have identified numerous potential engineering and environmental applications. Because of the fact, that mineralogical circumstances have also a strong effect on the IP response, it is possible to use IP as a tool for lithological discrimination and hydraulic characterization. In the last decades environmental applications of the IP method were carried out for the detection of organic or chemical contaminants. Many advances have been made in the science and practice of the IP method both in the theory of the phenomena and in the instrumentation, including multichannel receivers. Improved signal processing and inversion methods were developed making it possible to obtain 2D or 3D images of the polarizability of the subsurface (Park and Van 1991;Li and Oldenburg 2000).
The method of time domain IP was used for the laboratory measurement of numerous core and rock samples in the Geophysics Department of the University of Miskolc (Takács 1983). One group of the specimens was collected from an area of polymetallic ore deposits, another one originated from coal seam complexes, while the third included samples from shaly-sand beds. The research comprised the analogue measurement of the overall length of the IP decay curve. 120 samples were taken from the section of decay as a combined linear-logarithmic time series. This sampling rate allowed for the stable spectral and time constant distribution analyses of IP signals. The time constant spectrum w(τ) was used for the characterization of the decay curve the determination of which was made by the TAU transformation of the observed apparent polarizability curve η a (t) introduced by Turai (1981).
The strongly monotonously decreasing apparent polarizability curve observable in the time domain can be described by the following first kind Fredholm type integral equation where t is the time passed by after the switch-off of the electric current excitation, τ is time constant of the polarization, w(τ) is time constant spectrum to which the normalization assumption is demanded.
The time constant spectrum contains all the spectral information on polarization observable by the IP measurement (Turai 1985), hence its determination can be useful even when petrophysical model for the description of the IP phenomenon is not available.
The inverse of Eq. (1) was named TAU transformation by Turai (1981) which allows for the determination of the time constant spectra correspondent with the observed apparent polarizability curves. There is no analytic solution for the TAU transformation. Turai (1985) elaborated polynomial interpolation and Fourier transformationbased methods for the approximate solution for Eq. (3). General solution can be given by a series expansion-based inversion method (Turai and Dobróka 2001). The series expansionbased inversion method was previously applied to solve the non-linear well-logging inverse problem by Simulated Annealing as a global optimization method (Szabó 2004). An easier way for the computer implementation of the inversion method is the use of the Monte Carlo procedure (Turai et al. 2010). (1) 1 3

TAU transformation solved by series expansion-based inversion
The time constant spectrum w(τ) is a continuous real-valued function, which can be estimated only with a finite accuracy from finite number of measurement data. To describe it with finite number of data, in the first step, it must be discretized. Let the time constant spectrum be written in the form of series expansion where Φ q is the q-th basis function, B q is the q-th expansion coefficient and Q is the number of the expansion coefficients. Since the basis functions are a priori given, the extraction of the time constant spectrum reduces to the determination of unknown expansion coefficients.
If the TAU transformation is defined as an inverse problem, the vector of series expansion coefficients ⇀ B is the unknown model vector and Eq. (1) is the response function used for solving the forward problem, which, combined with Eq. (4), provides a connection at measured time t k By introducing the following notation the calculated data can be generated by the expression which in matrix form is Consider the introduction of the deviation vector of the measured ( ⃗ obs ) and calculated apparent polarizability data TAU transformation can be solved as an inverse problem. If the discretization of the time constant spectrum is performed properly, it can be achieved that the number of measurement data will be more than that of the expansion coefficients. In that case, the inverse problem is overdetermined the solution can be derived by the minimization of calc = .
the L2-norm of Eq. (9). This is the least squares method, which leads to the classical normal equation (in case of linear inversion) By solving the above equation the time constant spectrum w( ) can be calculated by substituting the expansion coefficients into Eq. (4). As an analogous example, the series expansion-based inversion methodology was used for estimating the lateral variations of layer-boundaries and resistivities of shallow formations in geo-environmental site characterization (Gyulai and Szabó 2014;Gyulai et al. 2014).
The choice of the type of basis function Φ q depends on the given problem. In the general case (from the numerical point of view), series expansion using Legendre polynomials is favorable. Interval-wise constant functions have also an important role in environmental geophysical applications, because the domain of time constants is relatively narrow, which can be efficiently split into four intervals, only. In this case, the basis function in the range of 2Δ around q can be chosen as which allows for the analytic calculation of matrix elements in Eq. (6). (In a typical case of equal intervals Δ = max ∕8 .) The advantage of this type of discretization is that the expansion coefficient directly represents the constant value of time constant spectrum valid in the given interval. In that case, Eq. (5) can be calculated easily as In order to describe line spectra of the time-constant Dirac-delta functions can be used. In this special case, Eq. (4) takes the form as where q is the time constant belonging to the q-th polarization mechanism. The delta function ( − q ) is sometimes thought of as an infinitely high, infinitely thin spike at q , with total area one under the spike. From a purely mathematical viewpoint, the Dirac delta is not strictly a function, it makes sense as a mathematical object only when it appears inside an integral. In our case Eq. (5) expressing the forward problem takes the form which because of the integration property of Dirac-delta simplifies to The previous equation clearly shows that the meaning of expansion coefficient B q is the amplitude belonging to time constant q . The deviation vector of observed and predicted data and its L2-norm calculated similarly as earlier led to Eqs. (10, 11). The inverse problem can be solved by several linear and global optimization techniques. By choosing the latter, the application of Monte Carlo method is presented in the paper.

Inversion-based solution of TAU transformation using Monte Carlo procedure
The curve η(t) is measured in discrete moments of time by the induced polarization method. The spectral amplitude belonging to time constant τ contains information on the rock environment. The basic relationship (2) for the time constant spectrum after discretization will be as follows: where w n is the amplitude of the n-th discrete component and N is the number of components. The number of components, i.e. types of polarization effects, should be determined before the Monte Carlo simulation of parameters (τ, w). In this study, maximum ten components are chosen for the inversion procedure (10 ≥ N), in which the forward problem takes the form of Dirac-delta function-based series expansion (expansion coefficients correspond to the spectral amplitudes). After the generation of pair N (τ, w), the goodness of fit between the calculated and observed polarizability curves is measured by the L 2 -norm of individual deviations normalized to the measured data where η obs is the measured polarizability curve, η calc is that of calculated by the Monte Carlo method, N t is the number of moments of time.
The measurement in Telkibánya (Turai et al. 2010) and Ráckeve (Turai and Nádasi 2020) area was carried out by an IRIS SYSCAL instrument using Wenner-array with 72 electrodes, which produced IP data measured in 744 reference points at N t = 20 moments of time. The fitting should be studied in each spatial coordinate. The algorithm checks 100,000 randomly generated curves to seek the calculated curve, which fits the best to the measured one. If no one is found under the error level of 1%, then the acceptance error is increased by 1% and 100,000 new curves are generated with randomly chosen parameters. This loop is continued until a decay curve is calculated, which meets the above requirements. At the end of the procedure, we accept the values w and τ of the fittest calculated curve and assign them to every reference points. The acceptance error level was 1% in both cases of synthetic and field data. Using the estimated time constants and their amplitudes, geoelectric profiles are plotted showing the effects of polarization components. The repetition number of 100,000 can be changed arbitrarily. The further increase of this number will not give much better result. On the other hand, smaller number of trying may degrade the result of inversion. It is confirmed also by the experience that the average error calculated on the data set is always approximately the same, though the profiles constructed from the calculated data may show some difference. It was assumed that the measured decay curve η(t) was strongly monotonic (decreasing) function. Additional to it, the curve was not allowed to slip over into the negative domain. Both assumptions are rather rigorous, which may decrease the noisiness of the data set, but at the same time, it can also cause the significant loss of data by eliminating those measuring points which are disagreeable with the above conditions. The procedure was proved to be feasible to practical applications. However, later we improved the numeric Monte Carlo software as follows. First, random numbers were generated I times in the interval (τ max − τ min ), where I is the number of the random generated components. Let these random numbers be the values of the time constants (τ i,rnd ). Secondly for each τ i,rnd value, a M number of random numbers were generated in the interval (w max − 0). The M was usually five thousand. These random numbers were considered as the w i,rnd amplitudes of the τ i,rnd time-constant component. Third, the previous τ i,rnd and w i,rnd random generations were repeated N times. The N was five thousand too. In this way, the algorithm performed (I × N × M) random generation. For each pair (τ i,rnd , w i,rnd ), the polarizability curve was obtained by random generation. The one with the smallest data distance (minimal L-norm) from the measured polarizability curve was selected. With the improved algorithm, the same or higher accuracy can be achieved with a smaller number of random generations. We have already used the advanced Monte Carlo algorithm for the inversion of Mongolian pole-dipole IP measurements.

Estimation of degree and type of contamination
The method of TAU transformation was developed for the laboratory test of core samples with dispersed ore mineralization (Turai 1981;Takács 1983). The first field test using IP data was conducted for delineating a municipal waste site near Offheim (Germany) and the assessment of the degree of contamination (Turai et al. 1992). This application was followed by several Hungarian case studies of ore exploration, environmental protection and protection of the water base.
Based on time constant spectra reconstructed by the processing of field data different types of polarizations can be distinguished, i.e. filtration, membrane, redox and metallic polarizations (in the order of increased values of time constants). The distribution of time constant values is used to estimate the type of polarization. Their geological reasons are summarized in Turai (2011) and Table 1. Filtration polarization is encountered in porous soils/rocks with highly conductive pore fluid, membrane polarization is caused by the presence of dispersed clay minerals. Redox polarization is developed by oxidation and reduction processes induced by contaminations, while metallic polarization appears because of metallic soil constituents.
Several case studies (more than fifty) revealed that the components represented by small time constants (τ < 1 s) can be associated with filtration and membrane polarizations, which are not risky from the point of view of environmental pollution. Redox and metallic polarizations (in metal solutions with free electrons) made by dangerous electrochemical and metallic effects and dielectric polarizations caused by non-conducting dispersed materials (finegrained plastics, dispersed aromatic hydrocarbon derivatives, insulating organic materials) are recorded in the domain of greater time constants (τ > 1 s). Metallic and dielectric polarizations give the longest remaining and highest amplitude signal. In all cases when the conductivity between both sides of a boundary in a sequence of strata shows big difference the given surface charges up. The rate of charging (i.e. storage capacity) is directly proportional to the size of the surface, thus the longest remaining and highest amplitude IP signals with large time constants can be measured in case of dispersed (fine-grained, steamy, blobby) contaminants and ore formations. The conductivity of dispersed metallic contaminants (metals, metallic salts, coal and graphite, conductive organic materials) affecting metallic polarization can be one or two orders of magnitude higher than that of the surrounding deposits containing dissociation ions in their pore spaces. On the contrary, in case of dielectric polarization, the conductivity of dispersed non-conductive contaminants (fine-grained plastics, dispersed aromatic hydrocarbon derivatives, insulating organic materials) is orders of magnitude lower than that of the surrounding deposits with ionic conduction.
Based on the above experiences, Turai (2004) introduced the amplitude spectrum weighted by the time constant (WAV-Weighted Amplitude Value) to quantify the degree of contamination and ore formation, which increases the weight of dangerous contaminations (connected with redox, metallic and dielectric polarizations), while it reduces the effect of less hazardous contaminations (connected with filtration and membrane polarizations). The WAV parameter can also be called polarization contamination using a simplified name.
The reprocessing of earlier field measurements (Turai 2012a, b) and the results of new surveys confirmed that the use of the WAV value was feasible in every measurement area. The spatial distributions of dangerous contaminations can be delineated by the WAV values. According to our knowledge, the WAV values indicates the following types of contaminations or ore formations: WAV > 20 percent-very highly contaminated area or ore formation is very high, 20 percent > WAV > 10 percent-highly contaminated area or ore formation is high, 10 percent > WAV > 5 percent-medium contaminated area or ore formation is medium, 5 percent > WAV > 2 percent-weakly contaminated area or ore formation is small, WAV < 2 percent-clean (uncontaminated or no ore formation) area. Turai et al. (2008) introduced the notion of corrected conductivity as the multiplication of time constant spectrum and specific conductivity (σ) Bad or non-conductive materials (disperse plastic materials, aromatic disperse hydrocarbons, electrically non-conductive organic materials) in porous rocks with conductive fluid 1 3 The corrected conductivity range above 100 mS/m derived from field measurements allows for the separation of areas with strong ionic contamination.

Results of field data processing
In the framework of the "clean drinking water" project, we used IP measurements in the examination of the water base of the Danube near Ráckeve (Turai and Nádasi 2020). The aim of the studies was to designate areas suitable for the abstraction of water for human consumption and to exclude contaminated parts of the water base. The IP measurements performed in the area of Ráckeve were processed by Monte Carlo inversion, assuming a single average polarization component. After inversion of the multielectrode resistivity and IP measurements in the Wenner electrode array, corrected apparent conductivity profiles were calculated. Figure 1 shows the corrected apparent conductivity distribution calculated under the Rac 1 profile. The corrected apparent conductivity distribution under the Rac 2 profile is shown in Fig. 2. The values of the time constants in this area were less than 0.8 s, which showed that the polarization was mainly caused by membrane and filtration effects.
In both figures, the blue contour line of 10 mS/m indicates the boundary between the ionically contaminated and uncontaminated pore water. Below this value, the water base is recommended to open for water abstraction. Due to the ion concentration-increasing effect of the dispersed clay content in the corrected conductivity values, clay formation also appears as a polarization impurity-increasing effect in the range greater than 40 mS/m (red contour line). In these parts, the dispersed clay content reduces porosity and permeability, as well as binds ionic impurities, so opening the water base for water extraction is not recommended here. In the areas between the blue and red contour lines, water extraction is only recommended with continuous hydrochemical cleaning.  Turai and Nádasi 2020, in Hungarian) IP surveys were made near Rudabánya in 2011 for the assessment of rock dump. IP soundings using Schlumberger array were carried out in the dumping yard of Felsőtelekes. The sounding points were located at the corner points of a 100 × 100 m grid. The grid center was at point FT103, and its coordinates were EOVX = 340,468 and EOVY = 767,513. The inversion processing was performed by the Monte Carlo method, using a single average polarization component. The vertical distribution of WAV values observed at point FT103 against half the current electrode spacing is illustrated in Fig. 3. It must be mentioned that the approximate true depth of reference equals AB/5, so the approximate depth range is between 0.8 and 20 m. It is concluded from the figure that the shallowest part of the section is weakly contaminated, which is followed by a strongly contaminated interval between 1 and 4 m, then the degree of contamination varies from medium to high levels.
The lateral variation of contamination for depth level of 2 m under the colliery dump is plotted in Fig. 4. It is seen from the figure that the highest contamination occurs in the middle of the area, then the degree of contamination decreases gradually towards to the sides of the area. The vertical variation of the contamination level can be inferred from Figs. 4 and 5. The spatial distribution of polarization contamination in depth of 2 m can be seen in Fig. 4, while Fig. 5 refers to the depth of 4 m. It is concluded from the figures that the center of area of the strong contamination shifts growingly to the northern boundary of the investigated area with increasing depth and the extension of contamination increases, too. The strong contamination can be associated with the effect of redox polarization in the depth of 10 m (Fig. 6), which implies that the rock material of the dump has significant amount of metal and metallic salt content even in this depth.  By Monte Carlo inversion of IP data measured in a pole-dipole electrode array in a Mongolian gold mine, we have shown that ore formation is multi-component (polymetallic). The measurement area is located in the region of the Yamaat gold mine, 110 km from Ulaanbaatar, where detailed geological and geophysical surveys were carried out between 2011 and 2017. The number of polarization (optimal) components (I opt ) was determined. The number of optimal components (I opt ) can be determined by performing a random search with one-component, two-component, three-component,…, and tencomponent curves at some measurement points of the measured section and accepting the optimal number of components for the smallest data distance (D). The determination of the optimal component number for the Mongolian IP measurement is shown in Fig. 7, where the minimum of data distance is at 3 components. From the time constants (τ i ) and amplitudes (w i ) of the three components, an average WAV parameter was calculated according to the following equation: One vertical average WAV section calculated by Monte Carlo inversion is shown in Fig. 8, where the polymetallic ore deposition is mainly below the 200 m depth level, but can also be found near the surface at lower concentrations.

Conclusions
The results of the presented case studies confirm that the extension, degree and type of polarization contamination can be estimated by the time constants and the new parameters are introduced for the characterization of polarization contamination of the media. The method was successfully applied in the investigation of the ionic contamination of the Ráckeve water base. IP measurements made in the area of the Felsőtelekes rock dump are applicable to the classification of the material of debris. By means of the WAV values, the lateral and vertical distributions of contaminations can be detected and the reason of the polarization contamination can be given for the site. The strong level of polarization contamination in the investigated waste-rock pile can be explained by the metal (metallic polarization) and metallic salt (redox polarization) contents of rock debris. IP measurements in the area of gold ore mine in Mongolia by Monte Carlo inversion have shown that ore formation has three components and the concentration of ore formation increases below the 200 m depth level.