Constraining the lithostratigraphic architecture of a buried bedrock valley using surface electrical resistivity and seismic refraction tomography

Buried bedrock valleys play an important role in regional groundwater flow systems (e.g., Russell et al. 2004; Cummings et al. 2012; Sharpe et al. 2013, 2018; Bajc et al. 2018). These features represent an erosional structure commonly found in northern mid-latitude environments including across Canada (e.g., Butler et al. 2004; Al et al. 2005; Davies et al. 2008; Cummings et al. 2012; Bajc et al. 2018; Broster et al. 2013; Steelman et al. 2018), the northern United States (e.g., Kehew and Boettger 1986; Shaver and Pusc 1992; Seyoum and Eckstein 2014; Korus et al. 2017), and northern Europe (e.g., Sandersen and Jørgensen 2003; Steuer et al. 2009). These regions have been subjected to glaciofluvial and paleoalluvial processes resulting in localized bedrock incision and subsequent infilling with glacial sediments.

Valley formation can take many millennia and be polygenetic as a result of both physical and chemical weathering processes (e.g., Kehew and Boettger 1986; Jørgensen and Sandersen 2006; Gao 2011; Priebe et al. 2019). Valley incision is influenced by regional tectonics, structural features such as faults and fractures, glacial history including mechanical influences on the bedrock fracture conditions from loading and unloading, and pre- and postglacial paleoenvironments (Cummings et al. 2012). While buried bedrock valleys are associated with glaciated terrains, there is debate on the relative importance of preglacial, subglacial, glaciofluvial, and proglacial erosion in their formation (e.g., Shaver and Pusc 1992; Jørgensen and Sandersen 2006; Meyer and Eyles 2007; Gao 2011; Steelman et al. 2018). Detailed characterization of the geology in and around buried bedrock valleys is required to further unravel how these features form and ultimately the spatial variability of hydrogeologic characteristics in and around these features.

Buried bedrock valleys preserve sediments that can be used to reconstruct regional glacial histories and understand preglacial drainage patterns (Greenhouse and Karrow 1994). Such valleys can influence local and regional groundwater flow systems and play a role in recharge of bedrock aquifers (Cole et al. 2009; Cummings et al. 2012). They can also act as preferential pathways, enhancing a deeper aquifer’s susceptibility to contamination (e.g., Shaver and Pusc 1992; Sharpe et al. 2013; Seyoum and Eckstein 2014; Burt 2018). Many municipal, private, industrial, and agricultural wells are completed in aquifers within buried bedrock valley sediments (Russell et al. 2004; Burt 2018), as well as targeting the fractured bedrock aquifers nearby to buried bedrock valleys where the bedrock may possess enhanced transmissivity (e.g., Cole et al. 2009). Russell et al. (2018) identified buried bedrock valleys as under-represented targets in groundwater resource studies, noting poor data support in regional studies and the need for more cost-effective methods for locating and understanding these features. Strategies that lead to more efficient and detailed characterizations of valley morphology and infill sediment type will improve our understanding of regional groundwater flow and municipal water supply strategies (e.g., Cummings et al. 2012; Bajc et al. 2018; Burt 2018).

A variety of geophysical methods have been used to investigate the lithostratigraphy in and around buried bedrock valleys (e.g., Greenhouse and Karrow 1994; Pugin et al. 1996; Butler et al. 2004; Bajc et al. 2012; Oldenborger et al. 2013; Stumpf and Ismail 2013; Oldenborger et al. 2016; Korus et al. 2017; Sharpe et al. 2018; Steelman et al. 2018; Morgan et al. 2019). Compared to conventional coring, surface geophysical surveys such as provided by electrical, seismic, and gravity-based instruments, are relatively inexpensive and offer the potential for more efficient data collection over a larger volume or area with sufficient spatial resolution suitable for delineating the presence and spatial continuity in two- or three-dimensions of local-scale hydrogeologic features (e.g., Kanta et al. 2013; Sirieix et al. 2014; Ciruzzi and Lowry 2017; Steelman et al. 2018). Conventional geophysical techniques do not directly measure hydrogeological properties, which means these data are most effective when constrained or interpreted (e.g., calibrated or confirmed) alongside core and (or) hydraulic measurements (Boaga 2017).

In this study, we demonstrate the value of combining co-located, high-resolution, ground-based electrical resistivity and seismic refraction tomography surveys with existing continuous core and an electrofacies model based on downhole geophysical logs along two transects centred over a buried bedrock valley in Elora, Ontario, Canada. This study follows a multiple lines of evidence approach using geophysical methods in combination with existing data to assess the valley architecture and refine the conceptual site model for erosional and depositional characteristics of this feature. The insights are relevant to understanding the sediment variability influencing 3D flow system conditions in and around the Elora municipal wellfield, but the fundamental insights and experience with data sets will also support similar studies elsewhere. This work provides for a more robust understanding of buried bedrock valley evolution and ultimately the hydrogeological significance and characteristics of these features.

The bedrock geology of the upper ∼150 m of the study area near Elora, Ontario comprised Middle Silurian dolostone, and shale sequences (Brunton and Brintnell 2020). These rocks were deposited in a shallow water carbonate ramp that draped over the Algonquin Arch, a tectonic high dividing the subsiding Michigan and Appalachian basins (Brunton and Brintnell 2020). This sedimentary system comprised a diverse range of sediments from paleoenvironments including quiet backwater reef lagoons and high-energy reefal environments recording multiple sea level regressions and transgressions over time, resulting in varying unit thicknesses (Brunton and Brintnell 2020). Differential erosion of younger units of the dipping strata formed a Paleozoic bedrock cuesta known as the Niagara Escarpment, which forms the eastern boundary of an 80–140 m thick groundwater flow system of fractured carbonate sedimentary rock successions. The Silurian dolostone serves as one of the best freshwater bedrock aquifers in Ontario, with numerous private wells and major cities including Hamilton, Cambridge, Guelph, Elora, and Fergus reliant on these aquifers for water supply (Singer et al. 2003). Key bedrock aquifers include the Guelph, Goat Island, and Gasport formations with regional flow generally occurring from north to south (Brunton and Brintnell 2020). The Eramosa and Cabot Head formations are generally considered regional aquitards, but the former has been demonstrated to be discontinuous and is absent from the study area (e.g., Golder Associates 2013; Matrix Solutions Inc. 2017). While the bedrock units can be highly fractured, flow typically occurs along dissolution-enhanced sequence stratigraphic breaks, bedding plane partings, or lithologic contrasts (Priebe et al. 2019; Munn et al. 2020b). The regional Quaternary geology consists of glacial and nonglacial sediments deposited since the late Pleistocene through multiple cycles of glaciations across Ontario (Karrow 1967; Chapman and Putnam 1984; Burt 2018).

An extensive network of buried bedrock valleys exists throughout southern Ontario with the Quaternary sediments infilling these valleys frequently distinct from and potentially older than the unconsolidated glacial deposits outside the valleys, making aquifer–aquitard sequences within the valley fills difficult to predict (e.g., Bajc et al. 2012, 2018; Burt 2018). Many of these valleys follow the regional jointing and fracturing pattern and parallel modern-day river systems (Eyles et al. 1997).

Valley incision is known to influence the physical properties of the surrounding bedrock (Matheson and Thomson 1973; Ferguson and Hamel 1981; Hamel 2018). Erosion of bedrock along incised channels reduces confining pressures and results in stress-relief fracturing, which may increase bedrock transmissivity (Gross and Engelder 1991; Hamel 2018). Cole et al. (2009) observed a slightly positive correlation (r² = 0.4675) between the rated capacity of municipal wells in Guelph, Ontario, and their proximity to the Guelph buried bedrock valley, with wells proximal to the buried bedrock valley producing more water than those located farther away. The hydraulic influence of the Guelph buried bedrock valley was described as extending up to 4 km (Cole et al. 2009). Priebe et al. (2019) analysed pumping test results collected in a 600 km² area around the City of Guelph in southern Ontario, but their results did not show a strong statistical correlation (i.e., r² = 0.3787–0.4281) between hydraulic conductivity within the Gasport formation and valley proximity. Nevertheless, they noted the highest hydraulic conductivity values within their study area were associated with karstic sequences that formed where vertical conduits such as bedrock valleys penetrate most deeply, facilitating recharge flux into, and dissolution enhancement within, the Gasport formation. Cole et al. (2009) and Priebe et al. (2019) hypothesized that bedrock transmissivity may be associated with the extent of valley incision into the Gasport formation, implying that the Gasport formation would possess more discrete flow features including dissolution enhanced channels and karst proximal to the valley if incision were able to advance through the Guelph, Eramosa, and Goat Island formations. Therefore, evolution of a valley is an important aspect to consider in assessments of groundwater flow system complexity with consideration of hydraulic conductivity distributions between fractured bedrock and unconsolidated sediments.

The study area is located approximately 2 km southeast of Elora, Ontario, and the present-day Grand River valley within the Township of Centre Wellington (Fig. 1a), approximately 20 km equidistant between the Orangeville Moraine to the northeast, Waterloo Moraine to the southwest, and Paris Moraine to the east. During past glaciation, the study area was likely part of the Georgian Bay–Ontario interlobate zone. In this zone, both the Lake Ontario ice lobe and the Georgian Bay ice lobe would have covered the area, with these advancing from the northwest and southeast, respectively (Karrow 1967; Greenhouse and Karrow 1994). The local area is draped by sandy and silty to clayey tills, glaciofluvial deposits, and ice-contact stratified deposits including sand and gravel (Burt 2018). These glacial drift deposits can exhibit extensive lateral variation over tens of metres and range in thickness from nonexistent (i.e., areas of outcropping bedrock) near the Elora Gorge to approximately 150 m thick over buried bedrock valleys or across the Orangeville Moraine (Burt and Dodge 2016; Steelman et al. 2018).

The rock depression now known as the Elora buried bedrock valley was first officially identified in the late 1960s from water well records (Greenhouse and Karrow 1994). The Elora buried bedrock valley was subsequently mapped using gravity measurements between Belwood Lake and Inverhaugh (Greenhouse and Karrow 1994). The only direct measurements of the Quaternary sediments available in the immediate study area come from two locations along the axis of the Elora buried bedrock valley (UW1 and UW2/UWBH34-78 in Fig. 1a). The southwest location consists of a geophysical log, UW1, and the northeast location consists of a geophysical log, UW2, and an adjacent continuous core, UWBH34-78 (Greenhouse and Karrow 1994). Geophysical logs reported by Greenhouse and Karrow (1994) include formation resistivity, natural γ, and neutron. Greenhouse and Karrow (1994) assessed the downhole geophysical logs for UW2 and divided the valley infill into three different “electrofacies” (Fig. 2a). These electrofacies were assumed to represent distinct depositional environments and sediment properties largely based on their clay content and total porosity (or pore size) as interpreted from the combined geophysical logs and correlated with the lithology observed in the continuous core record at UWBH34-78 (refer to Greenhouse and Karrow (1994) for details on electrofacies model development). According to Greenhouse and Karrow’s (1994) division of UW2, electrofacies A extends from approximately 0 to 18 m and includes Port Stanley, Maryhill, Catfish Creek, and Pre-Catfish Creek I, II, and III tills; electrofacies B extends from 18 to 35 m and consists of Pre-Catfish Creek Till IV; electrofacies C extends from 35 to 64 m and consists of a fine to medium sand. While these electrofacies are dominated by the main lithologic units they contain, interbeds of sand, clay, and gravel within and between the till units are expected within these electrical divisions. With no continuous core at the southwest location, Greenhouse and Karrow (1994) relied on a comparison of the downhole geophysical logs from UW1 with UW2 to interpret the infill in this section of the valley. Based on this geophysical comparison, Greenhouse and Karrow (1994) inferred that UW1 consists only of electrofacies C extending from approximately 50 to 15 m and B extending from approximately 15 m to surface, with electrofacies A absent at this location (Fig. 2b). A continuous core completed by the Ontario Geological Survey (OGS), BH11-OF-2008, approximately 3 km north of the study area (Fig. 1a), indicates the Quaternary sediment outside the Elora buried bedrock valley consists of approximately 8 m of Catfish Creek Till overlain by approximately 2 m of Maryhill and Port Stanley Till, along with interbeds of glaciolacustrine sediments (Burt and Dodge 2016). Bedrock depth control of the valley comes from UW1 and UW2 and two additional holes, Holes 3 and 4 (represented as H3 and H4; Fig. 1a) from Hilton (1978). Additionally, continuous bedrock core and geophysical logs from a nearby monitoring well (MW2; Fig. 1a) completed by Golder Associates (2013) provides an indication of the bedrock lithostratigraphic contact elevations near the transects (Table 1).

A regional-scale 3D Quaternary model was recently developed by the OGS between Orangeville and Fergus, Ontario (Burt and Dodge 2016; Burt 2018). This model was constructed using 43 continuous core logs, together with numerous low- and medium-quality borehole records (mainly from water well records) and a regional gravity survey. The OGS model shows the bedrock interface with the general trend of a bedrock depression from Belwood Lake toward Inverhaugh south of Elora (Figs. 1b and 1c), and the expected stratigraphy of the infill sediments. This model predicts that the Elora buried bedrock valley infill contains five main lithologic units, each classified as an aquitard or aquifer based on their sediment type from grain size (Table 2). However, the model remains to be verified in the Elora buried bedrock valley on a local scale.

The Elora buried bedrock valley and its Quaternary infill has been further investigated using a combination of co-located surface electrical resistivity and seismic refraction measurements, together with previously reported continuous core records and an electrofacies model of the valley infill. Our surface geophysical investigations focus on the lateral characteristics of the Quaternary deposits and bedrock along two high-resolution transects proximal to the historical data from Hilton (1978) and Greenhouse and Karrow (1994) (Fig. 1a). Based on the regional bedrock and Quaternary maps (Figs. 1b and 1c), valley incision along the northeastern “Road” transect is expected to be relatively deep and narrow, whereas the southwestern “Trail” transect is shallow and broad. These two transects, 3.5 km apart, provide contrasting conditions of the valley morphology and Quaternary infill architecture.

Electrical resistivity tomography (ERT) surveys were collected using a Syscal Pro Switch (Iris Instruments, Orleans, France) with a 96 electrode multicore cable. Multicore cables included 5 and 10 m electrode spacings that allowed for maximum electrode separations of 475 and 950 m, respectively. Electrodes consisted of 30 cm stainless steel spikes inserted approximately 20 cm into the ground. Electrodes were paced approximately 1 m from the shoulder of the Road and Trail with ground cover being primarily grassy vegetation. A Wenner–Schlumberger array was selected for its more balanced sensitivity of vertical and horizontal structures, as well as its higher signal-to-noise ratio compared to other commonly used arrays like dipole–dipole and Wenner (Everett 2013). Surveys were consecutively collected using 5 and 10 m electrode spacings along each transect. The 5 m spacing provided improved resolution of spatial variability within the upper 15 or 30 m for the Road and Trail transects, respectively, while the 10 m spacing provided greater depth of investigation along both transects. Arrays were designed with an expected maximum depth of investigation of 70 and 150 m for the 5 and 10 m arrays, respectively. Electrode arrays were rolled along in “leap-frog” fashion to cover the full length of the profile.

For each transect location, the 5 and 10 m data sets were combined into a single file for inversion. Pre-processing of apparent resistivity data involved data merging and removal of erroneous potentials based on the following criteria: extremely high and localized resistivity values; anomalously low measured potentials suggesting low signal-to-noise ratio; high standard deviation of repeat measurements. Apparent resistivity measurements were corrected for topography based on an interpolation of Differential GPS elevations collected along the transect. The elevation data were collected at a 30 m station interval, which was deemed reasonable given the limited topographic relief along the transects.

ERT inversions were completed using RES2DINV ver. 4.10.14 (Geotomo Software, Malaysia), which uses the Gauss–Newton least-squares method (Loke and Dahlin 2002). A robust inversion scheme (L₁-norm) was used with a homogenous initial model. A moderate damping factor of 0.15 was applied to both transects. The width of the model cells was set to 5 m. All other model parameters were selected to minimize the absolute error within four to six model iterations (Loke 2004).

Seismic refraction surveys were completed using a Geode Seismograph (Geometrics, San Jose, CA) with 48 channels and vertical component geophones with a natural frequency of 100 Hz spaced 5 m apart. A total of 85 and 136 shots were completed for the Road and Trail transects, respectively, using a 40 kg accelerated drop weight or a 10 kg sledgehammer (due to equipment failure, the drop weight was only used for the first half of the Road transect) with a 2 m lateral offset. Shots were stacked two to three times to suppress low-frequency noise. All data were recorded over a 235 ms window with a 0.125 ms sampling rate. Shots were gathered at the first geophone and every third geophone, thereafter, resulting in a 15 m shot interval along the spread. A final shot gather was collected at the end of the geophone spread. To achieve the full coverage of the targeted valley, lines were rolled along in “leap-frog” fashion, which resulted in repeat shots along the overlapping portions of the cable spread.

Seismic refraction data processing was completed in ReflexWin ver. 8.5 (Sandmeier Geophysical Research, Karlsruhe, Germany). Initial processing consisted of a time window cut of 130 and 140 ms for the Road and Trail transects, respectively, and a bandpass frequency filter with a lower cut of 5 Hz and an upper cut of 150 Hz. A background removal filter was applied for the Road transect to suppress low-frequency background noise associated with wind. This was achieved by calculating an averaged trace from the impacted area of the shot file and subtracting this trace from each individual trace within the impacted area. The first break arrival times of the refracted wave was manually selected at each geophone for each shot gather resulting in 2380 and 4252 detectable first arrivals for the Road and Trail transects, respectively. First arrivals were somewhat more difficult to delineate for the Road transect than the Trail transect, with first breaks selected for an average of 56% and 60% of geophone traces per shot for the two transects, respectively. Inversion of first break picks were completed using a seismic refraction tomography algorithm available in the ReflexWin software package. This inversion method automatically fits synthetic travel time databased on 2D velocity model to the selected travel times using the simultaneous iterative reconstruction technique and a curved ray travel time calculation with a finite-difference approximation of the eikonal equation (Sandmeier 2012).

Models were discretized at a 2.5 × 2.5 m (x, z) interval, which represents half the geophone spacing. Topographic variations along the lines were neglected due to the limited topographic relief along a single 48 geophone spread. However, a topographic shift was applied post inversion based on ground surface elevations. The decision to neglect topographic effects during the modelling process was reasonable given the relatively low slopes (less than 10%) along the transects (Sandmeier 2012). An initial starting model was specified with a velocity of 500 m/s at surface and a vertical velocity gradient of 100 m/s per model cell. The model domain was set for a maximum depth of 30 and 50 m below ground surface for the Road and Trail transects, respectively. Inversion parameters were varied to assess their overall impact on model results. Final parameter settings yielded a model that was consistent with the anticipated velocities of unconsolidated and consolidated materials. The quality of the tomographic results was further assessed through a forward ray tracing calculation, which provided a measure of the root mean squared error (RMSE) between measured arrival times and calculated arrival times. Tomographic inversions were repeated with different starting models and statistical conversion criterion to achieve a relatively low RMSE with a minimal number of iterations to reduce noise fitting within the model.

Electrical resistivity models reveal a zone of lower resistivity (<350 Ω m) underlain by higher resistivity (>550 Ω m) along the Road and Trail transects (Fig. 3). The steepest resistivity gradient correlates to the 350 Ω m contour interval, with this corresponding with available depth control points within the valley (e.g., UW2 and UW1). Based on the available core logs, these contour intervals likely represent the transition from unconsolidated sediments to consolidated bedrock.

The upper 10 m of unconsolidated sediment are characterized by resistivities between 25 and 125 Ω m. These lower resistivity sediments increase in thickness over the valley at the Road transect (Fig. 3a) and increase in thickness outside the valley along the Trail transect (Fig. 3b). A zone of higher resistivity sediments ranging from 125 to 325 Ω m appears within the uppermost section of the valley infill along the Road transect. However, these sediments are much thicker and more broadly distributed along the Trail transect, and also appear outside the valley. These electrical resistivity trends indicate the infill sediments along the Trail transect are more variable compared to the Road transect.

These resistivity data show that the valley incision is deeper and narrower along the Road transect with a slightly asymmetric U-shaped bedrock profile with steep sidewalls at 280 and 480 m lateral distances along the transect. These data also indicate minimal internal variability within the bedrock flanking the valley which decrease below the valley floor. In comparison, the Trail transect exhibits a shallower and wider valley profile with more internal variability in the bedrock zone. However, the valley profile is relatively less defined. The transition from unconsolidated to consolidated material along the Trail transect is characterized by an undulating interface depicting a wider valley profile with gently sloping sidewalls at 600 and 1000 m (Fig. 3b). Morphology of the buried valley at the two locations is shown most clearly using the same elevation and distance scales (Fig. 4a), which shows the elevation change in thalweg position at the two transect locations (335–320 masl). Overall, the bedrock along the Trail transect is much less resistive (<3000 Ω m) than that of the Road transect (>6000 Ω m). These contrasting conditions suggest that the rock surface along the Road could be more competent (i.e., less weathered or porous) than the Trail.

Seismic refraction measurements revealed multiple refracting events including the water table at ∼1500 m/s, as well as a second and third refractor at ∼2500 and ∼3500 m/s, respectively, based on visual inspection of the time-offset plots. Inversion of first arrival picks based on a 2D velocity model revealed a zone of low velocity between 300–460 m and 680–940 m distance along the Road and Trail transects, respectively (Fig. 5). These low-velocity zones correspond to the position of the buried bedrock valley feature identified in the resistivity surveys. The sharpest velocity gradient was encountered at approximately 3500 m/s, which likely represents the transition from unconsolidated to consolidated material. This value is within the range reported by previous workers locally (e.g., Lee 1975; Hilton 1978; Steelman et al. 2018). It should also be noted that the bedrock interface may be obscured by relatively high velocities associated with the Catfish Creek Till as well as a potentially heavily weathered or fractured bedrock interface. The valley walls and floor were not well constrained due to the steep slopes and limited depth of investigation.

These velocity models show minimal lateral variability within the uppermost low velocity regions (i.e., unconsolidated materials). The velocity of unconsolidated sediments was lower (560–3500 m/s) over the Road transect with a steeper velocity increase with depth (i.e., ∼200 m/s per m). The Trail transect exhibited a higher velocity (940–3500 m/s) with a more gradual velocity increase with depth (i.e., ∼100 m/s per m). Within the valley feature, the Road transect shows a lower velocity relative to directly above flanks of the valley, whereas along the Trail transect the velocity within the valley feature is similar to along the valley flanks.

For the 10 m below the 3500 m/s velocity contour, bedrock velocities are similar (both range to approximately 4000 m/s) along the two transects, providing minimal insight into the relative bedrock competency along the bedrock interfaces; the exact relationship between bedrock velocity and degree of weathering is unknown, and any velocity variations may be within the error expected within this region of the velocity models (i.e., RMSE ∼5% corresponding to ±200 m/s in the upper 10 m of bedrock). The velocity model for the Road transect does not extend deeper into the bedrock. For the seismic refraction survey method, a velocity increase with depth is required to generate critically refracted waves that can be recorded as first arrivals at increasing geophone offsets on the surface (i.e., beyond the “crossover distance” for that unit) (Everett 2013). The difference in depth of investigation could, therefore, result from either uniform or decreasing velocity bedrock at depth along the Road transect (i.e., critically refracted waves not generated) and increasing bedrock velocity along the Trail transect (i.e., critically refracted waves generated). Higher velocity could be associated with higher bedrock competency, potentially indicating a differing trend to that inferred from the resistivity measurements in this deepest region of the bedrock. There are, however, factors not related to the condition of the bedrock that may also contribute to the difference in depth of investigation for refraction surveys. These include effect of velocities and thicknesses of overlying units on the offset distance required for a critically refracted wave from a deeper unit to be recorded as a first arrival, with greater crossover distances increasing the likelihood that these arrivals are obscured by noise; geology of overlying units (e.g., sand vs. diamict) effects the degree of signal attenuation and the amount of energy remaining to resolve deeper units; signal-to-noise ratio of the data set with the shot-ground and geophone-ground coupling, and background noise conditions (i.e., wind and traffic) all influencing the ability to select late first breaks from deeper bedrock refractors. While assessing the relative contribution of these factors on the survey data sets is difficult, it can be noted that the noise level—resulting from windier conditions during data collection—on the Road transect was greater than along the Trail transect, likely limiting the ability to identify first arrivals from deeper refractors and resulting in a relatively lower depth of investigation.

Bedrock river morphology is often influenced by bedrock physical properties (e.g., Tinkler and Wohl 1998; Whipple et al. 2000). Locally, the valley thalweg appears to follow an erosion-resistant boundary that is associated with a lithostratigraphic unit of the Silurian Paleozoic sequence. Surface resistivity data show the thalweg is at an elevation of approximately 330 and 325 masl along the Road and Trail transects, respectively (Fig. 4a). The stratigraphy from a nearby borehole (MW2; Fig. 1a) was reinterpreted (Fig. 4b) based on photographs of the continuous core record and previously collected downhole geophysical logs (Golder Associates 2013), with three lithostratigraphic contacts identified at elevations comparable to the thalweg of the buried valley: Guelph Formation and Ancaster member of the Goat Island formation at 329.5 masl, Ancaster member and Niagara Falls member of the Goat Island formation at 326.6 masl, and Niagara Falls member and Gasport formation at 324.6 masl. Given the position of MW2 relative to the Road and Trail transects (Fig. 1a), together with knowledge of the bedrock strata in southern Ontario (Brunton and Brintnell 2020), a single erosion-resistant unit may underlie the entire bedrock reach of the Elora buried bedrock valley. The Ancaster member forms the caprock of much of the cuesta that makes up the Niagara Escarpment including the lip of the Niagara Falls (Brunton and Brintnell 2020), as well as forming the riverbed along sections of the Eramosa River, a modern-day bedrock river in the Guelph area (Steelman et al. 2015a; Steelman et al. 2015b; Capes et al. 2018). The unit is described as finely crystalline, siliceous dolostone with abundant chert nodules resulting from the predominance of siliceous sponges in the Goat Island Facies (e.g., Brunton and Brintnell 2020). In rock core, the dolostones of the Ancaster member are notably harder and more resistant to scratching than the overlying Eramosa Formation and underlying Niagara Falls member. Moreover, hydrogeological studies utilizing multilevel systems with many depth-discrete ports and (or) deployment of numerous pressure sensors with depth in and around the Guelph area (e.g., Nunes 2015; Capes 2017; Skinner 2019; Johnson 2020; Munn et al. 2020b; Nunes et al. 2021; Munn et al. 2020a) show hydraulic head profiles with a significant hydraulic head loss across the Ancaster member, indicating that it acts as a competent aquitard in many locations. Based on this evidence, the Guelph Formation could have supported valley incision until the thalweg reached the mechanically resistive Ancaster member. With the two transects showing similar cross-sectional areas for the valley (∼10 000 m²; Fig. 4a), and the Ancaster member impeding further vertical incision at the Trail transect, increased lateral erosion resulting in a widening of the channel profile may have occurred to accommodate the same volumetric flows. Given this hypothesis, these erosion-resistant layers, which often act as confining units or aquitards in groundwater flow systems, are likely important hydrologic features influencing both present-day and paleorivers.

Lithologic variability between the Road and Trail transects may have also contributed to the contrasts in local channel morphology. Surface geophysics revealed a narrower, deeper, and smoother valley profile along the Road transect compared to a wider and more complex profile along the Trail transect. Plucking of loose joint blocks may be one of the most significant erosional processes in bedrock rivers (Whipple et al. 2000), with more weathered and fractured rock likely offering lower resistance to lateral erosion. Along the Trail transect, the Elora buried bedrock valley has a complex cross-sectional profile, possibly indicating where lithological variability and bedding planes bound vertical joints influencing spalling surfaces that offered varying degrees of resistance to erosion. Surface resistivity measurements also revealed a marked decrease in bedrock resistivity in the upper bedrock units (i.e., Guelph Formation) that supports this hypothesized increased porosity and decreased competency change from northeast to southwest. This variability corresponds with the change from a relatively straight channel in the northeast to a more meandering channel in the southwest (Greenhouse and Karrow 1994; Conway-White et al. 2022) as well as the observed wider and shallower incision in the channel. The seismic velocity models for the bedrock interface shows similar velocities along both transects. While a lower seismic velocity would typically be anticipated for a more weathered bedrock, the exact nature of this relationship is unknown and the order of magnitude of velocity change may be within the associated error in the upper bedrock region of the two velocity models. A relatively higher bedrock velocity could be interpreted for depths >10 m below the interpreted bedrock interface along the Trail transect based on the relative depth of investigation of the two transects (see Section 4.1.2). This may demonstrate an opposing bedrock competency trend for the deeper bedrock region than that inferred from the resistivity measurements. However, given other factors could also account for the more limited depth of investigation along the Road transect, significant uncertainty associated with this interpretation of the deeper bedrock remains. Continuous coring into the bedrock below the Quaternary sediment interface at these two transect positions will be required to evaluate these bedrock trends and their relationships with the resolved valley morphologies.

Valley infills are often variable with complex successions of sediments with varying permeabilities (e.g., Kehew and Boettger 1986; Jørgensen and Sandersen 2006; Cummings et al. 2012; Steelman et al. 2018). Greenhouse and Karrow (1994) first described the Elora buried bedrock valley infill and, using geophysical logs (i.e., UW1 and UW2), reduced these geologic layers to a series of “electrofacies” A, B, and C, which they correlated with lithology from a nearby cored hole (UWBH34-78) (Fig. 2; see Section 2.2). While UW2 is broadly consistent with the ERT measurements along the Road transect, Greenhouse and Karrow’s (1994) electrofacies interpretation for UW1 (Fig. 2; Table 2) is not consistent with the trends observed along the Trail transect (Fig. 5) and the expected surficial geology outside of the Elora buried bedrock valley (e.g., BH11-OF-2008). The resistivity model resolves two distinct zones above the valley within the elevation range for electrofacies B at UW1. Here, a zone of lower resistivity (herein referred to as electrofacies A1) overlies a zone of higher resistivity (herein referred to as electrofacies A2). While electrofacies A2 has a similar resistivity as electrofacies B (i.e., 125–270 Ω m versus 130–190 Ω m, respectively), this unit is shown to be continuous across the entire Trail transect, suggesting it may represent a more widespread unit that extends beyond the confines of the buried bedrock valley. Given electrofacies B is inferred to be Pre-Catfish Creek Till IV (Fig. 2; Greenhouse and Karrow 1994)—a unit not expected outside the valley—the observed continuity suggests electrofacies A2 may instead represent the regionally extensive Catfish Creek Till that is expected as the basal unit above the bedrock outside the valley (e.g., a 8 m thick Catfish Creek Till is identified above the bedrock in BH11-OF-2008, approximately 3 km north of the study area). The continuity of this unit across the Trail transect contrasts with the Road transect where the higher resistivity unit is confined to the valley at elevations consistent with electrofacies B in UW2. Given this evidence, we interpret the upper tills (i.e., from bottom to top, Catfish Creek, Maryhill, and Port Stanley Till) as being present above the valley along the Trail transect, not Pre-Catfish Creek IV as previously proposed (Greenhouse and Karrow 1994). A transitional resistivity zone (herein referred to as electrofacies D) can also be added to represent a thin basal “gravel” layer or rubble top of rock layer located along the valley floor as described in core logs and drilling notes (Hilton 1978; Greenhouse and Karrow 1994) (Table 2).

These data were used to construct a 2D electrofacies model of the Quaternary deposits over the valley along the Road and Trail transects (Fig. 6) to help assess the infilling of the Elora buried bedrock valley. It is important to note that multiple cycles of deposition and erosion are likely, resulting in unconformities that may not be resolvable using surface geophysics. Further, it is also likely that interbedding of sand, clay, and (or) gravel typical of water-lain sedimentary facies (e.g., UWBH34-78 in Greenhouse and Karrow 1994) occurred within and between the electrofacies, although these would be difficult to resolve given inherent limitations in surface geophysical measurement resolution.

Electrofacies “D” represents the oldest infill unit with the ERT models showing resistivities ranging from 200–350 Ω m directly above the bedrock (Fig. 6). The presence of this electrofacies is inferred from drilling notes from Hilton (1978) that document a marked loss in drilling fluids near the valley floor, suggestive of a transmissive sand/gravel layer or weathered bedrock zone. This unit would be consistent with a contact aquifer and represents a transitional zone between electrofacies C and the competent bedrock. This unit has been described as coarse boulder gravel and was thought to be fluvial in origin (Hilton 1978; Greenhouse and Karrow 1994). Given the limitations of the ERT method at depth, the thickness of this unit cannot be determined here, but it is anticipated to be 3–5 m and likely limited to the floor of the valley.

The basal gravel rubble layer is overlain by electrofacies C, a relatively thick sequence comprising water-lain fine to medium sand with resistivities ranging from 80 to 130 Ω m. This supposed aquifer unit is described by Hilton (1978) as a clayey sand. Given the uniform nature of this sediment, Greenhouse and Karrow (1994) concluded that it had a lacustrine origin with pollen analyses indicating deposition during an interstadial period. The date of this sediment remains unknown but based on many ice fluctuations depositing several different till units above, the unit is proposed to be at least as old as the last interglacial period (>100 000 years ago) (Greenhouse and Karrow 1994). Based on the ERT results for the two transects, this unit may be more homogeneous at the northeast end of the valley (Fig. 3a) but becomes progressively more heterogeneous and laterally discontinuous toward the southwest (Fig. 3b). The origin of this lateral heterogeneity is not clear; however, the areas of lower and higher resistivity within the sediment may be interpreted as variations in clay and sand content. In addition, while the valley bottom is relatively smooth along the northeastern Road transect (Figs. 3a and 6a), the valley cross-sectional profile becomes more complex with undulations and bedrock mounds at the more southwestern Trail transect, coincident with a decrease in bedrock valley height and increased width (Figs. 3a and 6b). The thickness of the unit may vary as a result of the uneven bedrock topography, with thicker deposits infilling between the bedrock mounds.

The sandy unit (electrofacies C) is overlain by electrofacies B (130–190 Ω m) or A2 for the northeastern Road transect or southwestern Trail transect (125–270 Ω m), respectively. These units are both expected to be stiff, sandy tills with the Pre-Catfish Creek IV Till contained within the valley along the Road transect and the Catfish Creek Till overlying the shallower valley along the Trail transect. Since Pre-Catfish Creek Till is proposed to have formed from a major ice advance (Greenhouse and Karrow 1994), it is possible this unit was initially deposited across the entire study area with subsequent erosion and (or) incorporation into more recent sediments. This would explain its absence farther southwest along the valley; the northeastern portion of the valley is deeper and narrower than the southwestern portion, offering better preservation of Pre-Catfish Creek Till. Lateral variations in these sediments can be seen along the transects with resistivity varying over distances of several hundred metres.

The uppermost sediments over the valley axis and flanks consist of the Catfish Creek, Maryhill, and Port Stanley Till that are defined by electrofacies A on the Road transect and A1 and A2 on the Trail transect with a resistivity range of 25–270 Ω m (Fig. 6). While the electrical properties of the Maryhill and Port Stanley Till are uniform along both transects, the resistivity of the Catfish Creek Till may vary (Table 2). For the southwestern portion of the study area, the Catfish Creek Till can be distinguished from the overlying Maryhill and Port Stanley Till by the relatively higher resistivity (i.e., >100 Ω m) anticipated for a sandier or overconsolidated diamicton. However, for the northeastern section of the study area, the Catfish Creek Till cannot be differentiated using ERT from the other two upper tills despite a relatively thick 5–10 m sequence of sediments likely to be present above the bedrock outside the valley (e.g., BH11-OF-2008; Burt and Dodge 2016). Such resistivity trends are consistent with the expected properties of the Catfish Creek Till, which are known to vary (Karrow 1967; Burt and Dodge 2016). For unconsolidated sediments, an increase in resistivity could be caused by multiple lithological factors including a decrease in clay content, a decrease in porosity (e.g., an increase in cementation or overconsolidation of the sediment), or a decrease in porewater conductivity. While direct analysis of the sediments would be required to verify the cause of the observed resistivity contrast, it can be noted that Karrow (1967) identified a typically fine-grained Catfish Creek Till in the Guelph area, as well as noting similar occurrences near Kitchener, Ontario. This textural variability may result from the pre-existing sediments and (or) bedrock in the area that were incorporated into the Catfish Creek Till and the timing of deposition, with earlier Catfish Creek Till tending to be more clay rich (Karrow 1967; Burt and Dodge 2016). Therefore, it can be hypothesized that the Catfish Creek Till in the northeastern part of the study area may be more clay rich, becoming sandier farther southwest along the valley axis.

As observed in the resistivity models, the contact between the Catfish Creek Till (electrofacies A2; Fig. 6b) and the overlying Maryhill Till (A1; Fig. 6b) can be quite complex within the study area. Undulations along this interface are interpreted to indicate variable erosion by glacial meltwater that incised into the Catfish Creek Till. Interbedded sands described in the cores by Greenhouse and Karrow (1994) also support such an erosional contact between the different till units. The contacts between the different units are more complex for the southwestern sections of the study area than in the northeast. A series of discontinuous higher resistivity (∼300 Ω m) zones can be seen along the top of the Catfish Creek Till (A2) above the bedrock valley (Fig. 6b). This feature could be interpreted as a braided river system that would have followed the valley. This suggests the continued importance of fluvial processes as the valley infill continued to evolve over time.

Bedrock and valley infill interpretations based on surface geophysical measurements and historical core records were used to generate an electrofacies and lithologically informed lithostratigraphic conceptualization of the Quaternary deposits and Silurian bedrock of the northeastern and southwestern extent of the Elora buried bedrock valley south of Elora, Ontario (Fig. 7).

Although the Guelph, Goat Island, Gasport, and Cabot Head formations are not readily differentiable within these surface geophysical measurements, resistivity measurements did provide insight on relative bedrock competency around the incised valley, particularly near the bedrock–Quaternary interface. Variations in bedrock competency (e.g., weathering, fracturing, and enhanced porosity) will control channel morphology and hydrogeological properties. A more competent bedrock—containing fewer fractures and lower porosity—would be expected to have relatively higher electrical resistivity and seismic velocity (e.g., Powers et al. 1999; Rønning et al. 2014). While the seismic velocity is similar in the upper ∼10 m of bedrock along the two transects, and therefore provides minimal insight here, the bedrock does exhibit decreasing resistivity from northeast to southwest along the valley. This may indicate an increase in bulk porosity (fractures and matrix) and possibly fracture interconnectivity down the axis of the valley thalweg; consequently, zones of higher transmissivity may be more developed at the bedrock contact along the southwestern portion of the valley. For the bedrock region >10 m below the bedrock interface, the greater depth of investigation of the seismic refraction survey along the Trail transect may be interpreted as a higher velocity. However, given the relatively noisier data set along the Road transect, it is unclear if the difference in depth of investigation can be attributed to actual differences in bedrock conditions, particularly when considered next to the clear but opposing trends observed in the resistivity models. Along each of the two transects, the lowest and most variable resistivity is observed closest to the valley, with higher resistivities farther away from the valley (Fig. 3). This is hypothesized to represent a zone of more fractured or weathered rock, which would be consistent with the observations made by Cole et al. (2009) who conceptualized enhanced dissolution and weathering of the bedrock in proximity to buried bedrock valleys to occur as a result of changes in hydraulic pressures associated with valley incision. Depending on the nature of the sediments infilling the paleovalley, this increased bedrock transmissivity may lead to enhanced bedrock recharge from overlying Quaternary sediments.

Construction of a hydrostratigraphic model of the Quaternary deposits based on the geophysical measurements requires an approach to classify observed petrophysical variability (i.e., electrofacies) as either aquifers or aquitards. In general, variations in the resistivity of sediment are related to pore volume, pore surface area, and sediment mineralogy, with these properties also controlling the hydraulic conductivity of the sediment (Slater 2007). Clay minerals decrease the resistivity of sediments, while coarser sediments such as sand and gravel increase resistivity (Everett 2013). In the absence of direct hydraulic measurements, the clay fraction within a layer is often used to infer whether it will act as an aquifer or aquitard as a first-order hydrostratigraphic interpretation. Identifying aquifer and aquitard sequences from resistivity surveys using “thresholds” or defined resistivity ranges has been used by Morgan et al. (2019) with some success for developing a lithostratigraphic model of a buried bedrock valley network. However, geophysical core logs from Greenhouse and Karrow (1994) and Hilton (1978) indicate that the resistivity of the aquifer unit within the valley also contains some clay (i.e., it is not a clean sand), and thus, may display a lower resistivity than would be anticipated for a typical aquifer. With this context, we correlate electrofacies units from this study to existing regional lithostratigraphic units used by (Burt and Dodge 2016) to infer hydrostratigraphic units, but caution that the latter deserves hydraulic data for confirmation to generate an appropriate level of confidence in the hydrogeologic conceptual model.

The OGS lithostratigraphic model represents the most complete conceptualization of the Quaternary deposits around and within the Elora buried bedrock valley to date, with five different Quaternary units represented within the study area. OGS Quaternary units (database release from Burt and Dodge 2016) were correlated to electrofacies based on the existing core and geophysical logs (Greenhouse and Karow 1994) and the surface geophysics of this study. Overall, the OGS lithostratigraphic classifications are consistent with the expected lithological properties of the electrofacies from Greenhouse and Karrow (1994) and Hilton (1978) with two notable exceptions. First, electrofacies D is thought to be a highly transmissive basal gravel rubble unit (Hilton 1978); this type of sediment would be significantly more permeable than the “Older Till” (ATG1) assigned to the base of the valley and described as a valley till with aquitard properties in the OGS model (Burt 2018). This apparent discrepancy cannot be resolved with the available surface geophysical data as these two lithologies could have similar electrical properties. The interpretation of electrofacies D is more consistent with the Contact Aquifer (AFH1) mapped elsewhere in the OGS model. Second, electrofacies B is described as a Pre-Catfish Creek IV Till by Greenhouse and Karrow (1994), while the OGS model indicates this layer to be Catfish Creek Till (ATC1). Interpretation of the Pre-Catfish Creek Till as part of the Catfish Creek Till by the OGS follows a tendency to divide tills into fewer distinct units than was historically common (Burt (personal communication, 2020)), what was previously considered distinct tills created by separate ice advances is now often recognized as belonging to a single unit characterized by internal variability. In addition, the scale of the OGS model necessitated simplification of the number of modeled units, and the purpose of the model to represent lithologically derived hydrostratigraphy did not require differentiation of the Pre-Catfish Creek Till within the model. Therefore, a pragmatic approach that combined the Pre-Catfish Creek Till with the Catfish Creek Till unit was taken by the OGS (Burt (personal communication, 2020)).

The level of detail identified in the geophysical transects over the Elora buried bedrock valley builds upon the existing regional-scale OGS model in several ways. First, the Catfish Creek Till, while previously known to vary, is suggested from the geophysical results to be more clay-rich in the northwest along the Road transect, becoming sandier in the southeast at the Trail transect. Second, the Catfish Creek Till above the valley in the southwest is observed to have lenses of coarse-grained sediments; these lenses were not observed farther to the northeast and are expected to impact local permeability and potential contaminant pathways. Finally, the more permeable sandy unit infilling the valley is seen to be more variable in terms of clay and sand content between the two transects and likely more discontinuous in the southwest than the northeast where the unit is more uniform.

This study demonstrates the application of a multidimensional interpretive framework using surface electrical resistivity and seismic refraction tomography together with a continuous core and electrofacies model based on geophysical logs in the context of a regional lithostratigraphic model to improve our understanding of bedrock morphology and infill architecture of the Elora buried bedrock valley in Centre Wellington, Ontario. The co-located surface geophysical surveys completed along two transects provided spatial information needed to extrapolate interpretation of conventional core and downhole geophysical datapoints, resulting in a more robust conceptual model of the Elora buried bedrock valley morphology and infill history.

Surface geophysics revealed a linear U-shaped valley in the northeast that deeply incised bedrock but becomes progressively shallower and broader downslope towards the southwest. The valley thalweg elevation may be controlled by bedrock mechanical properties, with the erosion-resistant Ancaster member of the Goat Island formation hypothesized to impede downward vertical erosion. The hydrogeological condition imposed by the Ancaster member during incision may have further influenced valley morphology, with the river channel widening as the river reached this unit and was forced to eroded laterally to accommodate volumetric flows from further upstream. Notably, this hypothesis would support the growing evidence for this lithostratigraphic unit being an important aquitard surface throughout the region. Surface geophysical measurements revealed variability in the physical properties of the incised bedrock, also likely influencing valley morphology. A reduction in bedrock competency along the valley thalweg as inferred from the resistivity may have also contributed to the broader valley profile with a more pronounced weathered zone toward the southwest. As suggested by Cole et al. (2009) and others (e.g., Shaver and Pusc 1992; Kunert et al. 1998; Smerdon et al. 2005; Steelman et al. 2018), enhanced understanding of bedrock valley morphology and associated spatial variability of bedrock hydraulic properties proximal to buried bedrock valleys could have important implications for future groundwater resource assessments in this region.

Electrical resistivity tomography measurements of the infill sediment revealed a discontinuous low resistivity unit and the presence of coarse-grained sediment within units that would otherwise be conceptualized as continuous aquitards. The regionally widespread Catfish Creek Till exhibited increasing resistivity from northeast to southwest across the study area, indicating variability in clay content or degree of consolidation within this unit. Further, the presence of potential Pre-Catfish Creek sediments in the northeast, not included in Burt’s (2018) conceptualizations of the Elora buried bedrock valley, was confirmed through these geophysical measurements. The absence of Pre-Catfish Creek sediments in the southwest, a unit previously suggested within and over the Elora buried bedrock valley instead of Catfish Creek, Maryhill, and Port Stanley Till, was also identified. These examples demonstrate the value of incorporating surface geophysics with conventional core and downhole geophysical measurements to improve resolution and accuracy of hydraulically relevant lithostratigraphic units within a buried bedrock valley as well as below the valley floor.

Although the regional hydrogeological implications of valley incision and infilled sediments remain to be assessed, the surface geophysical methods demonstrated in this study offer an approach to better resolve lithologic and physical property variability that are relevant for 3D flow system and well-field hydrogeological investigations. Geophysical data, in combination with available core logs and regional geologic knowledge, can be used to establish a more robust lithostratigraphic model that, once hydrologically calibrated, provides an essential framework for a hydrostratigraphic model of an area. Our proposed conceptualization of the Elora buried bedrock valley will be valuable for groundwater practitioners interested in understanding the influence of bedrock valleys on groundwater flow systems within well-informed lithostratigraphic frameworks. Geophysical measurements may be most useful where a detailed conceptualization of the valley morphology and geometry is required to assist with the development of groundwater flow models recognizing the importance of preferential zones of recharge, flow pathways that support municipal water budgets, sustainable yield estimates, and local area risk assessments. The style of physical property variability observed electrophysically, yet calibrated to lithostratigraphy, will provide 2- and 3D frameworks for hydrologic calibration and validation of regional-scale hydrostratigraphic models with improved representation of site conditions and prediction.

Please contact the authors for research data availability.

Conceptualization: OCW, CMS, EA, HU, JDM, BLP

Data curation: OCW, CMS

Formal analysis: OCW, CMS

Funding acquisition: CMS, HU, BLP

Investigation: OCW, CMS, EA

Methodology: OCW, CMS, EA, HU, BLP

Project administration: CMS

Resources: CMS, BLP

Supervision: CMS, EA, BLP

Validation: OCW, CMS

Visualization: OCW

Writing – original draft: OCW

Writing – review & editing: CMS, EA, HU, JDM