Imaging an Enigmatic Paleovalley with Passive Seismic Methods (Unaweep Canyon, Colorado, United States)

Ambient noise seismic interferometry (SI) is a well‐established technique for shear‐wave velocity imaging at continental scale (e.g., Shapiro and Campillo, 2004; Lin et al., 2008), basin scale (e.g., Cheng et al., 2021; Jia and Clayton, 2021), and subkilometer scale (e.g., Hannemann et al., 2014; Sergeant et al., 2020). To directly image seismic impedance contrasts (reflectors) from ambient noise, autocorrelation (AC) has been applied for illuminating deep crustal structures (e.g., Tibuleac and von Seggern, 2012; Kennett, 2015). However, due to the lack of well‐defined and steeply incident body waves, only a few studies have used ambient noise AC for shallow targets (e.g., Saygin et al., 2017). Other studies have applied preprocessing of ambient noise to select body waves from local seismicity (Dangwal and Behm, 2021) and local microearthquakes (Polychronopoulou et al., 2020) to improve shallow reflectivity retrieval from AC. An alternate approach uses teleseismic signals with near‐vertical incidence to image deep crustal (Abe et al., 2007) and shallow structures (Ph?m and Tkalčić, 2018).

In this study, we apply passive seismic methods to image the sediment fill and the basement geometry along a 5 km long longitudinal profile of an upland valley. We take advantage of steeply incident teleseismic waveforms to extract shallow (<500 m depth) basement reflections using AC, resulting in zero‐offset P‐wave reflectivity section along the profile. Furthermore, we use ambient noise SI to retrieve surface‐wave propagation in the same depth range. We then invert the frequency–velocity dispersion trends of these surface waves to obtain a 2D shear‐wave velocity model. These results are interpreted in context with the existing well data and the intersecting active seismic section.

Unaweep Canyon is a ∼70 km long northeast–southwest‐trending gorge that bisects the Uncompahgre plateau of Colorado. It is an enigmatic landform, because two creeks too small to have carved the canyon emanate from a nearly imperceptible divide (at 2142 m elevation) within the middle of the canyon–East Creek flowing northeast (toward Whitewater), where it joins the Gunnison River, and West Creek flowing southwest (toward Gateway) where it joins the Dolores River (Fig. 1a). The canyon cuts through the Mesozoic strata and into Precambrian crystalline basement, and hosts sedimentary fill locally at least 330 m thick (Massey well; Soreghan et al., 2007) of Quaternary and possible pre‐Quaternary age. The competing hypotheses for the formation of Unaweep Canyon are: (1) Cenozoic fluvial incision, (2) Cenozoic glacial incision, and (3) late Paleozoic glacial incision and burial followed by Cenozoic fluvial exhumation. The widely accepted fluvial hypothesis proposes incision by the ancestral Gunnison River and/or ancestral Colorado River (e.g., Gannett, 1882; Lohman, 1961; Cater, 1966), but many geomorphological attributes of the canyon are difficult to explain by an entirely fluvial genesis. For example, the U‐shaped, amphitheater‐like side valleys, hanging valleys, and apparent truncated spurs have been cited as evidence for the recent (Quaternary) glacial origin (Cole and Young, 1983). However, the absence of Quaternary glacial deposits in the canyon and low elevations confound this claim. The third hypothesis posits formation by Paleozoic glaciation followed by later burial and then partial exhumation by the ancestral Gunnison River (Soreghan et al., 2007, 2008, 2014, 2015). Yet, the Paleozoic glacial hypothesis remains controversial owing in part to the equatorial setting of the Uncompahgre uplift (a paleohighland) during the late Paleozoic ice age.

Given the importance of determining the shape of the canyon for understanding its origin, and thus paleoclimatic implications, several studies have been attempted to assess the basement geometry using geophysical methods. Soreghan et al. (2008) demonstrated that a transverse U‐shaped basement structure best matched the gravity observations relative to a V‐shaped basement. In addition, prestack time migration (PSTM) imaging along a 2.5 km long transverse seismic reflection profile has revealed a U‐shape buried valley floor (Patterson et al., 2021; Fig. S1, available in the supplemental material to this article). More recently, a drilling expedition in February 2022 (well UDR‐1A; see Data and Resources) struck basement at 368 m depth (elevation 1553 m).

A 5 km long section of highway CO‐141 traversing Unaweep Canyon was chosen to deploy the passive seismic profile (Fig. 1). This section lies in the western canyon, ∼15 km east of Gateway (near the western mouth of the canyon). The section was chosen because (1) it intersects the previously mentioned active seismic profile acquired in 2017, (2) three drilled and partly cored wells are located nearby, and (3) a gravel quarry near the western end of the section as well as the highway itself act as a source of surface‐wave energy. Along this section 120 Fairfield MagSeis ZLand 3C nodes (5 Hz corner frequency) were deployed with 40 m spacing. The nodes recorded continuously for a period of 35 days from 21 August to 25 September 2020 with a sampling rate of 250 Hz. The recorded ambient noise is dominated by high‐frequency (1–10 Hz) surface waves from traffic and quarry operations. Although this area and the Colorado plateau, in general, are seismically quiet, induced seismicity caused by brine disposal injection operations at the Bureau of Reclamation’s Paradox Valley Unit well (Block et al., 2021) is a potential source of body waves. However, most of the induced events detected by the Paradox Valley Seismic Network are low magnitude (⁠$Mw<1$⁠) with insufficient signal strength to be detected in our data (∼54 km away). In addition to the ambient noise, 10 teleseismic earthquakes (magnitude >6; Fig. S2) with high signal‐to‐noise ratio are also present in our dataset.

The objective of this study has two complementary aspects: (1) delineating the geometry of the buried basement interface and (2) imaging the depth and velocity of the sedimentary fill. We use teleseismic coda‐wave AC for the first aspect. For the second aspect, we use ambient noise SI to reconstruct dispersive surface waves traveling between the receivers, which are subsequently inverted for a 2D shear‐wave velocity profile. The methods described in this section are applied to the vertical‐component data.

AC of the transmission response of a 1D layered earth recorded at a surface receiver can be used to estimate the zero‐offset reflection response (Claerbout, 1968). This has been established for vertically incident energy propagating from a deep source below the receiver. However, in the absence of a local source of vertically incident seismic energy, teleseismic waves from distant earthquakes can be used to achieve the same goal (Abe et al., 2007; Ph?m and Tkalčić, 2018). To minimize the influence of source‐side multiples and complex source signatures of these strong events, we use the P‐wave coda for AC processing. The presence of multiple receiver‐side reflections in the scattered P‐wave coda favors the retrieval of shallow reflectivity structure (Ph?m and Tkalčić, 2018).

To ensure validity of the near‐vertical incidence assumption, we shortlist 10 teleseismic events in the epicentral distance range of 60°–95° from a global catalog of large‐magnitude events (magnitude >6) that occurred during the study period (August 2020–September 2020; Fig. S2). The earthquake catalog used in this study was accessed through the Incorporated Research Institutions for Seismology (IRIS) Wilber 3 system (see Data and Resources). The P‐wave first arrival times of these events at the Unaweep array are predicted using the IASP91 reference model (Kennett and Engdahl, 1991). We then define the P‐wave coda for these events as a 50 s long time window starting 100 s after the first P arrival estimated by the reference model.

Preprocessing of the teleseismic coda waves includes: (1) downsampling from 250 Hz recording rate to 50 Hz, (2) detrend and demean, (3) the 5 Hz low‐pass filtering to emphasize the frequency band dominated by the teleseismic coda while suppressing the local high‐frequency surface waves, and (4) automatic gain control using a 0.1 s window length to normalize amplitudes from varying source strength. We calculate the average reflectivity response at a station by first autocorrelating each teleseismic coda‐wave event recorded at that station and then stacking the AC responses of each event (Fig. 2). Similarity between the AC responses of individual teleseismic coda‐wave events suggests that the method retrieves the receiver‐side structural reflections and is not related to the source signature. Moreover, any source‐side multiples that may be present in the teleseismic coda waves are suppressed by the stacking process. The retrieved reflections have a negative phase due to the negative reflectivity coefficient (−1) for P wave at the free surface. We, therefore, reverse the polarity of the AC response before interpreting it as a reflection.

To strengthen the confidence in our AC‐derived shallow basement reflectivity response, we also attempt to retrieve the Moho reflection. For this purpose, AC is applied to the ballistic (first arrival) teleseismic phases. The results show a laterally coherent, strong phase at a two‐way travel time (TWT) of about 15.5 s, which we interpret as the Moho reflection (Fig. S3). This hypothesis is supported by analysis of Consortium for Continental Reflection Profiling (COCORP) active source reflection data in this part of the Colorado plateau where the Moho reflection and corresponding depth were interpreted at 16 s TWT and 50 km, respectively (Hauser and Lundy, 1989).

The SI entails the cross correlation of the seismic recordings at two receivers to estimate the virtual response recorded at one of the receivers due to a virtual source at the other receiver location (Wapenaar, Draganov, et al., 2010). The retrieved virtual response (or interferogram) is thus an approximation of the medium’s response to wave propagation between the receivers.

The continuously recorded data are segmented into 1 hour long windows before applying the preprocessing steps listed in the previous section. It is noted that here the surface waves in the <5 Hz band are desirable for imaging the sedimentary fill down to ∼500 m depth. We, therefore, low pass the ambient noise segments at 5 Hz. The cross‐coherence method (Wapenaar, Slob, et al., 2010) is used to calculate the interferograms. Because this method already includes spectral whitening, separate prewhitening is not required. For a station pair, an interferogram is calculated for each 1 hr time segment. These interferograms are then phase‐weighted stacked to obtain an average interferogram for the station pair. Finally, all the station pair interferograms are sorted into 120 virtual source gathers (VSGs; Fig. S4). A VSG emulates the seismic shot gather obtained from an active (explosive) source experiment conducted at a receiver location.

Owing to strong contribution of ambient noise from local surface sources (traffic, machinery, etc.) and the application of a 5 Hz low‐pass filter, surface waves dominate the VSGs in the 1–5 Hz frequency band. Visual analysis of the VSGs indicates high signal‐to‐noise ratio in offset ranges up to1000 m. We therefore limit the VSGs to 1000 m offset for further processing. Next, the VSGs are f‐k‐filtered to suppress the zero‐lag peak and the apparent low‐velocity surface‐wave reflections from the canyon walls. To derive a 2D S‐wave velocity model along the profile, the preprocessed VSGs (Fig. 3a–c) are then subjected to a workflow similar to Socco et al. (2009). This approach is based on the MASW (multichannel analysis of surface waves) technique and takes advantage of 2D multifold data, for example, for which the source–receiver geometry has a high degree of overlap. This is usually the case for ambient noise data collected along profiles for which each receiver is turned into a virtual source. The processing steps include receiver‐sorting and binning, dispersion image calculation via the Slant frequency–wavenumber method (Serdyukov et al., 2019), stacking of all dispersion images within each bin (Fig. 3d–f), semiautomatic extraction of Rayleigh‐wave dispersion curves for each bin, and a joint and laterally constrained inversion of all dispersion curves. This more sophisticated procedure provides structurally similar but more detailed and robust results when compared to conventional methods such as phase‐shift dispersion image calculation of individual VSGs and their spatially nonconstrained linearized inversion (Xia et al., 1999). The final result is a 2D shear‐wave velocity model along the profile (Fig. 4a).

The teleseismic coda‐wave AC stack (Fig. 2d) retrieves a laterally coherent phase between 0.2 and 0.6 s, which is interpreted as zero‐offset reflection from the high‐impedance contrast at the sedimentary–basement interface in the subsurface. With 3 Hz dominant frequency, the AC reflection has a lateral resolution (radius of first Fresnel zone) of ∼250 m at 0.2 s TWT and ∼400 m at 0.5 s TWT. For a complex 3D subsurface geometry, this implies superposition of reflections from different subsurface points. This effect likely plays a role at the eastern end where the profile approaches the basement outcrop (Fig. 1b) and later reflections appear. However, the clear and unambiguous appearance of the reflection along all other parts of the profile suggest a wide transverse extent of the valley floor, as also interpreted in the PSTM image of the active source line (Patterson et al., 2021; Fig. S1).

The study of Patterson et al. (2021) was conducted prior to the recently drilled well UDR‐1A, and the PSTM was depth‐converted using the smoothed PSTM velocity model in the absence of well information. However, the sonic log in the well showed significantly lower interval velocities and subsequently a shallower basement. We, therefore, re‐evaluate the PSTM time‐to‐depth conversion using the sonic log from the well (Fig. S1) and also use these velocities to depth‐convert the AC stack. This yields a good match in the basement depths (∼360 m) estimated at the intersection point of the AC stack and PSTM stack (Fig. 4). Along the 5 km long passive seismic profile, the depth‐converted AC reveals that the deepest basement interface is approximately 500 m (at ∼2600 m profile distance) below the modern surface and shallows to a depth of nearly 200 m near the eastern end. We note that the potential lateral velocity variations in the sediments are not considered in our depth conversion. However, given the limited range of $VS$ velocities (see subsequently) and the overall similarity of the sonic data and the sedimentary types in the UDR‐1A and Massey wells, we assume that the significant depth changes of the AC stack are not to be expected.

Next, we analyze the 2D $VS$ inversion along the passive profile (Fig. 4a; referred to as “passive $VS$ model”). The model indicates a layer of $VS<650m/s$ extending from the modern surface down to ∼200 m. In the same depth range, a moderate velocity increase (650–700 m/s) is observed between profile distances 1700–2500 and 4000–4700 m. These velocity variations could be the result of basement rock debris in the shallow sediments. Such debris are also exposed at that location. The overall low sediment velocities result in $VP/VS$ ratios ranging from 3 to 4 when compared to active source travel‐time tomography (Behm et al., 2019) and well log data. Similar $VP/VS$ ratios were also reported by Behm et al. (2019) for the sediment fill on the active line.

At depths below 200–250 m—the low‐velocity zone, $VS$ transitions from 800 to >1250 m/s within a depth interval of ∼200 m. We refrain from interpreting this transition zone as a geological unit, because this is most likely an artifact of the $VS$ inversion process when applied to an abrupt velocity change at the sedimentary–basement interface. This assumption is supported by the well logs, which show a rather homogenous lithology (lacustrine sediments) for the deeper section.

Combined interpretation of the passive $VS$ model and the depth‐converted AC stack (Fig. 4a) shows that the 1250 m/s $VS$ contour coincides with the basement interface reflection in the left and central parts of the profile. The low $VS$ for the gneissic basement may suggest strong weathering and fracturing, which is also observed in the well cores. However, the velocity is still significantly low when compared to the basement sonic log (⁠$VP∼4500m/s$⁠) and would imply an unrealistic high basement $VP/VS$ ratio > 3. Further, the velocity increases as the basement shallows near the eastern end (profile distance 4000–4700 m). We, therefore, attribute the low basement velocity in the western and central part to the reduced sensitivity of the 1–5 Hz surface waves to the deep layers.

Behm et al. (2019) applied surface‐wave analysis of the active source data with a traditional MASW approach (Xia et al., 1999) using phase‐shift dispersion imaging. To compare to the passive $VS$ model, we reprocess this data set with the workflow and identical parameters as outlined in Ambient noise SI and 2D V_S inversion section. The result indicates realistic $VS$ (upto 2000 m/s) for the shallow (<200 m) gneiss basement. But it lacks to correctly image the deep part where the active and passive profiles intersect. We again attribute this to the lack of low frequencies in the data, reducing the sensitivity to the deeper layers.

The proximity of the eastern part of the passive profile to the northern basement outcrop suggests that these stations image the northern (shallow) slope of the buried valley. Stations in the western part (profile distances about 1000–1800 m) were deployed close to the southern slope. The geometry thus implies that the passive profile crosses the canyon diagonally (Fig. 1b). Consequently, the central part of the profile, where the basement appears the deepest (profile distance 2200–3200 m), should cross the axis of the valley floor. The deepest point of the sedimentary–basement interface on the passive profile (elevation ∼1470 m at profile distance ∼2600 m; Fig. 4a) is deeper than the lowest‐elevation outcrops of crystalline basement near the western outlet of the canyon (elevation ∼1552 m, Fig. 1a), meaning that the basement surface slopes toward the east here. Given that rivers cannot flow uphill, this implies that the longitudinal profile is, at least locally, overdeepened. Overdeepening refers to a valley floor that is eroded below the fluvial base level, which is impossible for a river, but common in alpine glacial settings (Preusser et al., 2010; Fig. 4c). A structural explanation is unlikely because of the lack of evidence for neotectonics, and lack of evidence for major vertical offset in the active and passive seismic data. Glaciation in the Quaternary is untenable, given the low elevation relative to Pleistocene glaciation in the region as well as the burial by nonglacial Pleistocene sediments (Soreghan et al., 2007). These results are the most consistent with the hypothesis that the buried valley of Unaweep canyon formed during the late Paleozoic glaciation, although more data are needed to unambiguously confirm or reject this hypothesis.

We have successfully applied teleseismic P‐wave coda AC to image the geometry of the buried basement interface along a 5 km longitudinal section of a paleovalley from 10 large‐magnitude teleseismic events. We also estimate the longitudinal shear‐wave structure of the sedimentary fill in the valley using SI applied to high‐frequency (1–8 Hz) surface waves from highway traffic and local quarry operations. In our case, passive seismic imaging from short‐duration deployments (1 month) allows imaging of shallow targets (<500 m depth) with sufficient resolution and significantly reduced effort and costs compared to active source acquisition.

Our results show an undulating valley floor in depths ranging from 200 to 500 m and are validated by recently acquired active seismic and well data. In the context of the debate on the evolution of Unaweep Canyon, our findings support glacial carving of the bedrock. An integrated interpretation and discussion of further geodynamical implications will benefit from the UDR‐1A core analysis, which is currently ongoing.

Archiving of the continuously recorded passive seismic data used in this study is underway at Incorporated Research Institutions for Seismology (IRIS). For the teleseismic catalog, the IRIS Wilber 3 system was searched using https://ds.iris.edu/wilber3/ (last accessed September 2022) The teleseismic events used for autocorrelation (AC) and the processed virtual source gathers (VSGs) generated by interferometric processing are publicly available at https://osf.io/ckdxp/ under the project “Unaweep Passive Imaging.” The information from UDR‐1A well is available at https://osf.io/5uhck/ under the project “UDR_public.” The supplemental material contains intermediate processing products (e.g., raw teleseismic waveforms, raw VSGs) and products that complement our passive seismic workflow (active‐source prestack time migration [PSTM] stack, AC‐derived Moho reflection).

The authors declare that there are no conflicts of interest recorded.

The authors acknowledge the dedicated efforts of the field team (R. Ng, P. Ratre, Y. Wang, and Z. Wang) during the data acquisition. The authors are grateful to the Moores family for allowing access to their property and the Gateway community for their hospitality. The project was funded through the National Science Foundation Grants (NSF‐EAR‐1849623). Additional funding was provided by the Eberly Family Chair in the University of Oklahoma School of Geosciences.