# Characterizing Near-Surface Fractured-Rock Aquifers: Insights Provided by the Numerical Analysis of Electrical Resistivity Experiments

Serdar Demirel, Delphine Roubinet, James Irving, Emily Voytek

► To cite this version:

Serdar Demirel, Delphine Roubinet, James Irving, Emily Voytek. Characterizing Near-Surface Fractured-Rock Aquifers: Insights Provided by the Numerical Analysis of Electrical Resistivity Experiments. Water, 2018, 10 (9), pp.1117. 10.3390/w10091117. [https://doi.org/10.3390/w10091117]

HAL Id: hal-01976565

https://hal.science/hal-01976565

Submitted on 10 Jan 2019

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
Article

# Characterizing Near-Surface Fractured-Rock Aquifers: Insights Provided by the Numerical Analysis of Electrical Resistivity Experiments

Serdar Demirel <fnref n="1" />,†, Delphine Roubinet <fnref n="2" />,† ID, James Irving <fnref n="1" />,* ID and Emily Voytek <fnref n="1" /> ID

Received: 16 July 2018; Accepted: 16 August 2018; Published: 23 August 2018

<figure><img src="image_2.png" /></figure>

check for
updates

**Abstract:** Fractured-rock aquifers represent an important part of the groundwater that is used for domestic, agricultural, and industrial purposes. In these natural systems, the presence and properties of fractures control both the quantity and quality of water extracted, meaning that knowledge about the fractures is critical for effective water resource management. Here, we explore through numerical modeling whether electrical resistivity (ER) geophysical measurements, acquired from the Earth's surface, may potentially be used to identify and provide information about shallow bedrock fractures. To this end, we conduct a systematic numerical modeling study whereby we evaluate the effect of a single buried fracture on ER-profiling data, examining how the corresponding anomaly changes as a function of the fracture and domain characteristics. Two standard electrode configurations, the Wenner-Schlumberger (WS) and dipole-dipole (DD) arrays, are considered in our analysis, with three different spacing factors. Depending on the considered electrode array, we find that the fracture dip angle and length will impact the resistivity anomaly curves differently, with the WS array being better adapted for distinguishing between sub-horizontal and sub-vertical fractures, but the DD array leading to larger overall anomaly magnitudes. We also find that, unsurprisingly, the magnitude of the resistivity anomaly, and thus fracture detectability, is strongly affected by the depth of overburden and its electrical resistivity, as well as the fracture aperture and contrast between the fracture and bedrock resistivities. Further research into the electrical properties of fractures, both above and below the water table, is deemed necessary.

**Keywords:** fractured rock; geophysics; electrical resistivity; numerical simulations; sensitivity

## 1. Introduction

Groundwater, which represents about 98% of all liquid freshwater on Earth, is the most extracted raw material from the subsurface in the world. It is used for irrigation, domestic, and industrial purposes, and it provides drinking water for approximately half of the world's population [1,2]. With the exception of sand and gravel aquifers, groundwater is stored in hard-rock geological environments that are characterized by discontinuities such as fissures, voids, and fractures. Although some of these discontinuities may act as barriers to flow, they are typically highly permeable structures that have a strong control on both the quantity and quality of groundwater extracted [3]. In fractured-rock aquifers, for example, fractures act as preferential flow paths that can greatly increase the risk of rapid contaminant propagation over large distances. This concerns not only surface-originating pollutants coming from industry and agriculture, but also subsurface-originating

<footnote n="1">¹ Institute of Earth Sciences, University of Lausanne, 1015 Lausanne, Switzerland; serdar.demirel@unil.ch (S.D.); emily.voytek@unil.ch (E.V.)</footnote>

<footnote n="2">² Geosciences Montpellier (CNRS), University of Montpellier, 34090 Montpellier, France; delphine.roubinet@umontpellier.fr</footnote>

<footnote>* Correspondence: james.irving@unil.ch</footnote>

<footnote>† Authors contributed equally to this work.</footnote>
contamination related to increased resource extraction, energy production, and geological storage activities [4–6]. As a result, identifying the presence of fractures, as well as determining their properties and connectivity, are key requirements for the sustainable management and protection of groundwater.

Numerous studies have focused on the development of methods to identify and characterize subsurface fractures and fracture networks, in order to improve conceptual and numerical models of groundwater flow and contaminant transport [7–10]. In this context, applied geophysical methods are promising, because many of the physical properties to which these methods are sensitive are strongly influenced by the presence of fractures, and the corresponding measurements can be acquired rapidly and non-invasively from the ground surface and/or from boreholes. Examples include the use of seismic [11,12], ground-penetrating radar [13,14], electrical resistivity [15,16], self potential [17,18], induced polarization [19,20], and electromagnetic methods [21,22].

The electrical resistivity (ER) method is of particular interest in fractured-rock environments because (i) numerous field, laboratory, and theoretical studies have shown that the electrical conductivity (inverse of resistivity), and therefore ER data, are sensitive to the presence and characteristics of fractures, such as their density, orientation, and aperture [23–25]; (ii) ER measurements can be obtained over a range of spatial scales [26,27]; and (iii) analogies between fluid and electric current flow in fractured media suggest that ER data may be used to quantify subsurface hydrogeological properties [28,29]. For these reasons, ER studies in fractured rock have been numerous, and include the development and application of anisotropic tomographic inverse methods [30–32]; the use of azimuthal resistivity surveys to estimate predominant fracture orientations [15,25,33] and hydraulic anisotropy [34–37]; and the acquisition of electrical resistivity tomography (ERT) data for the identification and characterization of fracture zones [38–41].

Despite the large number of studies where ER methods have been applied to fractured rock, quantitative knowledge of how individual fractures influence ER data appears to be missing from the existing literature. This knowledge is critical for understanding what information ER data may contain regarding discrete fractures and their properties. Indeed, all of the studies cited above have focused on either the impact of multiple fractures or fracture networks on the bulk resistivity characteristics of the subsurface, or the influence of fracture zones on ERT data. Recent numerical modeling results, however, suggest that, despite the well-known limits in resolution of ERT methods from a tomographic inverse standpoint [42–44], even a single millimeter-scale fracture in the near-surface can significantly influence the distribution of subsurface electric potential [45], which in some cases can produce a strong anomaly in ER profiling data [46]. A systematic study of the magnitude of changes in such data, corresponding to different fracture configurations, is therefore necessary to evaluate whether, and under what conditions, ER data might be used to characterize individual fractures.

In this work, we explore whether surface-based ER measurements may be used to recover information about fractures in near-surface bedrock aquifers. To this end, we conduct a detailed numerical analysis using a highly efficient and accurate 2.5D modeling approach that is based on an explicit representation of the fractures, with no limit to the complexity of the considered fracture distribution [45,46]. After presenting the methodological background for our study, we show the results of ER profiling over a single buried fracture considering various values for the fracture length, dip angle, aperture, and electrical resistivity, as well as different values for the overburden thickness and resistivity. Two commonly utilized and complementary electrode configurations, the Wenner-Schlumberger and dipole-dipole arrays, are considered. The presented results and their detailed interpretation lead to improved understanding of (i) the signature of a buried fracture in ER-profiling data and how it changes as a function of the fracture and overburden characteristics; (ii) the importance of this signature in comparison with the inherent noise in ER measurements; and (iii) the impact of the electrode array on fracture detectability and resolution with ER methods.
## 2. Methodological Background

To numerically investigate the impact of the presence and properties of shallow bedrock fractures on ER measurements, we consider a 50 × 50-m simulation domain (Figure 1a). The domain consists of homogeneous bedrock having electrical resistivity $ρ_b$, which is covered by a layer of overburden having thickness $d$ and resistivity $ρ_o$. A single fracture having length $ℓ$, aperture $b$, dip angle $α$, and resistivity $ρ_f$ intersects the bedrock surface 22 m from the left-hand side of the domain.

<figure><img src="image_3.png" /><figcaption>**Figure 1.** (a) Single-fracture configuration considered in our simulations with the fracture, bedrock, and overburden represented using blue, white, and brown, respectively. (b) Diagram showing the Wenner-Schlumberger and (c) dipole-dipole electrode arrays, where *a* is the minimum considered electrode spacing and *n* is the spacing factor (*n* = 2 shown). A and B denote the current electrodes, whereas M and N denote the potential measurement electrodes. The symbols below the arrays indicate the center point of each measurement along with the relative depths of the three considered spacing-factor values, and correspond to the symbols used in Figures 5–8.</figcaption></figure>

The acquisition of ER-profiling data over the domain in Figure 1a is simulated using the recently developed 2.5-D electrical modeling code of [46], which builds on the 2-D discrete-dual-porosity modeling approach developed in [45]. This methodology results in efficient and accurate simulations of resistivity experiments, whereby fractures are represented explicitly and the effects of the electrically conductive matrix are properly taken into account. The subsurface geological structures and material properties are assumed to extend infinitely in the out-of-profile direction. In contrast to 2-D models, this 2.5-D approach correctly represents the point-source nature of the resistivity electrodes.

To conduct the ER simulations, the domain in Figure 1a is discretized into square 0.5-m cells. Resistivity electrodes are considered every 1 m along the ground surface from 2 m to 48 m. Two commonly used and complementary electrode configurations are considered for the measurements. The Wenner-Schlumberger array (Figure 1b) involves measurement of the difference in potential between two electrodes located between the current electrodes, and is known to have greater sensitivity to vertical changes in subsurface resistivity. In contrast, with the dipole-dipole array (Figure 1c), the potential measurement is made adjacent to the current electrodes, which results in greater sensitivity to lateral resistivity changes [27]. For both arrays, we consider a minimum electrode spacing of *a* = 2 m and spacing-factor values of *n* = 1, 2, 3 (Figure 1b,c). Increasing the latter parameter increases the length of the array, and, thus, the effective depth of investigation [27]. Measurements
are simulated for each array configuration at 1-m intervals across the top of the domain. For each measurement, a 100 mA electric current is passed between electrodes A and B, and the potential difference between electrodes M and N is recorded (Figure 1b,c).

The raw data acquired in an ER experiment are the voltages between the potential electrodes M and N. As these voltages depend on the amount of current passed between electrodes A and B, it is customary to transform the results into apparent resistivity $\rho_a$, which is defined as the resistivity of the isotropic, homogeneous half-space that would produce the equivalent voltage [27]. The apparent resistivity is calculated using

$$ \rho_a = K \cdot \frac{\phi_{MN}}{I_{AB}}, \qquad (1) $$

where $\phi_{MN} = \phi_M - \phi_N$ is the difference in potential between electrodes M and N, $I_{AB}$ is the current passed between electrodes A and B, and K is the electrode-array-dependent geometric factor. For the Wenner-Schlumberger array (Figure 1b), the geometric factor is given by

$$ K = \pi n (n + 1) a. \qquad (2) $$

The geometric factor for the dipole-dipole array is [27]

$$ K = -\pi n (n + 1) (n + 2) a. \qquad (3) $$

To evaluate the sensitivity of the simulated ER profiling data to the presence of the fracture in Figure 1a, we define the apparent resistivity anomaly, $\rho_a^*$, as the relative difference (in %) between the apparent resistivity of the considered fractured domain, $\rho_a$, and that of the background unfractured domain, $\rho_a^u$. That is,

$$ \rho_a^* = \frac{\rho_a - \rho_a^u}{\rho_a^u} \times 100. \qquad (4) $$

Consideration of $\rho_a^*$ in our analysis, as opposed to $\rho_a$, facilitates the comparison of the observed anomaly against an anticipated level of measurement error. For ER experiments, this is commonly found to be on the order of 1% of the measured value based on reciprocal analysis [39,47–49]. For each fracture configuration, electrode array, and n value considered in our study, we plot $\rho_a^*$ versus the midpoint of the electrode array along the profile. For comparison between scenarios, we define the *anomaly magnitude* as the absolute difference between the maximum and minimum $\rho_a^*$ values, and the *anomaly width* as the distance between the first and last points where $|\rho_a^*|$ is equal to 5% of the anomaly magnitude (Figure 2). Note that the anomaly magnitude and width are calculated based on values observed at discrete electrode locations. As a result, the values do not necessarily correspond to the true magnitude and width if the anomaly were to be evaluated continuously. Nonetheless they provide important quantitative information regarding the impact of changes in the different fracture scenarios on the observed resistivity response.

In this study, we systematically examine the effect of the different fracture and overburden properties described in Figure 1a on the simulated ER-profiling data. In all simulations, the electrical resistivity of the bedrock is kept at a constant value of 10,000 Ω·m. Table 1 summarizes the fracture, overburden, and bedrock properties for the 16 test cases that are considered in our analysis. In Cases 1–6, we study the impact of the fracture angle and length on the apparent resistivity anomaly curves. We consider in these cases a fracture electrical resistivity that is three orders of magnitude less than the bedrock resistivity [45,46], whereas the overburden is assumed to be one order of magnitude less resistive than the bedrock. Fracture dip angles of 10°, 45°, and 80° are considered, as well as fracture lengths of 5 m and 10 m. In the subsequent configurations (Cases 7–16), we use Case 2 as a base configuration and examine how the apparent resistivity anomaly magnitude and width are affected by
changes in overburden thickness (Cases 7–9), overburden resistivity (Cases 10,11), fracture aperture (Cases 12–14), and fracture resistivity (Cases 15,16). For these cases, the fracture dip angle and length are held constant at 45° and 5 m, respectively.

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Diagram showing how the apparent resistivity anomaly, $\rho_a^*$, is plotted as a function of the horizontal measurement position, $x$, and how the anomaly magnitude and width are determined. The anomaly magnitude is defined as the absolute difference between the maximum and minimum observed anomaly values, in %. The anomaly width is defined as the distance between the first and last points where the absolute anomaly value is equal to 5% of the anomaly magnitude, the latter of which is indicated by a red band in the figure. Note that the measurement position corresponds to the midpoint of the considered electrode array.</figcaption></figure>

<table><caption>**Table 1.** Parameter values corresponding to the 16 test cases considered in our analysis. For a description of the variables, see Figure 1.</caption><thead><tr><th>Cases</th><th>ρ<sub>b</sub> [Ω·m]</th><th>ρ<sub>f</sub> [Ω·m]</th><th>ρ<sub>o</sub> [Ω·m]</th><th>d [m]</th><th>α [°]</th><th>b [mm]</th><th>ℓ [m]</th></tr></thead><tbody><tr><td>1, 2, 3</td><td>10,000</td><td>10</td><td>1000</td><td>1</td><td>10, 45, 80</td><td>1</td><td>5</td></tr><tr><td>4, 5, 6</td><td>10,000</td><td>10</td><td>1000</td><td>1</td><td>10, 45, 80</td><td>1</td><td>10</td></tr><tr><td>7, 8, 9</td><td>10,000</td><td>10</td><td>1000</td><td>0.5, 1.5, 2</td><td>45</td><td>1</td><td>5</td></tr><tr><td>10, 11</td><td>10,000</td><td>10</td><td>100, 10,000</td><td>1</td><td>45</td><td>1</td><td>5</td></tr><tr><td>12, 13, 14</td><td>10,000</td><td>10</td><td>1000</td><td>1</td><td>45</td><td>2, 3, 4</td><td>5</td></tr><tr><td>15, 16</td><td>10,000</td><td>1, 100</td><td>1000</td><td>1</td><td>45</td><td>1</td><td>5</td></tr></tbody></table>

### 3. Results and Discussion

We present below the results of our analysis of the various test cases described in Table 1. In Section 3.1, we focus on the interpretation of the resistivity anomaly curves corresponding to Cases 1–6, in order to assess (i) how the curves are influenced by the particular details of the fracture configuration; and (ii) whether the sensitivity of the ER-profiling data to fracture dip angle and length is high enough to suggest that we might be able to eventually invert for these parameters. In Section 3.2, we focus on the effect of changing overburden thickness and electrical resistivity, as well as fracture aperture and resistivity, on the apparent resistivity anomaly magnitude and width. Here, the aim is to understand how these overburden and fracture properties will influence the potential detectability of the fracture. For all of the considered cases, Tables 2 and 3 contain the anomaly magnitude and width data, respectively.
<table><caption>**Table 2.** Apparent resistivity anomaly magnitudes (in %) for the test cases described in Table 1. Measurement configurations are denoted by WS (Wenner-Schlumberger) or DD (dipole-dipole) followed by the spacing factor (n = 1, 2, 3). Values greater than 2.0% are in bold. See Figure 2 for how these values are determined.</caption><thead><tr><th>Case</th><th>WS-1</th><th>WS-2</th><th>WS-3</th><th>DD-1</th><th>DD-2</th><th>DD-3</th></tr></thead><tbody><tr><td>1</td><td><strong>3.3</strong></td><td><strong>4.2</strong></td><td><strong>4.8</strong></td><td><strong>5.0</strong></td><td><strong>5.6</strong></td><td><strong>6.0</strong></td></tr><tr><td>2</td><td>1.4</td><td>1.7</td><td>1.8</td><td><strong>3.3</strong></td><td><strong>4.5</strong></td><td><strong>4.7</strong></td></tr><tr><td>3</td><td>1.0</td><td>1.2</td><td>1.2</td><td><strong>3.4</strong></td><td><strong>4.5</strong></td><td><strong>4.6</strong></td></tr><tr><td>4</td><td><strong>3.6</strong></td><td><strong>4.5</strong></td><td><strong>4.4</strong></td><td><strong>5.0</strong></td><td><strong>7.4</strong></td><td><strong>7.8</strong></td></tr><tr><td>5</td><td>1.4</td><td>1.8</td><td>1.9</td><td><strong>3.3</strong></td><td><strong>4.7</strong></td><td><strong>5.0</strong></td></tr><tr><td>6</td><td>1.1</td><td>1.3</td><td>1.3</td><td><strong>3.5</strong></td><td><strong>4.8</strong></td><td><strong>5.0</strong></td></tr><tr><td>7</td><td><strong>2.6</strong></td><td><strong>3.0</strong></td><td><strong>3.4</strong></td><td><strong>6.5</strong></td><td><strong>7.9</strong></td><td><strong>8.2</strong></td></tr><tr><td>8</td><td>0.8</td><td>1.1</td><td>1.2</td><td>1.7</td><td><strong>2.9</strong></td><td><strong>3.2</strong></td></tr><tr><td>9</td><td>0.5</td><td>0.8</td><td>0.9</td><td>0.9</td><td><strong>2.0</strong></td><td><strong>2.4</strong></td></tr><tr><td>10</td><td>0.2</td><td>0.2</td><td>0.2</td><td>0.4</td><td>0.6</td><td>0.7</td></tr><tr><td>11</td><td><strong>4.5</strong></td><td><strong>6.0</strong></td><td><strong>7.1</strong></td><td><strong>9.2</strong></td><td><strong>13.2</strong></td><td><strong>14.4</strong></td></tr><tr><td>12</td><td><strong>2.2</strong></td><td><strong>2.7</strong></td><td><strong>2.9</strong></td><td><strong>5.2</strong></td><td><strong>7.5</strong></td><td><strong>8.1</strong></td></tr><tr><td>13</td><td><strong>2.7</strong></td><td><strong>3.4</strong></td><td><strong>3.6</strong></td><td><strong>6.6</strong></td><td><strong>9.8</strong></td><td><strong>10.7</strong></td></tr><tr><td>14</td><td><strong>3.1</strong></td><td><strong>4.0</strong></td><td><strong>4.2</strong></td><td><strong>7.7</strong></td><td><strong>11.7</strong></td><td><strong>12.8</strong></td></tr><tr><td>15</td><td><strong>4.5</strong></td><td><strong>6.1</strong></td><td><strong>6.3</strong></td><td><strong>11.4</strong></td><td><strong>18.7</strong></td><td><strong>21.0</strong></td></tr><tr><td>16</td><td>0.3</td><td>0.4</td><td>0.5</td><td>0.5</td><td>0.7</td><td>0.7</td></tr></tbody></table>

<table><caption>**Table 3.** Apparent resistivity anomaly widths (in m) for the test cases described in Table 1. Measurement configurations are denoted by WS (Wenner-Schlumberger) or DD (dipole-dipole), followed by the spacing factor (n = 1, 2, 3). See Figure 2 for how these values are determined.</caption><thead><tr><th>Case</th><th>WS-1</th><th>WS-2</th><th>WS-3</th><th>DD-1</th><th>DD-2</th><th>DD-3</th></tr></thead><tbody><tr><td>1</td><td>15.8</td><td>20.8</td><td>24.7</td><td>14.9</td><td>19.1</td><td>21.7</td></tr><tr><td>2</td><td>17.0</td><td>22.4</td><td>26.8</td><td>13.9</td><td>17.0</td><td>19.6</td></tr><tr><td>3</td><td>17.8</td><td>23.4</td><td>27.1</td><td>13.4</td><td>16.7</td><td>19.6</td></tr><tr><td>4</td><td>19.6</td><td>25.3</td><td>29.1</td><td>13.6</td><td>18.2</td><td>21.1</td></tr><tr><td>5</td><td>19.2</td><td>25.2</td><td>29.2</td><td>15.7</td><td>19.0</td><td>21.7</td></tr><tr><td>6</td><td>19.4</td><td>25.5</td><td>28.8</td><td>14.4</td><td>17.9</td><td>21.0</td></tr><tr><td>7</td><td>15.5</td><td>21.0</td><td>25.3</td><td>12.0</td><td>15.4</td><td>17.9</td></tr><tr><td>8</td><td>18.4</td><td>23.5</td><td>27.3</td><td>15.7</td><td>18.0</td><td>20.6</td></tr><tr><td>9</td><td>19.8</td><td>24.3</td><td>27.8</td><td>17.5</td><td>19.2</td><td>21.9</td></tr><tr><td>10</td><td>19.0</td><td>26.0</td><td>28.5</td><td>16.1</td><td>19.7</td><td>22.6</td></tr><tr><td>11</td><td>13.7</td><td>19.3</td><td>23.4</td><td>7.8</td><td>11.4</td><td>14.5</td></tr><tr><td>12</td><td>17.4</td><td>22.8</td><td>27.0</td><td>14.1</td><td>17.1</td><td>19.5</td></tr><tr><td>13</td><td>17.8</td><td>23.1</td><td>27.2</td><td>14.5</td><td>17.3</td><td>19.6</td></tr><tr><td>14</td><td>18.1</td><td>23.4</td><td>27.4</td><td>14.8</td><td>17.6</td><td>19.8</td></tr><tr><td>15</td><td>19.1</td><td>24.2</td><td>28.1</td><td>15.6</td><td>18.0</td><td>20.4</td></tr><tr><td>16</td><td>16.4</td><td>21.3</td><td>25.8</td><td>14.5</td><td>18.2</td><td>21.0</td></tr></tbody></table>

### 3.1. Sensitivity to Fracture Geometrical Characteristics

We begin our analysis by studying how the fracture geometrical properties, i.e., fracture position, dip angle, and length, affect the apparent resistivity anomaly, $\rho_a^*$. Figures 3 and 4 show $\rho_a^*$ plotted as a function of the midpoint of the electrode array along the profile line, x, for Cases 1–3 and 4–6 from Table 1, respectively. Curves for the Wenner-Schlumberger (WS) and dipole-dipole (DD) arrays for spacing-factor values n = 1, 2, 3 are shown. Videos showing the distribution of electric potential in the subsurface as the electrode array is moved along the line are also included as supplementary material. The description of these videos for Cases 1–3 with n = 1 (Appendix A) provides insight into the behavior of the difference in electric potential between the fractured and corresponding unfractured domains, $\phi^*(x, z)$. This in turn leads to a detailed understanding of the measured potential difference between electrodes M and N relative to background, $\phi_{MN}^*$, along the profile. The conclusions drawn in
Appendix A are used to help interpret the behavior of $\rho_a^*(x)$, which is related to the measured potential difference between M and N through Equations (1)–(4). Note that, as the geometric factor for the dipole-dipole array is negative (Equation (3)), we expect the opposite trends between the apparent resistivity anomaly and the measured difference in potential for this array.

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Apparent resistivity anomaly $\rho_a^*$ (in %) as a function of horizontal position for a single, 1-mm-aperture fracture having length $l = 5$ m, electrical resistivity $\rho_f = 10$ Ω·m, and hosted in bedrock having resistivity $\rho_b = 10,000$ Ω·m. The fracture is covered by a 1 m thick layer of overburden, having resistivity $\rho_o = 1000$ Ω·m. Anomalies for fracture angles of 10°, 45°, and 80° are shown, corresponding to Cases 1, 2, and 3 from Table 1, respectively. Measurement configurations are denoted by WS (Wenner-Schlumberger) or DD (dipole-dipole) followed by the spacing factor ($n = 1, 2, 3$).</figcaption></figure>

### 3.1.1. Impact of the Fracture Position and Dip Angle

From the ER profiling movies described in Appendix A for Cases 1–3, we know that the impact of the fracture geometrical properties on $\rho_a^*(x)$ will differ significantly depending on the electrode array considered. For the WS array with $n = 1$, for example, we learned that the top extremity of the fracture determines the position at which the minimum value of $\phi_{MN}^*$ (and thus $\rho_a^*$) is observed, whereas the fracture dip angle controls the overall shape of the anomaly curves. In Figure 3 (WS-1), the minimum value of $\rho_a^*$ is seen to occur at $x = 24$ m for all fracture angles ($\alpha = 10°, 45°, \text{and } 80°$). For the sub-horizontal fracture configuration ($\alpha = 10°$), $\rho_a^*(x)$ displays a monotonic decrease and increase before and after the minimum value, respectively. A more complex behavior is observed before the minimum value as the fracture dip angle increases, leading to a peak at the position of the top extremity of the fracture ($x = 22$ m) for the sub-vertical fracture configuration ($\alpha = 80°$). These results suggest that it may be possible to distinguish between sub-horizontal and sub-vertical fractures from the shape of the WS-1 anomaly. In contrast, this distinction is not possible using the DD array for $n = 1$. Indeed, we see in Figure 3 (DD-1) that the fracture dip angle impacts the measurement location at which the
minimum value is observed, as well as the anomaly width; however, there is nothing in the overall form of the anomaly that would allow us to distinguish between sub-horizontal and sub-vertical fractures. The latter behavior occurs because the measurement electrodes are, in comparison with WS configuration, further away from the injection current electrode, implying that the measured differences in potential are less sensitive to vertical changes in properties. Note also that, with few exceptions, the anomaly magnitudes for the WS and DD arrays with *n* = 1 decrease with increasing dip angle (Table 2), with the greatest change being seen between Case 1 (α = 10°) and Case 2 (α = 45°). This reflects the fact that, for lesser dip angles, more of the fracture lies closer to the surface where the ER data are most sensitive to variations in resistivity. The anomaly magnitudes are also larger for the DD array compared to the WS array, most notably for α = 45° and α = 80°.

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Apparent resistivity anomaly $\rho_a^*$ (in %) as a function of horizontal position for a single, 1-mm-aperture fracture having length *l* = 10 m, electrical resistivity $\rho_f$ = 10 $\Omega \cdot m$, and hosted in bedrock having resistivity $\rho_b$ = 10,000 $\Omega \cdot m$. The fracture is covered by a 1 m thick layer of overburden, having resistivity $\rho_o$ = 1000 $\Omega \cdot m$. Anomalies for fracture angles of 10°, 45°, and 80° are shown, corresponding to Cases 4, 5, and 6 from Table 1, respectively. Measurement configurations are denoted by WS (Wenner-Schlumberger) or DD (dipole-dipole) followed by the spacing factor (n = 1, 2, 3).</figcaption></figure>

The ER profiling movies, provided as supplementary material, also help us to understand the impact of increasing the spacing factor *n* for both the WS and DD arrays. In the movies, we observe that the overall amplitude of $\phi_{MN}^*$ decreases as *n* increases, which is related to a greater proportion of electric current traveling deeper into the subsurface and below the fracture, meaning that larger areas of the unfractured bedrock are investigated. This leads to less current flow through the fracture, and thus to smaller values of the potential difference measured between electrodes M and N relative to the background value. Quite importantly, a different behavior is observed for $\rho_a^*$, in the sense that the anomaly magnitude tends to increase as *n* increases (Figure 3, Table 2). As increasing *n* results in locating the potential electrodes M and N further away from the current electrodes, smaller absolute values of $\phi_{MN}$ (and thus $\rho_a$) tend to be measured, leading to larger absolute values of the relative
anomaly as defined in Equation (4). Interestingly, these results suggest that increasing the spacing factor may help us to detect fractures when considering typical levels of measurement error based on a percentage of the measured value, simply because lower values of the potential are recorded as *n* increases. However, it must be remembered that, at some point, the measured potential differences will be of such low amplitude that they will become indistinguishable from other sources of error not considered in our analysis, for example background noise related to telluric currents or self-potential effects. Also note that, although increasing *n* may serve some benefit in terms of increasing the ER anomaly magnitude, it results in a broader, "stretched-out" anomaly, with less ability to resolve individual fractures (Figure 3, Table 3) due to the corresponding increase in the array length.

Noticeable changes are also observed in the shape of the apparent resistivity anomaly for the sub-horizontal fracture configuration ($\alpha = 10°$) as *n* increases. For the other fracture dip angles, changes in the anomaly shape are relatively minor (Figure 3). For the WS array, we see that for $\alpha = 10°$, the position at which the minimum value of $\rho_a^*$ is observed does not change as *n* increases. However, the simple decrease and increase that occur around $x = 24$ m when *n* = 1 are replaced with a series of decreases and increases when *n* = 2 and *n* = 3. In the latter cases, the larger spacing between the current and potential electrodes allows us to distinguish between the effects of each of the current electrodes. In particular, the ER profiling movies indicate that (i) the decrease and increase that are observed when *n* = 2 and *n* = 3 for $15 ≤ x ≤ 20$ m correspond to the impact of moving electrode B closer to, and then further away from, the leftmost (top) extremity of the fracture; and (ii) the decrease and increase that are observed when *n* = 3 for $30 ≤ x ≤ 35$ m correspond to the impact of moving electrode A above the fracture, and then away from the fracture. For the DD array, on the other hand, we see that the single negative peak in $\rho_a^*(x)$ when *n* = 1 is replaced by two negative peaks when *n* = 2 and *n* = 3. These two peaks correspond to measurement positions where the potential- and current-electrode pairs are located above the fracture.

### 3.1.2. Sensitivity to the Fracture Length

The impact of the fracture length on the apparent resistivity anomaly $\rho_a^*(x)$ is studied by increasing the length parameter $\ell$ from 5 m (Cases 1–3) to 10 m (Cases 4–6), and comparing the corresponding curves (Figures 3 and 4, respectively), as well as their magnitudes (Table 2) and widths (Table 3). From these results, we observe that varying the fracture length does not lead to major changes in the apparent resistivity anomaly when $\alpha = 45°$ and $\alpha = 80°$. However, significant differences are observed when $\alpha = 10°$. In particular, extending the fracture length by 5 m leads to an increase in the range of measurement positions over which $\rho_a^*$ varies from its minimum (peak) value to zero when moving towards the right-hand side of the domain.

Focusing on the sub-horizontal fracture ($\alpha = 10°$) case with the WS array, we see that the minimum value of $\rho_a^*$ in Figures 3 and 4 is consistently observed at $x = 24$ m, where the potential electrodes are at their closest position to the top extremity of the fracture. However, the specific behavior that is observed for *n* = 3 and $\ell$ = 5 m, i.e., an increase, decrease, and then increase in $\rho_a^*$ when $x > 24$ m, is not reproduced when $\ell$ = 10 m. In the former case ($\ell$ = 5 m), the shape of the curve after $x = 24$ m results from (i) moving potential electrodes M and N across and away from the fracture, resulting in an increase in $\rho_a^*$ towards zero; (ii) moving current electrode A across the fracture, thereby decreasing $\rho_a^*$; and (iii) moving electrode A away from the fracture, resulting in another increase in $\rho_a^*$ towards zero. This behavior corresponds to two distinct regimes in which the values of $\rho_a^*$ are influenced by the position of the potential electrodes when $24 ≤ x ≤ 29$ m, and the position of electrode A when $x > 29$ m. The transition between these regimes occurs when $\rho_a^* = 0$ at $x = 29$ m, which corresponds to a configuration for which the positions of both current and potential electrodes do not impact the apparent resistivity anomaly. In contrast, when $\ell$ = 10 m, this behavior is no longer observed, because there is no clear transition between the effects of the positions of the potential and current electrodes. That is, when electrode A starts to be located above the fracture, thus impacting the values of $\rho_a^*$, these values are still influenced by the position of the potential electrodes, which are also located
above the fracture. Consequently, there is no measurement position for which the potential and injection electrodes are located on either side of the fracture, and there is not a clear transition between the two regimes where the current and potential electrodes independently affect the values of $\rho_a^*$.

Considering now the results obtained with the DD array for the sub-horizontal fracture case, we observe a nearly symmetric behavior about $x = 24$ m for $l = 5$ m when $n = 2$ and $n = 3$ (Figure 3), which becomes strongly asymmetric for $l = 10$ m (Figure 4). The ER profiling movies demonstrate that the symmetry in the former case arises from the the current and potential electrodes being located on either side of the fracture. The first negative peak in $\rho_a^*$ around $x = 22$ m comes from the potential electrodes moving across the fracture, whereas the second negative peak around $x = 26$ m comes from moving the current electrodes across the fracture. When $l = 10$ m, on the other hand, the behavior of the curves is similar until $x = 22$ m. At this point, the potential electrodes are still located above the fracture when the current electrodes approach the fracture. This implies that there is no measurement position for which the potential and current electrodes are located on either side of the fracture, meaning that we do not observe the strong increase towards zero between the negative peaks that is observed when $l = 5$ m. Instead, $\rho_a^*$ tends to decrease to its minimum value, which is reached around $x = 26$ m and corresponds to the cumulative impact of having both the current and potential electrodes located above the fracture.

## 3.2. Impact of Changes in Material Properties

We saw in our analysis of Cases 1–6 from Table 1 that some single-fracture configurations, under the considered subsurface conditions and electrical properties, can lead to rather substantial apparent resistivity anomalies having magnitudes of up to 7.8% (Table 2). We now wish to explore the impact of the overburden and fracture properties on the anomaly characteristics, with the aim of understanding the effects of these properties on fracture detectability and resolution. To this end, we vary the overburden thickness and electrical resistivity values, as well as fracture aperture and resistivity, for the Case 2 configuration which is considered as a base. We examine the impact of these changes on the resistivity anomaly magnitude and width, whose calculation are again described in Figure 2 with results presented in Tables 2 and 3, respectively. The magnitude, expressed in percent, provides us with a basic measure that allows us to determine whether the fracture would have any chance of being detected in the field. Here, we argue that values upward of 2%, which are in bold in Table 2, represent the bare minimum requirements for detection considering that the standard deviation of typical ER measurement noise is widely assumed to be on the order of 1% of the measured value [39,47–49]. The anomaly width, on the other hand, provides quantitative information on how well the anomaly corresponding to a single fracture might be separated from that of another fracture. Longer electrode arrays, for example, may result in larger anomaly magnitudes (Figures 3 and 4), but will be less useful for resolving multiple fractures because of the greater associated anomaly width.

### 3.2.1. Overburden Properties

We first examine the effect of changing overburden thickness (Cases 2, 7, 8, and 9 from Table 1). Figure 5 shows the apparent resistivity anomaly magnitude and width, plotted for both the WS and DD arrays and for different *n* values, as a function of the overburden thickness *d*, whose value ranges between 0.5 m and 2 m. We see in the figure that doubling the overburden thickness results in an approximate halving of the anomaly magnitude in all cases, meaning that *d* has a strong control on whether the fracture may be detected with ER profiling. For this example, the magnitudes obtained with the DD array are, on average, roughly 2.5 times those obtained with the WS array, suggesting that the former configuration may be better suited to fracture detection. For a fracture whose top extremity is situated at 2 m depth, there is little hope in being able to detect it using the WS array because the anomaly magnitude is well under 2%. Detection using the DD array will also be difficult, as the maximum change across the profile (with *n* = 3) is only 2.4% (Table 2).
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Apparent resistivity anomaly magnitude and width as a function of overburden thickness for the case of a single, 1-mm-aperture fracture having length *l* = 5 m, electrical resistivity ρf = 10 Ω·m, and dipping at 45°. The fracture is hosted in bedrock having resistivity ρb = 10,000 Ω·m, whereas the overburden resistivity is ρo = 1000 Ω·m. The measurement configuration is denoted by WS (Wenner-Schlumberger) or DD (dipole-dipole) with spacing factor *n*. Cases 2, 7, 8, and 9 from Table 1 are considered.</figcaption></figure>

As the spacing factor *n* increases, we see an increase in the ER anomaly magnitude, which is consistent with the results presented in Section 3.1. Note, however, that the factor by which the magnitude changes with *n* increases substantially as the overburden depth increases. Consider, for example, the changes in magnitude for the dipole-dipole array for Cases 7 and 9, which correspond to overburden depths of 0.5 m and 2 m, respectively, between *n* = 1 and *n* = 3. For Case 7, the magnitude changes by a factor of 1.3, whereas for Case 9, it changes by a factor of 2.7 (Table 2). This indicates that there is more benefit derived from increasing *n* when the fracture is located deeper in the subsurface, which we believe results because much of the current travels through in the fracture for *n* = 3 when *d* = 2 m, whereas it will propagate largely below the fracture when *d* = 0.5 m.

Regarding the anomaly width, trends seen in Figure 5 between different *n* values largely reflect the relative differences in array length, as observed previously. This is also the case between the WS and DD arrays, the latter of which is longer for *n* > 1. The WS anomaly width, for example, is consistently higher than that for the DD array, and anomaly widths are observed to consistently increase with increasing *n* value. With increasing overburden thickness (or equivalently fracture depth), we see that the anomaly width increases approximately linearly for the same *n* value and electrode array type. This finding is consistent with geophysical potential field theory [27] which predicts an increasing anomaly width with increasing target depth, and has implications for the ability to resolve multiple fractures, which will evidently decrease as the fracture depth increases.

We next examine the effect of changing overburden resistivity (Cases 2, 10, and 11 from Table 1). Figure 6 shows the calculated anomaly magnitude and width for the three resistivity values considered in our analysis, which were varied logarithmically between 100 Ω·m and 10,000 Ω·m to reflect the wide natural range of this parameter in geological materials [50]. Note that the highest considered value of ρo is equal to the bedrock resistivity ρb (Case 11). As a result, this case may also be viewed as
a situation in which there is no overburden and the fracture is embedded in the bedrock with its top extremity located at depth *d*.

We observe in Figure 6 that, as the overburden resistivity value is increased towards that of the bedrock, there is a marked increase in the apparent resistivity anomaly magnitude, with values ranging from 0.2% to 0.7% when $ρ_o = 100 Ω·m$, to between 4.5% and 14.4% when $ρ_o = 10,000 Ω·m$ (Table 2). The factor by which the magnitude increases is significantly greater in going from $ρ_o = 100 Ω·m$ to $ρ_o = 1000 Ω·m$ (8.5 for the WS configuration with $n = 1$), than when going from $ρ_o = 1000 Ω·m$ to $ρ_o = 10,000 Ω·m$ (3.5 for the same configuration). The strong sensitivity of anomaly magnitude to overburden resistivity is not surprising, as the presence of a low-resistivity surface layer will have the effect of concentrating electric current flow in this layer, thereby prohibiting current from effectively traveling in the underlying bedrock. Practically, this means that ER surveys conducted directly on the bedrock surface, or in highly resistive overburden, will lead to a much easier detection of fractures.

With respect to the anomaly width, Figure 6 confirms that, because of the differences in the lengths of the various electrode arrays considered, the DD array will yield more narrow anomalies than the WS array for the same *n* value, and anomaly widths will increase as *n* increases. With increasing overburden resistivity towards that of the bedrock, however, we also see a decrease in the anomaly width. For the DD configuration with *n* = 1, for example, the calculated width value goes from 16.1 m in Case 10 to 7.8 m in Case 11 (Table 3). The change in anomaly width with changing overburden resistivity value is a consequence of the distortion of current flow patterns that occurs when there is a contrast in resistivity between the overburden and the bedrock [27]. That is, for $ρ_o = 100 Ω·m$, the potential field and current flow lines between electrodes A and B are significantly distorted compared to the homogeneous case ($ρ_o = 10,000 Ω·m$), with the effect that the difference in electric potential measured to the sides of the fracture represents a greater portion of the maximum potential difference value measured over the fracture.

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Apparent resistivity anomaly magnitude and width as a function of overburden resistivity for the case of a single, 1-mm-aperture fracture having length *l* = 5 m, electrical resistivity $ρ_f = 10 Ω·m$, and dipping at 45°. The fracture is hosted in bedrock having resistivity $ρ_b = 10,000 Ω·m$. The overburden thickness is 1 m. The measurement configuration is denoted by WS (Wenner-Schlumberger) or DD (dipole-dipole) with spacing factor *n*. Cases 2, 10, and 11 from Table 1 are considered.</figcaption></figure>
### 3.2.2. Fracture Properties

We now explore how the properties of the fracture impact the characteristics of the observed apparent resistivity anomalies. Figure 7 shows the trends in anomaly magnitude and width as the aperture of the fracture is increased from 1 mm to 4 mm (Cases 2, 12, 13, and 14 from Table 1). As the aperture of the fracture is doubled, we see that the anomaly magnitude is multiplied by a factor of approximately 1.5 for all electrode arrays. Indeed, for the WS-1 configuration, changing the aperture from 1 mm to 2 mm results in a change in anomaly magnitude from 1.4% to 2.2%, whereas changing the aperture from 2 mm to 4 mm results in a further change to 3.1%. This behavior indicates that the aperture of the fracture has a strong control on its detectability using ER methods, and that these methods will detect only subsurface fractures that are wide enough to generate a strong anomaly, while perhaps missing many others that may be important in the context of contaminant migration. In contrast, changing the fracture aperture has no significant effect on the anomaly width (Figure 7). Although more electric current will evidently flow through a larger-aperture fracture, the overall form of the anomaly does not change.

<figure><img src="image_9.png" /></figure>

**Figure 7.** Apparent resistivity anomaly magnitude and width as a function of fracture aperture for the case of a single fracture having length $l = 5$ m, electrical resistivity $ρ_f = 10 \ \Omega \cdot m$, and dipping at 45°. The fracture is hosted in bedrock having resistivity $ρ_b = 10,000 \ \Omega \cdot m$, whereas the overburden resistivity is $ρ_o = 1000 \ \Omega \cdot m$. The overburden thickness is 1 m. The measurement configuration is denoted by WS (Wenner-Schlumberger) or DD (dipole-dipole) with spacing factor $n$. Cases 2, 12, 13, and 14 from Table 1 are considered.

Finally, Figure 8 plots the observed anomaly magnitude and width as a function of the fracture resistivity, the latter of which was varied logarithmically between values of 1 Ω·m and 100 Ω·m (Cases 2, 15, and 16 from Table 1). The fracture resistivity is arguably one of the most important properties in our analysis, as its value can vary dramatically depending on the context of the study. Assuming here that this parameter depends on the filling of the fracture and the effects of surface conductivity, both of which are largely unknown, we consider a wide range of values for which the corresponding ratio between the matrix and fracture resistivities varies between two and four orders of magnitude. We observe in Figure 8 that, as could be expected, the fracture resistivity
has a tremendous impact on whether the fracture might be seen with ER methods. In the case of a four-order-of-magnitude contrast between the matrix and fracture resistivity (Case 15), the ER anomaly magnitudes vary between 4.5% and 22.0% (Table 2), which suggests that the fracture may be detected in consideration of typical levels of measurement error using any of the considered electrode arrays. Conversely, in the case of only a two-order-of-magnitude difference in resistivity (Case 16), the fracture becomes virtually undetectable with anomaly magnitude values between 0.3% and 0.7%. An order-of-magnitude increase in $ρ_f$ has the effect of decreasing the ER anomaly magnitude by a factor that varies from 3.2 to 6.4. With respect to the anomaly width, no significant trend in ER anomaly width with fracture resistivity can be identified.

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Apparent resistivity anomaly magnitude and width as a function of fracture resistivity for the case of a single, 1-mm-aperture fracture having length *l* = 5 m and dipping at 45°. The fracture is hosted in bedrock having electrical resistivity $ρ_b$ = 10,000 Ω·m, whereas the overburden resistivity is $ρ_o$ = 1000 Ω·m. The overburden thickness is 1 m. The measurement configuration is denoted by WS (Wenner-Schlumberger) or DD (dipole-dipole) with spacing factor *n*. Cases 2, 15, and 16 from Table 1 are considered.</figcaption></figure>

## 4. Conclusions

We have sought to quantify in this work the influence of near-surface fractures on ER-profiling experiments through a systematic analysis of simulated data. In contrast to analytical solutions for ER methods that have been developed for simple configurations related to fractures, for example an infinite low-resistivity vertical dyke embedded in a homogeneous half-space [27], the use of numerical modeling in our paper allows for consideration of more complex and realistic scenarios involving various fracture angles and lengths, as well as different background properties including the presence of a layer of overburden. This provides an extended and more comprehensive understanding of the shape of apparent resistivity anomaly curves acquired over fractures, their relationship to fracture and medium properties, and their dependence on the ER survey configuration utilized.

We find in our work that fracture detectability increases as (i) the overburden resistivity approaches that of the bedrock; (ii) the overburden thickness decreases; (iii) the fracture aperture
increases; and (iv) the bedrock-fracture resistivity contrast increases. In this regard, one of the greatest limitations to the detectability of fractures is the presence of a layer of overburden. Any material located above the top extremity of a fracture will reduce its detectability, but low-resistivity material presents a particular challenge because electric current flow becomes concentrated in this upper layer. Another key parameter controlling fracture detection is the fracture resistivity. Whereas much previous work has focused on understanding the electrical properties of rocks and soils, there is little information on what represent realistic values for the electrical resistivity of a fracture, above and below the water table. Consequently, future work should include further understanding of the electrical properties of fractures in order to place our results in the context of real-world applications.

Fracture detectability is also dependent upon the chosen ER survey configuration. However, no single survey type can be identified as best for fracture detection, because results depend strongly on the fracture scenario considered, along with the desired property of interest. For instance, the WS array appears to be more sensitive to the dip angle of a fracture, and should be used if the primary goal is to distinguish between low- and high-angle fractures. However, in many cases the DD array results in greater resistivity anomaly magnitudes, thereby increasing overall fracture detectability. We also find that increasing the spacing factor may increase the detectability of a fracture, as long as total signal strength remains above noise levels. In the examples presented, our analysis is limited to $n \le 3$, but it is possible that greater electrode spacings may improve even further fracture detection.

Finally, although our work analyzes the impact of a large number of parameters over a wide range of values, it is important to note that the presented results and their interpretation are limited to the considered fracture model, i.e., a single thin fracture embedded in homogeneous bedrock and covered by a homogeneous layer of overburden. Future work will investigate the impact of multiple fractures, as well as realistic heterogeneities in the subsurface resistivity distribution, on the observed anomaly curves. Indeed, the latter may result in variability along the ER profile that is significantly greater than the 2% metric that we consider here for fracture detection. Follow-up research will also examine the information brought by other geophysical methods, such as streaming-potential measurements, which may help to detect fractures and estimate their properties.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4441/10/9/1117/s1, Video supplementary: 36 videos that are referenced in Appendix A.

**Author Contributions:** Conceptualization, J.I. and D.R.; Methodology, J.I. and D.R.; Software, S.D. and D.R.; Validation, S.D. and D.R.; Formal Analysis, S.D., J.I., and D.R.; Investigation, S.D.; Resources, J.I.; Writing—Original Draft Preparation, J.I., D.R., and E.V.; Writing Review and Editing, J.I., D.R., and E.V.; Visualization, S.D., J.I., and E.V.; Supervision, J.I. and D.R.; Project Administration, J.I.; Funding Acquisition, J.I.

**Funding:** This research was funded by the Swiss National Science Foundation, grant number 200021-143758.

**Acknowledgments:** We thank Ludovic Baron at the University of Lausanne for many helpful discussions.

**Conflicts of Interest:** The authors declare no conflict of interest. The founding sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; and in the decision to publish the results.

# Appendix A. Description and Interpretation of ER Profiling Movies Supplementary Materials

## Appendix A.1. Overall Description of the Supplementary Material

In order to acquire a full understanding of the results presented in Section 3, movies of the simulated ER experiments are provided as supplementary material for Cases 1–3 from Table 1, for both the Wenner-Schlumberger and dipole-dipole electrode arrays with spacing factor $n = 1$. These movies show, as the considered electrode array is displaced along the profile line, (a) the distribution of electric potential $\phi(x, z)$ throughout the fractured domain; (b) the distribution of the difference in electric potential $\phi^*(x, z)$ between the fractured and corresponding unfractured (or background) domains; and (c) the difference in $\phi^*$ between the measurement-electrode locations, which we denote by $\phi_{MN}^*$ and which represents the potential difference between electrodes M and N relative to the background
value. In (a) and (b), the current and measurement electrodes are shown as red stars and black circles, respectively. In (c), the midpoint of the electrode array, which represents the measurement location, is shown as a black star. Contours corresponding to the electric potential distribution in the unfractured background domain are shown in (a) as dotted white lines for reference.

Overall, the ER profiling movies show changes in $\phi^*(x, z)$ that are localized around the fracture, and whose characteristics vary as the current electrodes are displaced from the left to the right-hand side of the model domain. These changes correspond to a preferential propagation of electric current through the fracture, where the lowest and highest values of $\phi^*(x, z)$, represented in blue and yellow, correspond to the positions at which the electric current enters and exits the fracture, respectively. The impact of these changes in potential on $\phi_{MN}^*$ depends on (i) how the electric current circulates in the fracture, which in turn depends on the position of the current electrodes (A and B) in relation to the fracture position; and (ii) how the potential measurement electrodes (M and N) are located in relation to the fracture position and the current electrodes. The impact of the electrode and fracture positions on $\phi_{MN}^*$ is described in detail below, using the notation $x_J (J = A, B, M, N)$ for the position of electrode J along the x-axis, and $x_{f_i}$ ($i = 1, 2$) for the position along the same axis of the top and bottom extremities of the fracture, respectively.

## Appendix A.2. General Behavior of the Electric Current Flow

In the movies, we observe that when the current electrodes A and B are both located to the left of the fracture (i.e., when $x_A$ and $x_B$ are less than $x_{f_1}$), the electric current will preferentially propagate within the fracture from the bottom to the top extremity. As the array is displaced to the right and electrode B becomes situated above the fracture, however, the exit position of the electric current circulating in the fracture is modified. More precisely, when $x_{f_1} \le x_B \le x_{f_2}$, the current exits the fracture at the position that is closest to electrode B. At the same time, the position of electrode A determines where the electric current will enter the fracture. As $x_A$ approaches $x_{f_1}$, the current will begin to enter the fracture at its top extremity rather than at the bottom. Moving electrode A further towards the right-hand side of the domain and past $x_{f_1}$ results in shifting of this entrance position along the fracture when the dip angle is not too steep (i.e., for $\alpha = 10^\circ$ and $45^\circ$).

For all configurations, the propagation of electric current in the fracture implies that strong negative and positive values of $\phi^*$ will occur around electrodes A and B, respectively, when these electrodes are close to the fracture. The impact of these anomalies on $\phi_{MN}^* = \phi_M^* - \phi_N^*$ depends on the position of electrodes M and N in relation to electrodes A and B. For the Wenner-Schlumberger array with $n = 1$, for example, electrodes M and N are inside and relatively close to electrodes A and B, which implies that $\phi_{MN}^*$ will usually be negative. In contrast, for the dipole-dipole array, electrodes M and N are located outside and to the right of electrodes A and B, and are mostly impacted by the high positive values of $\phi^*$ around electrode B, because it is closest. As electrode M is closer to electrode B than electrode N, $\phi_{MN}^*$ will usually be positive in this case.

## Appendix A.3. Detailed Description for the Wenner-Schlumberger Array

The movie provided for Case 1 ($\alpha = 10^\circ$) for the Wenner-Schlumberger array with $n = 1$ shows that moving the electrodes across the domain from left to right leads to a decrease and then an increase in $\phi_{MN}^*$, with the minimum value being reached at a measurement location of $x = 24$ m. In the beginning of the movie, the electrode array is located close to the left-hand side of the domain where the influence of the fracture is minimal, and, thus, $\phi_{MN}^*$ is close to zero. As the array is moved closer to the fracture, however, $\phi_{MN}^*$ becomes increasingly negative, with $x = 24$ m representing the location where electrodes A and B are close to the left and right fracture extremities, respectively. This allows for a significant proportion of electric current to travel through the fracture, and results in a strong anomaly. As the array moves to the right away from this location, the influence of the fracture decreases and $\phi_{MN}^*$ moves again towards zero. The fact that the fracture is not completely horizontal means that the measured anomaly in (c) is not perfectly symmetric around $x = 24$ m.
From the movies for Cases 2 and 3 where the fracture dip angle is increased (α = 45° and 80°), we observe that the behavior of $\phi_{MN}^*$ is different for $x < 24$ m in comparison to Case 1 described above. Here, $\phi_{MN}^*$ is seen to decrease from the left-hand side of the domain to position $x = 20$ m, increase between $x = 20$ m and $x \approx 22$ m, and then decrease again until $x = 24$ m. The electric current behavior explaining this series of decreases and increases is easily identified for α = 80°, where we observe that $\phi_{MN}^*$ tends to zero when $x = 22$ m. At this position, current electrodes A and B are located at roughly equal distances on opposite sides of the sub-vertical fracture, whose presence only slightly modifies the electric current propagation because the flow is predominantly orthogonal to the fracture. That is, the presence of the sub-vertical fracture for this measurement location does not allow for any preferential propagation of electric current.

In contrast, we observe that the behavior of $\phi_{MN}^*$ for $x \ge 24$ m is similar for α = 10°, 45° and 80°. When the electrode array is centered on $x = 24$ m, electrodes A and M are located on either side of the top extremity of the fracture. The value of $\phi_{MN}^*$ is lowest here because: (i) The small distance between electrode A and the top (and leftmost) extremity of the fracture means that a significant proportion of electric current will enter the fracture and travel downwards and to the right towards electrode B; and (ii) the small distance between electrode M and the top extremity of the fracture implies that the strong negative values of $\phi^*$ will be sensed. These results importantly demonstrate that the minimum value of $\phi_{MN}^*$ observed along the ER profile line with the Wenner-Schlumberger array will indicate the position of the top extremity of the fracture, independent of the dip angle or position of the bottom extremity of the fracture.

### Appendix A.4. Detailed Description for the Dipole-Dipole Array

With the dipole-dipole array, we see from the ER-profiling movies for Cases 1–3 that, as the array is moved from left to right across the model domain, $\phi_{MN}^*$ tends from zero to a weak negative value before the fracture, then a strong positive value over the fracture, followed by another weak negative value after the fracture before returning to zero. As explained previously, when the current electrodes A and B are located to the left of the top extremity of the fracture (i.e., for $x_A$ and $x_B$ less than $x_{f1}$), the electric current preferentially circulates from the bottom to the top of the fracture. When the potential electrodes M and N are also located to the left of this position (i.e., for $x_M$ and $x_N$ less than $x_{f1}$), slightly larger values of $\phi^*$ occur in N than in M, resulting in weak negative values for $\phi_{MN}^*$. These negative values peak in the three movies at $x = 19$ m, because this is the closest location to the fracture where all of the electrodes are still to the left of its top extremity.

As the dipole-dipole array is moved further to the right, electrodes M and N pass over the top extremity of the fracture, and the positive values of $\phi^*$ around this location, due to the exit of electric current, now result in a positive value for $\phi_{MN}^*$. This positive value becomes larger as the current electrodes approach (and move over) position $x_{f1}$, and eventually peaks at $x = 24$ m for the 10° fracture, and $x = 22$ m for the 45° and 80° fractures, indicating that the fracture dip angle has an influence on the position at which the anomaly is maximized. As electrode A passes over the top extremity of the fracture, preferential current flow in the fracture is reversed, now traveling from top to bottom.

Moving the dipole-dipole array even further to the right of the fracture, weak negative values of $\phi_{MN}^*$ are again observed. For the 10° fracture, this occurs after $x = 30$ m when all of the electrodes are located to the right of the bottom fracture extremity (i.e., for $x_A, x_B, x_M$ and $x_N$ greater than $x_{f2}$). In this case, there is minimal preferential circulation of electric current in the fracture, and the presence of the fracture only slightly modifies the electric potential distribution. For the 45° and 80° fractures, on the other hand, the negative values are observed after approximately $x = 25$ m, as the bottom extremity of the fracture has less of an impact because of its deeper position.
References

1. European Commission. *Groundwater Protection in Europe. The New Groundwater Directive—Consolidating the EU Regulatory Framework*; Office for Official Publications of the European Communities: Luxembourg, 2008.

2. Margat, J.; van der Gun, J. *Groundwater Around the World: A Geographic Synopsis*, 1st ed.; CRC Press: London, UK, 2013.

3. Singhal, B.; Gupta, R.; Gupta, R. *Applied Hydrogeology of Fractured Rocks*, 2nd ed.; Kluwer Academic Publishers: Roorkee, India, 1999.

4. Gallegos, T.J.; Varela, B.A.; Haines, S.S.; Engle, M.A. Hydraulic fracturing water use variability in the United States and potential environmental implications. *Water Resour. Res.* **2015**, *51*, 5839–5845. [CrossRef] [PubMed]

5. Harvey, O.R.; Qafoku, N.P.; Cantrell, K.J.; Lee, G.; Amonette, J.E.; Brown, C.F. Geochemical implications of gas leakage associated with geologic CO₂ Storage—A qualitative review. *Environ. Sci. Technol.* **2013**, *47*, 23–36. [CrossRef] [PubMed]

6. Vidic, R.D.; Brantley, S.L.; Vandenbossche, J.M.; Yoxtheimer, D.; Abad, J.D. Impact of shale gas development on regional water quality. *Science* **2013**, *340*, 1235009. [PubMed]

7. Berkowitz, B. Characterizing flow and transport in fractured geological media: A review. *Adv. Water Resour.* **2002**, *25*, 861–884. [CrossRef]

8. Bonnet, E.; Bour, O.; Odling, N.E.; Davy, P.; Main, I.; Cowie, P.; Berkowitz, B. Scaling of fracture systems in geological media. *Rev. Geophys.* **2001**, *39*, 347–383. [CrossRef]

9. Davy, P.; Le Goc, R.; Darcel, C.; Bour, O.; De Dreuzy, J.R.; Munier, R. A likely universal model of fracture scaling and its consequence for crustal hydromechanics. *J. Geophys. Res. Solid Earth* **2010**, *115*. [CrossRef]

10. National Academies of Sciences, Engineering, and Medicine. *Characterization, Modeling, Monitoring, and Remediation of Fractured Rock*; The National Academies Press: Washington, DC, USA, 2015. [CrossRef]

11. Herwanger, J.; Worthington, M.; Lubbe, R.; Binley, A.; Khazanehdari, J. A comparison of cross-hole electrical and seismic data in fractured rock. *Geophys. Prospect.* **2004**, *52*, 109–121. [CrossRef]

12. Pytharouli, S.; Lunn, R.; Shipton, Z.; Kirkpatrick, J.; do Nascimento, A. Microseismicity illuminates open fractures in the shallow crust. *Geophys. Res. Lett.* **2011**, *38*, L02402. [CrossRef]

13. Dorn, C.; Linde, N.; Doetsch, J.; Borgne, T.L.; Bour, O. Fracture imaging within a granitic rock aquifer using multiple-offset single-hole and cross-hole GPR reflection data. *J. Appl. Geophys.* **2012**, *78*, 123–132. [CrossRef]

14. Tsoflias, G.P.; Halihan, T.; Sharp, J.M. Monitoring pumping test response in a fractured aquifer using ground-penetrating radar. *Water Resour. Res.* **2001**, *37*, 1221–1229. [CrossRef]

15. Lane, J.; Haeni, F.; Watson, W. Use of a square-array direct-current resistivity method to detect fractures in crystalline bedrock in New-Hampshire. *Ground Water* **1995**, *33*, 476–485. [CrossRef]

16. Robert, T.; Caterina, D.; Deceuster, J.; Kaufmann, O.; Nguyen, F. A salt tracer test monitored with surface ERT to detect preferential flow and transport paths in fractured/karstified limestones. *Geophysics* **2012**, *77*, B55–B67. [CrossRef]

17. Roubinet, D.; Linde, N.; Jougnot, D.; Irving, J. Streaming potential modeling in fractured rock: Insights into the identification of hydraulically active fractures. *Geophys. Res. Lett.* **2016**, *43*, 4937–4944. [CrossRef]

18. Wishart, D.N.; Slater, L.D.; Gates, A.E. Fracture anisotropy characterization in crystalline bedrock using field-scale azimuthal self potential gradient. *J. Hydrol.* **2008**, *358*, 35–45. [CrossRef]

19. Marescot, L.; Monnet, R.; Chapellier, D. Resistivity and induced polarization surveys for slope instability studies in the Swiss Alps. *Eng. Geol.* **2008**, *98*, 18–28. [CrossRef]

20. Schmutz, M.; Ghorbani, A.; Vaudelet, P.; Revil, A. Spectral induced polarization detects cracks and distinguishes between open- and clay-filled fractures. *J. Env. Eng. Geophys.* **2011**, *16*, 85–91. [CrossRef]

21. Donadille, J.M.; Al-Ofi, S. Crosswell electromagnetic response in a fractured medium. *Geophysics* **2012**, *77*, D53–D61. [CrossRef]

22. Steelman, C.; Parker, B.; Kennedy, C. Evaluating local-scale anisotropy and heterogeneity along a fractured sedimentary bedrock river using EM azimuthal resistivity and ground-penetrating radar. *J. Appl. Geophys.* **2015**, *116*, 156–166. [CrossRef]

23. Berryman, J.G.; Hoversten, G.M. Modelling electrical conductivity for earth media with macroscopic fluid-filled fractures. *Geophys. Prospect.* **2013**, *61*, 471–493. [CrossRef]
24. Jinsong, S.; Benyu, S.; Naichuan, G. Anisotropic characteristics of electrical responses of fractured reservoir with multiple sets of fractures. *Pet. Sci.* **2009**, *6*, 127–138. [CrossRef]

25. Taylor, R.W.; Fleming, A.H. Characterizing jointed systems by azimuthal resistivity surveys. *Ground Water* **1988**, *26*, 464–474. [CrossRef]

26. Everett, M.E. *Near-Surface Applied Geophysics*, 1st ed.; Cambridge University Press: Cambridge, UK, 2013.

27. Telford, W.M.; Geldart, L.P.; Sheriff, R.E. *Applied Geophysics*; Cambridge University Press: Cambridge, UK, 1990.

28. Brown, S.R. Transport of fluid and electric current through a single fracture. *J. Geophys. Res. Solid Earth* **1989**, *94*, 9429–9438. [CrossRef]

29. Van Siclen, C.D. Equivalent channel network model for permeability and electrical conductivity of fracture networks. *J. Geophys. Res. Solid Earth* **2002**, *107*. [CrossRef]

30. Greenhalgh, S.; Zhou, B.; Greenhalgh, M.; Marescot, L.; Wiese, T. Explicit expressions for the Fréchet derivatives in 3D anisotropic resistivity inversion. *Geophysics* **2009**, *74*, F31–F43. [CrossRef]

31. Herwanger, J.V.; Pain, C.C.; Binley, A.; De Oliveira, C.R.E.; Worthington, M.H. Anisotropic resistivity tomography. *Geophys. J. Int.* **2004**, *158*, 409–425. [CrossRef]

32. Pain, C.C.; Herwanger, J.V.; Saunders, J.H.; Worthington, M.H.; de Oliveira, C.R. Anisotropic resistivity inversion. *Inverse Probl.* **2003**, *19*, 1081. [CrossRef]

33. Busby, J. The effectiveness of azimuthal apparent-resistivity measurements as a method for determining fracture strike orientations. *Geophys. Prospect.* **2000**, *48*, 677–695. [CrossRef]

34. Boadu, F.; Gyamfi, J.; Owusu, E. Determining subsurface fracture characteristics from azimuthal resistivity surveys: A case study at Nsawam, Ghana. *Geophysics* **2005**, *70*, B35–B42. [CrossRef]

35. Ritzi, R.W.; Andolsek, R.H. Relation between anisotropic transmissivity and azimuthal resistivity surveys in shallow, fractured, carbonate flow systems. *Ground Water* **1992**, *30*, 774–780. [CrossRef]

36. Skinner, D.; Heinson, G. A comparison of electrical and electromagnetic methods for the detection of hydraulic pathways in a fractured rock aquifer, Clare Valley, South Australia. *Hydrog. J.* **2004**, *12*, 576–590. [CrossRef]

37. Steinich, B.; Marin, L.E. Hydrogeological investigations in northwestern Yucatan, Mexico, using resistivity surveys. *Groundwater* **1996**, *34*, 640–646. [CrossRef]

38. Lesparre, N.; Boyle, A.; Grychtol, B.; Cabrera, J.; Marteau, J.; Adler, A. Electrical resistivity imaging in transmission between surface and underground tunnel for fault characterization. *J. Appl. Geophys.* **2016**, *128*, 163–178. [CrossRef]

39. Robert, T.; Dassargues, A.; Brouyère, S.; Kaufmann, O.; Hallet, V.; Nguyen, F. Assessing the contribution of electrical resistivity tomography (ERT) and self-potential (SP) methods for a water well drilling program in fractured/karstified limestones. *J. Appl. Geophys.* **2011**, *75*, 42–53. [CrossRef]

40. Robinson, J.; Johnson, T.; Slater, L. Evaluation of known-boundary and resistivity constraints for improving cross-borehole DC electrical resistivity imaging of discrete fractures. *Geophysics* **2013**, *78*, D115–D127. [CrossRef]

41. Sharma, S.; Baranwal, V. Delineation of groundwater-bearing fracture zones in a hard rock area integrating very low frequency electromagnetic and resistivity data. *J. Appl. Geophys.* **2005**, *57*, 155–166. [CrossRef]

42. Torleif, D.; Bing, Z. A numerical comparison of 2D resistivity imaging with 10 electrode arrays. *Geophys. Prospect.* **2004**, *52*, 379–398. [CrossRef]

43. Friedel, S. Resolution, stability and efficiency of resistivity tomography estimated from a generalized inverse approach. *J. Appl. Geophys.* **2003**, *153*, 305–316. [CrossRef]

44. Loke, M. *Tutorial: 2-D and 3-D Electrical Imaging Surveys*; Geotomo Software: Penang, Malaysia, 1996–2004.

45. Roubinet, D.; Irving, J. Discrete-dual-porosity model for electric current flow in fractured rock. *J. Geophys. Res. Solid Earth* **2014**, *119*, 767–786. [CrossRef]

46. Caballero Sanz, V.; Roubinet, D.; Demirel, S.; Irving, J. 2.5-D discrete-dual-porosity model for simulating geoelectrical experiments in fractured rock. *Geophys. J. Int.* **2017**, *209*, 1099–1110. [CrossRef]

47. Koestel, J.; Kemna, A.; Javaux, M.; Binley, A.; Vereecken, H. Quantitative imaging of solute transport in an unsaturated and undisturbed soil monolith with 3-D ERT and TDR. *Water Resour. Res.* **2008**, *44*. [CrossRef]
48. Slater, L.; Binley, A.; Daily, W.; Johnson, R. Cross-hole electrical imaging of a controlled saline tracer injection. *J. Appl. Geophys.* **2000**, *44*, 85–102. [CrossRef]

49. Zhou, B.; Dahlin, T. Properties and effects of measurement errors on 2D resistivity imaging surveying. *Near Surf. Geophys.* **2003**, *1*, 105–117. [CrossRef]

50. Schön, J.H. *Physical Properties of Rocks: Fundamentals and Principles of Petrophysics*, 2nd ed.; Elsevier: New York, NY, USA, September 2015.

<figure><img src="image_11.png" /></figure>

© 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).