Evaluating local-scale anisotropy and heterogeneity along a fractured sedimentary bedrock river using EM azimuthal resistivity and ground-penetrating radar

Journal of Applied Geophysics 116 (2015) 156–166

Contents lists available at ScienceDirect

Journal of Applied Geophysics journal homepage: www.elsevier.com/locate/jappgeo

Evaluating local-scale anisotropy and heterogeneity along a fractured sedimentary bedrock river using EM azimuthal resistivity and ground-penetrating radar C.M. Steelman a,b,⁎, B.L. Parker a,b,c, C.S. Kennedy a,c a b c

G360 Centre for Applied Groundwater Research, University of Guelph, Canada School of Engineering, University of Guelph, Canada School of Environmental Sciences, University of Guelph, Canada

a r t i c l e

i n f o

Article history: Received 30 September 2014 Received in revised form 10 February 2015 Accepted 2 March 2015 Available online 4 March 2015 Keywords: EM azimuthal resistivity Ground-penetrating radar Fractured rock Bedrock rivers Anisotropy

a b s t r a c t Fractured sedimentary bedrock rivers exhibit complicated flow patterns controlled by the geometry, extent and connectivity of fractures and dissolution-enhanced conduits. In a bedrock river environment these features variably connect surface water to groundwater. Given the nature of discrete fracture and conduit networks, these flow systems can be anisotropic and heterogeneous over a wide range of scales. Portable and non-invasive geophysical methods are ideal for the initial characterization of shallow hydrogeologic systems in ecologically sensitive environments. Here, we evaluate the utility of the electromagnetic (EM) azimuthal resistivity method for characterizing shallow vertical and subvertical fracture set orientations in the presence of localized dissolution-enhanced features near a sedimentary bedrock river located in Southern Ontario, Canada. Multiple EM coil spacings ranging from 1 to 10 m in vertical and horizontal dipole modes were applied at two locations within a 10 × 25 m study plot; azimuthal rotations were conducted using symmetric and asymmetric geometries. The observed anisotropic response was evaluated using 100 MHz ground-penetrating radar (GPR) measurements collected across the study plot area. A joint analysis of EM and GPR data revealed that fracture-based anisotropy can be identified in the presence of local heterogeneities (i.e., dissolution-enhanced features) at scales considerably smaller than those of previous studies. These non-invasive geophysical data sets can be used to optimize the design of multi-borehole groundwater monitoring stations along bedrock river channels, such that angled boreholes could be emplaced and orientated to intercept major fracture networks and local dissolution features. These data could also improve 3-D fracture network conceptualizations utilizing discrete information from angled and vertical boreholes. © 2015 Elsevier B.V. All rights reserved.

1. Introduction Portable and non-invasive geophysical methods can provide valuable insight into hydrologic processes in ecologically-sensitive environments. Given the dynamic, complex and interdependent nature of these ecohydrological systems, critical interfaces where groundwater interacts with surface water, such as the hyporheic zone, are more likely to be impacted by groundwater pumping, surface water withdrawal, and contamination (Winter et al., 1998). While the importance of hyporheic processes has been extensively examined (e.g., Vannote et al., 1980; Meyer, 1997; Woessner, 2000), recent advances in geophysical sensors, portable drills, and depth-discrete multilevel monitoring systems have enabled the collection of higher-resolution information with reduced impact to surrounding ecosystems (e.g., Conant et al., 2004; Hatch ⁎ Corresponding author at: University of Guelph, 50 Stone Road East, Guelph, Ontario N1G 2W1, Canada. E-mail address: [email protected] (C.M. Steelman).

http://dx.doi.org/10.1016/j.jappgeo.2015.03.003 0926-9851/© 2015 Elsevier B.V. All rights reserved.

et al., 2006; Schmidt et al., 2007; Gabrielli and McDonnell, 2011; Ward et al., 2012). These developments have improved our capacity to study dynamic processes. A river's inherent sensitivity to anthropogenic impacts makes minimally-invasive geophysical methods well-suited for conceptualization of complex hydrogeologic systems. While the application of geophysical methods for enhanced hydrologic process characterization in alluvial rivers has received reasonable attention (e.g., Birkhead et al., 1996; Crook et al., 2008; Ward et al., 2010; Musgrave and Binley, 2011; Wojnar et al., 2013), there has been very limited studies of rivers that flow directly along a fractured sedimentary bedrock surface (e.g., Oxtobee and Novakowski, 2002). By definition, a bedrock river flows directly on an exposed bedrock streambed, without the attenuative effects of sand or sediment overburden. Because bedrock rivers are more challenging to instrument (Tinker and Wohl, 1998b), our hydrogeologic understanding of fractured rock riverbeds has become biased toward alluvial river conceptual models (Woessner, 1998; Winter et al., 1998).

C.M. Steelman et al. / Journal of Applied Geophysics 116 (2015) 156–166

The majority of bedrock river studies focus on fluvial geomorphological processes in the bedrock channel (e.g., Tinker and Wohl, 1998a; Hodge et al., 2011) and comparisons between alluvial and bedrock constrained multi-channel river networks (Meshkova and Carling, 2013). Oxtobee and Novakowski (2002) conducted a comprehensive field investigation of a bedrock creek using geographical, hydrological and geochemical data to evaluate groundwater–surface water interaction between Twenty Mile Creek and the local aquifer near Smithville, ON, Canada. They concluded that the interaction between the river and aquifer was limited with N95% of the groundwater underflowing the creek during baseflow condition. Within their study area, the absence of strong vertical fracturing limited groundwater exchange in the vicinity of the creek. Given the importance of bedrock rivers to groundwater supply and the potential threats to surrounding ecological environments, there is an intrinsic need to develop integrated investigative field applications capable of providing high-resolution spatial information about the fracture and dissolution-enhanced conduit networks beneath these riverbeds (Hancock, 2002). Non-invasive geophysical data could inform localized hydrogeologic observations, and collectively advance our static and dynamic conceptual understanding of flow and transport along discrete fracture networks which variably connect groundwater to surface water. The Eramosa River, in Ontario, Canada, exhibits a bedrock river environment where the utility of portable geophysical methodologies can be tested, along with the capacity to detect shallow fracture networks. This river flows along a variably exposed, highly-fractured and fragmented Silurian dolostone floodplain with abundant karst features visible at the surface (Kunert et al., 1998). Groundwater and surface water interaction is hypothesized to occur through discrete vertical fractures and dissolution-enhanced conduits (e.g., Steelman et al., 2015); the distribution and extent of these flow paths is expected to vary with depth and along the reach of the river. Effective conceptualization of dynamic hyporheic zone processes requires direct measurement of groundwater and surface water conditions; however, optimized emplacement of groundwater monitoring stations in a fractured sedimentary bedrock river requires prior knowledge of shallow fracture networks and suspected karst features at depth. Non-invasive geophysical information could be used to locate and design a distribution of groundwater monitoring stations. The objective of this study is to assess the capacity of the azimuthal electromagnetic (EM) resistivity method for the detection of local-scale vertical and subvertical fractures and features possibly affecting groundwater–surface water exchange along a sedimentary bedrock river. To achieve this, a series of azimuthal resistivity surveys were conducted at two stations located within a 10 × 25 m intensive study plot. The azimuthal measurements were assessed using full-resolution groundpenetrating radar (GPR) measurements using 100 MHz antennas to support interpretation of fracture sets and possible dissolution features or heterogeneities impacting the anisotropic response.

157

flow patterns, which may be influenced by dissolution processes. Integrated geophysical applications yielding high-resolution spatial information about the fracture geometry, distribution, frequency and intensity can improve conceptualizations of hydrogeologic flow systems across larger spatial scales, thereby advancing our understanding of shallow hydrologic processes along bedrock riverbeds. The primary interaction of groundwater and surface water along a sedimentary rock riverbed will be controlled by fractures and conduits connecting the rockbed interface to subsurface hydrostratigraphic units. Geophysical research on shallow bedrock flow systems has focused on the characterization of major dissolution features and structures, such as epikarst, caverns, and sinkholes (e.g., Batayneh et al., 2002; Al-fares et al., 2002; Kruse et al., 2006), and local-to-regional scale fracture characteristics (e.g., Taylor and Fleming, 1988; Lane et al., 1995; Busby, 2000). More recently, Skinner and Heinson (2004) developed a set of borehole-to-surface direct current and electromagnetic induction methods to delineate laterally extensive preferential hydraulic pathways in a shallow fractured bedrock environment. Using the azimuthal resistivity method, which detects directional changes in apparent horizontal electrical conductivity, they were able to identify the azimuth of maximum hydraulic connectivity along vertical and sub-vertical conduits. Although the authors noted the presence of inhomogeneities in the rock, and expressed the potential impact on the interpretation of homogeneous fracture anisotropy, they did not examine the nature of the inhomogeneities and their relationship to the local-scale fracture network. Although regional-scale fracture properties may be useful for largescale conceptual models (Taylor and Fleming, 1988; Busby, 2000), these data sets are of limited value in the evaluation of local-scale processes, such as water interactions at the hyporheic interface. The highresolution detail of field-based spatial representations of vertical fracture networks, bedding plane geometries and localized dissolutionenhanced conduits is fundamental to the conceptualization of fractureflow in close proximity to the bedrock riverbed interface. The identification of heterogeneous anisotropy, caused by the presence of isolated flow features in a fractured system, can provide valuable insights about potential preferential pathways. Given the success of previous electrical investigations on the detection of anisotropy in heterogeneous environments (e.g., Slater et al., 1998; Watson and Barker, 1999; Skinner and Heinson, 2004; Watson and Barker, 2005; Watson and Barker, 2010; Yeboah-Forson and Whitman, 2014) we have applied the EM azimuthal resistivity method to the conceptualization of vertical fracture networks in a shallow bedrock river environment. However, unlike previous studies we have applied multiple coil spacings, magnetic dipole orientations, and survey configurations over an intensive 250 m2 study area to characterize horizontal fracture anisotropy, while also evaluating the influence of localscale inhomogeneities (e.g., cavities, dissolution-enhanced features) on the anisotropic signal. 2.2. Utility of azimuthal resistivity method

2. Background 2.1. Fracture networks and dissolution features along sedimentary bedrock rivers Surface and subsurface connectivity along rivers that flow directly on sedimentary bedrock surfaces with exposed vertical fracture networks will exhibit very complicated flow patterns relative to their alluvial counterparts largely due to the nature of their effective porosities (i.e., interconnected pores vs. discrete fracture networks) (Tinker and Wohl, 1998b). Consequently, bedrock rivers can be anisotropic and heterogeneous over a wide range of spatial scales (Fig. 1). Primary fracture orientation, spacing, aperture, length and connectivity can vary over tens (local) to thousands (regional) of meters. This extreme variation in contributing factors often results in anisotropic and heterogeneous

The electrical resistivity of fractured bedrock composed of vertical or subvertical joint sets with a preferred orientation will vary depending on the direction of the applied electrical current (Habberjam, 1975; Taylor and Fleming, 1988). Measured apparent electrical resistivity will be highest when the current flows parallel to the dominant fracture sets and lowest when the current flows orthogonal to the fracture sets. This so-called paradox of anisotropy can be explained by an increase in the density of the injected or induced electric current that is forced to flow along the strike of the anisotropy (Wasscher, 1961; Al-Garni and Everett, 2003). For decades researchers have been advancing the application of the azimuthal resistivity method to jointed bedrock systems to better understand hydraulic transmissivity (e.g., Ritzi and Andolsek, 1992), hydraulic connectivity (e.g., Skinner and Heinson, 2004; Abdullahi and

158

C.M. Steelman et al. / Journal of Applied Geophysics 116 (2015) 156–166

Fig. 1. Conceptualization of vertical fracture network underlying a bedrock river where groundwater (GW) and surface water (SW) interact along discrete fractures. Laterally extensive sub-vertical to horizontal bedding planes and fractures have been removed for clarity. Different effective sampling volumes (marked by sampling area) are shown to illustrate the spatial dependency of the geophysical response to fracture and dissolution features.

Osazuwa, 2011), geological attributes (e.g., Lane et al., 1995; Busby, 2000; Matias, 2002), and hazard assessment (e.g., Busby and Jackson, 2006). Recently, azimuthal self-potential gradient measurements have been applied to characterize preferential groundwater flow paths (e.g., Sandberg et al., 2002; Wishart et al., 2008). Much of this work has focused on the use of galvanic methods which require insertion of electrodes into the ground to inject a direct current. While the direct current method has been highly effective for delineating fracture anisotropy, due to its high signal-to-noise ratio, applications have been limited to environments where: (1) surface bedrock makes for difficult electrode implant; and (2) restricted accessibility or workable space presents further installation challenges, given the necessary extent of collinear lines (N30 m). Frequency-domain electromagnetic (EM) induction methods provide a measure of apparent electrical conductivity using separate transmitting and receiving coils. These methods have much smaller footprints relative to conventional galvanic electrical methods and do not require direct contact with the subsurface; therefore the method is much more suitable for application in ecologically sensitive and spatially restricted areas with an exposed bedrock surface. Although this method is typically accompanied by higher signal noise relative to galvanic methods (Slater et al., 1998), it is much more efficient and simpler to apply over a wide-range of environments with variable ground conditions and workable space.

layer of Lower Goat Island Formation (1.6–1.8 m) which is underlain by a vuggy Gasport Formation. The Gasport Formation becomes highly weathered and karstic at a depth of 12 m (Steelman et al., 2015). These cores exhibited abundant bedding planes with rust-to-black colored staining indicating groundwater flow. Fig. 2a shows the position of a 10 × 25 m study plot (250 m2) located adjacent to a section of the Eramosa River. Geophysical measurements were conducted during low-stage conditions, as depicted in Fig. 2, which enabled imaging into a representative portion of the floodplain which is submerged during high-stage conditions. Topography over the study plot is reasonably flat with small (relative to the GPR footprint) localized depressions in the north end which show evidence of rock dissolution at surface. A portion of the bedrock surface consisted of a thin veneer of weathered rock fragments overlying competent

2.3. Study site description The exposed bedrock surface along the Eramosa River has been described as a medium to thick bedded, massive, bioclastic, and biohermal dolomite of the Amabel formation (Kunert et al., 1998); recent nomenclature introduced by Brunton (2009) describe these units as the Lower Goat Island and underlying Gasport Formations. These dolostone units are susceptible to the development of epikarst, as well as vuggy and cavernous porosity, which have resulted in a highly-transmissive regional bedrock aquifer approximately 15 m below the Eramosa River in the immediate study area. Many of the fractures have been enhanced through dissolution, which significantly increased the bulk porosity and permeability of these rocks. Core logs collected approximately 250 m south of the immediate study area indicate that the river is flowing along a thin

Fig. 2. a) View of study plot over a section of exposed bedrock along the east bank of the Eramosa River. Geophysical surveys were conducted during low-stage conditions as depicted. The topography is nearly-flat lying across the test plot, with a slight dip toward the water in the south-west quadrant. b) Close-up of the bedrock surface within the northwest quadrant of the test plot. Note the vertical dissolution-enhanced fractures roughly orientated north–northwest and east–west. The study plot is located at 43°34′53″N, 80°08′49″W.

C.M. Steelman et al. / Journal of Applied Geophysics 116 (2015) 156–166

rock with exposed fracture joints, many of which show signs of dissolution (Fig. 2b). Steelman et al. (2015) presented evidence of numerous depth-discrete vertical and sub-vertical fractures, multiple horizontal to slightly dipping bedding planes, and karst features in the upper 20 m along this portion of the river. 3. Methodology 3.1. Processing and interpretation of azimuthal resistivity data While a comprehensive review of the EM azimuthal resistivity method can be found in Slater et al. (1998), a basic summary of the elements applied in this study has been provided below. The detection of electrical anisotropy is achieved by rotating the transmitter (TX) and receiver (RX) coil at equal increments through 360° with respect to a point at the surface defined as the survey position. Plotting the apparent resistivity in polar plot form yields a qualitative map of directional changes in apparent electrical conductivity. These changes can be interpreted to represent dominant fracture sets, where the highest resistivity (i.e., long axis of the ellipse) will define the strike orientation (Taylor and Fleming, 1988). A single fracture set of preferred orientation will produce a symmetric ellipse-like pattern over 180° (Al-Garni and Everett, 2003); additional fracture sets with different preferred orientations can produce multiple superimposed ellipses. Azimuthal resistivity plots can be interpreted using spectral analysis to quantify the level of noise and assess potential impacts of local inhomogeneities on the signal (e.g., Bolshakov et al., 1997; Slater et al., 1998). First, the frequency content of the EM azimuthal resistivity dataset is assessed through the Discrete Fourier Transform:

Χ ðkÞ ¼

N −1 X

−ið2π=NÞkn

xðnÞ  e

ð1Þ

n¼0

where Χ(k) is the Discrete Fourier Transform of the finite sample series pffiffiffiffiffiffiffiffi x(n), N is the period of the signal, k and n are integers and i = −1. Here, the EM azimuthal resistivity dataset, x(n) is the apparent resistivity at azimuthal position n, which ranges from 0 to 360°/A, where A is the angular spacing in degrees between successive measurements. A more complete discussion of this analysis can be found in Slater et al. (1998). The spectral analysis approach requires coupled symmetric and asymmetric azimuthal array geometries. Here, a symmetric azimuthal resistivity array was performed with the point of rotation located at the midpoint between the TX and RX coils. Rotation through 360° resulted in the relocation of each position with the TX and RX interchanged at 180°. With the asymmetric array the point of rotation was not located at the midpoint; rather, the RX coil represented the point of rotation (i.e., the array was shifted by half coil spacing with respect to the symmetric array). In this case, rotation through 360° resulted in a different TX location for each position. Spectral analysis of symmetric and asymmetric data sets at a single measurement point provides quantitative information on the degree of symmetry (i.e., strength of the anisotropic response) contained in an azimuthal resistivity polar plot. In a homogeneously anisotropic environment, symmetric arrays rotated over 360° should yield symmetry over 180° due to the principle of reciprocity, which states that the measurement should be the same if the TX and RX coils are interchanged; therefore, any deviation from this would indicate the presence of measurement noise. While asymmetric arrays in the same environment will also exhibit symmetry over 180° when rotated through 360°, the presence of inhomogeneities can result in a non-symmetric response over 180°, depending on the magnitude of the inhomogeneity and the strength of the anisotropy. Given that measurement noise may also contribute to deviations over 360°, asymmetric arrays require rotation through 720° to assess corresponding measurement noise, and thus,

159

confirm an apparent heterogeneous response. A theoretical description of the spectral analysis technique supported by controlled laboratory experiments using symmetric and asymmetric rotations can be found in Slater et al. (1998). In this study, anisotropy and heterogeneity are evaluated using spectral analysis of symmetric arrays rotated over 360° and asymmetric arrays rotated over 720°. Here, we simply calculate the ratio (V) of the total energy contained in the power spectrum relative to the energy in the odd harmonics: N −1 X

V¼

k¼0 N=2−1 X 2

2

Χ ðkÞ :

ð2Þ

Χ ð2  k þ 1Þ

k¼0

Slater et al. (1998) showed that for perfect symmetry over 180° (i.e., a purely anisotropic environment) all odd harmonics of the power spectrum will be zero. Thus, this calculated parameter represents the signal-to-noise ratio for symmetric arrays, and the total energy relative to a combination of noise and heterogeneity for asymmetric arrays. A joint interpretation of these data sets enables an assessment of strength and representative spatial scale of the anisotropic response in the presence of measurement noise and local inhomogeneities. 3.2. EM azimuthal resistivity measurements An EM induction measurement collected over a layered earth represents a nonlinear average of the individual layer electrical conductivities within the effective sampling volume of the instrument. In general, EM induction measurements collected in horizontal dipole (vertical coplanar) and vertical dipole (horizontal coplanar) modes will be sensitive to depths of approximately 0.75 and 1.5 times the coil spacing, respectively. These values represent approximate estimates since the electrical conductivity distribution and applied frequency also affects the sampling volume. The earth response to different EM induction configurations can be estimated using an exact EM forward model that incorporates local sensitivities (McNeill, 1980); reliable layer electrical conductivity models of the subsurface hinges on accurate instrument calibrations (e.g., von Hebel et al., 2014). In this study, relative changes in EM induction response are interpreted to identify directional changes in resistivity. EM azimuthal resistivity surveys were collected at stations A and B within the 10 × 25 m plot (Fig. 3). Symmetric and asymmetric arrays were collected at 22.5° increments relative to north, and rotated through 360° and 720°. All measurements were collected using commercially available ground conductivity meters manufactured by Geonics Limited (Mississauga, Canada), and included the EM-38, EM31 and EM-34 devices, operated in horizontal and vertical dipole configurations. The EM-38 and -31 are manufactured with fixed coil spacings of 1 m and 3.66 m, respectively. Although the EM-34 is capable of 10, 20 and 40 m coil spacings, only the 10 m spacing was considered for this study due to spatial restrictions. Approximate sensing depths for these instruments in horizontal and vertical dipole modes in a homogeneous halfspace are as follows: EM-38 — 0.75 and 1.5 m; EM-31 — 2.75 and 5.5 m; EM-34 (10 m) — 7.5 and 15 m. Measurements collected with the EM-38 and -31 were averaged over 100 readings at each orientation using an automatic data logger, while readings were manually recorded with the EM-34. A manufacturer recommended calibration of the EM-38 instrument was performed prior to data collection, while no field calibrations were applied to the EM-31 and -34. Prior to data acquisition the devices were turned on for a period to allow instrumentation to stabilize to atmospheric temperature. Both the EM-38 and -31 symmetric and asymmetric arrays were rotated through 720° while the symmetric EM-34 arrays were only rotated through 360° (i.e., an asymmetric array was not performed

160

C.M. Steelman et al. / Journal of Applied Geophysics 116 (2015) 156–166 Table 1 GPR acquisition and signal processing parameters. Acquisition parameters

Signal processing

System: PulseEKKO 100 Nominal frequency: 100 MHz Step size: 0.1 m Antenna separation: 1 m Sampling interval: 0.8 ns Time window: 500 ns

Trace geometry editing Subtract-mean (dewow): 10 ns Bandpass frequency filter: 0–40–160–200 MHz Move start time: file header Time window cut: 350 ns Gain function: start = 50 ns; linear = 1 (pulse width−1); exponent = 1 db/m; max = 500 Migration: Kirchhoff 2D-velocity (mean 1D interval velocity model) Time-depth conversion: 0.086 m/ns Complex trace analysis: envelope Arithmetic function: square

Transmitter voltage: 1000 V Stacking: 64 traces Recoding: manual trigger Reflection line interval: 0.25 m CMP offset range = 1–10 m

using a dewow filter to remove dc shifts, bandpass frequency filtering to enhance signal, a gain function to compensate for attenuation, and migration analysis to collapse diffraction hyperbolas to their source. A complex trace analysis (i.e., trace envelope) was used to isolate strong reflection events over the 2-D section; these enveloped traces were subsequently squared to enhance high amplitude events relative to background. This arithmetic processing step increased the amplitude of collapsed diffraction hyperbolas and primary horizontal reflection events, improving visualization of high energy features. All postacquisition GPR data processing and visualization was performed using ReflexW software (Sandmeier Software, Germany), which included one-dimensional vertical velocity model generation using the semblance statistic within the CMP velocity analysis module. Finally, fullresolution 3-D models were generated using the parallel 2-D reflection profiles with a linear interpolation between lines. 4. Results and interpretation Fig. 3. Schematic of the 10 × 25 m study plot showing the position of two EM azimuthal resistivity stations (station A and station B). Azimuthal resistivity measurements were collected at 22.5° increments relative to north. GPR reflection profiles were collected across the study plot along parallel 2D lines from south to north spaced every 0.25 m (x-axis). GPR traces were collected every 0.1 m (y-axis) with the antennas orientated perpendicular to the line. CMP soundings collected at each azimuthal station at 0°, 45°, 90° and 135°, resulting in four velocity models at each station.

with EM-34 due to spatial restrictions). All instruments were placed directly on the ground surface during data collection. 3.3. Enhanced visualization of bedrock features using GPR GPR imaging provided optimal depth of investigation relative to the EM induction instruments. Although the dominant signal wavelength will limit the instrument's ability to resolve closely spaced fractures (Davis and Annan, 1989), major fracture features and dissolutionenhanced conduits should be apparent with 100 MHz antennas in the form of diffraction hyperbolas (e.g., Grasmueck et al., 2005, 2013). Reflection profiles and common-midpoint (CMP) soundings were collected using a Sensors & Software PulseEKKO 100 (Mississauga, Canada) equipped with 100 MHz antennas. Reflection profiles were collected along parallel lines from south to north originating from the south-west corner and spaced every 0.25 m along x-axis (west to east) (Fig. 3). GPR traces were collected along the y-axis (south to north) every 0.1 m using a manual trigger with the antennas orientated perpendicular to the line. Vertical velocity profiles were collected at stations A and B through a series of CMP soundings centered at each station; four soundings were collected at each station in a north–south, northeast–southwest, east–west, and southeast–northwest orientation. A detailed summary of the acquisition and signal processing parameters used in this study are provided in Table 1. GPR data was processed

4.1. Detection of electrical anisotropy using EM azimuthal resistivity EM azimuthal resistivity measurements at stations A and B are shown in Fig. 4 along with regional vertical fracture strike information of the underlying Amabel Formation collected along outcrops a few kilometers north of the study (Kunert et al., 1998). These polar plot data identify anisotropic trends over varying effective sampling volumes of the bedrock riverbed. The long-axis of the apparent resistivity ellipse marks the azimuth of the electrical anisotropy. It is evident that the anisotropic responses are variably dependent on coil spacing, magnetic dipole orientation and array type (symmetric vs. asymmetric). These data also exhibit considerable variability in the symmetry of the anisotropic ellipse. The high reproducibility of the individual measurements over 720° indicates low measurement noise. Regional fracture orientations from Kunert et al. (1998) indicate a dominant vertical and subvertical fracture trend of 30–40°. The azimuthal measurements collected at our study plot indicate greater spatial variability than that of the regional data. Further, EM azimuthal plots show multiple lobes (i.e., anisotropic responses) that vary with space and depth. Although there is evidence of a regionally trending first-order anisotropic response at 20–40° at both stations which roughly agrees with the regional fracture trend of 30–40°, additional secondorder trends were identified at 350–10° (north–south) and 80–100° (east–west). However, these secondary responses are variable over the study plot. A few third-order anisotropic responses were identified at marked peaks along the individual polar plots (e.g., station A: EM38 HD at 22.5° and EM-34 VD at 135°). The spatial variability in anisotropy becomes apparent through a comparison of EM-34 surveys at each station. Horizontal dipole modes at stations A and B show a consistent northeast–southwest trend, while the relatively deeper sampling vertical dipole modes identify an

C.M. Steelman et al. / Journal of Applied Geophysics 116 (2015) 156–166

161

Fig. 4. EM azimuthal resistivity measurements collected at station A and station B using an EM-38, EM-31, and an EM-34 with a 10 m coil spacing with vertical and horizontal magnetic dipole orientations. Measurements were collected at an equal rotational step of 22.5° through 360° (single rotation) to 720° (double rotation). Solid lines for EM-38 and EM-31 represent the average of two full rotations (720°). The EM-34 was only rotated through 360°. For symmetric arrays the TX and RX coils were set an equal distance from the point of rotation. For the asymmetric arrays the RX was positioned at the point of rotation, while the TX was rotated. The rose diagram represents regional joint orientations from the underlying Amabel Formation by Kunert et al. (1998) a few kilometers north of the study area.

east–west and north–south trend, respectively. It is apparent that the long-axis of the ellipses at station A are much broader than the corresponding ellipses at station B, particularly for the EM-38 and -31. Azimuthal polar plots were processed using spectral analysis to assess the strength of the anisotropic response in the presence of measurement noise and suspected spatial inhomogeneities. This approach

provides a semi-quantitative assessment of the relative strength of the anisotropy. Table 2 summarizes the spectral analysis results of the data presented in Fig. 4. Spectral ratios (V) were used to assess the nature of the apparent anisotropic response corresponding to the effective sampling volumes of the individual instruments. Vertical fracture anisotropy was defined as weakly, moderately or strongly homogeneous

Table 2 Spectral ratios (V) based on the EM azimuthal resistivity data for vertical dipole (VD) and horizontal dipole (HD) orientation. The symmetric quantities represent the signal-to-noise ratio of the apparent anisotropic response measured over 360°. The quantity for the asymmetric array collected over 360° represents the total energy relative to a combination of noise and spatial heterogeneities; the corresponding 720° quantity represents the signal-to-noise ratio. Using these values the anisotropic responses of the fractured bedrock are described in terms of whether they were weakly (w), moderately (m), or strongly (s) homogeneous or heterogeneous with respect to their effective sampling volume. Station A

EM-38 EM-31 EM-34

VD HD VD HD VD HD

Station B

Symmetric

Asymmetric

Anisotropic response

360°

360°

720°

21.6 20.2 3.9 24.1 6.0 28.2

3.2 11.1 3.8 2.2 – –

14.7 26.1 5.3 13.7 – –

s-heterogeneous s-homogeneous w-homogeneous s-heterogeneous m-homogeneous s-homogeneous

Symmetric

Asymmetric

360°

360°

720°

Anisotropic response

5.4 4.5 2.2 21.2 10.8 8.0

3.4 4.1 2.4 2.8 – –

9.4 4.7 11.6 61.3 – –

m-heterogeneous m-homogeneous w-homogeneous s-heterogeneous m-homogeneous m-homogeneous

162

C.M. Steelman et al. / Journal of Applied Geophysics 116 (2015) 156–166

or heterogeneous at each azimuthal station on a relative scale based on the following criteria: a relatively high signal-to-noise ratio for a symmetric array suggested a strongly anisotropic media, while a relatively low ratio for the corresponding asymmetric array over 360° at the same location indicated that the anisotropic response was likely being influenced by local inhomogeneities. Because an asymmetric array rotated over 360° incorporates the effects of measurement noise and inhomogeneities, a second rotation through 720° was used to assess the reproducibility of the measurements. If a relatively higher ratio was encountered with the 720° rotation, the initial heterogeneous response observed through 360° would likely be caused by local inhomogeneities. Spectral analysis results (Table 2) indicate a high-spatial dependency in depth and position over the study plot. While an anisotropic response can be observed in each data set, the nature of the response (i.e., homogeneous or heterogeneous), and its strength (i.e., weak, moderate, strong) with respect to noise and local inhomogeneities, varies considerably. For instance, station A exhibited a stronger anisotropic response in the near surface (e.g., EM-38 VD and HD; EM-31 HD; EM-34 HD) compared to station B. The strength of the anisotropic response was generally greater in the deeper bedrock sections for the EM-31 and EM-38 (i.e., vertical dipole resistivity was greater than the horizontal dipole) except for EM-38 responses collected at station B. Conversely, the horizontal dipole resistivity was greater than the vertical dipole for the EM-34 response. Overall, the influence of inhomogeneities on the anisotropic response appears to be moderate to strong, becoming progressively more dominant toward the rockbed interface, especially in the north end of the study plot.

porosity) or a geomechanical boundaries (i.e., change in fracture density or pattern). Depth-interval velocity models generated across four azimuths at each station are shown in Fig. 5; interval velocities exhibit some variability at each station. Station B shows greater velocity variability in the lower two layers, while station A shows relatively higher variability in the uppermost layer. Average interval velocity models obtained at each station and across the study plot are shown in Table 3. Fig. 6 provides a series of cube model representations of the study plot area from four perspectives. A quadrant of each block was removed to show the GPR responses near the azimuthal stations. These models highlight varying degrees of reflectivity with depth and multiple reflection features, such as strong hyperbolic events, signal scattering patterns, and highly reflective horizontal interfaces. These models indicate that interface 1:2 (i.e., interface between layers 1 and 2) is horizontal and reasonably continuous across the study plot. Interface 2:3 shows strong continuity and exhibits a slight 4° dip to the west, while interface 3:4 is much weaker and is highly discontinuous with an apparent 3° dip of to the east. The majority of the GPR energy over the study plot arises from layer 2 in the form of discontinuous horizontal reflections and scattered energy; the strength of the response increases toward the west. Large diffraction hyperbolas are evident from the base of layer 1 through the bottom layer 4, with highest accumulation along interface 1:2. Reflectivity is highest in the northeast quadrant (i.e., east of station B). Fig. 7 presents the 2-D reflection profile extracted from the model at position x = 5 m from south to north through stations A and B. The processed unmigrated section (Fig. 7a) shows clear diffractions along the bottom of layer 1 and top of layer 2. These events were collapsed using a 2D Kirchhoff migration (Fig. 7b) with a 1-D mean interval velocity. The three reflection boundaries identified in the CMP soundings,

4.2. Interpretation of electrical anisotropic response Azimuthal resistivity measurements indicate moderately to highly variable electrical anisotropy with strong indications of heterogeneous features affecting the apparent anisotropic response over a small 250 m2 area. These anisotropic responses are interpreted to represent the orientation of preferred fracture sets in proximity to the EM coils (i.e., effective sample volumes of the individual instruments and their coil configurations). Symmetric arrays exhibited reasonable symmetry over 180° consistent with a vertical or subvertical fracture set; however, moderate to high deviations were encountered with corresponding asymmetric arrays which suggest that local inhomogeneities may be influencing or causing the apparent anisotropic response. Therefore, vertical fracture networks and local inhomogeneities are believed to be variably distributed across the study plot area. Azimuthal data also suggests that there are relative changes in the strength of the anisotropy across the study plot. Given that these EM conductivity instruments respond to the vertical magnetic field resulting from induced horizontal eddy currents in the ground, a reduction in the electrical current density due to an increase in conductive area along the strike of the anisotropy would effectively reduce the magnitude of the measured anisotropic response (Al-Garni and Everett, 2003). Based on this idea, station B could be positioned over a more developed fracture network or dissolution-enhanced zone relative to station A. Field observations at station B indicate visible bedrock dissolution at surface. The relative differences in apparent resistivity between corresponding azimuthal measurements at stations A and B indicate that the intensity or connectivity of the vertical fracture network varies over small distances. These resistivity trends support the presence of a spatially complex vertical fracture network with varying degrees of connectivity (i.e., interconnection of fractures, intensity, aperture, and dissolution) between the river and aquifer. 4.3. Visualization of fracture and dissolution-enhanced features using GPR CMP soundings identified three strong reflection events interpreted to represent major lithostratigraphic boundaries (i.e., changes in matrix

Fig. 5. CMP interval velocity models at station A and station B within the upper 10 m. Semblance analysis was used to construct an interval velocity depth model at four orientations (0°, 45°, 90°, 135°) for each station. The average interval velocity models are summarized in Table 3.

C.M. Steelman et al. / Journal of Applied Geophysics 116 (2015) 156–166 Table 3 Average CMP interval-velocity model results shown in Fig. 4. Interval

Layer 1 Layer 2 Layer 3

Station A

Station B

Mean

m

m/ns

m

m/ns

m

m/ns

1.95 4.05 3.20

0.082 0.085 0.091

2.13 3.56 3.40

0.080 0.081 0.095

2.04 3.80 3.30

0.081 0.083 0.093

and also denoted in the block model, can be identified in the migrated section. The squared envelope of these traces (Fig. 7c) indicate that reflection 1 is largely demarcated by a series of high-energy zones corresponding to the source of diffraction hyperbolas, while reflection 2 is much more uniform across the plot. The lowermost reflection event exhibits markedly lower energy relative to the overlying events, which weakens toward the south. Isolated high-energy pockets associated with collapsed diffraction events and short discontinuous reflections can be observed throughout the profile across the upper 15 m. To evaluate the spatial distribution of these high-energy zones across the study area, reflection energy was computed over the data cube by summing the squared enveloped amplitudes integrated over 1 m depth slices from 0 to 15 m. The resulting relative amplitude energy plots (Fig. 8) identify numerous variably distributed linear features, as well as broader zones of higher reflectivity. The density, orientation and distribution of linear features vary laterally and with depth, which supports a variably connected vertical fracture network. Clear zones of higher reflectivity are evident to the south from 2 to 3 m depth and the north from 11 to 15 m depth. Fig. 9 presents rose plot models

163

based on azimuthal picks of the major linear features identified over the study plot (i.e., dotted lines superimposed over interpreted linear features in Fig. 8). Plots were generated for depth intervals of 0–2 m, 2–6 m, 6–9 m, and 9–15 m, roughly corresponding to the layers defined by the major GPR reflecting interfaces. Layer 1 exhibits a reasonably strong linear set orientated at 20–30°. Layer 2 exhibits two major sets orientated at 50–60° and 120–140°, while layer 3 is marked by sets at 30–50° and 160–170°. Finally, the underlying bedrock is marked by a single set at 0–10°. 4.4. Geophysical conceptualization of fracture networks and dissolution enhanced features Full-resolution GPR improved visualization of suspected features influencing the electrical anisotropic response, including: hyperbolic events likely associated with vertical fractures or local discontinuities; signal scattering patterns resulting from vuggy or dissolution enhanced bedrock; and highly reflective horizontal interfaces associated with larger water filled fractures and conduits. The spatial distribution of reflected energy across the study plot supports the existence of linear features of preferred orientation. Although these linear high-energy events may not explicitly represent fracture sets, they do support the existence of features with high dielectric contrast and preferred orientation. Given their high reflectivity, these zones likely represent water filled (i.e., electrically conductive and low radar velocity) features in the bedrock matrix. Observed changes in layer velocity may be indicative of varying degrees of heterogeneity in the form of total porosity. For instance, measured interval velocities at station A were slightly

Fig. 6. Block model representations of GPR reflectivity and major reflection features across the study plot. The dot represents the southwest corner and origin of the reflection profile grid with north directed along the y-axis. Layers 1–3 correspond to the CMP interval velocity models, while layer 4 represents the remaining depth section. (i) Strong hyperbolic reflections are observed throughout; (ii) zones exhibiting signal scattering associated with dissolution or vuggy sections of rock; and (iii) highly reflective interfaces with strong diffraction hyperbolas. A time-to-depth conversion was applied using a velocity of 0.086 m/ns.

164

C.M. Steelman et al. / Journal of Applied Geophysics 116 (2015) 156–166

Fig. 7. (a) Processed unmigrated 100 MHz reflection profile at x = 5 m. Profile exhibits readily identifiable diffraction hyperbolas (pink arrows) along the top of layer 2. Intermittent reflection events are evident to depths of 15 m. (b) Migrated reflection profile using 2D Kirchhoff migration with mean 1D interval velocity model from station A and station B. The three major reflecting boundaries (black arrows) used in the CMP interval velocity models are identified. (c) The squared envelope amplitude of the migrated section was calculated to accentuate position and depth of high-energy events, including collapsed diffraction hyperbolas shown in (a). A time-to-depth conversion was applied to all sections using a velocity of 0.086 m/ns.

higher which indicates lower porosity in the south end of the study plot. This is consistent with the weaker degree of anisotropy observed at station B which is thought to be associated with increased fracture density and/or concentrated dissolution-enhanced conduits, resulting in relatively higher porosity and hydraulic connectivity.

5. Discussion The application of EM azimuthal resistivity supported by highresolution GPR measurements improved our conceptualization of vertical and sub-vertical fracture networks across a complex, shallow bedrock flow system. A thorough evaluation of the azimuthal responses at this site was achieved using a range of coil spacings and magnetic dipole orientations at two locations, providing high-resolution information

over a range of representative volumes considerably smaller than those of previous studies (e.g., Lane et al., 1995; Busby, 2000; Skinner and Heinson, 2004; Obiadi et al., 2013). These azimuthal data were augmented by GPR reflectivity measurements over the study plot. Together, these data indicated a high degree of spatial variability in vertical fractures and dissolution-enhanced features to depths of 15 m over a very small 250 m2 area. While previous work focused on strategies for removing heterogeneous effects on the anisotropic signal (e.g., Slater et al., 1998), we found that the identification of inhomogeneities within a variably anisotropic environment provided valuable local-scale information considering the discrete nature of hydrologic flow paths along fractured rockbed interfaces. Although the presence of inhomogeneities within such a small sampling volume led to more challenging interpretations of azimuthal responses, these issues were partially rectified through an integrated analysis of EM azimuthal resistivity measurements using multiple coil spacings, magnetic dipole orientations at multiple locations, coupled with high-resolution GPR measurements. This approach enabled a more complete conceptualization of spatially complex fracture and conduit flow features over a representative volume of rock. It is important to recognize the value of decoupling spatial inhomogeneities from the local anisotropic signal over relatively small sample volumes. In practice, large sampling volumes are used to reduce potentially negative effects of inhomogeneities on the anisotropic signal (e.g., Busby, 2000). Although our data indicates that vertical fracture networks and conduits varied considerably over very small spatial scales, alternative larger-scale measurements would provide limited insight into important local-scale processes given their diluted representation of the local-scale environment. For instance, the interpretation of hydraulic response data collected in a fractured environment relies upon an accurate conceptual understanding of fracture networks (Ritzi and Andolsek, 1993); our local-scale geophysical methodology would enable non-invasive characterization of a volume of rock consistent with that needed to understand radial flow around shallow wells. Specifically, spatial changes in anisotropy could be used to qualitatively assess variations in the intensity and effective apertures of hydraulically connected fracture zones. This application was demonstrated in this study through contrasting anisotropic responses and GPR reflectivity at station A and station B. Conducting EM azimuthal measurements using multiple coil spacing and configurations at the same location should lead to more effective conceptualizations of homogeneous and heterogeneous anisotropy, thus providing spatially representative information on lateral and vertical connectivity of fractures networks underlying bedrock river environments. The design and construction of multi-borehole groundwater monitoring stations along bedrock river systems could be optimized through prior non-invasive geophysical characterizations of the local fracture-conduit flow system, such that angled boreholes could be emplaced and orientated to intercept major fracture networks and local dissolution features. Conversely, existing angled and vertical borehole data can be placed into a 3-D fracture network context using these complementary geophysical data sets. It is recognized that many challenges remain in the effective characterization of depth-discrete azimuthal resistivity measurements obtained from an EM induction data set. While a recent study by von Hebel et al. (2014) has shown that EM measurements employing multiple coil configurations can be jointly inverted using advanced algorithms to resolve individual subsurface layers, the strong horizontal anisotropy of fractured bedrock would likely require more sophisticated 2 and 3 dimensional analyses utilizing the azimuthal response, to fully constrain a multilayered fractured bedrock model. 6. Conclusions Our integrated analysis of azimuthal resistivity and GPR reflectivity over the study area indicated spatial variations in the density,

C.M. Steelman et al. / Journal of Applied Geophysics 116 (2015) 156–166

165

Fig. 8. Cumulative amplitude plots of squared envelope over discrete depth slices. Depth-slices correspond to the full 10 × 25 m plot area with position of azimuthal stations A and B (north is pointing up). These relative energy plots identify numerous linear features interpreted to represent major vertical to subvertical fractures or dissolution-enhanced features (water filled fractures or conduits). Broader areas of high amplitude may indicate strong dissolution and along major lithologic bedding planes. Dotted lines mark interpreted azimuth of a high-energy reflecting feature.

orientation and distribution of linear features in the upper 15 m of bedrock over a 250 m2. Specific comparisons between EM azimuthal responses and the distribution of GPR enveloped amplitude with depth indicated that fracture anisotropy will likely exhibit a strong spatial dependency along the river given the irregular length, distribution and magnitude of local-scale trends. The varying heterogeneous anisotropic responses observed through spectral analysis of EM azimuthal resistivity measurements provided valuable insights about the connectivity and continuity of vertical fracture networks within a naturally heterogeneous bedrock system. Therefore, the utility of EM azimuthal resistivity and GPR for the characterization of discrete vertical fracture networks at scales representative of hydrogeologic processes along a fractured bedrock river, i.e., a scale consistent with hypothesized groundwater and surface water exchange patterns following discrete fracture networks across the rockbed surface, was found to be reasonably strong. These results confirmed the existence of a complex and variably distributed vertical fracture network with dissolution-enhanced features that may readily connect surface water to groundwater. Given the spatial limitations and the inherent ecological sensitivity of river environments, it was important to develop a non-invasive

methodology capable of providing relevant information about the local flow system without unnecessary impacts to surrounding environments. Fracture anisotropy and the presence of inhomogeneities within the shallow bedrock flow system were successfully characterized without the need for large survey geometries and invasive methods. This integrated approach could be used to assess lateral variations in groundwater–surface water connectivity along rivers with exposed bedrock or relatively thin sediment cover. Our geophysical methodology might prove useful in the design and emplacement of hydrogeological and hydrogeophysical monitoring networks aimed at advancing our understanding of complex groundwater–surface water flow patterns in shallow bedrock flow systems. Acknowledgments This research was made possible through funding by the Ontario Research Fund (ORF-RE 03-061: Sustainable Bedrock Water Supplies for Ontario Communities) secured by Dr. Parker and an NSERC PDF and Banting Postdoctoral Fellowship to Dr. Steelman. The authors would like to thank Charles and Anna Simon, Richard Lay and the Eden Mills

Fig. 9. Rose plot representation of interpreted high-energy reflecting features over the study plot between (a) layer 1 (0–2 m), (b) layer 2 (2–6 m), (c) layer 3 (6–9 m), and (d) material beneath interface 3 (9–15 m) based on apparent trends in Fig. 8.

166

C.M. Steelman et al. / Journal of Applied Geophysics 116 (2015) 156–166

Mill Pond Conservation Association for fostering this research, as well as the many students and staff at G360 who assisted with data acquisition in the field. References Abdullahi, N.K., Osazuwa, I.B., 2011. Geophysical imaging of municipal solid waste contaminant pathways. Environ. Earth Sci. 62, 1173–1181. Al-fares, W., Bakalowicz, M., Guérin, R., Dukhan, M., 2002. Analysis of the karst aquifer structure of the Lamalou area (Hérault, France) with ground penetrating radar. J. Appl. Geophys. 51, 97–106. Al-Garni, M., Everett, M.E., 2003. The paradox of anisotropy in electromagnetic loop–loop responses over a uniaxial half-space. Geophysics 68 (3), 892–899. Batayneh, A.T., Abueladas, A.A., Moumani, K.A., 2002. Use of ground-penetrating radar for assessment of potential sinkhole conditions: an example from Ghor al Haditha area, Jordan. Environ. Geol. 41, 977–983. http://dx.doi.org/10.1007/s00254-001-0477-8. Birkhead, A.L., Heritage, G.L., White, H., van Niekerk, A.W., 1996. Ground-penetrating radar as a tool for mapping the phreatic surface, bedrock profile, and alluvial stratigraphy in the Sabie River, Kruger National Park. J. Soil Water Conserv. 51, 234–241. Bolshakov, D.K., Modin, I.N., Pervago, E.V., Shevnin, V.A., 1997. Separation of anisotropy and inhomogeneity influence by the spectral analysis of azimuthal resistivity diagrams. Proceedings of the 3rd EEGS — European Section Meeting, Aarhus, Denmark, pp. 147–150. Brunton, F.R., 2009. Update of revisions to the Early Silurian stratigraphy of the Niagara Escarpment: integration of sequence stratigraphy, sedimentology and hydrogeology to delineate hydrogeologic units. Summary of Field Work and Other Activities 2009, Ontario Geologic Survey, Open File Report 6240, pp. 25-1–25–20. Busby, J.P., 2000. The effectiveness of azimuthal apparent-resistivity measurements as a method for determining fracture strike orientations. Geophys. Prospect. 48, 677–695. Busby, J., Jackson, P., 2006. The application of time-lapse azimuthal apparent resistivity measurements for the prediction of coastal cliff failure. J. Appl. Geophys. 59, 261–272. Conant Jr., B., Cherry, J.A., Gillham, R.W., 2004. A PCE groundwater plume discharging to a river: influence of the streambed and near-river zone on contaminant distributions. J. Contam. Hydrol. 73, 238–279. Crook, N., Binley, A., Night, R., Robinson, D.A., Zarnetske, J., Haggerty, R., 2008. Electrical resistivity imaging of the architecture of substream sediments. Water Resour. Res. 44, W00D13. http://dx.doi.org/10.1027/2008WR006968. Davis, J.L., Annan, A.P., 1989. Ground-penetrating radar for high-resolution mapping of soil and rock stratigraphy. Geophys. Prospect. 37, 531–551. Gabrielli, C.P., McDonnell, J.J., 2011. An inexpensive and portable drill rig for bedrock groundwater studies in headwater catchments. Hydrol. Process. 26, 622–632. Grasmueck, M., Weger, R., Horstmeyer, H., 2005. Full-resolution 3D GPR imaging. Geophysics 70, K12–K19. http://dx.doi.org/10.1190/1.1852780. Grasmueck, M., Quinta, Pomar, K., Eberli, G.P., 2013. Diffraction imaging of sub-vertical fractures and karst with full-resolution 3D ground-penetrating radar. Geophys. Prospect. 61, 907–918. http://dx.doi.org/10.1111/1365-2478.12004. Habberjam, G.M., 1975. Apparent resistivity, anisotropy and strike measurements. Geophys. Prospect. 23, 211–247. Hancock, P.J., 2002. Human impacts on the stream — groundwater exchange zone. Environ. Manag. 29, 763–781. http://dx.doi.org/10.1007/s00267-001-0064-5. Hatch, C.E., Fisher, A.T., Revenaugh, J.S., Constantz, J., Ruehl, C., 2006. Quantifying surface water — groundwater interactions using time series analysis of streambed thermal records: method development. Water Resour. Res. 42, W10410. http://dx.doi.org/ 10.1029/2005WR004787. Hodge, R.A., Hoey, T.B., Sklar, L.S., 2011. Bed load transport in bedrock rivers: the role of sediment cover in grain entrainment, translation and deposition. J. Geophys. Res. Earth Surf. 116, F04028. http://dx.doi.org/10.1029/2011F002032. Kruse, S., Grasmueck, M., Weiss, M., Viggiano, D., 2006. Sinkhole structure imaging in covered Karst terrain. Geophys. Res. Lett. 33, L16405. http://dx.doi.org/10.1029/ 2006GL026975. Kunert, M., Coniglio, M., Jowett, E.C., 1998. Controls and age of cavernous porosity in Middle Silurian dolomite, southern Ontario. Can. J. Earth Sci. 35, 1044–1053. Lane Jr., J.W., Haeni, F.P., Watson, W.M., 1995. Use of a square-array direct-current resistivity method to detect fractures in crystalline bedrock in New Hampshire. Ground Water 33, 476–485. Matias, M.J.S., 2002. Square array anisotropy measurements and resistivity sounding interpretation. J. Appl. Geophys. 49, 185–194. McNeill, J.D., 1980. Electromagnetic terrain conductivity measurement at low indiction numbers. Technical Note, TN-6. Geonics Ltd. Meshkova, L.V., Carling, P.A., 2013. Discrimination of alluvial and mixed bedrock-alluvial multichannel river networks. Earth Surf. Process. Landf. 38, 1299–1316. Meyer, J.L., 1997. Stream health: incorporating the human dimension to advance stream ecology. J. N. Am. Benthol. Soc. 16, 439–447.

Musgrave, H., Binley, A., 2011. Revealing the temporal dynamics of subsurface temperature in a wetland using time-lapse geophysics. J. Hydrol. 396, 258–266. Obiadi, I.I., Onwuemesi, A.G., Anike, O.L., Ajaegwu, N.E., Anakwuba, E.K., Nwosu, C.M., Akpunonu, E.O., Onuigbo, E.N., Onuba, O.L., 2013. Determining subsurface fracture characteristics from the azimuthal square array resistivity survey at Igarra, Nigeria. Pure Appl. Geophys. 170, 907–916. Oxtobee, J.P.A., Novakowski, K., 2002. A field investigation of groundwater/surface water interactions in a fractured bedrock environment. J. Hydrol. 269, 169–193. Ritzi Jr., R.W., Andolsek, R.H., 1992. Relation between anisotropic transmissivity and azimuthal resistivity surveys in shallow, fractured, carbonate flow systems. Ground Water 30, 774–780. Sandberg, S.K., Slater, L.D., Versteeg, R., 2002. An integrated geophysical investigation of the hydrogeology of an anisotropic unconfined aquifer. J. Hydrol. 267, 227–243. Schmidt, C., Conant Jr., B., Bayer-Raich, M., Schirmer, M., 2007. Evaluation and field-scale application of an analytical method to quantify groundwater discharge using mapped streambed temperatures. J. Hydrol. 347, 292–307. Skinner, D., Heinson, G., 2004. A comparison of electrical and electromagnetic methods for the detection of hydraulic pathways in a fractured rock aquifer, Clare Valley, South Australia. Hydrogeol. J. 12, 576–590. Slater, L.D., Sandberg, S.K., Jankowski, M., 1998. Survey design procedures and data processing techniques applied to the EM azimuthal resistivity method. J. Environ. Eng. Geophys. 3, 167–177. Steelman, C.M., Kennedy, C., Parker, B.L., 2015. Geophysical conceptualization of a fractured sedimentary bedrock riverbed using ground-penetrating radar and induced electrical conductivity. J. Hydrol. 521, 433–446. Taylor, R.W., Fleming, A.H., 1988. Characterizing jointed systems by azimuthal resistivity surveys. Ground Water 26, 464–474. Tinker, K.J., Wohl, E.E. (Eds.), 1998a. Rivers over rock: fluvial processes in bedrock channels. Geophysical Monograph 17. American Geophysical Union, Washington, D.C. Tinker, K.J., Wohl, E.E., 1998b. A primer on bedrock channels, in Rivers over rock: fluvial processes in bedrock channels. In: Tinker, K.J., Wohl, E.E. (Eds.), Geophysical Monograph 17. American Geophysical Union, Washington, D.C., pp. 1–18. Vannote, R.L., Minshall, G.W., Cummins, K.W., Sedell, J.R., Cushing, C.E., 1980. The river continuum concept. Can. J. Fish. Aquat. Sci. 37, 130–137. Von Hebel, C., Rudolph, S., Mester, A., Huisman, J.A., Kumbhar, P., Vereecken, H., van der Kruk, J., 2014. Three-dimensional imaging of subsurface structural patterns using quantitative large-scale multiconfiguration electromagnetic induction data. Water Resour. Res. 50. http://dx.doi.org/10.1002/2013WR014864. Ward, A.S., Gooseff, M.N., Singha, K., 2010. Characterizing hyporetic transport processes — interpretation of electrical geophysical data in coupled stream–hyporheic zone systems during solute tracer studies. Adv. Water Resour. 33, 1320–1330. Ward, A.S., Fitzgerald, M., Gooseff, M.N., Voltz, T.J., Binley, A.M., Singha, K., 2012. Hydrologic and geomorphic controls on hyporheic exchange during base flow recession in a headwater mountain stream. Water Resour. Res. 48, W04513. http://dx.doi.org/ 10.1029/2011WR0114161. Wasscher, J.D., 1961. On four-point resistivity measurements on anisotropic conductors. Philips Res. Rep. 16, 301–306. Watson, K.A., Barker, R.D., 1999. Differentiating anisotropy and lateral effects using azimuthal resistivity offset Wenner soundings. Geophysics 64, 739–745. http://dx.doi. org/10.1190/1.1444583. Watson, K.A., Barker, R.D., 2005. Modelling azimuthal resistivity sounding over a laterally changing resistivity subsurface. Near Surf. Geophys. 3, 3–11. Watson, K.A., Barker, R.D., 2010. Tank modelling of azimuthal resistivity syrveys over anisotropic bedrock with dipping overburden. Near Surf. Geophys. 8, 297–309. http:// dx.doi.org/10.3997/1873-0604, 2010019. Winter, T.C., Harvey, J.W., Franke, O.L., Alley, W.M., 1998. Ground water and surface water: a single resource. USGS Circular 1139, Denver, Colorado. Wishart, D.N., Slater, L.D., Gates, A.E., 2008. Fracture anisotropy characterization in crystaline bedrock using field-scale azimutahl self potential gradient. J. Hydrol. 358, 35–45. http://dx.doi.org/10.1016/j.jhydrol.2008.05.017. Woessner, W.W., 1998. Changing views of stream-groundwater interaction. In: Van Brahana, J., Eckstein, U., Ongley, L.K., Schneider, R., Moore, J.E. (Eds.), Proceedings of the Joint Meeting of the XXVIII Congress of the International Association of Hydrogeologists and the Annual Meeting of the American Institute of Hydrology: Gambling with Groundwater, Las Vegas, Nevada, Sept.28 – Oct.2, 1998. AIH, St.Paul, MN, pp. 1–6. Woessner, W.W., 2000. Stream and fluvial plain ground water interactions: rescaling hydrogeologic thought. Ground Water 38, 423–429. Wojnar, A.J., Mutiti, S., Levy, J., 2013. Assessment of geophysical surveys as a tool to estimate riverbed hydraulic conductivity. J. Hydrol. 482, 40–56. Yeboah-Forson, A., Whitman, D., 2014. Electrical resistivity characterization of anisotropy in the Biscayne aquifer. Ground Water 52, 728–736.