Reflection and Love-wave imaging of a buried valley using 2D 3C land-streamer seismic data Open Access

We conduct reflection and Love-wave imaging on three 2D, 3C land-streamer seismic lines to characterize a buried bedrock valley underlying glacial deposits situated within southern Ontario, Canada. Understanding the valley’s features is important for designing above-ground structures, investigating groundwater and surface water interactions, and understanding erosion processes driven by glacial cycles. These characteristics include width, extension, overburden thickness, spatial facies distribution, and geometric configuration of the overburden-bedrock contact. Reflection imaging on the west-southwest–east-northeast trending profiles reveals the bedrock geometry and quaternary sediment stratigraphy, including a shallow reflection associated with a lithostratigraphic transition zone. We use the dispersive properties of Love waves to derive pseudo-2D S-wave velocity profiles along the seismic transects. The inversion process is constrained by the P-wave refraction velocities from reflection processing. We emphasize the significance of data preconditioning prior to surface-wave analysis and develop an optimal processing sequence for moderately to highly noisy gathers. Likewise, we construct bedrock-reaching velocity profiles for depth conversion of reflectivity sections by using S-wave velocity from Love-wave inversion, regression analysis, and water well data. Our signal preconditioning and integrated velocity modeling approach can be implemented in other areas where land-streamer seismic surveys are available, enhancing quantitative imaging and depth conversions. Blind borehole records demonstrate average depth-conversion errors of less than 2% following this approach. The Love-wave velocity imaging and the body-wave reflection imaging correlate well with each other and with lithologic changes observed in water wells. This is demonstrated by the alignment of velocity contrasts with the geometry of the shallow reflection event at the lithostratigraphic transition zone. We find that the bedrock valley is approximately 60 m deep with a southward thinning width, possibly due to the valley’s geometry shifting from a north–south to an east–west orientation.

Canada is evaluating a site in South Bruce, southwestern Ontario, as a potential location for a deep geological repository for spent nuclear fuel. Although the project’s primary focus is on deep subsurface characterization, near-surface studies are being conducted concurrently. This includes the study of a known buried valley consisting of a heterogeneous combination of fluid-saturated, porous, and permeable sand and gravels intermixed with silts and clay layers. Understanding the structural and stratigraphic characteristics of the valley is important for the engineering and architectural planning of above-ground structures, as well as for developing hydrogeologic models of the transient dynamics of groundwater and surface water interactions. Likewise, knowledge of past glacial erosion in the study area is essential for anticipating the effects of erosion processes induced by glacial cycles at repository level, such as valley overdeepening and paleochannels. The 3D and 2D seismic data from land streamers have proven to be very successful in near-surface high-resolution imaging (Dolena et al., 2008; Andre et al., 2009; Pugin et al., 2013; Maries et al., 2017; Brodic et al., 2021; Papadopoulou et al., 2023). The documented success of land-streamer data motivates us to acquire multiple 3C 2D seismic land-streamer profiles in the South Bruce site using a horizontal transverse seismic source. The acquisition of a 3D seismic array does not offer substantial advantages over a 2D configuration in our investigation area. This is because there is minimal heterogeneity along the strike direction of the valley, based on available borehole data. In addition, the shooting lines are parallel to the direction of maximum inclination, along which lithologic changes predominantly occur. These conditions result in minimal cross-dip effects and negligible cross-profile velocity variations, which in a predominantly horizontally layered medium validate the use of conventional common midpoint (CMP) reflection processing and surface-wave analysis using 2D geometries (Foti et al., 2018; Villamizar et al., 2022).

We demonstrate the land-streamer potential for high-resolution S-wave reflection imaging and velocity (⁠$VS$⁠) inversion of Love waves. We implement multicomponent reflection processing along three land-streamer seismic profiles. We find that the shear transverse (ST) component (with transverse source orientation) produces higher-quality reflectivity images, exhibiting enhanced continuity and coherence of reflection events, compared with those derived from the shear radial (SR) or vertical component. Similarly, Love waves on the ST component are observed and used for $VS$ inversion. Although there is extensive research on the use of Rayleigh waves for S-wave velocity inversion (Brodic et al., 2015; Ikeda and Tsuji, 2020; Wang et al., 2023; Niu et al., 2024), research demonstrating the effectiveness of Love waves in quantitative velocity imaging is comparatively limited (Pan et al., 2016; Foti et al., 2018; Li et al., 2024).

However, there are documented advantages to using Love waves for S-wave inversion. For instance, the inversion of their dispersion data tends to be more stable than that of Rayleigh waves (Safani et al., 2005). They are less reliant on initial models and produce dispersion images with a higher signal-to-noise ratio (S/N) (Xia et al., 2012; Dal Moro and Mazanec, 2024). In addition, they are not as affected by variations in P-wave velocities as Rayleigh waves are, and Love waves are less likely to generate higher dispersion modes (Foti et al., 2018). Although uncorrelated seismic gathers in their raw form are generally adequate to produce dispersion images (Park et al., 1999), we find that shot gathers exhibiting moderate to strong presence of noise (e.g., reflected waves, guided waves, multiples, converted modes, and anthropogenic noise) require additional treatment. We propose a signal processing workflow that proves essential in obtaining a coherent distribution of phase velocities with frequency for noisy shot gathers. This workflow is grounded in traditional surface-wave processing techniques, such as the elimination of near-offset data, but it also incorporates velocity-designed surgical mutes and frequency-wavenumber filtering for the removal of high-wavenumber components. By following this processing sequence and adopting the surface-wave analysis method introduced by Park et al. (1999), we demonstrate the effectiveness of Love waves in retrieving high-quality $VS$ images of shallow subsurface layers.

The primary focus of this case study lies in the integration of body-wave (shear) reflection and surface-wave imaging techniques to characterize the buried valley within the South Bruce site. In the body-wave reflection analysis, we conduct structural and stratigraphic interpretation by examining reflectivity seismic sections in the time domain. In the surface-wave analysis, we use the quantitative imaging of Love waves to interpret geologically meaningful velocity contrasts in the first few meters of the subsurface, and we establish a velocity modeling workflow for time-to-depth conversion of reflectivity seismic sections. Given that seismic data is collected in the time domain, its conversion to depth introduces uncertainty into the interpreted lithostratigraphic framework. Depth conversion is particularly challenging when targeting shallow layers (less than 100 m) due to limited borehole data with velocity information. Common practice involves using normal moveout (NMO) velocities or depth-to-bedrock markers from shallow boreholes, such as water wells, to construct velocity functions (Miller et al., 1999; Keho et al., 2002; Oldenborger et al., 2016; Dietiker et al., 2024). However, the accuracy of this approach is heavily reliant on the number of available boreholes and their proximity to the seismic transects. In this case study, we propose a velocity modeling workflow whereby the quantitative imaging of Love waves is integrated with regression analysis and information from water wells to produce bedrock-reaching pseudo-2D velocity models for depth conversion. Our approach better accounts for the true depth positioning of geologic features, enabling our ability to describe the bedrock top geometry and delineate the stratigraphy and facies distribution within the overburden.

The study area is located 3 km west of the town of Teeswater, within the South Bruce township, in Southern Ontario, Canada. Southern Ontario is underlain by thick successions of undeformed Paleozoic sedimentary rocks, ranging in age from Cambrian to late Devonian and Mississippian in some regions (Carter et al., 2021). Bedrock underlying the site consists of the Devonian-aged Lucas Formation and is unconformably overlain by variable thickness of quaternary glacial sediments. The Lucas Formation is regionally extensive comprising limestone, dolostone, and anhydritic beds with local sandy limestone (Armstrong and Carter, 2010). Most glacial sediments in the area have been deposited by the Huron and Georgian Bay lobes of the Laurentide Ice Sheet during episodic advances and retreats during the Late Wisconsinan glacial episode, dated between 23,000 to 10,000 years ago (Barnett, 1992). These glacial episodes are responsible for shaping the current landscape of the study area.

Figure 1 shows the spatial distributions of glacial deposits within the study area consisting of glacial fluvial outwash, lacustrine sediments, and ice-contact drift, predominantly comprising sand and gravel, with some silt content leading to elevated permeabilities (Ontario Geological Survey, 2012). Glacial tills identified in the study area include the Elma and Rannoch tills. Cowen et al. (1978) and Sharpe and Edwards (1979) suggest the Elma till signifies continuous ice coverage from 25,000 to 14,000 years before present, characterized by a stony to sandy silt till. In contrast, the younger Rannoch till consists mainly of silt till deposited as ground moraine (Cowen and Pinch, 1986). Approximately 70% of the study area is composed of sands and gravels originating from glaciofluvial or glaciolacustrine processes, formed during the retreat of the late Wisconsinan glacial event. Glaciofluvial deposits, associated with substantial meltwater discharge from glaciers, typically transport coarse-grained sediment. Both ice-contact stratified drift deposited proximal to glaciers, and outwash deposits further downstream through channels or river deltas, are recognized in the area. Glaciolacustrine deposits represent sediments carried by glacier meltwater and subsequently settled in glacier-fed lakes. The remaining deposits within the study area are associated with more recent Holocene events, resulting in deposition within various bogs, swamps, and river flood plains.

Overburden thickness mapped by Gao et al. (2006) across southern Ontario show thicknesses exceeding 200 m with average values ranging from 30 to 60 m. Based on water well records within the study area (water well information system [WWIS]), glacial deposits have an average thickness of 21 m, which, however, exceeds 50 m thick along a north-trending corridor associated with a segmented northern extension of the Hutton Heights bedrock valley. Bedrock valleys exist beneath the quaternary deposits throughout much of southern Ontario (Russell et al., 2007, 2018; Gao, 2011; Oldenborger et al., 2016). Based on Gao (2011), bedrock valleys result from erosion from subglacial meltwater beneath the Huron ice lobe, probably during or shortly after the Late Wisconsinan glacial maximum. Near the town of Teeswater, overburden is thin and Paleozoic bedrock is exposed.

The database of water wells (WWIS) provides critical information on the type of subsurface material. This network was augmented by the addition of seven newly drilled vertical monitoring wells through overburden extending into the underlying Paleozoic bedrock (see Figure 1) (Geofirma, 2023). Drilling was completed using a 0.203 m (8 in.) diameter hollow stem auger with split spoon sampler to obtain a continuous log of overburden (Figure 2). Overburden logs are overwhelmingly rich in sand and sandy gravel, accounting for approximately 80% of the well intervals, with the remaining 20% comprised of thin intermixed silt and clay layers. This is consistent with mapped sand and gravel-rich glaciofluvial and glaciolacustrine deposits (Figure 1). Sand and sandy gravel intervals were logged with mean thickness of approximately 7 m, separated by thinner silt and clay lenses (approximately 3–4 m thick). A suite of geophysical logs was also acquired with natural gamma ray measured for overburden through a polyvinyl chloride casing. Elevated natural gamma logs are locally consistent with silt and clay content, within the predominantly sand and gravel deposits, and are used for stratigraphic interpretation.

The land-streamer seismic surveys were acquired along west-southwest to east-northeast trending profiles adjacent to Bruce Road 6, Concession Road 8, and Concession Road 10 (Figure 1). The surveys were strategically positioned to transect the bedrock valley buried beneath the overburden of the investigation area. A vibroseis system, equipped with a horizontally transverse oscillating plate, serves as a shear source across all 1942 shotpoints. We use a stream of 72 3C 10 Hz receivers, each separated by a uniform distance of 1.5 m. In seismic imaging, especially in the context of surface-wave analysis for velocity inversion, the use of receivers with a low natural frequency is preferred (below 5 Hz). These receivers facilitate the probing of deeper layers and allow greater coverage of spatial bandwidth by retrieving lower wavenumber features. Nonetheless, research has shown that 10 Hz receivers can record Rayleigh and Love waves with high S/N at frequencies below their fundamental vibrational frequency (Ivanov et al., 2006; Comina et al., 2017; Mahvelati et al., 2017). The nearest source-receiver offset of the streamer is set to be 4.5 m, establishing an array ranging from 4.5 to 111 m. This configuration mitigates the influence of near-field effects, reinforcing the validity of the plane-wave assumption inherent in surface-wave analysis methods (Park et al., 2002), while preserving the shorter wavelengths needed for optimal migration of reflected arrivals (Mora, 1989).

The vibroseis unit was operated to vibrate at intervals of 3 m along each of the three profiles shown in Figure 1. After multiple sweep tests with varying vibration lengths and sweep frequencies, we determine 1.2 s as the optimal length, with usable frequencies between 8 and 120 Hz. Production sweeps consisted of two horizontal sweeps, transverse to the shooting orientation and 180° out of phase. Each shotpoint was thus excited by two sweeps, initiated in opposite directions. For horizontal-component data sets, the two out-of-phase records are correlated and subtracted, whereas for vertical-component data sets, they are correlated and added together. This double-sweep acquisition approach enhances the S/N in all data sets, reducing P-wave arrivals in the horizontal components and mitigating random noise in the vertical components. Table 1 provides the acquisition parameters for the three land-streamer profiles.

Despite using 3C geophones, we observe that the quality of S-wave data within the transverse component (ST) generally surpasses that of the radial (SR) and vertical component (see also Pugin and Yilmaz, 2019). This is shown in Figure 3a and 3b, which presents raw shot gathers for each component at shotpoints located around the central segment of line 2 and eastern segment of line 1, respectively, along with their corresponding amplitude spectra. Note that the time axis for all seismic data presented in this study represents two-way traveltimes (TWT). The SR component in Figure 3a exhibits contamination of converted modes within the 200 to 500 ms window, particularly at mid- to long-offsets (red arrows). These modes adversely impact reflection and Love-wave imaging by masking reflected arrivals and distorting the dispersion properties of Love waves. In contrast, the ST component shows a clear distinction between the refracted S wave and the lower velocity, notably more coherent package of surface waves, with virtually no converted mode interference. Likewise, the ST shot gather in Figure 3b shows a set of reflection events at approximately 300–500 ms (red arrow), likely originating from either the top of bedrock or slightly shallower lithologic units within the overburden. These events are only faintly discernible on the SR and vertical component. These observations highlight the inherent differences in how each component captures seismic energy. Although seismic events are generally expected to be traceable across all components, certain phases may be more prominently recorded in one component while appearing attenuated or even absent in others due to factors such as propagation effects and component-specific sensitivities. The normalized amplitude spectra for all shot gathers exhibit a consistent pattern, where at least 20% of the spectral strength spans a frequency range from 8–10 Hz up to 120 Hz. The ST component tends to have a somewhat bimodal frequency distribution, with elevated power concentrations observed in low- and high-frequency bands.

We process the 2D 3C seismic land-streamer data through two complementary workflows: reflection imaging and quantitative imaging. The reflection imaging workflow focuses on body wave (S wave) processing, emphasizing the analysis of reflected arrivals to produce high-resolution migrated seismic images of the subsurface. In contrast, the quantitative imaging workflow uses surface-wave (Love wave) analysis to derive S-wave velocity profiles of the shallow subsurface. These velocity models, constrained by borehole data, are essential for constructing bedrock-reaching velocity models and facilitating the time-to-depth conversion of the migrated seismic images. Figure 4 shows this integrated workflow for S-wave imaging, highlighting the complementary roles of reflection and quantitative imaging in subsurface interpretation.

The final depth-domain seismic sections enable the interpretation of the bedrock-overburden architecture. This includes delineating the geometry of the bedrock-overburden contact, identifying lithostratigraphic variations within the quaternary overburden, and characterizing subsurface features of geologic significance.

We use a conventional imaging strategy for the reflection processing of the land-streamer data set, similar to that used by Hunter et al. (2022) and Villamizar et al. (2022). Table 2 provides the main processing steps. We begin by applying cross-correlation of the vibroseis pilot sweep with the recorded signal to derive a meaningful representation of the subsurface impulse response. This step is followed by data reformatting and the assignment of geometry information to each trace within the data sets. Subsequently, the different components are separated into transverse, radial, and vertical. To remove direct and refracted arrivals, muting is applied. Frequency filtering is then implemented to mitigate the presence of converted modes, air waves, and surface waves, as well as to enhance S/N of reflected arrivals. Static corrections are not applied due to the shallow survey and nearly flat-lying stratigraphy, which minimized the impact of long-wavelength static effects. Trim statics were applied within common depth point (CDP) gathers during velocity analyses to mitigate short-wavelength static effects. Velocity analyses are conducted on a semblance map every 50 CMP gathers to facilitate subsequent NMO correction. The flattened data sets are stacked to produce a 2D zero-offset image beneath the processing line. Finally, this image is migrated using a finite-difference poststack time migration algorithm. Although we produce stack sections using the vertical and SR component, their quality is suboptimal compared with the ST component in imaging the bedrock top. The vertical and SR components exhibit significant energy defocusing and amplitude smearing. Consequently, we have chosen not to include them here.

The final ST-component poststack time-migrated seismic profiles are shown in Figure 5 for all lines. No horizon interpretation is overlain on the images to aid the reader’s visualization. Two prominent reflections are identified: the shallow reflection (R1) and the overburden-bedrock reflection (R2). Figure 5 also highlights additional features, such as bedrock multiples (M), coherently stacked noise (N), areas with significant amplitude smearing (S), and local topographic changes (T1 and T2) of R1. The local topographic changes (T1 and T2) are particularly relevant for correlation with the quantitative imaging of surface waves. The reflection identified as a multiple (M) of the overburden-bedrock reflection consistently appears at a time delay that reflects the additional travel time due to multiple interactions with subsurface interfaces. This multiple generally has a weaker amplitude compared with the primary reflection. Moreover, borehole drilling data from locations between the survey lines, extending down to 100 m depth, show no significant variations in lithology or rock type beneath the bedrock contact. This strongly supports the interpretation that no distinct reflector exists below the overburden-bedrock contact, confirming that the observed reflection is a multiple rather than a separate geologic interface.

We find that the shallow reflection R1 is the most consistent and continuous reflection event across all three seismic lines, aside from the overburden-bedrock reflection. R1 is generally flat lying; though it appears smeared in some sections, particularly in the middle segments of lines 2 and 3, R1 remains recognizable across the entire length of the transects. Similarly, R2 clearly depicts the geometry of the bedrock top, which has the greatest impedance contrast and is readily identifiable along the full extent of the profiles. We indicate several areas suffering from amplitude smearing (S), likely stemming from poor signal alignment and energy defocusing during stacking. This phenomenon is largely associated with lateral compositional and lithologic heterogeneities that cause higher amplitude attenuation, resulting in smeared reflection events. For instance, shallower intervals around and above R1 exhibit more amplitude distortions than areas at greater depths, likely due to increased porous space and water content, as discussed subsequently in the paper. Multiples (M) and coherently stacked noise (N), both fundamentally introduced by the overburden-bedrock lithologic contact, are identified on the migrated images, predominantly affecting lines 2 and 3.

The geometric description of the buried valley, such as its width, thickness, and extent, as well as the intralayer characterization of the overburden material, is discussed in subsequent sections after the depth conversion of the time-domain reflectivity images shown in Figure 5.

Thus far, we have conducted conventional CMP seismic imaging focusing on the reflected arrivals within the data set. This enables us to visualize, in the time domain, the structure of the shallow earth in terms of seismic reflectivity contrasts. However, it provides limited information of the physical properties of the subsurface. Therefore, we incorporate the analysis of another type of seismic arrivals present in the recorded waveforms. We use Love waves to retrieve quantitative pseudo-2D velocity profiles along the three land-streamer transects. These velocity profiles not only facilitate quantitative interpretation within the first few meters of the subsurface but also contribute to mitigating the uncertainty associated with time-to-depth conversion of reflectivity images.

We use the dispersive nature of surface-wave propagation to understand velocity-frequency dependency (Miller et al., 1999; Foti et al., 2018; Yang et al., 2024). Once their dispersion properties are established, an inverse problem is formulated to retrieve velocity information. The inversion of S-wave velocity entails iteratively realizing theoretical dispersion data, generated by numerical forward modeling. The iterative process continues until theoretical dispersion data aligns with experimental dispersion data, updating the underlying parameter model (Wathelet, 2005).

We assess the geometric dispersion properties of Love waves by transforming each shot gather from their space-time domain into their phase velocity-frequency domain via frequency-wavenumber transforms (Yilmaz, 2001; Wathelet, 2005; Naskar and Kumar, 2022). A fundamental step prior to the analysis of dispersion data is the conditioning of the input shot gathers. Seismic data conditioning for the analysis of surface waves commonly involves the removal of near-offset traces to meet the plane-wave assumption, and temporal frequency filtering to mitigate the presence of other arrivals, such as body waves. Nonetheless, we find that 30% of our data set, which is approximately 600 shots acquired across marshlands, wetlands, and highly unconsolidated materials, requires the implementation of additional processing steps to render optimum dispersion data. This is because isolating the arrivals of interest while removing undesired events is challenging in these shot gathers. This challenge arises because the distinct characteristics of Love waves — low velocity, low frequency, and high amplitude — overlap with those of other seismic phases such as low-frequency reflected arrivals, guided waves confined within the overburden, and late-arriving dispersive phases.

To mitigate these challenges, we suggest a simple but robust five-step seismic conditioning sequence: (1) eliminating spurious traces, (2) removing near-offset traces, (3) applying frequency filtering, (4) implementing surgical velocity-based muting, and (5) eliminating residual high-wavenumber components. Identifying and removing nonsource-generated, high-amplitude noisy traces is the most fundamental step within any conditioning sequence because it removes portions of the data that notably impact the dispersion image. Likewise, near-field effects, generated by the interference of body waves in the vicinity of the source, can introduce biases resulting in erroneous estimation of phase velocity. We minimize these effects by processing most shot gathers with a nearest source-receiver offset of at least one wavelength. Although some authors propose a half-wavelength criterion (Xia et al., 2000; Park et al., 2002), we restrict near-source offsets to the range between 12 and 20 m, which exceeds the half-wavelength for our data set and complies with the plane wave assumption. Upon analyzing the spectral characteristics of shot gathers, we determine that the frequency bandwidth of Love wave phase velocities lies within the 8–32 Hz range. As a result, we implement band-pass frequency filters with corners 5–7–35–40 Hz. This filtering helps diminish the contribution of seismic arrivals such as refracted and reflected events, as well as high-frequency noise. Nonetheless, some events such as late-arriving air waves and low-frequency refracted phases persist. To address this, we use velocity-based muting defined by the slopes at the upper and lower boundaries of the surface-wave package. We find that apparent velocity ranges between 50 and 600 m/s are optimal for mitigating seismic coda and earlier refracted arrivals. Lastly, remaining high-wavenumber components that are not identifiable in the space-time (T-X) domain are eliminated in the frequency-wavenumber (f-k) domain. In this step, unlike in reflection processing where we identify and remove the specific slopes of surface-wave energy, our goal is simply to preserve coherent information from Love waves. We achieve this by tapering off high-wavenumber information across all frequency bands.

Figure 6 shows the implementation of the optimal conditioning sequence for a shotpoint located in the westernmost segment of line 3. Figure 6a and 6b shows the shot gather in the T-X domain (left) and in the phase velocity-frequency domain (right), presented on logarithmic scales. Figure 6a shows raw data characterized by interference from reflected waves (red arrows), refracted phases (green arrow), and noisy traces, hindering the identification of the fundamental mode of Love waves on the dispersion image. In Figure 6b, the same shot gather is presented after trace editing and frequency filtering. The improvement of the dispersion image is evident; however, low-frequency components appear distorted. As this low-frequency information lies below the natural frequency of the receivers, they display higher sensitivity to other low-frequency arrivals, such as remnant reflected phases (blue arrow). To mitigate semblance distortions and enhance energy focusing, we introduce velocity-based muting, using apparent velocities ranging from 80 to 600 m/s, along with f-k filtering.

This is demonstrated in Figure 6c, where the dashed black lines on the T-X graph represent the velocity-designed pass cone. The corresponding dispersion image exhibits a more coherent semblance content. In the low-frequency range, the maximum semblance of the fundamental mode shows a smooth continuous variation of phase velocity with frequency, effectively improving the interpretability of the dispersion image. The high phase velocity observed in Figure 6c, approximately 530 m/s, correlates with the shallow bedrock found along the west-southwest section of line 3. These conditioning steps are important to correctly interpret the fundamental mode of Love waves in some of our shot gathers. As shown in Figure 6a and 6b, without these steps, the interpretation of the maximum semblance would have been unreliable due to contamination from unwanted phases.

After data conditioning, Love-wave dispersion images are generated for each shot gather. Dispersion curves are extracted by tracing the maximum semblance associated with the fundamental mode. Although we acknowledge the presence of higher modes in some of our shot gathers, we opt not to include them in our workflow due to their ambiguous interpretation.

To invert our experimental dispersion data, the underlying model parameters should first be established. In our model parameterization, we define the P-wave velocity and density using a configuration of two layers plus a half-space. The S-wave velocity is defined within a configuration of nine layers plus a half-space, and the Poisson’s ratio is defined within a single layer. The locations of the land-streamer transects coincide with a 3D seismic reflection survey in our study area. This overlap allows us to incorporate the refraction velocities obtained from tomographic inversion to parameterize P-wave velocity and density variations, as well as their layer thicknesses within our $VS$ inversion workflow. P-wave velocities range between 500 and 1000 m/s in the first layer, 1000 and 1500 m/s in the second, and 2000 and 3000 m/s within the half-space. Density is estimated following Gardner’s relationship, $ρ(g/cc)=1.74Vp0.25$ (Gardner et al., 1974). We validate that the calculated densities correlate with literature values corresponding to the sediment types present within the overburden, primarily silt, sand, and gravel (Telford et al., 1976; Brocher, 2008). Poisson’s ratio is set to fluctuate between 0.2 and 0.5. To mitigate the generation of unrealistic velocity models in the S-wave velocity parameterization, while ensuring that the inverse solution is not excessively constrained, we establish a shear velocity variation range of 50 to 800 m/s. The determination of layer thickness is governed by a minimum and maximum quarter-wavelength criterion. These wavelengths are characterized by the corresponding phase velocity-frequency pairs within dispersion images. For all models in our study, the layer thickness varied between approximately 1 and 10 m. This range allows the algorithm to probe a wide variety of potential layer thicknesses within set limits.

The inversion uses a global search optimization approach based on the neighborhood algorithm and a Love-wave forward solution, as mathematically described in the work of Wathelet (2005). This method mitigates the risk of convergence to local minima, which can occur with local-based methods, particularly when the initial model is not within the basin of attraction of the global solution. Given the 1D nature of the inverse problem, the size of the data set, and the available computational resources, the global optimization approach allows for a more reliable and comprehensive exploration of the model parameter space. Few conditioning parameters are set during the inversion stage: 200 initial randomly generated models, 50,000 models searched per iteration, and 100 models considered for resampling. To reduce solution nonuniqueness, at least three trials are conducted per run, ensuring minimal variability among trials, and permitting an average root-mean-square (rms) error not exceeding 1%. This error is calculated using a misfit function designed to measure the distance between a theoretical dispersion curve and an observed dispersion curve. Figure 7a–7c shows examples of the fit between experimental and theoretical dispersion curves, and Figure 7d and 7e shows their corresponding normalized rms residuals for the last inversion. The experimental dispersion data shown in Figure 7a–7c correspond to shot gathers in the central region of line 1, the west-southwest segment of line 2, and the east-northeast segment of line 3, respectively. Higher frequency content is observed in Figure 7b compared with Figure 7c, in agreement with local geology. The frequency spectrum experiences more attenuation in the east-northeast segment of line 3 due to the greater dissipation of seismic energy within a section of thicker overburden. Most inversion runs converge to their global solution within 20,000 model updates, as shown by representative examples in Figure 7d and in the magnified view window in Figure 7e. In these instances, the misfit decreases to within 1% after approximately 5000 model updates. We opt to parameterize the inversion with 50,000 model updates, as other profiles exhibit varying convergence trends.

The inverted model retrieves velocity variations to an approximate depth of 12 m, enabling the description of shallow layers within the overburden, such as the strong reflection observed above the bedrock top in Figure 5 (R1). After the inversion phase, we assign geometry information to the inverted profiles, constructing a pseudo-2D velocity model for each land-streamer transect. Data are gridded to a depth sampling of 0.1 m. To mitigate the vertical fabric that is intrinsic to the 1D nature of the method, we apply a 2D Gaussian function to the pseudo-2D profiles. Figure 8 shows the final velocity profiles in the time domain, overlain on the seismic reflectivity images corresponding to the three land-streamer seismic transects. The velocity profiles are shown in the time domain to facilitate a comparative analysis between the quantitative imaging derived from Love waves and the reflectivity imaging derived from reflected waves. Accordingly, we stretch the final $VS$ models from depth to time using $t(x,z)=∫z0z(dz′/v(x,z′))≈∑in(Δzi/v(x,zi))$⁠, where $t(x,z)$ is the traveltime, $v(x,z)$ is the final velocity profile, $z0$ is the starting depth from datum, and $Δzi$ in the summation is the depth difference between consecutive samples. The numerical approximation of the integral is achieved by the second portion of the equation, discretization into a finite sum.

The black lines shown in Figure 8a–8c correspond to the shallow reflection, R1, shown in Figure 5. The gray-shaded areas indicate regions where equivalent depths exceed 12 m and only reflection data are available. S-wave velocities range approximately from 100 to 650 m/s. In most profiles, however, we observe that vertical velocity changes are predominantly limited to the first few meters of the subsurface (6 m or less), becoming nearly constant or exhibiting a steeply dipping depth-velocity gradient at greater depths. Lateral velocity variations appear to be predominantly influenced by the stratification of the overburden, which in turn is largely controlled by the bedrock geometry. This velocity relationship is demonstrated by the remarkable correlation between the shallow reflection geometry interpreted on reflection data (dashed black line) and the velocity contrast derived from Love-wave inversion. For instance, a notable topographic low is observed toward the east-northeast portion of line 1 in Figure 8a, as indicated by the geometry of the shallow reflection. This topographic depression is also identified by the velocity model, which displays lower $VS$ values compared with the surrounding areas, quantitatively exhibiting a similar basin-like behavior that correlates with a thicker overburden column. The inverted velocity profile not only captures large-scale variations but also detects subtle patterns. This is evident from velocity variations associated with local topographic changes, where the velocity behavior corresponds with observed local morphological changes in reflected-wave data. For instance, the topographic depression in Figure 8a (T1) shows a local thickening of low-velocity near-surface overburden materials, with a relative velocity difference of approximately 22%. This aligns with the shallow reflection geometry observed around the west-southwest segment in the migrated image of line 1, locally referred to as T1 (Figure 5a). Similarly, the topographic high observed in Figure 8c, indicated by T2, exhibits a velocity contrast that mimics the geometry of the shallow reflection (see T2 in Figure 5c). The apparent correlation between the reflectivity patterns in the migrated images and the velocity contrasts in the $VS$ models from Love-wave inversion confirms the effectiveness of these techniques in providing reliable quantitative descriptions of the overburden’s near-surface configuration.

Depth conversion of seismic data is essential for accurately mapping facies distribution and bedrock geometry. Typically, shallow boreholes are used to determine velocity functions through downhole surveys and/or by generating bedrock velocity maps from lithologic markers. However, the accuracy of this approach significantly depends on the number and spatial distribution of the available boreholes. The water wells in the study area indicate a median bedrock depth of approximately 22 m, as exemplified in Figure 2. This suggest that velocity predictions over a short depth are adequate to account for most of the overburden-bedrock contact. We propose integrating the velocity data obtained from Love-wave inversion with water well information to build bedrock-reaching velocity profiles.

Due to our inverted velocity profiles extending only to a depth of approximately 12 m, we use regression analysis to estimate $VS$ at greater depths. We investigate the degree of linearity of the inverted $VS$ profiles to assess the applicability of linear regression. Selected 1D profiles are shown for line 1 (Figure 9a and 9b), line 2 (Figure 9c and 9d), and line 3 (Figure 9e and 9f). We fit a linear regression model to each inverted $VS$ profile, creating a line of best fit (dashed red line) through the original data points (solid blue line). The first 3 m are ignored due to their low-velocity, highly attenuating near-surface sediments, which lead to increased velocity nonlinearity. The gray shading represents velocity errors of ±10%, reflecting the typical range of velocity uncertainty (10%–15%) related to S-wave velocity inversion of surface waves (Ivanov et al., 2006; Boaga et al., 2011; Foti et al., 2018). Figure 9a and 9c shows representative profiles with a good linear fit, whereas 9d–9f shows representative profiles with a poor linear fit. We find that most profiles exhibit a linear depth-velocity relationship, as shown in Figure 9a, 9c, and 9e, where we identify a consistent steeply dipping depth-velocity trend in the inverted S-wave profiles. Even for profiles with reduced linearity, such as those in Figure 9b, 9d, and 9f, applying a 2D Gaussian function mitigates these deviations in the context of velocity modeling for time-to-depth conversion. This adjustment brings the predictions within the method’s inherent uncertainty, allowing the assumption of linearity to hold more than short depth intervals. The validity of this assumption relies on the absence of significant compositional and lithologic differences that could cause pronounced depth-velocity nonlinearities. Our analysis of the water well data reveals no evidence of such scenarios.

The velocity prediction is further constrained by a velocity map along the overburden-bedrock contact, derived from water well data. This map is generated by creating a two-way traveltime map from the interpreted bedrock on reflectivity sections and a depth-to-bedrock map based on water well lithologic markers. Although this serves as a proxy of velocity variations only along the overburden-bedrock contact, it is incorporated into the workflow by setting upper velocity bounds for $VS$ predictions, preventing unrealistic extrapolation at greater depths. This is achieved by setting a rate of change based on the Love-wave inverted $VS$ samples and using it to adjust predicted values exceeding the upper bounds defined by the bedrock velocity map. We use the integrated velocity model for time-to-depth conversion of the reflectivity images, as shown in Figure 10. Figure 10a, 10c, and 10e shows the final time-domain bedrock-reaching velocity models for lines 1, 2, and 3, respectively, superimposed on the time-domain reflectivity images. Because the assumption of linearity does not apply beneath the bedrock, a gray-shaded region has been incorporated into the images, following the bedrock’s morphology. Figure 10b, 10d, and 10f shows show the final depth-converted sections for lines 1, 2, and 3, respectively. We find that the overburden reaches its maximum thickness between lines 1 and 2, with an apparent maximum thickness of approximately 60 and 50 m, respectively. This correlates with water well data; for instance, borehole 1405385 near line 2 registers the bedrock depth at approximately 51 m. There is no clear indication of the valley’s presence around line 3. This absence may be due to the valley becoming shallower or changing orientation, potentially aligning with the current orientation of the Teeswater River (alluvial deposit in Figure 1). No interpretation has been added to Figure 10 to aid readers in visualizing the reflectivity contrasts.

Because borehole velocity measurements are unavailable, the velocity model used for depth conversion cannot be directly validated. Instead, the accuracy of the time-to-depth conversion results is assessed to indirectly validate the underlying velocity model. Boreholes 1408389, 1404645, 1405385, and 7245533 serve as “blind” boreholes and are excluded from all analyses conducted in our study. The appraisal approach involves comparing lithologic markers associated with two primary seismic reflections: the bedrock reflection (R2) and the shallow reflection (R1). Table 3 provides the comparative exercise organized by borehole-seismic distances, ranging from closest to farthest. We compute residuals by subtracting the actual depth of lithologic markers from their corresponding depths on the depth-converted reflection images. The residuals exhibit remarkably low errors, with an average value of less than 2%. Borehole 7245533, being shallower, does not penetrate the overburden-bedrock contact. By comparing prominent reflections on the depth-converted reflectivity sections with ground truth data, we validate the accuracy and reliability of the underlying velocity model.

The overburden deposits, as shown in Figure 2, predominantly comprise sand and gravel, with some clay and silt layers. Due to the poor sensitivity of seismic reflectivity to subtle compositional and lithologic differences between various material types, we use the depth-converted sections to compute volume-amplitude technique (TECVA) and dissimilarity attributes (Chopra and Marfurt, 2005; Kumar and Sain, 2018; Ercoli et al., 2020), as shown in Figure 11. By highlighting variations in texture and continuity of the seismic response, these attributes can make small-scale heterogeneities more apparent. Figure 11 shows the interpretation of the overburden-bedrock contact, depicting the geometry of the bedrock valley, and the interpretation of the shallow reflection within the overburden. Surficial geology is shown above the profiles, and boreholes within a 250 m radius are displayed on the interpreted TECVA images. Figure 11a shows the TECVA attribute, Figure 11b shows the dissimilarity attribute, and Figure 11c shows a masked version of the TECVA attribute incorporating borehole information.

Lines 1 and 2 present wells that penetrate a thicker overburden associated with the bedrock valley, terminating in the Paleozoic bedrock. Overburden logs indicate that the fill of the valley consists primarily of gravel near the base and sand throughout most of the channel with clay and silts interbeds. This is consistent with results from overburden logs SB_MW01 and SB_MW06 (Figure 2), which show significant thickness of overburden material located within the broader part of the valley.

The interval between the ground surface and the shallow reflection is designated as unit A, while the underlying interval, extending down to the overburden-bedrock contact, and is designated as unit B. In all lines, particularly in lines 1 and 2, distinct reflection patterns are evident below the prominent shallow reflection (R1) within the thickest section of the overburden. These patterns result from stratigraphic heterogeneities, as indicated by boreholes SB-MW01, SB-MW05, and SB-MW06 at depths greater than 15 m (Figure 2). However, we refrain from interpreting the geologic significance of these reflection patterns because they are somewhat smeared and locally clustered, making it challenging to correlate them with the sparse borehole coverage near the seismic transects.

The shallow reflection, R1, is recognizable in all three seismic profiles and in the attribute images. This prominent reflection serves as the boundary separating units A and B, which is suspected to be a transition zone involving clay to silty clay and sand to silty sand materials. The transition between these materials and their corresponding seismic responses is not always distinct. For instance, boreholes in line 1 show a transition from clay to sand, whereas borehole 1405385 in line 2 shows a transition from sand to clay (Figure 11c). This variability is attributed to the presence of a mixture of materials with different properties, such as varying material type, permeability, and porosity, which increases lateral heterogeneity. In fact, continuous areas of high dissimilarity in line 1 correlate well with locations of mapped organics, clay, and peat, which also correspond to wetlands. However, this observation is not always consistent due to lateral heterogeneity. Elsewhere, high dissimilarity in unit A is associated with sands and gravels of alluvial deposits, such as in line 2, or with glaciofluvial outwash and ice-contact drift deposits in lines 2 and 3. Along the east-northeast segment of line 2, a thin high dissimilarity response is observed along the base of unit A. Areas of high dissimilarity appear to be associated with near-surface low-velocity zones, mainly within unit A. For instance, the eastern section of line 1 exhibits lower S-wave velocities (as shown in Figure 8a), which corresponds to the high dissimilarity interval noted in Figure 11b within unit A. We believe that these differences in dissimilarity may arise from either water-saturated intervals or from a mixture of heterogeneous materials, such as poorly sorted sands and gravels, leading to increased signal attenuation and relatively low velocities.

The bedrock horizon, overlain on the TECVA image, follows a clear, coherent, and high-amplitude reflection. Locally, within the deeper parts of the bedrock valley, the signal is lost but only over short distances. This allows a full geometric description of the buried valley along the seismic profiles. Seismic lines 1 and 2 in Figure 11 present the vertical and lateral extents of the bedrock valley. The valley generally shows an asymmetric geometry in line 1, gradually becoming shallower toward the east-northeast. Based on this geometry, the east-northeast side of the valley may represent floodplain-type deposits. The width of the valley is approximately 850 m as observed in line 1, with approximately 200 m accounting for the possible floodplain. Line 2 shows the valley with a more symmetric geometry, characterized by shallower dipping sides and an approximate width of 600 m. There is no geometric evidence in line 3 of a deeply incised bedrock valley, as observed in lines 1 and 2. This likely indicates that the valley becomes shallower or changes orientation, as previously stated.

The bedrock interpretation seems not to align with the bedrock marker in borehole 7330479 (Figure 11c, line 3). This discrepancy arises from the projection of the borehole onto the seismic line, which does not account for local variations in bedrock geometry. The edge of a local depression in the bedrock, resembling a channel and running in an north-northwest–south-southeast direction, is indicated by the red arrow located just west of borehole 7330479. Because this borehole is situated off the line, the reported bedrock depth corresponds to the edge of the channel-like bedrock feature rather than its bottom (green arrow). Hence, when the borehole data is projected onto the line, it appears to mismatch the seismic interpretation because the projection does not follow the strike of the channel-like structure.

In this study focused on a bedrock valley in southern Ontario, we use seismic imaging techniques to characterize its structural and stratigraphic configuration. Through common-midpoint reflection processing and Love-wave inversion using land-streamer seismic data, we successfully achieve high-resolution imaging of the overburden-bedrock contact and the horizontal sedimentary layers above it. Although multicomponent reflection processing is attempted, we find that satisfactory imaging results are primarily obtained from the transverse component of the 3C seismic data sets. We demonstrate the necessity for, and suggest, a more comprehensive signal preconditioning sequence in the generation of Love-wave dispersion images from noisy gathers. Our imaging approach not only delineates the geometry and extent of the valley but also provides insights into the distribution of quaternary sediments sitting above the limestone bedrock. Through the integration of reflection data for structural and stratigraphic interpretation and Love-wave inversion for velocity imaging, we identify and characterize the spatial distribution of a prominent lithostratigraphic transition zone within the overburden. Likewise, using reflection and surface-wave data, we establish a methodology for depth conversion, crucial for accurately positioning geologic features. Our study contributes to refining the understanding of near-surface geologic conditions at the South Bruce site, providing a framework applicable to analogous environments where detailed characterization of the shallow subsurface is necessary for engineering and numerical modeling assessments.

The authors gratefully acknowledge A. Pugin, Research Scientist at Natural Resources Canada, for his invaluable insights during the revision of this work. We also thank A. Parmenter, Director of the Geoscience Division at NWMO, for his support of this project and for granting permission for publication. Finally, we appreciate the constructive feedback from the editors and the three anonymous reviewers, which greatly contributed to the improvement of this paper.

Data associated with this research are confidential and cannot be released.

Brian Villamizar is a geophysicist with more than 10 years of experience in subsurface imaging and characterization in sedimentary and crystalline environments. He received a B.E. from Simón Bolívar University (Venezuela), an M.Sc. from China University of Petroleum (China), and a Ph.D. from Western University (Canada), all in geophysics. He began his career at PDVSA in seismic processing, interpretation, and structural modeling, and subsequently served as technical lead at Novamera Inc. Currently, he serves as senior geophysicist lead at the Nuclear Waste Management Organization, where he oversees the acquisition, processing, and integration of borehole and seismic data for site characterization. He is a member of the Association of Professional Geoscientists of Ontario (APGO) and SEG.

Aaron DesRoches received a Ph.D. (2017) in geophysics and a B.Sc. (2004) in geology from the University of New Brunswick, as well as an M.Sc. (2008) in earth science from Simon Fraser University. He joined the Nuclear Waste Management Organization in 2011 and he has held roles as senior geoscientist, section manager, and currently as manager of geoscientific systems integration. His research focuses on integrating multidisciplinary geoscientific data — such as geology, geophysics, hydrogeology, geochemistry, and geomechanics — to develop comprehensive 3D geologic models. He is an active member of professional organizations, such as APGO.