P-wave velocities presented in this manuscript were acquired by two methods: shallow seismic refraction (SSR) and sonic logging (SL). Both methods rely on the refraction of elastic energy at elastic boundaries. The elastic properties are dependent on the density and bulk and shear moduli of a given material. Refraction is advantageous because the velocity of a material can be determined simply by the first-arrival time of the seismic energy—the phase and amplitude of the arrival can be ignored.

During a SSR survey, the travel times between a seismic source to an array of geophones are measured and inverted to generate a 2D tomogram of seismic velocity (Woodward, 1992; Sheehan et al., 2005; Gance et al., 2012; Zelt et al., 2013). This type of SSR surveying is a robust method that has been used on crustal scales (Mooney and Brocher, 1987; Korenaga et al., 2000; Hayes et al., 2013), hillslope scales (Whiteley and Eccleston, 2006; McClymont et al., 2012; Thayer et al., 2018; Leone et al., 2020), and over outcrops (Holbrook et al., 2014; Flinchum et al., 2018b). SSR sources utilize a low-frequency (10s of Hz) source that produces relatively large wavelengths (10s−100s of meters). In an SSR survey seismic energy can travel 10s−100s of meters both vertically and laterally before being recorded at a geophone.

SL is less common in the CZ because it requires an uncased, large-diameter (2” or greater) borehole. In addition, sonic logs only work beneath the water table, since water couples the sonic source and receivers to the surrounding rock. Due to these constraints, only a limited number of studies have used sonic velocities to study the top 100 m of Earth's surface (Watson et al., 2005; Martí et al., 2006; Holbrook et al., 2014; Raef et al., 2015; Flinchum et al., 2018a; Gu et al., 2020a). The acquisition of a sonic velocity via SL works on the same principles as the SSR survey; that is, velocity is obtained by measuring the time taken to travel between a source and multiple receivers. In the SL case, however, most of the energy is traveling vertically, at right angles to the majority of travel paths in SSR.

There are a few important differences between the SSR and SL. First, during a SL survey, the maximum offset between the source and receivers will be limited by the tool length, often to 10′s of centimeters to a few meters, in contrast to the 10′s−100's of meters that might separate a source and geophone during SSR. The tool is moved up and down the hole in centimeter increments, providing estimates of velocity every few centimeters. Second, the sonic tool uses a source frequency that is 2–3 orders of magnitude higher than those used in a seismic refraction survey (our study is comparing ~80–15,000 Hz). The higher frequencies produce wavelengths 2–3 orders smaller than refraction. Finally, given that water couples the source and receiver to the borehole wall, the first arrival from energy will be critically refracted along the water/rock interface making it sensitive to a thin region along the borehole wall.

We reprocessed two SSR refraction lines that were part of a larger study presented by Flinchum et al. (2018a,b). L29, which extends from the top of the gently undulating ridge into the valley to the east, intersects boreholes BW-1 and BW-4 and comes within 10 m of holes BW-6 and BW-7 (Figure 3). A sledgehammer source was used by striking a stainless steel plate; hammer swings were stacked at least 8 times to improve the signal-to-noise. L29 was collected with 240 channels at 1 m spacing with shots every 12 m. L25 was collected with 96 channels at 2.5 m spacing with shots every 10 m. On L25, three off-end shots extending 30 m before the first geophone and 30 m past the last geophone were also collected. There were no off-end shots collected on L29.

For the SSR surveys, arrival times were picked manually. To pick the arrivals, shot gathers were trace normalized. We made as many picks as possible on the trace-normalized data without aid of any filtering. After picking on the unfiltered data, we applied a band-pass filter between 10 and 150 Hz to help identify far-offset arrivals. After picking, L29 contained 3935 first arrival picks and L25 contained 2039 picks (Supplementary Figure S1). In cases where the first arrival was noisy or could not be determined, no pick was made, to avoid adding errors into the final models.

To invert the first-arrival times and produce a 2D image of seismic velocities we use the Python Geophysical Inversion Modeling Library (PyGIMLi) (Rücker et al., 2017). PyGIMLi utilizes a shortest path algorithm (Dijkstra, 1959; Moser, 1991; Moser et al., 1992), which models seismic energy propagation as rays, not as waves. The shortest path algorithm is used in many other travel-time tomography codes (e.g., Zelt et al., 2013). The shortest path algorithm calculates a path through the velocity model—referred to as a ray path. The ray paths are used to mask areas in the final velocity that are not constrained by the data. In other words, if a ray path does not pass through a given cell in the model it will not be plotted. Ray paths tend to be curved in appearance because they refract toward the surface when velocity increases with depth. As a result, the resolved depths of the models vary spatially as driven by the data.

PyGIMLi utilizes a deterministic Gauss-Newton inversion scheme (Gance et al., 2012) which requires a data weight for each travel-time pick. These data weights are used to drive the error weight matrix and ultimately influence the final model fits. Here, instead of assigning a constant error value to each pick we used a linear equation as a function of offset (distance from source) line length, (1)error=0.003*offsetline length+0.0005 where offset is the distance [L] between the source and the geophone. Using Equation 1 means that the errors start at 0.5 ms and increase linearly to a maximum of 3 ms at the spread length. We chose a maximum error of 3 ms because of the high quality of the data. Each associated error is read into the inversion and used to weight the travel-times, such that travel times with the smallest errors are given the most weight.

To compare P-wave velocities derived from SSR with those from SL, it is important to quantify uncertainty on the velocities produced by the tomography. To quantify uncertainty, we randomly resampled the travel-time data, to utilize the large amount of input travel-times and randomly remove influences for each one. To achieve this, both L29 and L25 were inverted 20 times on the same mesh with the same inversion parameters, with 40% of the travel-time picks randomly removed before each inversion. Velocities at each pixel in the mesh were averaged together from all inversion realizations to produce an average 2D velocity model as well-standard deviations expressed as a percentage of velocity (Supplementary Figures S2, S3). We then use the average velocity profile to trace all the picked arrival times, so that the final fitting criteria for the tomograms shown here use all of the first arrival picks. We also compute the vertical velocity gradient (dv/dz) for each averaged tomographic model to highlight vertical changes in velocity.

For the SL acquisition, we used a Mount Sopris 2AFWS (Colorado, USA). The downhole instrument is ~ 3.0 m long with an insulator separating the source from the three receivers. The distance between the source and receivers is fixed at 0.91, 1.22, and 1.52 m. The source has a center frequency of 15 kHz. To process the downhole data, we recorded the elastic wavefield response at each of the three receivers. The elastic wavefield recording allows us to determine P-wave velocities by picking the first arrival times. As the tool travels through the borehole a source is generated at 5 cm intervals. To aid in the picking of arrivals, a semblance plot is created by summing the amplitudes of the waveforms along fifteen different slowness values (1/velocity) using a 50 μs window for every depth sampled. The best-fitting velocity from each source location (every 5 cm) was auto-picked based on the highest semblance value. After the automatic picks were made across the entire borehole depth, the picks are visually inspected, and erroneous picks are manually adjusted or removed. Because water in the borehole is required to couple the seismic source and geophones, sonic velocity logs do not exist for the unsaturated zone near the surface. Our data set consists of P-wave velocity profiles below the casing and below the water table for all five boreholes in this study (Figure 3).

The CZ structure in the weathered and fractured granite at the study site is laterally heterogeneous. A Mount Sopris QL40-OBI-2G Optical Televiewer was deployed to provide a 360° unwrapped optical borehole image (OBI) at each of the boreholes. We converted the images to black and white and boosted the contrast to emphasize fractures and heterogeneity in the rock. The OBI logs and sonic logs based on both seismic velocities from the four deep (> 50 m) boreholes, all located within 70 m of each other, show strong lateral heterogeneity in the fractured rock directly below the casing (Figure 5). Downhole hydraulic conductivity logging revealed that hydraulic conductivity varies over two orders of magnitude but generally decreases with depth (Ren et al., 2019). Other geophysical studies have shown that the influence of heterogeneous remnant fabric from fractures extends upward into the saprolite (Novitsky et al., 2018) and that the depth to the intact bedrock is also laterally heterogeneous on the 100s−1,000s of m scale, extending as deep as 60 m below the surface (Flinchum et al., 2018b; Keifer et al., 2019; Wang et al., 2019a). In addition, nuclear magnetic resonance (NMR) logging was completed at all of the boreholes in the study area and show low but variable water contents (<0.15 m³/m³) in the fractured rock region of the CZ; the water is focused in visibly fractured regions (Flinchum et al., 2019; Ren et al., 2019). The fracture density in the fractured rock layer inferred from NMR is highly heterogeneous, with porosity varying from 0 to 0.15 m³/m³.

At L29, the SSR profile is generally divided into two domains. The first domain is defined by velocities < ~1,800 m/s (warmer colors in Figure 6A). The second domain is defined by velocities > ~4,000 m/s (cooler colors in Figure 6A). These two domains are separated by a thin transitional region where velocity increases rapidly over a short vertical distance (~5–10 m) with gradients of 100–300 m/s/m (Figure 6B). The vertical velocity gradient image (Figure 6B) highlights this transition or boundary separating these two domains. This boundary is deepest (> 50 m) under the ridge (between 25 and 100 m along the profile) (Figures 6A,B) and comes within 10 m of the surface under the valley (between 210 and 240 m along the profile) (Figures 6A,B). The strong increase in velocity can also be observed in the travel-time picks and represents the most significant seismological boundary revealed by the surface refraction data (Supplementary Figure S1).

Sonic velocity logs from BW1, BW-4, BW-6, and BW-7, which are on or near seismic profile L29 (Figure 3), show significant variability above 40 meters depth (Figure 6C). This variability is surprising, given that BW-4, BW-6, and BW-7 lie within 30 m of each other and BW-1 is only 70 m from BW-4 (Figure 3). Between 20 and 25 meters depth, for example, the velocity at BW-4 is ~3,000 m/s, while at BW-6 it is ~5,000 m/s. This variability disappears below about 40 m depth, where the sonic velocities converge to a velocity of ~4,500 m/s (Figure 6C).

Extracted 1D velocity profiles from the surface refraction data are smooth and lack both the detailed vertical structure and the lateral heterogeneity that is observed in the sonic velocities (Figure 6C). Above 40 m depth, where strong variability in sonic velocities occurs, refraction velocities are slower by as much as 4,000 m/s (Figure 6C). Despite the large differences in the shallow parts of the borehole, the tomographic velocity and sonic velocities converge to a common velocity between 4,000 and 5,000 m/s (Figure 6C), close to the 4,000–4,500 m/s measured by surface refraction on granite outcrops in the Southern Sierras (Holbrook et al., 2014) and nearby in the Laramie Range (Flinchum et al., 2018b). The discrepancy between the sonic velocities and travel-time velocities above 40 m is not a processing or inversion artifact (Figure 6C). We conducted a forward model study using the sonic velocities to predict SSR travel times for L29. The forward modeling results show clearly that the sonic velocities are too fast to explain the observed travel times (see Supplementary Materials S2, S3). Furthermore, similar discrepancies have been observed in the Piedmont of South Carolina (Figure 2) (Holbrook et al., 2019) and at a minimum of three other sites worldwide (Watson et al., 2005; Martí et al., 2006; Raef et al., 2015).

Just 500 m northwest of profile L29, the sonic and refraction velocities on line L25 show a much better match at all depths with overlapping observations (Figure 7). As with L29, the SSR profile can be divided into two domains: one with velocities < ~1,800 m/s, and one with velocities > ~4,000 m/s (Figure 7A). In this second refraction profile (Figure 7A), the seismic refraction velocity shows a high velocity zone that comes to the surface and extends 100s of meters laterally (Figure 7A).

The L25 profile crossed a depression in the landscape, so the water table is at ~7 m depth. This shallow water table allowed us to acquire sonic velocities closer to the surface than in any of the other boreholes (Figure 3). At this location, the sonic velocities and the seismic refraction velocities are in much better agreement (Figure 7C). The match is not perfect, however: the SSR velocity at L25 is about ~10% higher (~4,300 compared to ~3,900 m/s) than the sonic velocities observed in the borehole (Figure 7C). This is most likely caused by the strong vertical velocity gradient near the surface, which is a difficult structure for travel-time tomography to resolve. This area of the model has our highest uncertainty estimates (Supplementary Figure S3e). The uncertainty in the strong vertical velocity gradient near the surface would have the tendency to make the velocities near the surface slightly slower and the velocities at depth slightly faster to fit the observed travel times. Thus, we believe that the 1D velocities at L25 (Figure 7C) agree within error. Additionally, the surface refraction velocity profile still lacks the velocity detail observed in the sonic profile, and both velocities converge to an average velocity of ~4,000–4,500 m/s, consistent with velocities collected over outcrop in the Southern Sierras (Holbrook et al., 2014) and the Laramie Range (Flinchum et al., 2018b).

One of the challenges of studying CZ processes is integrating a wide variety of observational data across different spatial scales (e.g., mineralogical vs. hillslope vs. landscape scales). The discrepancy between refraction and sonic velocities (Figures 2B, 6C) might appear to beg an obvious question: which velocity is “correct”? Here we argue that the answer is that both velocities are correct, at the wavelengths over which they are measured. Neither velocity is in error; rather, the velocities differ because they sense fracture densities that are scale-dependent. Before presenting our argument, we must first address whether the observed sonic/refraction velocity discrepancy is due to vertical anisotropy. Anisotropy is the ratio between seismic velocities measured in different directions. In general, seismic velocities travel slower across foliation or fractures (Novitsky et al., 2018; Eppinger et al., 2021). The seismic waves produced by SSR travel largely in the horizontal direction, whereas the seismic waves produced by SL largely travel in the vertical direction, raising the idea that pervasive, near-vertical fractures would cause slower refraction velocities than sonic velocities. However, vertical anisotropy cannot explain the velocity discrepancy we observe for several reasons. First, the anisotropy values required to explain the SSR/SL discrepancy would need to be as high as ~100%, which is far higher than observed anisotropy due to fracturing elsewhere (e.g., Crampin et al., 1980; Maultzsch et al., 2003). Second, faster velocities in the vertical direction would imply pervasive near-vertical fracturing. While fracturing exists at our site, the amount of anisotropy observed in the saprolite (Novitsky et al., 2018) is modest (maximum 12%), and we are unaware of processes that would explain 100% anisotropy in the fractured bedrock with only ~10% anisotropy in the saprolite. Third, numerous near-horizontal fractures are observed in the borehole optical logs (Figures 5A–D)—a pattern that should create faster velocities in the horizontal direction. Finally, decreases in sonic velocity between the boreholes at Blair correlate strongly with fracture density observed in the boreholes (Figures 5A–D). This suggests that sonic velocities are controlled by the number of fractures in a given volume (i.e., fracture density), not the orientation of the fractures.

In the case of our sonic and seismic waves passing through the fractured bedrock layer, the wavelengths differ by two orders of magnitude (v = fλ): a seismic P-wave with a dominant frequency of 40 Hz traveling at 2,000 m/s would have a dominant wavelength of 50 m, while a sonic P-wave at 6 kHz and 3,000 m/s has a wavelength of 50 cm (Figure 8). P-waves propagate by compressing and dilating the material at the scale of the seismic wavelength. This, in turn, means that the seismic wave senses all the fractures present in a volume ~50 m in diameter, while the sonic wave only senses fractures within a 0.5 m diameter. The volume that is deformed by this compression is related to the wavelength and is known as the representative elementary volume. The volume can change with scale depending on the heterogeneity of the material. The representative elementary volume is an important concept that can used to explain the variation of a physical property with upscaling (Matonti et al., 2015; Bailly et al., 2019).

To use a hydrogeology analogy, we postulate that comparing sonic and refraction velocities is similar to comparing hydraulic conductivity measurements of a fractured aquifer obtained with a pump (i.e., transmissivity) test with measurements made in the lab on a small sample. Apparent hydraulic conductivity is well-known to scale with the measurement, where small scale measuremetns tend to show lower hydraulic conducitvites (Neuman, 1988; Clauser, 1992; Rovey and Cherkauer, 1995; Schulze-Makuch and Cherkauer, 1998). In most cases, depending on how the sample was obtained, the lab sample will have a lower hydraulic conductivity because the pump test is an average of high and low hydraulic conductivity fractures along the entire length of the borehole (Clauser, 1992) (equivalent to the large yellow area in Figure 8). It is not uncommon in fractured rock aquifers to have one or two highly productive fractures while others do not contribute to flow (Sausse and Genter, 2005; Lee and Kim, 2015).

The analogy of hydraulic conductivity is useful but imperfect, in one respect: in the case of hydraulic conductivity, the length of the screened area or the location of the sample are known, while in the case of seismic refraction, the wavefield is complex (Supplementary Figure S5), and wavelength varies as a function of velocity and source frequency (Supplementary Figures S4, S5). The velocity is affected by the porosity distribution of the CZ, which is, by definition, site specific. (Liu et al., 1976; Xia et al., 2003; Tisato and Marelli, 2013). Fluid in pore space can cause frequency dependent attenuation in the form of Wave Induced Fluid Flow (WIFF) (Müller et al., 2010). In seismic refraction, the elastic energy can travel 100s of meters from source to receiver (Figure 8). The combination of these physical phenomena makes it difficult to explicitly define the sample region represented by the travel-time observation.

We hypothesize that velocities measured across such vastly different scales in a porous medium would only be equal if the fraction of pore space was equal across those scales—essentially, if the distribution of pore space is uniform. Refraction and sonic velocities would match if we collected a seismic refraction profile over a borehole with sonic velocities in a location where the underlying rock has a homogenous porosity distribution—that is, if porosity is equal regardless of the sample volume. While such a fairly uniform distribution of porosity can occur in a sedimentary rock (e.g., a clean sandstone or a gravel aquifer), it is unlikely to be true when porosity is due to fractures. Fractured materials can have highly non-uniform fracture distributions across a wide variety of scales. We illustrate this point with a simple example: a layer of bedrock pervasively broken by evenly spaced fractures (Figure 8).

If fractures are 1 m apart and each have an aperture of 4 cm, then within the 50-m wavelength of a seismic P-wave, there will be 50 × 0.04 = 2 m of open fracture space, or a fracture porosity of 4% (2/50). Given the strong sensitivity of P-velocities to fracture porosity, adding 4% porosity to the bedrock could easily reduce seismic velocities from, say, 4,500–2,000 m/s (Figure 4). The seismic wave travels at this relatively slow velocity throughout the fractured bedrock layer, and its travel times in the refraction data reflect that slow propagation. The sonic tool, in contrast, with its wavelength of around 50 cm, would measure non-weathered and unfractured rock (velocity 4,000–4,500 m/s) except when crossing an individual fracture (Figure 8), when the velocity would drop precipitously—exactly as observed in the sonic logs of Figures 6B, 7B.

For this reason, we hypothesize that the persistent discrepancy between seismic and sonic velocities in the fractured bedrock layer of the CZ in South Carolina (Holbrook et al., 2019) (Figure 2) and the Laramie Range (Figure 6) should be viewed not as an error, but rather as valuable information. Seismic velocities that consistently underestimate sonic velocities are a sign of pervasive fracturing at scales between the respective seismic and sonic wavelengths. Where the two velocities agree, in contrast, then porosity is likely uniformly distributed. This is the case in Figures 2, 6B at depths >40 m, where the relatively high velocities suggest that the meter-scale fractures are largely closed due to ambient pressure (e.g., St. Clair et al., 2015).

Under our proposed hypothesis, velocity models produced by seismic refraction tomography will be blind to small-scale (10s of meters) heterogeneity but will still be able to quantify the thickness of the CZ across landscapes over 100s of m. Seismic refraction still remains a particularly powerful geophysical tool to explore the deeper regions of the CZ, especially if the CZ boundaries of interest are controlled by changes in porosity or fracture density. The velocity profiles derived from seismic refraction represent large-scale averages that reflect the dominant seismic wavelength. It remains difficult to quantify the volume that this average represents, since quantifying the volumes would require a fully 3D seismic waveform solution. Full waveform inversions, even in 2D, are still in their infancy when working in the strong velocity gradients associated with near-surface materials (Gao et al., 2006, 2007; Virieux and Operto, 2009; Romdhane et al., 2011).

Despite the lack of sensitivity to small-scale structures, seismic refraction remains one of the only geophysical tools that can quantify the general location of boundaries defined by changes in porosity in the CZ over landscape to watershed scales. The seismic refraction method suffers from the classic tradeoff between velocity detail and spatial coverage, and users need to be aware of the relevant wavelength scales to interpret the images correctly. Although high-frequency sonic velocity logs in the CZ remain rare, our data show that P-wave seismic velocity is scale-dependent. Given the interdisciplinary nature of CZ science, it is likely that the 2D profiles of velocity will be used to inform or test conceptual models of CZ evolution (West et al., 2019; Leone et al., 2020). It is important to realize that the velocities in seismic refraction tomography represent large-scale averages over the scale of their dominant wavelength. While smaller-scale heterogeneity certainly exists beneath the lateral and vertical resolution of the refraction inversion, the velocity discrepancy itself is not a consequence of the inversion, but rather of the physics of wave propagation at different wavelengths. Understanding this large-scale averaging is particularly important at sites where velocities are used quantitatively to estimate porosity (Mota and Santos, 2010; Flinchum et al., 2018a; Callahan et al., 2020; Gu et al., 2020b) or strain (Hayes et al., 2019).

Our results show that the refraction and sonic velocities converge at depths where the fracture density is uniform, consistent with the idea of a protolith that is defined by a small, uniform porosity. The sample volume is not straightforward, as wave propagation physics is complicated, but a good rule of thumb is that the wavelengths, which depend on the velocity and frequency content of the source, provide a rough estimate of the averaging scale that the velocities represent.