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

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.


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

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 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 ρ 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 where φ MN = φ M − φ 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 The geometric factor for the dipole-dipole array is [27] K = −π n (n + 1) (n + 2) a.
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, ρ * a , as the relative difference (in %) between the apparent resistivity of the considered fractured domain, ρ a , and that of the background unfractured domain, ρ u a . That is, Consideration of ρ * a in our analysis, as opposed to ρ 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][48][49]. For each fracture configuration, electrode array, and n value considered in our study, we plot ρ * 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 ρ * a values, and the anomaly width as the distance between the first and last points where |ρ * 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][13][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.  Diagram showing how the apparent resistivity anomaly, ρ * 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.

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 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.

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, ρ * a . Figures 3 and 4 show ρ * 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, φ * (x, z). This in turn leads to a detailed understanding of the measured potential difference between electrodes M and N relative to background, φ * MN , along the profile. The conclusions drawn in Appendix A are used to help interpret the behavior of ρ * 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. . Apparent resistivity anomaly ρ * a (in %) as a function of horizontal position for a single, 1-mm-aperture fracture having length = 5 m, electrical resistivity ρ f = 10 Ω·m, and hosted in bedrock having resistivity ρ b = 10, 000 Ω·m. The fracture is covered by a 1 m thick layer of overburden, having resistivity ρ 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).

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 ρ * 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 φ * MN (and thus ρ * 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 ρ * a is seen to occur at x = 24 m for all fracture angles (α = 10 • , 45 • , and 80 • ). For the sub-horizontal fracture configuration (α = 10 • ), ρ * 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 (α = 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 • . . Apparent resistivity anomaly ρ * a (in %) as a function of horizontal position for a single, 1-mm-aperture fracture having length = 10 m, electrical resistivity ρ f = 10 Ω·m, and hosted in bedrock having resistivity ρ b = 10, 000 Ω·m. The fracture is covered by a 1 m thick layer of overburden, having resistivity ρ o = 1000 Ω·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).
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 φ * 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 ρ * 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 φ MN (and thus ρ 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 (α = 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 α = 10 • , the position at which the minimum value of ρ * 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 ρ * 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.

Sensitivity to the Fracture Length
The impact of the fracture length on the apparent resistivity anomaly ρ * a (x) is studied by increasing the length parameter 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 α = 45 • and α = 80 • . However, significant differences are observed when α = 10 • . In particular, extending the fracture length by 5 m leads to an increase in the range of measurement positions over which ρ * 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 (α = 10 • ) case with the WS array, we see that the minimum value of ρ * 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 = 5 m, i.e., an increase, decrease, and then increase in ρ * a when x > 24 m, is not reproduced when = 10 m. In the former case ( = 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 ρ * a towards zero; (ii) moving current electrode A across the fracture, thereby decreasing ρ * a ; and (iii) moving electrode A away from the fracture, resulting in another increase in ρ * a towards zero. This behavior corresponds to two distinct regimes in which the values of ρ * 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 ρ * 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 = 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 ρ * 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 ρ * 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 = 5 m when n = 2 and n = 3 (Figure 3), which becomes strongly asymmetric for = 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 ρ * 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 = 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 = 5 m. Instead, ρ * 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.

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][48][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.

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).  Table 1 are considered.
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.  Table 1 are considered.

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.  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.  Table 1 are considered.

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 ≤ 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. 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 φ * (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 φ * (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 φ * 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 φ * 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 ≤ x B ≤ 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 α = 10 • and 45 • ).
For all configurations, the propagation of electric current in the fracture implies that strong negative and positive values of φ * will occur around electrodes A and B, respectively, when these electrodes are close to the fracture. The impact of these anomalies on φ * MN = φ * M − φ * 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 φ * 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 φ * around electrode B, because it is closest. As electrode M is closer to electrode B than electrode N, φ * MN will usually be positive in this case.

Appendix A.3. Detailed Description for the Wenner-Schlumberger Array
The movie provided for Case 1 (α = 10 • ) 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 φ * 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, φ * MN is close to zero. As the array is moved closer to the fracture, however, φ * 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 φ * 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 φ * MN is different for x < 24 m in comparison to Case 1 described above. Here, φ * 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 ≈ 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 φ * 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 φ * MN for x ≥ 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 φ * 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 φ * will be sensed. These results importantly demonstrate that the minimum value of φ * 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, φ * 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 f 1 ), 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 f 1 ), slightly larger values of φ * occur in N than in M, resulting in weak negative values for φ * 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 φ * around this location, due to the exit of electric current, now result in a positive value for φ * MN . This positive value becomes larger as the current electrodes approach (and move over) position x f 1 , 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 φ * 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 f 2 ). 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.