Recovery of the Alberta Basin Sedimentary Structures in Southern‐Central Alberta Using Nonlinear Inversions of Receiver Functions

Secondary converted waves from receiver functions are highly sensitive to physical properties below the Earth's surface. When modeled properly, the waveforms of converted waves offer direct constraints on the impedance contrast, depth, and P‐to‐S velocity ratio pertaining to sedimentary, crustal and mantle interfaces. In this study we introduce a nonlinear waveform inversion algorithm that matches the first 5 s of receiver functions recorded in the Alberta Basin within the Western Canada Sedimentary Basin (WCSB). Our algorithm searches for the optimal thickness of the sedimentary cover and shear velocities of appropriately selected layers within and below it. Combining inversions with forward simulations, we determine the supracrustal stratigraphy from 80 regional broadband seismic stations in the WCSB. The inverted models show east tapering sedimentary layers with their thicknesses ranging from ∼6 km beneath the foothills of the Canadian Rocky Mountains to 3–4 km beneath the Alberta Basin. This finding is consistent with the sedimentary strata determined from regional well‐logging data. The sedimentary layer contains low velocity zones of variable thicknesses and amplitudes that, depending on the locations, may be caused by mechanisms involving deposition, composition or deformation history. Our shear velocity models near the top of the basement complement the existing sonic‐logs or single component seismic data and offer new constraints on the subsidence history of the WCSB. The resolved range of depths (0–14 km) effectively bridges the gap between the vertical scales of well logging (0–6 km) and those of traditional broadband analysis (>10 km) involving receiver functions and surface waves.

Beneath the sedimentary succession of the WCSB lies a crystalline basement at the western margin of the Canadian Shield (Burwash et al., 1994) that underwent more than 2 billion years of tectonic evolution. Early knowledge of the Alberta Basin structure and history benefited from petrography, as well as geochemical and isotropic aging analyses (Burwash & Kupricka, 1969, 1970Burwash & Culbert, 1976;Collerson et al., 1988;McGregor et al., 1990;Wanless, 1970) of drill core samples from a number of exploratory wells. Since the early 1990s (Ross et al., 1991(Ross et al., , 1994, a series of regional aeromagnetic, gravity and seismic surveys (Clowes, 2010;Eaton, Ross, & Hope, 1999;Welford & Clowes, 2006) provided unprecedented opportunities to systematically explore the crust beneath the WCSB. A complex network of basement domains (known as terranes) was proposed (Ross et al., 1991), which laid the groundwork for an improved understanding of the overall tectonic framework of the cratonic crystalline basement. A recent study conducted by Alberta Geological Survey (AGS) on sedimentary strata sheds new light on the Precambrian basement by interpolating between wellbore logging information (Branscombe et al., 2018). Still, limited numbers of deep wells and rising costs of drilling have hampered early initiatives in basement mapping, as many deep geological boundaries in and around the Alberta Basin remain inconclusive due to the extensive use of interpolations and extrapolations.
This study aims to provide new seismic constraints on the sedimentary strata of the Alberta Basin based on continuous waveform data collected from nearly two decades of regional passive seismic deployments. Our goal is to develop integrated shear velocity models of the sediment-basement transition beneath each receiver. By 10.1029/2022JB024813 3 of 19 accurately imaging the transition depth between exploration (<5 km) and passive (>10 km) scales using a common dataset, the resulting self-consistent shear velocity models shed new light on both the sedimentary and tectonic processes beneath the Alberta Basin.

Data and Method
Our dataset consists of 80 broadband seismic stations from seven seismic networks in southern-central Alberta and along the Cordillera (Figure 1). The recorded broadband seismic signals were generated by global earthquakes larger than Mw 5.5 from 2001 to 2018, in the epicentral distance range between 30° and 90° (Rondenay, 2009). Individual station back-azimuths differ slightly from the averages (Figure 1b), but generally offer relatively unbiased data coverage. Our test on the effect of azimuth suggests that the inverted models based on arrival angles in different azimuth bins do not differ significantly from the main findings (see Figure S2 in Supporting Information S1). To compute receiver functions (from here on, RFs), we first isolate the P wave coda using a 150-s time window (i.e., 30 s prior to and 120 s after the direct P wave), and then apply a Butterworth bandpass filter with corner frequencies of 0.03 and 5 Hz with an order of 4 to remove incoherent noise. To ensure sufficient data quality for our analysis, we further define a signal-to-noise ratio (SNR) criterion as the ratio between the standard deviation of the vertical component of an RF in the time window of 1-25 s after the onset of P (Dziewonski & Anderson, 1981) and that of a noise window of 105-5 s prior to P. For seismograms with an SNR greater than 2, we compute the RFs by deconvolving vertical component seismogram from the radial component with a water-level deconvolution method (Aki & Richards, 2002;Ammon, 1991) and apply a Gaussian filter to minimize high frequency noise. We incorporate two different frequency ranges by setting the Gaussian filter width parameters to 2.5 and 5.0, which correspond to the center frequencies of 1.2 and 2.4 Hz, respectively (Ammon, 1991). As the last step of quality control, we manually select more than 14,000 high-quality RFs for the subsequent processing, with an approximately 100 RFs retained at each station.
To recover the sedimentary structures, we perform non-linear joint inversion based on the genetic algorithm (Chen et al., 2015;Dokht et al., 2016) to match the first 5 s of the continuous waveforms after the onset of P phase. Our inversion is based on the forward modeling approach improved by Svenningsen and Jacobsen (2007). The cost function of the inversion is where W is the weighting matrix; is the forward operator (for either lowor high-frequency RFs), and m is the ground truth model. In this equation, D is a discrete approximation to the first-order derivative matrix, m 0 is the initial model or the model from the previous iteration, and μ is a regularization parameter that governs the trade-off between misfit and model norm; this parameter is determined by the trade-off curve through trial and error (see Figure  S1 in Supporting Information S1). Assuming the average S-wave velocity and P-wave velocity of 2.0 km/s and 4.0 km/s, respectively, for the sedimentary cover, the 5-s inversion window after the direct P phase is sensitive to S-wave velocities and P-wave velocities as deep as 10 and 20 km, respectively, considerably greater than the reported basement depth from regional well-logs. Necessitated by the complexity of the sedimentary structures, the genetic algorithm (from here on, GA, Chen et al., 2015;Davis, 1990;Dokht et al., 2016;Haupt & Haupt, 2004;Stoffa & Sen, 1991) seeks the optimal solutions through a global search. By introducing random perturbations to solutions in each generation, GA is resistant to local minima associated with improper choices of initial models. Furthermore, by performing a joint inversion of RFs from two distinct frequency bands at each station, our model values are sensitive to layers of variable thicknesses.
To set the initial model for the GA inversion, we adopt the CRUST 1.0 model (Laske et al., 2013) and modify its layer thicknesses and Poisson's ratios to meet the reported sedimentary conditions in the WCSB. In addition, to improve the accuracy of our inversion results, we divide the original depth increments in the published model into fixed, approximately 200 m wide segments. To improve the Poisson's ratios from CRUST 1.0, a relatively coarse global model, we also incorporate the lithology and the distribution of pore and crack shapes (Tatham, 1982) by integrating with reported values from well-logging data (Tobon, 2012) to limit the range of V p /V s ratios. The reported regional V p /V s ratios vary from 1.7 at the top of the basement to ∼5 near the free surface in a typical sedimentary basin. According to the following conversion between Poisson's ratio and V p /V s , where represents the Poisson's ratio, we adopt a linear gradient model in which Poisson's ratio varies from 0.25 at the base of the sedimentary strata to 0.48 near the free surface in correspondence with the V p /V s ratios. Within each layer, the velocity perturbation is permitted to deviate from the initial value by up to 15% during inversions. Our final inversion results generally recover over 90% of the amplitude in the observed waveforms ( Figure 2).
We adopt H-κ stacking, which was first introduced to study the depth of Moho (Gu et al., 2018;Zhu & Kanamori, 2000), to provide additional constraints on the properties (i.e., thickness and Vp/Vs ratio) of the sedimentary strata. In this approach, H represents the depth of the basement and k is the V p /V s ratio of the sedimentary layer. This algorithm performs a grid search for the most energetic stack of the direct Ps phase as well as the trailing multiple-converted phases on the basis of predicted delays relative to the direct P phase. The stacked amplitude should reach a maximum value when all phases in consideration are stacked coherently with the optimal thickness (H) and V p /V s ratio (κ).

Robustness of Inversion Results
Based on the inverted models, shear velocities in the sedimentary strata vary between 0.9 and 2.7 km/s ( Figure 3). The average value is 1.84 km/s, which agrees reasonably well with an earlier estimate of 1.9 km/s for Cretaceous sandstone in southeastern Alberta based on sonic well-logging data (Xie, 2014). Collectively, 45 out of 80 inverted velocity profiles contain a thin (∼0.2 km) layer near the surface with an anomalous high shear velocity of ∼2.0 km/s ( Figure 3a). To ascertain the existence and robustness of this layer, we perform forward simulations where a 0.2-km thick layer of high velocity (1 km/s increment) is inserted into the starting model at small depth increments within the top 1 km of the strata ( Figure 3b). The resulting waveforms from forward simulations are nearly identical for cases with and without this layer (Figure 3c), which suggests that the presence of near-surface artifacts and the limited vertical resolution is most likely responsible. Therefore, we fix the velocity of this abnormal layer as 0.8 km/s in later results.

Depth of the Crystalline Basement
The inverted waveforms, factoring the inherent noise embedded in the data, are sufficient to resolve horizontal stratifications with thicknesses greater than 200 m, especially the depths of the sediment and the shear velocity structures in the shallow basement. To assist the analysis of lateral changes in structures, we plot cross-sections of interpolated shear wave velocities from north to south and from west to east, respectively, with both profiles clearly revealing the boundaries between the sedimentary basin and the upper crust ( Figure 4a). Shear wave velocities increase with depth from 0.5 km/s to 4.2 km/s. Abnormally high (4.2 km/s) S-velocity in the upper crust may be caused by (a) a high velocity zone known as the Winagami reflection sequence (WRS) in the Alberta Basin between the depths of 6.5 and 16.5 km, caused by the presence of mafic sill complex intrusion (Ross & Eaton, 1997;Welford & Clowes, 2006); or (b) inversion artifacts induced by the mid-crustal LVZ between 11 and 28 km resulting from a relic of partial melting zone (Chen et al., 2015). The westward dipping sedimentary strata continue laterally, showing 4-6 layers with cumulative thicknesses of ∼6 km beneath the Cordillera and 3-4 km beneath the Alberta Basin. This coherent trend is consistent with that from the Alberta Geological Survey (AGS) based on well-logs (dots in Figure 4a, Branscombe et al., 2018).
Furthermore, the S-wave velocity at the basement is ∼2.7 km/s, which falls within the reported range of CRUST 1.0 (2.59-2.74 km/s in Alberta). To ascertain our inversion results, we compute the gradients of the velocity-depth curves in cross-sections AA' and BB' (Figure 4b). A strong positive gradient is evident on each cross-section, providing evidence for a sharp interface at variable depths. This interface roughly coincides with the S-velocity contour of ∼2.7 km/s (black lines in Figure 4b), which is interpreted as the top of the crystalline basement by AGS (red dots in Figure 4b).
To improve the accuracy of basement boundary depth estimations, we introduce a modified H-κ stacking method based on an earlier published algorithm of Gu et al. (2018) derived from Zhu and Kanamori (2000), which was routinely used in delineating the Moho interface (Yeck et al., 2013;Zhu & Kanamori, 2000). This modified method is more stable and flexible than the original algorithm due to a ray tracing-based approach (Niu et al., 2007) and a multilayer global crustal model. We then adopt the sedimentary structure model of CRUST 1.0 (Laske et al., 2013;Pasyanos et al., 2014) and calculate the most energetic stacks of both the direct Ps phase and the later multiple converted phases (i.e., PpPs and PpSs + PsPs). Because of a relatively weak converted phase resulting from the Precambrian crystalline boundary, as well as scattering and defocusing caused by sharply dipping interfaces and sedimentations, the H-κ stacking results are particularly accurate in resolving the depth and elasticity of relatively flat interfaces (Figure 4c).
From the computation, V p /V s ratio of the focus is ∼1.75 with uncertainties in the range between 0.01 and 0.06, which overlaps with the range provided by well-logging data (1.7-2.0; Tobon, 2012), and the uncertainties of basement depth are less than 0.4 km. Due to sedimentary structural complexities, a number of stations in our study show multiple energy foci in the H-κ domain (i.e., non-single H-κ maxima; e.g., stations TD001 and TD006; Figure 4c). In general, basement depths and V p /V s ratios computed from the H-κ stacking generally overlap with those of CRUST 1.0, showing large V p /V s ratios in the relatively shallow crystalline basement (e.g., at station REC in northern Alberta; Figure 4c). We then use the H-κ results to fix the basement depth from RF inversion, and for stations with multiple H-κ maxima we allow the inversion to determine the optimal depth. Our model of the Precambrian basement shows higher regional resolution than CRUST 1.0 (Figure 5a). In comparison with the available well-logging data (Figure 5b), our modeled basement depths at the majority of the stations are within 20% of the earlier reported values (Figure 5d). The uncertainty of the depth of the boundary between the sedimentary basin and upper crust increases toward the CDF however, which may be caused by the complex sedimentary structures along the fold-and-thrust belt.

Sedimentary Low Velocity Zone (LVZ) Revealed by the RF Inversion
In view of the various changes of the strata from west to east in the WCSB, we choose three clusters consisted of relatively concentrated stations with similar structures below in southwestern, central, and eastern Alberta (red rectangles in Figure 1). Average shear velocity maps of the three clusters ( Figure 6) are concordant with the west-to-east cross-section. We flag a boundary when the inverted velocity exhibits a sharp increase if (a) the gradient is within the depth range of the boundary of two layers from previous well-logging data, and (b) the S velocity is generally consistent with previous lithological records (Branscombe et al., 2018). Finally, we divide our inverted models into 5 to 6 layers in the sedimentary strata. Notable layers identified in the sedimentary strata are a thick (>1.5 km) Cretaceous layer and an apparent Devonian layer. In comparison, layers derived from other geological times such as the Jurassic, Triassic, and Permian are difficult to discern, which may reflect both lateral structural complexities and the resolution limit of RF imaging. Within each stratum, shear velocities exhibit significant similarities within each of the three clusters. A notable exception is the high shear velocities in the The percentage on the x-axis is calculated by the ratio of the absolute differences between our results and earlier data (normalized by the latter, i.e., ratio = |depth difference|/previous depth data). The majority of the stations show reasonable consistency (<25% difference).
upper Cretaceous layer in the western part of the study region, which is locally supported by the abundant shale and mudstone deposited during the middle of the late Cretaceous (AGS, 2019). As the upper Cretaceous layer thins toward central and eastern Alberta, and this high velocity structure diminishes as the RF inversion begins to reach its resolution limit.
In southwestern Alberta, we observe a thin, low-velocity layer at the depth of ∼2.7 km in the formation of the Jurassic/Triassic/Permian (gray shade in Figure 6a) in the average shear-velocity model. To verify the existence of this layer and quantify potential model biases, we perform several hypothesis tests based on forward waveform simulations using the reflectivity method (Bernasconi & Drufuca, 1990;Pirera & Zanzi, 1993). We also add random noise with the value of ∼40% of the signal, which is much larger than the noise level (<20%) from the actual observation, to further verify the capability of RF inversion to discern LVZ in sediments. In the first test, we generate synthetic RFs by a simple gradient model. The RFs are subsequently subjected to the same inversion process detailed earlier. The shear velocities of the recovered model are ∼4% slower than the initial model near to the boundary depth between sedimentary strata and upper crust. The extent of the recovery is reasonable, especially in view of the absence of notable artifacts surrounding the LVZ within the sedimentary layers of the inverted model (Figure 7a). In general, our inversions recover ∼90% of the peak energy of the input waveforms and the majority of the peaks and troughs of the simulated data are properly resolved (Figures 7b and 7c).  In the second test, we introduce a 400 m thick low-velocity layer at the depth of 2.8 km. The presumed shear velocity and spatial scale of the low-velocity layer resemble those of observed profiles beneath stations in the western cluster (see Figure 6a). The inversion result from a variable shear wave velocity initial model as the first test successfully recovers the low velocity layer (to ∼90% maximum amplitudes), showing nearly identical shear velocity and spatial scales as the input (Figure 7d) with an overall correlation coefficient of 0.90 (Figures 7e  and 7f). Since RF is more sensitive to velocity contrasts rather than absolute velocity (Ammon et al., 1990), we perform other inversions using two different initial models with constant S-velocity of 1.5 km/s and 2.5 km/s in the second test. The optimal solution models from 30 inversions (Figure 7d, blue and green lines) are capable of recovering the main features with acceptable differences (∼8% of maximum amplitudes and ∼7% of the correlation coefficient, respectively). To validate the robustness of the vertical resolution in our inversion, we also perform inversions based on different layer thicknesses (from 100 to 800 m, see Figure S6a in Supporting Information S1). The optimal inverted models from 20 inversions demonstrate that the model, especially the LVZ, can be fully recovered by using the selected layer thicknesses between 200 and 400 m (Figure 7g). According to the recorded well-logging data, the thickness of sedimentary strata is much thinner than 200 m. The finer details are beyond the vertical resolution limit of the RF inversions and may introduce artifacts to the outcomes. Further inversions of the stations in the southwestern cluster ( Figure 6a) using depth increments of 200 and 400 m (see Figure S6b and S6c in Supporting Information S1) confirm the existence of the LVZ, but the output model with 400-m depth interval exaggerates the LVZ thickness by ∼31%, thus suggesting that the optimal depth interval is ∼200 m.
Further analyses of waveform characteristics at broadband stations in our study reveal discrete thin low-velocity layers between 1.7 and 4.3 km in the sedimentary structures in the southwest corner of the WCSB, which is adjacent to the eastern edge of the cordillera (i.e., CDF, Figure 8a). The thicknesses of the low-velocity layers range between 200 and 600 m. In the average recovered model of stations underlain by LVZs, the depths of the shallow and deep sedimentary LVZs are between 1.7 and 2.8 km, and between 3.8 and 4.3 km, respectively ( Figure 8b). Inversion records also indicate the inverted LVZs reside within the area of better fitness (see Figure  S4 in Supporting Information S1).

Existence of Sedimentary LVZs
The RFs recorded by broadband seismic networks shed new light on the sedimentary strata beneath southern-central Alberta. Among the key findings, the RF data reveal a distinctive low velocity structure along the Cordilleran Deformation Front (CDF). It is tempting to link this observation to the Cordilleran orogeny and development of the thrust-and-fold belt, where the presence of extensive networks of large-scale cracks and faults (Price, 1994) would lead to reduced seismic velocities (Yang & Zhu, 2010). However, the earlier-reported large-scale faults (Price, 1994) are mostly confined within the fold-and-thrust belt and do not extend eastward beyond the CDF where the LVZ resides. Furthermore, these faults consistently dip westward at high angles, whereas the LVZ in our models exhibits a significantly shallower (by 0.5 deg) dip.
Data from the regional petroleum industry offers a more favorable explanation than collisional tectonics for the existence of the observed LVZs. According to the lithological data from AGS, Triassic, Permian, and Devonian formations in southern-central Alberta (Figure 6a) are mainly composited of sandstone, limestone, and dolostone, whereas middle Cretaceous and late Jurassic formations are dominated by shale and mudstone (Gradstein & Ogg, 2020;Panӑ & Elgr, 2013;Prior et al., 2013). These sequences are evidence for a system involving a cap over a reservoir, within which high total organic carbon content (TOC) in the depth range between 1.9 and 3.0 km (Passey et al., 1990) could be responsible for the observed LVZ (at 2-3 km depth) from RF inversions (e.g., Well 11-27-64-5W6). The inverted shear wave velocities of the LVZs range between 1.7 km/s and 2.6 km/s, which contrasts sharply with those of the overlying cap rock (2.3 km/s to 3.3 km/s). The abundance of shale (2.3-3.7 km/s) and sandstone (1.3-2.5 km/s) (Earle, 2015) could be responsible for the observed variations in velocity. Regional potential field observations lend further support for the inverted seismic velocities. For instance, the sedimentary LVZs revealed by RF inversions overlap with areas of known oil and natural gas depos its, which commonly exhibit reduced density (shown as lower gravity, Figure 9a) and greater magnetization while the subtle magnetic changes may be caused by the presence of extensive faults and cracks (Figure 9b).

Reconstruction of Sedimentary Structures
The Precambrian crystalline basement of the WCSB is characterized by a gentle southwestward dip from 3 to 4 km beneath central Alberta to ∼6 km beneath the foothills, the latter of which was caused by the compression and thickening of the near-surface part of the continental margin during the Cordilleran orogeny. The morphology of the sedimentary strata is direct evidence for the convergence between the accreted terranes and the North American craton. Specifically, the weight of compressed and thickened parts of the Cordillera and sediments during the orogeny induced the subsidence of the WCSB, gradually forming an east-tapering wedge ( Figure 10).
The RF imaging by broadband seismic networks provides further support for the recovery of sedimentary structures. Since the vertical dimension (defined by basement depth) is much smaller than the lateral extent of the West-to-East cross-section (Figure 4a), we make a simplifying assumption that changes near the surface (i.e., erosion and weathering) have only minor influences on the tectonic process. We then adopt a "stripping-away" method (Zheng et al., 2005, more details are shown in Figure S5 in Supporting Information S1) to reconstruct the geometry of the strata from Precambrian to Quaternary (Figure 11), as well as the historical basement subsidence (Figure 12a), by extracting the geological implications from our west-to-east shear velocity profile. When we Models of stations with and without LVZ can be found in Supporting Information S1 (see Figure S4 in Supporting Information S1).
focus on one stratum, we strip away the layers above it by applying an S-velocity threshold from previous models ( Figure 6) derived from well-logs and our inversions, and remove the portion of profile with velocity less than the threshold; this approach effectively takes a snapshot of the velocities at a certain geological time. Through multiple peeling processes, we are able to provide a quantitative, "frame-by-frame" reconstruction of tectonic evolution of the WCSB.
First, as the main driving force, the tectonic subsidence of the WCSB made major impact on the evolution of the Cordillera (Rohais et al., 2018). From the Precambrian to the Silurian, the shear velocity profile shows parallel deposits and a slight depression (Figures 11a and 11b, Figure 12aI), which could be attributed to Cambrian-Silurian rifting and the formation of a passive margin. Moving forward in time, from the late Cambrian to early Permian (∼401-∼270 Ma), the Devonian extension and the subsequent convergence from the Pennsylvanian to early Permian initiated eastward-dipping subduction (Monger & Price, 2002;Price, 1994), showing both a westward dip in the shear velocity cross-section ( Figure 11c) and a high average rate of the basement depth increase (Figure 12aII). Meanwhile, the retro-arc foreland basin to the east of the Cordillera began to develop. In the late Jurassic, a nearly flat curve on the historical basement subsidence map (Figure 12aIII) suggests extensive horizontal surficial deposition and an uplift of the basement (Figure 11d), mainly in response to the Farallon subduction and the ensuing isostatic rebound (Liu et al., 2008;Murphy et al., 2006;Rohais et al., 2018). From the late Jurassic to Eocene, the WCSB experienced the foreland basin stage. According to our reconstructed models (Figures 11e  and 11f), wedge-top and foredeep depozones initiated during this period and eventually formed the two significant depozones of the foreland basin system. The weight of the displaced and tectonically thickened structures above the crust induced the subsidence of the retro-foreland basin   (Beaumont, 1981) and the basement subsidence curve shows accelerated dipping rates (Figure 12aIV-V). Interestingly, deposits during the foreland basin stage, which took place less than 160 Ma ago, are notably absent from the interpreted geological time sequence in our reconstructed models. Possible causes are (a) detachment and (b) extrusion and erosion of older deposits, which eventually integrated into the younger deposits (Price, 1994).
Our RF-based reconstruction of the historical basement subsidence curve is consistent with the buried history curves of Well 8-17-53-21W4 in the Alberta Basin (Kalkreuth & McMechan, 1988;Nurkowski, 1985;Osadetz et al., 1990;Wright et al., 1994, Figure 12). To facilitate comparison, we compute the depths by averaging the cluster (Figure 2 round gray shade) closest to Well 8-17-53-21W4. Despite minor differences, likely affected by the resolution limit of the RFs, the subsidence rates, depths and trends of these two approaches reach an outstanding correlation coefficient of 0.92. The consistency between the basement subsidence curve derived from the RF inversion and that from well-logs shows a glimpse of the outstanding potential of multiple-station RF inversions in future basement evolution modeling (Figure 12c). While well-logging data remain "ground-truths", our new seismologically based approach is substantially more cost-effective in the investigation of the spatiotemporal evolution of sedimentary basins worldwide.

Conclusion
Understanding the formation history of the WCSB (concentrating on the Alberta Basin in this research) is crucial for both geophysical research and mineral resources exploration. In this study we introduce RF imaging, an independent, more cost-effective means than seismic exploration and well-logging efforts, to recover the geophysical properties of the sedimentary strata. Using earthquakes from hundreds of kilometers away, we are able to determine both the depth of the basement and the shear velocities within and below the sedimentary strata. Through imaging and inversions, we are able to discover low velocity sedimentary layers buried between 1.7 and 4.3 km depths near the wedge-top and southeastern foredeep depozones (in connection with tectonic deformation, crustal chemistry and lithology). Our findings enable us to reconstruct the tectonic history of the Precambrian basement, which sheds new light on the initiation and evolution of the WCSB. Figure 11. West-to-east cross-sections of shear wave velocity reformed to reconstruct the geometry features belonging to different strata from ∼600 Ma to the present. The X-axis of each plot represents the relative distance at the same spatial scales as those of Figure 4a.