Applying Electrical Resistivity Tomography in Ornamental Stone Mining: Challenges and Solutions

: In this study, the use of electrical resistivity tomography (ERT) as a tool to guide ornamental stone extraction is investigated. ERT is not conventionally used in highly resistive environments, such as on rock faces, due to the high contact resistances that can impede current injection. Here, the challenges of conducting ERT in such environments are discussed and possible solutions suggested. For this, an example of the application of ERT in a deep and narrow marble quarry is used. The marble deposit is affected by fracturing and karstiﬁcation. Due to the nature of these features, they present a signiﬁcant resistivity contrast to the background resistivity of the marble and thus excellent targets to test the application of ERT. Their location was mapped using ﬁeld observations and complementary ground penetrating radar data. By using an appropriate sensor deployment, a suitable resistivity meter, and advanced data processing routines, the derived 3D resistivity model is in good agreement with the independent observations. This shows that despite the challenges, ERT can be used as a non-invasive tool to obtain information on the stone properties prior to extraction. This will help in guiding quarry operations and will allow for a targeted, safe, and efﬁcient extraction of high quality stone, thereby increasing sustainability and economical competitiveness.


Introduction
Ornamental stone is a major industry with a value chain in excess of $47 billion per year [1]. Despite the value of the sector, exploration approaches are still heavily reliant on either surface observations or intrusive investigations. Intrusive approaches, such as drilling, despite providing a direct assessment, are limited in their volumetric representation of the rock mass, and assessing discontinuities and quality variations away from the drill hole remains difficult [2,3]. This is a serious issue because if stone quality variations cannot be adequately identified ahead of extraction, effective quarry planning cannot be undertaken, thereby impacting the economic viability of the mining operation. A related challenge faced by the industry is the large amount of waste that is produced by these operations. For example, the European mining and quarrying industry produces an annual waste volume of about 700 million tons, equal to 28% of the entire European waste production [4]. Waste production can be exacerbated by inadequate exploration in which areas of uneconomic material are not effectively characterised in advance of quarrying. To improve ornamental stone exploration and quarry design, a number of novel approaches are beginning to be considered. This includes the use of non-invasive geophysical imaging techniques [5], which have the potential to map subsurface property

Study Site
The study site, an operational narrow and deep open-pit dolomite quarry, is located in the Lombardy region of Northern Italy ( Figure 1). Geologically, it is situated within shallow marine sediments of the Corna formation of the Lower Jurassic, comprising grey, brown to whitish limestone and dolomite facies [28]. At the study site, the formation is heavily brecciated, giving it its commercial name "breccia aurora", which is marketed as a beige marble with unique and intense red to orange veins. At the time of the survey, the active extraction area of the quarry was divided into an upper and a lower floor.
The formation, and hence the quarry operation, are affected by faulting, which is opening potential moisture pathways into the dolomite formation, which in turn causes karstification ( Figure 2). This causes dissolution features of different sizes and shapes, which can be observed in the quarry faces and floor. The lower extraction floor showed a partly collapsed cave roof, which allowed visual inspection of the void. This void, located off-centre on the lower floor was air-filled in the upper 4.2 m and water-filled in the lower 2.8 m. The entrances to the cave were comparably small ( 1.8 × 0.5 m, and 0.5 × 0.5 m), but the cave seemed to expand significantly at about 4 m depth from the surface. A fault, developing into a clay-covered dissolution feature, was also observed east of the upper floor (outlined in red in Figure 2). Observations on the quarry faces and floors suggested that these features of the upper and lower floor were linked through a common fault. These dissolution features pose potentially hazardous features for the quarry operation, as they could cause collapse and thus fatalities and damage to excavators and other machinery. These features could also be economically detrimental. Quarried blocks affected by such features may show increased weathering, and thus lower stone quality, and a smaller exploitable area, thus only allowing for small scale tiles to be produced from them.

Methodology
The resistivity of rocks and soils depends on several parameters, such as porosity, saturation, pore water conductivity, temperature, and clay content [29,30]. In clay-free rocks, electrical current flows through electrolytic pathways in the pore fluid, as the rock matrix is usually assumed to be electrically insulating. Hence, increasing pore space and saturation lowers the electrical resistivity of a medium, while decreasing pore space and saturation increases its resistivity. An additional electrical conductance mechanism is added if clay minerals are present. Due to their electrical double layer, next to electrolytic conduction, surface conduction is also possible, thereby significantly lowering the electrical resistivity of a clayey rock. Ornamental stones, such as granites, limestones, or dolomites, usually have a very small porosity and small amounts of clay and thus tend to have resistivities well above 10 kΩ·m. However, alteration of this rock mass due to weathering, fracturing, or dissolution can change these properties significantly [31][32][33].
The technique considered here, electrical resistivity tomography (ERT), is used to reconstruct the subsurface electrical resistivity distribution through measurements of electrical potentials at discrete locations caused by the application of electrical currents. This is usually realized using linear arrays of electrodes, typically stainless-steel spikes driven into the subsurface, which are connected to a resistivity meter via multi-core cables. Modern resistivity meters are computer/micro-controller controlled to collect data autonomously at pre-defined pairs of current and potential electrodes. By using multiple receiver channels, a number of potentials can be measured at the same time, reducing the measurement time considerably. A subsurface resistivity model is then calculated through an inverse process, which aims at minimizing the misfit between modelled and measured data [12].
Different measurement configurations are commonly used. The most frequently applied are perhaps dipole-dipole and Wenner-Schlumberger configurations, which show different sensitivity to lateral and vertical changes in the subsurface electrical resistivity. A detailed discussion about the sensitivity of the most commonly applied measurement configurations can be found in [34]. In general, Wenner-type measurements are more sensitive to vertical changes, while dipole-dipole measurements provide better resolution of lateral changes. In addition, dipole-dipole configurations allow for an efficient use of multi-channel resistivity meters and therefore can be used to sample the subsurface more densely and in less time than Wenner-type configurations. In 3D, Chambers et al. (2002) [35] have shown that dipole-dipole surveys achieve higher resolution than comparable Wenner-type surveys. One common limitation of dipole-dipole surveys is the smaller signal-to-noise ratio, as receiver voltages are smaller compared to Wenner-type measurements. However, when applied to a highly resistive hardrock environment, voltages are high for both measurement types. Hence, in this study we employed dipole-dipole configurations.
Due to the quarry topography, the survey was split into two 3D ERT acquisitions; one on the upper floor and one on the lower floor ( Figure 3). Each part was acquired with an electrode and line spacing of 0.75 m. The upper part consisted of seven parallel lines employing 51 electrodes each, the lower part of 11 parallel lines employing 40 electrodes each. Two additional lines were acquired on the upper floor, employing the same electrode spacing, but with a 2.0 m line spacing. These additional lines used 23 electrodes each. The survey comprised 843 unique electrode locations, each of which was surveyed using high precision surveying equipment. Each line was measured independently using a Geolog2000 GeoTom resistivity meter, employing dipole lengths a of 0.75, 1.5, 2.25, 3.0, and 3.75 m, and dipole separations n of 1 to 8a. Other acquisition schemes exist that can further exploit the spatial distribution of electrodes, such as cross-line or equatorial measurement schemes [33,36]. However, Gharibi and Bentley (2005) [37] have shown that measurements along parallel lines can provide results comparable to full 3D acquisition schemes. A full set of reciprocal data was acquired for each line, where the reciprocal data set includes the same electrode positions as the forward measurement, but with swapped current and potential dipoles. The measured data are defined as the mean of these two measurements and the measurement error as the percentage standard error in the mean. During the course of the survey, 37,146 measurements were made, equaling 18,573 reciprocal pairs of measurements. This survey design was chosen in order to cover the majority of the active extraction area with a resolution sufficient to image fractures observed at the quarry walls and to investigate the void that was observed on the quarry floor. As these fractures and voids were clay and water covered, they were expected to have resistivities orders of magnitude smaller than the surrounding host rock. Thus, it was expected that features, which were significantly smaller than the electrode spacing, could be imaged.
Commonly, ERT employs stainless steel pins that are driven into the subsurface. In the case of a marble environment, this would require shallow drilling and destroy a significant volume of an economic resource, particularly considering the high density of electrode location. Hence, we followed an approach commonly used in archaeology for non-invasively examining ancient constructions, where instead of pins, steel plates are used [38]. As well as allowing non-destructive data acquisition, plates also have a larger surface area than pins, which is advantageous for reducing contact resistances. In this study, we used square, stainless steel plates with a size of 0.1 × 0.1 m, and improved their ground coupling by applying bentonite slurry to each electrode location. Given their separation of 0.75 m, a point-pole approximation for the data analysis is valid [39].
Several ground penetrating radar (GPR) datasets were also collected to provide an additional source of ground truth data for further validation of the ERT results. These data were acquired prior to the ERT survey (before the extraction of the lower floor). The data were acquired using a GSSI SIR 3000 instrument with a 200 MHz antenna along three profiles ( Figure 3). Processing of the data was limited and comprised system gain removal, background removal, band-pass filtering, and applying a gain function.

Contact Resistance
The electrically highly resistive environment makes electrical resistivity measurements challenging. This is because to inject a stable current flow into the ground, the ground (R Ground ) and contact resistances of the injection electrodes (R C1 , R C2 ) have to be overcome. This is best explained using Ohm's law, modelling the resistances as being in series In environmental or engineering geophysical problems, the ground resistance is usually significantly smaller than the contact resistance, and hence the current injection is mainly affected by the contact resistances R C1 and R C2 . In the case of a hard rock environment, the ground resistance is likely to exceed the contact resistance, and hence both ground resistance and contact resistances will determine the transmitter properties.
Conventional resistivity meters are limited in the voltage V tx they can apply. Hence, to overcome the resistance, the injection current I tx needs to be lowered. At the study site, contact resistances were measured for each of the 843 electrode locations. Despite using sensors with a large surface area and treating every single electrode location with bentonite slurry to lower contact resistance, they ranged from about 1 kΩ to >100 kΩ, with a median of 15.7 kΩ, and followed a lognormal distribution (Figure 4a). Split up into the two separate floors, the upper floor had a mean contact resistance of 35.27 kΩ, while the lower floor had a mean contact resistance of 12.37 kΩ. This was due to heavy rain and thus wetter conditions during data acquisition of the lower floor data set. Assuming a resistivity of the marble of about 50 kΩ·m [40], and a maximum transmitter voltage of 500 V, the median maximum transmitter current would be 6.1 mA. However, significantly larger contact resistances were measured and depending on its condition the marble could be of higher resistivity as well, maximum injection currents of 0.1 mA were expected. This has adverse effects on the measured data in two ways: (1) small injection currents are difficult to measure, and (2) the signal-to-noise ratio at the potential dipole deteriorates causing additional measurement errors due to unstable current injection [41,42].  [41] compared data of different commercially available resistivity meters. They focused on different makes and their performance regarding measurement error. More important for this study is a comparison of different systems employing different measurement principles, i.e., whether the source provides a constant current or a constant voltage. At a different site, but within a comparably resistive environment, we repeated dipole-dipole measurements using resistivity meters that employ the two different measurement approaches and analysed their corresponding reciprocal errors (Figure 4b). Note that both systems have comparable maximum transmitter voltages. In this highly resistive environment, we found that the constant current architecture provides superior data quality, as reciprocal errors were consistently smaller than for the constant voltage architecture. This is because the constant current architecture provides small enough injection currents to overcome the ground and contact resistances without the need to explicitly measure these currents. Hence, in such environments constant current resistivity meters are beneficial as they eliminate one contribution to the total measurement error.

Topography
The topography of the survey area and the exact electrode location is usually included in the data inversion, as they are known to have a significant impact on the measured data [43][44][45][46][47][48]. It has also been shown that topographic features outside the actual survey area can impact upon the measured resistivity data and can lead to misinterpretation, especially in the case of 2D data analysis [49,50]. In the study presented here, the survey area is located within a 30 × 40 m excavation that is 10-35 m deep ( Figure 3). Hence, in the immediate vicinity of the survey lines, rock walls of 10-35 m height are present ( Figure 2).
To account for the anticipated topographic effects, the entire quarry domain was modelled and extended by about 130 m in each direction to avoid boundary effects. This required sampling the digital elevation model (DEM, Figure 2) at 1 × 1 m resolution and using this point cloud to mesh the surface topography and underlying subsurface volume. Due to the complex surface topography, an unstructured mesh employing tetrahedral elements was used to discretize the quarry domain. The domain was discretized using TetGen [51] and comprised 1,508,934 elements (Figure 5a). To increase numerical accuracy, the mesh was refined around the electrode locations to have more than four elements between adjacent electrode locations (Figure 5b). The maximum cell volume in the actual imaging area was 0.5 m 3 , while in the extended model domain element volumes of up to 9 × 10 4 m 3 were allowed to reduce computational effort. The necessity to model the entire quarry domain is highlighted in Figure 6, where the potential resulting from a current injection close to the rock wall is shown. An injection current of 10 mA and a marble resistivity of 50 kΩ·m were assumed for this simulation. The potential iso-surfaces follow the topography and extent into the rock wall. The 0.5 V iso-surface propagates more than 15 m into the rock wall away from the point of current injection. Hence, subsurface resistivity variations that are away from the actual imaging area affect the potential distribution in the subsurface and have to be modelled, otherwise artifacts in the inverted resistivity model of the actual area of interest have to be expected.
This simulation highlights another challenge for the resistivity meter. Not only are small currents required to overcome the contact and ground resistance, small currents are also necessary to avoid overloading the receivers. In Figure 6, an injection current of only 10 mA was modelled, which results in potential drops of more than 10 V in the vicinity of the injection dipole. Usually, resistivity meters are designed to measure voltages in the mV range, and voltages that are orders of magnitude higher than that cause the receiver to overload, making it impossible to acquire accurate data. To circumvent receiver overloading, injection current is usually reduced. In the case of the Italian marble quarry, 6.99% of the data were measured using 10 mA, 40.29% using 1 mA, 48.93% using 0.1 mA, and 3.72% using 0.01 mA, highlighting the need to use a resistivity meter than can provide an accurate measure of very small injection currents [42].

Data Inversion
Accounting for the entire quarry topography increases the number of elements considerably. Hence, an inversion code needs to be chosen that can efficiently solve the inverse problem on a large mesh. This usually requires parallelization of the forward and inverse computations. Here, E4D is employed, an open-source, fully parallelized inversion code, that is based on Fortran and can be run on High Performance Computing (HPC) facilities [52]. This code was also chosen as it is designed to perform computations on unstructured grids, which are required to model the complex surface topography.
Limited pre-processing of the data is required before the inversion, which includes obtaining a data error model from the measured reciprocal errors and filtering the data. To develop an error model, an approach introduced by Koestel et al. (2008) [53] was followed, which requires binning the measured reciprocal errors based on their corresponding measured transfer resistance (in a logarithmic domain), estimating the standard deviation of each bin and fitting an equation to those points. During data acquisition, environmental parameters changed significantly between the measurements on the upper and the lower floor. While the upper floor was acquired on a comparably dry surface, heavy rain made the surface of the lower floor wet. Hence, separate error models were determined for the data of the two floors. Figure 7 shows the calculated data and fitted error models, where a linear relationship between measured transfer resistance and data error seemed to fit the data best (R 2 = 0.99 and 0.98 for the data of the upper and lower floor, respectively). The contact resistance is known to have an effect on the data quality, where increasing contact resistances usually leads to decreasing data quality [41]. Since lower contact resistances were measured for the lower floor, the relative resistance error is also smaller than for the data of the upper floor (0.05 and 0.07, respectively). The absolute resistance error was found to be 0.015 Ω for both floors. Before inversion, data were filtered on their measured reciprocal error. Data exceeding 10% reciprocal error were removed from the data set. Despite the high contact resistances, data quality was comparably good. Hence, this step removed only 3.8% of the data, reducing the data set to 17,870 unique data points.
The inversion uses an iteratively re-weighted least squares approach [52], where a L2-norm was employed on the model roughness, which favours spatially smooth resistivity changes. We also tested a L1-norm on the model roughness, which favours sharp resistivity contrasts and would perhaps describe the local settings better. The results were comparable, but the L2-norm resistivity model provided a superior data fit. The inversion converged well and fitted the data to their respective errors (χ 2 = 1.0). As the computational problem was large (>1.5 M cells, 17,870 data points; hence the Jacobian matrix had 2.7 × 10 10 float elements), the inversion was run on an HPC machine using 96 cores. The computation took about 4.5 h and required 18 iterations. Figure 8 shows the inverted resistivity model, for which the relative root-mean-squared (RMS) misfit between modelled and measured data was RMS = 2.1%. For visualizing and interpreting the resistivity models, all cells with sensitivities <5 × 10 −3 were removed from the model. Using this approach, the depth of investigation can be approximated to be about 10 m for the upper floor and about 8 m for the lower floor. Hence, it should be possible to resolve the full extent of the cave as measured in the field.

Results
The model shows resistivities ranging between 100 and 580,000 Ω·m, thus spanning four orders of magnitude. The location of the highest resistivity (ρ = 580 kΩ·m, lower floor at x = −16-−18 m, y = −21-−22 m) corresponds to the location of the collapsed roof of the cave and hence to an air-filled void. Other shallow high resistivity anomalies appear on the upper floor, but usually remain below 300 kΩ·m. Otherwise, the model shows reasonably homogeneous, highly resistive material (>20 kΩ·m, Figure 8c) in the majority of the model, indicating compact marble with a limited occurrence of clay-filled fractures. On the upper floor, this comparably homogeneous layer is intersected by two linear features of reduced resistivity; a shallow medium resistive feature (ρ < 10 kΩ·m) at y = −14 m, and a deep low resistive feature (ρ < 5 kΩ·m) at y = −24 m. Both features have a similar strike, and can be identified in a horizontal section through the model at z = 245.5 m (Figure 8c). The deeper reaching and lower resistive feature extends into the lower floor, where a similar linear low resistivity feature can be observed. With greater depth, both increase in size considerably and become less resistive (ρ < 1 kΩ·m).
At the surface, the different environmental conditions during data acquisition become apparent. Generally, the upper floor shows higher resistivities at the surface than those observed in the lower floor. This is due to heavy rainfall and hence increased surface moisture during the data acquisition on the lower floor. Some very low resistivity (ρ < 0.5 kΩ·m), shallow features can be observed on the surface (Figure 8a  slices; the model shows the resistivity distribution in the actual survey areas. Cells with sensitivities <5 × 10 −3 were removed. Note the resistivity color range, which spans 100 to 100,000 Ω·m. The maximum resistivity in the model exceeds this range at 580,000 Ω·m. Features representing surface water accumulation, the clay-lined fracture, and cave outline are apparent; field observations of faults are shown as dashed lines.

Discussion
Ground-truth information for the site based on direct visual observations (see 'Study Site' section) provide a means of assessing and validating the resistivity model. These observations of discontinuities recorded at the study site are shown in Figure 9, and show a remarkable agreement with low resistive features of the resistivity model. The observed fracture in the rock wall continues in the subsurface and extends in size with depth. In the lower floor, the location of the highest resistivities are located in the surveyed location of the largest void observed on the quarry floor, which was air-filled in the upper 4.2 m. The depth to the measured water table shows a very good agreement with the ρ = 2.5 kΩ·m iso-surface (which is shown fully opaque in Figure 9). The shape of it is confirmed by visual observations of the cave's shape, which seemed to extend significantly just above the water table. Similarly, the iso-surface gains elevation towards x = −10 m, which corresponds to a second, smaller void that was observed on the lower floor. The low resistive features observed beneath the upper and lower floor seem to be spatially connected, which is likely, as an expression of the fracture of the adjacent rock wall could be followed through the quarry on both floors. Figure 9b shows the resistivity profile at the location of the observed collapsed cave roof, and the indication of the air-and water filled parts of the cave. The field observations show good agreement with the model resistivity profile. The resistivity model shows very high values in the upper part of the air-filled cave, which reduce significantly towards the water-filled part. However, the low resistive feature, which represents the water-filled part of the cave, seems to be overestimated in its extent. This is likely to be an effect of the smoothness constraints of the inversion. They favour smooth resistivity contrasts, thereby underestimating the true resistivity of the caves fill, while overestimating its size. However, an indication of the lower boundary of the cave can be related to the change in the resistivity gradient, which approaches zero at about 8 m depth, corresponding well with the measured cave depth (7 m, Figure 9b). At 8.25 m depth, resistivities start to increase. However, the sensitivity becomes very small in those depths, and hence resistivities remain at lower values compared to the volumes of undisturbed marble found close to the surface. The results of the additional ground truthing using GPR are shown in Figure 10, together with the 2.5 kΩ·m iso-surface from the 3D ERT model and field observations. Generally, the data show some weak reflections and diffraction hyperbolas that can be related to fractures or voids in the otherwise homogeneous marble. At y = −20-−24 m the data show high amplitude hyperbolas that correspond in their location with the mapped fracture and voids. The hyperbolas of profiles P1 and P2 seem to originate from very shallow features, while on P3 the most pronounced hyperbola has its origin at ≈1.2 m depth (using a wave velocity of 0.117 m/ns for time-to-depth conversion, measured on an adjacent marble block). This hyperbola has a negative amplitude. Hence, it is likely that it corresponds to a contrast between marble and a water-filled void [54]. Note that this interface would be at a higher level than imaged in the ERT results. However, between the two surveys, extraction advanced significantly and the water table was lowered to allow for continuing extraction. Although the data quality in the study presented here was excellent, performing ERT measurements on rock faces remains challenging. This is mainly because of the very high contact resistances that are encountered in such environments. These are known to be detrimental to data quality [41]. Contact resistances can be lowered by employing electrodes with large surface areas and by treating electrode locations with a low resistive material to improve the galvanic ground coupling. Here, we used plate electrodes and a bentonite slurry to lower contact resistances. By doing so, contact resistances remained high, but sufficient coupling to the ground could be established to inject current. We have also shown that the resistivity meter has to be chosen carefully, as different measurement principles associated with different meter models are likely to have an impact on the data quality.
This application also requires an accurate DEM of a potentially complex quarry topography. Due to the high resistivities, potentials propagate far away from the current injection dipoles. Hence, resistivity measurements are sensitive to subsurface resistivity anomalies outside the actual survey area. While this has to be considered in the processing and interpretation of 3D data sets, it is even more important for 2D ERT data sets, where the intrinsic assumption of a 2.5D subsurface is likely being violated. As shown in the case study here, by modeling the entire quarry domain, these outer-space sensitivities [55] can be accounted for and artifacts in the resulting resistivity model avoided. However, this requires a fine discretization of the quarry domain, resulting in potentially large computational meshes, which require sufficient computational power and storage. Recent work aims to reduce the computational cost of ERT forward and inverse modelling [36,52,56], and will allow a wider application of this geophysical technique, allowing the requirement of HPC technology to be overcome.

Summary and Conclusions
Here, a case study is presented where electrical resistivity tomography (ERT) has been applied to characterise the volumetric electrical properties of an operational marble quarry. This technique was chosen as electrical resistivity is sensitive to changes in clay and water content, which are parameters linked to weathering, karstification and fracturing. These processes may alter rock properties that will impact the economic value of the quarried stone and its safe extraction.
Application of ERT within hard-rock environments is challenging. Contact resistances are very high and potentials propagate far from the injection dipoles, requiring a careful decision of the employed electrodes, measurement system, and data processing algorithms. We show how these challenges can be overcome, and that despite the unfavourable conditions for conducting electrical resistivity measurements, good quality data can be obtained. Comparing the derived resistivity model with field observations and complementary ground penetrating radar (GPR) data highlights the reliability of the resistivity data in imaging faults, and air-and water-filled voids. These results show that ERT can be used as a non-invasive exploration technique to image stone properties before excavation, and thus as a tool to aid the design of mining operations.
Future developments of the approach considered here could include the integrated interpretation of geophysical models with laboratory data to define stone quality classes. By using geo-statistical interpolation techniques, such as kriging or inverse distance interpolation, these integrated data analyses could be used to further refine resource estimates, and optimize and guide extraction progress. By targeting extraction of high quality material, quarries will be able to operate more economically, but also more sustainably, as this will also lower the amount of quarry waste that is produced. Although this paper focused on the application of ERT to ornamental stone quarrying, the challenges, solutions, and results shown will be transferable to other highly resistive environments, such as investigation of mountainous permafrost, rock slopes, or other mining problems.