2D Near‐Surface Full‐Waveform Tomography Reveals Bedrock Controls on Critical Zone Architecture

For decades, seismic imaging methods have been used to study the critical zone, Earth's thin, life‐supporting skin. The vast majority of critical zone seismic studies use traveltime tomography, which poorly resolves heterogeneity at many scales relevant to near‐surface processes, therefore limiting progress in critical zone science. Full‐waveform tomography can overcome this limitation by leveraging more seismic data and enhancing the resolution of geophysical imaging. In this study, we apply 2D full‐waveform tomography to match the phases of observed seismograms and elucidate previously undetected heterogeneity in the critical zone at a well‐studied catchment in the Laramie Range, Wyoming. In contrast to traveltime tomograms from the same data set, our results show variations in depth to bedrock ranging from 5 to 60 m over lateral scales of just tens of meters and image steep low‐velocity anomalies suggesting hydrologic pathways into the deep critical zone. Our results also show that areas with thick fractured bedrock layers correspond to zones of slightly lower velocities in the deep bedrock, while zones of high bedrock velocity correspond to sharp vertical transitions from bedrock to saprolite. By corroborating these findings with borehole imagery, we hypothesize that lateral changes in bedrock fracture density majorly impact critical zone architecture. Borehole data also show that our full‐waveform tomography results agree significantly better with velocity logs than previously published traveltime tomography models. Full‐waveform tomography thus appears unprecedentedly capable of imaging the spatially complex porosity structure crucial to critical zone hydrology and processes.


Introduction
Nearly all terrestrial life resides in the critical zone (CZ), the volume spanning the roof of vegetation down to the top of bedrock.Soil, saprolite, and weathered bedrock within the CZ support terrestrial life by supplying water and nutrients to vegetation (e.g., Brantley et al., 2007;Hahm et al., 2013;McCormick et al., 2021).Weathering processes are fundamental to CZ structure and function by creating porosity and permeability for groundwater and by releasing nutrients from bedrock for biological uptake (e.g., Dawson et al., 2020;Hahm et al., 2019;Klos et al., 2018;McCormick et al., 2021;Meunier et al., 2007;Navarre-Sitchler et al., 2015;Riebe et al., 2016).In eroding landscapes, weathering processes sculpt a CZ architecture that generally consists of, from top to bottom, soil, saprolite, weathered/fractured bedrock, and finally intact/unweathered bedrock.While this layered framework is a useful starting point, CZ structure varies greatly, both within and between sites (e.g., Bazilevskaya et al., 2013;St. Clair et al., 2015).Understanding the magnitude and scales of CZ heterogeneity requires improved knowledge of subsurface structure.
Because direct observations of the subsurface portion of the CZ are difficult, requiring trenches, soil pits, or boreholes, geophysical imaging is often used to study the shallow subsurface (e.g., Parsekian et al., 2015).Seismic imaging has the advantage of being primarily sensitive to porosity (e.g., Callahan et al., 2020;Flinchum et al., 2018;Hayes et al., 2019;Holbrook et al., 2014), which determines subsurface water storage capacity and reflects chemical and physical weathering in eroding landscapes.In the near-surface, the seismic methods most used are first-arrival traveltime tomography (FATT) and multichannel analysis of surface waves (MASW).These methods have been reliably applied in engineering and research contexts for decades (e.g., Pasquet et al., 2016;Xia et al., 1999).In particular, FATT has been used extensively to study variations in CZ architecture at scales of tens of meters (e.g., Befus et al., 2011;Callahan et al., 2022;Holbrook et al., 2014;Huang et al., 2021;St. Clair et al., 2015).
Despite the utility of the FATT and MASW methods, CZ outcrops and boreholes typically show much more compositional and structural heterogeneity than is captured in typical seismic images.For example, borehole logs and recovered core often show small-scale weathering zones, corestones, root systems, compositional variations, and fractures that are not visible in smooth FATT velocity models at the same location (e.g., Flinchum et al., 2022;Holbrook et al., 2019;Moravec et al., 2020).This implies that typical geophysical views of the subsurface elide details about CZ structure, and by extension, hydrological and weathering processes.Improved resolution would enable detection of smaller-scale heterogeneities relevant to the hydrology, biology, and geochemistry of the critical zone.Full-waveform tomography offers a means to accomplish this.
Full-waveform tomography (FWT) is a seismic imaging technique that can improve the flexibility and resolution of seismic tomography by modeling the phase and/or amplitude of seismic arrivals, rather than only the arrival time.FWT is more flexible than other methods because it can be applied to any part of the waveform (e.g., body waves, surface waves, reflections, etc.).Resolution is enhanced because FWT leverages more of the seismic waveform than other methods of seismic tomography (e.g., Fichtner, 2010;Schuster, 2017).As a result, FWT has been widely applied in global and exploration seismology (e.g., Choi & Alkhalifah, 2012;Lei et al., 2020;Mao et al., 2016;Pratt, 1999;Virieux & Operto, 2009).To date, however, FWT has only been used to study the near surface on a handful of occasions (e.g., Cai & Zelt, 2022;J. Chen et al., 2017;Köhn et al., 2018;X. Liu et al., 2022;Sheng et al., 2006;W. Wang, Chen, Keifer, et al., 2019;W. Wang, Chen, Lee, & Mu, 2019; Y. Wang, Miller, et al., 2019).
Applying FWT to near-surface seismic data faces numerous challenges: first is the considerable computational expense and technical overhead associated with FWT as compared to FATT and MASW.Another hurdle in applying FWT is the need for domain-specific inversion strategies.For example, workflows used for inverting global seismology data, land seismic data, and marine reflection data all vary greatly (e.g., Borisov et al., 2020;Lei et al., 2020;Mao et al., 2016).FWT in the near surface is also challenging because of the strong velocity contrasts and heterogeneity in elastic properties due to the rapid compaction of regolith (e.g., Köhn et al., 2018;X. Liu et al., 2022;Sheng et al., 2006).Given the limited application of FWT to study the CZ, best practices remain ambiguous.A major goal of this paper is to present an FWT workflow for CZ seismic data.
To date, applications of FWT to near-surface problems have been aimed at a wide diversity of targets and use a variety of approaches.Some studies have focused on exclusively using body waves to inform shallow subsurface structure (e.g., Cai & Zelt, 2022).For example, Sheng et al. (2006) used first arrivals to invert for p-wave velocity (Vp) while other researchers have used both P and S body-wave phases to constrain both Vp and shear-wave velocity (vs.) (e.g., J. Chen et al., 2017;X. Liu et al., 2022).Others have primarily inverted surface waves for archeological or engineering applications (e.g., Köhn et al., 2018;Pan et al., 2019;Smith et al., 2019;Wang et al., 2019c) or inverted surface waves extracted from ambient noise data to study CZ structure and weathering (W.Wang, Chen, Keifer, et al., 2019;W. Wang, Chen, Lee, & Mu, 2019).While this set of prior work informs how FWT can be applied in the CZ, there remain several areas in which major advancements can be made.First, with the exception of W. Wang, Chen, Keifer, et al. (2019) and W. Wang, Chen, Lee, and Mu (2019), all the aforementioned studies base their forward and adjoint modeling on finite difference methods, which are not well suited to areas with complex topography (e.g., Fichtner, 2010;Komatitsch & Vilotte, 1998).Second, nearly all previous applications of FWT in the CZ focused exclusively on inverting either body or surface waves, but not both, meaning much of the available data was left unused.Third, all these past works used proprietary codes not readily available to all.Finally, these prior studies were unable to ground truth their methods and results against borehole logs.Ground truthing FWT results against borehole logs is particularly important for near-surface case studies because of the complicated geophysical environment.In the CZ, dynamic ranges of velocity (e.g., St. Clair et al., 2015), porosity (e.g., Hayes et al., 2019;Holbrook et al., 2014), attenuation (Askan et al., 2007), and anisotropy (Eppinger et al., 2021;Novitsky et al., 2018) are present and can make the implementation of FWT challenging.As such, we think it immensely valuable for results to come with a ground truth (e.g., J. Chen et al., 2017).
Our work builds on previous applications of FWT in the CZ by creating a workflow that, for the first time, combines all the following features.First, our workflow enables full elastic wave propagation across complex topography via the 2D spectral element method (e.g., Komatitsch & Tromp, 1999;Komatitsch & Vilotte, 1998).Second, our FWT strategy is informed by sensitivity analyses (e.g., Tape et al., 2010) and tailored to elucidate CZ structure using both surface and body waves.Third, our method is implemented using readily available opensource packages, which enable optimized computation on standard workstations (Chow et al., 2020;Komatitsch & Tromp, 1999;Modrak et al., 2018).Finally, we selected a data set with which we can test our results against two borehole logs.
The resulting FWT workflow is broadly applicable to the seismic refraction data sets commonly acquired for CZ science.These results bring the CZ community one step closer to routinely imaging subsurface CZ heterogeneity, including weathering profiles, fracture zones, and corestones.In the following sections, we discuss our study site, describe our FWT method and workflow, benchmark the method against synthetic data, present FWT results that show remarkable heterogeneity corroborated by downhole data, and discuss the implications of this work for improving our understanding of CZ processes.

Study Site
The Blair Wallis (BW) area is located in the Medicine Bow National Forest, ∼21 km southeast of Laramie, WY.Over the past decade, the site was studied extensively by the Wyoming Center for Environmental Hydrology (WyCEHG).The topography of BW exhibits gently undulating hillslopes (Bradley, 1987;Chapin & Kelley, 1997;Eggler et al., 1969;Evanoff, 1990).BW receives about 620 mm of annual precipitation, most of which (>90%) is snow, and has a mean annual temperature of 5.4°C (Natural Resources Conservation Service, 2015).Along the ridges, sagebrush is the dominant vegetation, while aspens, lodgepole pines, and willows appear in topographic lows.

Methods
In July 2013, WyCEHG collected a 239 m, vertical component, seismic refraction line along a ridge in BW (Figure 1).First arrival traveltimes of these data were manually picked, inverted using FATT, and published first in Flinchum et al. (2018) and later in Flinchum et al. (2022).Our workflow followed these steps, each of which is described in more detail in Sections 3.3-3.6below.First, we used the Flinchum et al. (2022) Vp model as the starting Vp model to estimate source time functions.Second, we constructed an initial Vs model using wave equation dispersion inversion (e.g., Li et al., 2016).Third, we conducted sensitivity analysis of the phases we inverted using the adjoint state method (e.g., Tromp et al., 2005).Finally, we applied FWT to the data using a custom workflow tailored to the challenges of near-surface seismic data.

Earth and Space Science
10.1029/2023EA003248 EPPINGER ET AL.

Seismic Data
Data used in this study were acquired on a linear array of 240 verticalcomponent geophones spaced at 1 m intervals and sampled every 500 µs over a record length of 1 s.The frequency response of the geophones increases from 0-4.5 Hz and then is virtually flat up to 1,000 Hz.A total of 20 sledgehammer source points generated seismic energy every 12 m along the profile, with one missing source at 24 m.The recorded waveforms show multiple distinct phases, including first arrival p-waves, a p-wave coda phase directly behind the first arrival, a clear fundamental mode Rayleigh wave phase, and a higher-mode surface wave.The body wave phases display a broader bandwidth with more energy at higher frequencies (8-56 Hz), while the surface waves exhibit narrower bands concentrated around lower frequencies (6-22 Hz) (Figure 2).In Section 3.5, we perform sensitivity analysis on each of the four phases identified in the top right panel of Figure 2.

Forward and Adjoint Modeling
The crux of FWT involves computing forward and adjoint solutions to the elastic wave equation.Earth and Space Science 10.1029/2023EA003248 In Equation 1a, ρ is density and u ∈ R 2 is a 2D displacement field, although it should be noted that we only use the vertical component of u due to collecting exclusively vertical component data.The stress tensor, σ ∈ R 2x2 , is given by the constitutive relationship, σ = C:ε, where ε = (∇u + ∇u T )/2 is the strain tensor, : denotes tensor multiplication, and C is the stiffness tensor.In this study, we assume that the subsurface is isotropic, so C is parameterized by spatially varying Lamé moduli, μ(x,z) and λ(x,z), which relate to compressional wavespeed Vp = ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ (λ + 2μ)/ρ √ and shear wavespeed Vs = ̅̅̅̅̅̅̅ ̅ μ/ρ √ .Equation 1a can be written compactly as Equation 1b, where L is a differential operator and m represents parameters for an Earth model (e.g., Fichtner, 2010).
We also use specfem2d to compute adjoint solutions to the elastic wave equation, that is to say, to solve Equation 1c.The adjoint wavefield, u † , is excited by the adjoint source, f † , and also satisfies the elastic wave equation given in 1a, albeit with additional boundary conditions.Sensitivity kernels with respect to Vp and Vs (Κ Vp and Κ Vs ) are also computed with specfem2d by correlating the forward and adjoint wavefields via Equation 2.
The adjoint wavefield, u † , is excited by using windowed phases of a particular seismic trace for the adjoint source (Section 3.5).The gradient of the misfit function with respect to model parameters, ∇ m χ, is computed via the adjoint method by summing together kernels for each hammer swing, , where in this case, the adjoint source used is dependent on the choice of misfit functional, χ.We define both the misfit function and its corresponding adjoint source explicitly in Section 3.6.For a detailed review of the spectral element method, sensitivity analysis, and adjoint tomography, we refer the reader to the following: (a) Komatitsch & Tromp, 1999, (b) Tromp et al., 2005, (c) Q. Liu & Tromp, 2006, (d) Fichtner, 2010.

Source Estimation
The source time function (STF) at each hammer location may vary depending on local ground conditions, the individual generating the source, and potentially other factors (Figure S1 in Supporting Information S1).Therefore, it is necessary to estimate a unique STF at each source location.Estimating STFs via deconvolution also helps us account for unmodeled 3D effects (e.g., Köhn et al., 2018).Because our initial FATT Vp model accurately predicts first arrival times, it is suitable for estimating the STF for each shotpoint, using (e.g., Borisov et al., 2020;Pratt, 1999): In Equation 3, j is the source index, i is the index of a particular trace, u 0 is preprocessed observed data, g is preprocessed data modeled using a STF with a unit frequency spectrum, * denotes the complex conjugate, γ is a regularization parameter which also helps avoid division by zero, and s is the estimated STF.For the source estimation, the preprocessing steps include normalizing the waveforms such that the maximum amplitude of each trace is 1 and muting out all phases except the first arrivals by applying a time window behind the first arrival pick with a duration of f 1 0 , where f 0 = 30 Hz is roughly the dominant frequency of the data.We only use first arrivals to inform the STF estimates because these are the only phases reliably fit by the FATT model.We also use data within an optimal offset range to minimize near-source effects and avoid noisy traces far from the source.Hence, in Equation 3, i ∈ window refers to traces between 100 and 175 m away from the source.Using data at these offsets to inform the STF also helps account for some of the unmodeled effects of anelastic attenuation (e.g., Borisov et al., 2020).All 20 uniquely estimated STFs can be viewed in Figure S1 in Supporting Information S1.

Earth and Space Science
10.1029/2023EA003248 EPPINGER ET AL.

Initial Shear-Wave Velocity Model Building
With an estimate of the STFs in hand, next we build a suitable starting Versus model for FWT.Rather than using a scaled version of the starting Vp model for this purpose (e.g., X. Liu et al., 2022), we find that a more rigorous prior estimate of the Vs field is necessary to perform FWT on the surface waves.To do this, we leverage the wellknown phenomenon of surface wave dispersion using the wave equation dispersion inversion (WD) technique developed by Li et al. (2016).WD is a skeletonized data inversion strategy, meaning that it uses the same forward and adjoint modeling typically employed in FWT, but fits a significantly smaller portion of the data.While this results in a lower resolution inversion, unlike FWT, the convergence of WD is almost guaranteed.
The WD method minimizes the following functional where i is the source index, N s is the number of sources, ω is the angular frequency, and Δκ is the difference between the dispersion curves of observed and synthetic data.To compute Δκ, two Fourier transforms are performed on the preprocessed synthetic and observed shot gathers, U(t,x) and U 0 (t,x) to derive Ũ(ω,κ) and Ũ0 (ω,κ) respectively, transforming the shot gathers from the time-offset domain to the angular frequency-wavenumber domain.Then, for each ω, Δκ is calculated via cross-correlation such that with R{ .}taking the real part of a complex number.The preprocessing during the WD inversion involves normalizing all traces and muting data outside the 10-75 m offset range.The shear-wave velocity of the starting model for the WD inversion increases linearly with depth.Due to the limited depth sensitivity of the WD method, the lower portion of the final WD model is altered to be a scaled version of the FATT Vp model by a factor of two.

Sensitivity Analysis
With an initial model parameterized, we can compute the traveltime sensitivity kernels (also known as bananadoughnut kernels or Fréchet Derivatives) of the four readily identifiable phases: the first arrival, p-wave coda, higher mode surface wave, and fundamental mode surface wave (Figure 2).For demonstration purposes, we used the trace and windows shown in the upper right panel of Figure 2 to back project the time-reversed particle velocity at the receiver location through the initial model (e.g., Fichtner et al., 2008;Tape et al., 2010;Tromp et al., 2005) (Figure 3).Repeating this sensitivity analysis on other source-receiver pairs yielded similar results.The sensitivity analysis is an important step in our workflow as it provides insights into which waveforms are useful for updating certain model parameters.
For example, the traveltime kernel of the first arrival is characteristic of a diving wave, exhibiting the quintessential banana shape often associated with teleseismic body waves in vertical cross-section (top left panel of Figure 3).Along the ray path, the kernel is negative, meaning that a decrease in the wave speed will result in an increase in traveltime (e.g., Tromp et al., 2005).The traveltime kernel of the first arrival has a wide sensitivity zone, indicating a broad depth range determines the diving wave arrival time.The p-wave coda, in contrast, appears to travel primarily in the near-surface and is sensitive to a narrower depth range (second row of Figure 3).Regardless of the paths the energy takes, both types of body waves are primarily sensitive to Vp (Figure 3).The fundamental mode Rayleigh wave traveltime is primarily sensitive to the upper 15-20 m, while the higher mode Rayleigh wave shows a more complicated sensitivity kernel, with deeper and shallower sensitivity where energy focuses and defocuses respectively (third and last rows Figure 3 respectively).Although other clear phases exist in the vertical-component data (Figure 2), we were unable to use sensitivity analysis to confirm that any of these arrivals are shear or converted body waves.However, our method is applicable to such phases if observed.

Full-Waveform Tomography Workflow
In this study, we used a fork of the FWT workflow manager SeisFlows published by Modrak et al. (2018), which primarily serves as a wrapper for the forward and adjoint (an)elastic wavefield solver, specfem2d (Komatitsch & Tromp, 1999).The use of the spectral element method in this study is particularly important given the topographic variation in the ground surface of our model.Specifically, the free surface boundary condition at the ground surface is rigorously fulfilled by the spectral element method, unlike in other modeling strategies, such as gridbased finite difference methods (e.g., Fichtner, 2010).
During the FWT portion of our workflow, we defined the functional to be minimized, χ, using the normalized correlative (NC) misfit norm, where ûi,j = u i,j /‖u i,j ‖ 2 and û0 ⃦ ⃦ 2 are synthetic and observed traces normalized to their L2 norms, N s is the number of sources, and N r is the number of receivers.∫ T ûi,j • û0 i,j dt is called the correlation coefficient and measures the similarity of two time series.We chose the NC norm because it is both noise-resistant and fits phase information while being insensitive to amplitude versus offset phenomena (e.g., Borisov et al., 2020;Choi & Alkhalifah, 2012;X. Liu et al., 2022).X. Liu et al. (2022) used the NC norm to apply FWT in the CZ despite working with relatively noisy data, while Borisov et al. (2020) found that using the NC norm helped contend with unmodeled anelastic and 3D effects.In choosing the NC norm, the corresponding adjoint source for each shot at each receiver location is given by f † i,j = ûi,j ∫ T ûi,j • û0 i,j dt) -û0 i,j (Choi & Alkhalifah, 2012), where i and j are the shot and receiver indices respectively.To minimize the NC norm, we iteratively update the velocity model, m, according to In the above equation, ∇ m χ, the gradient with respect to the misfit functional, χ, is computed via the adjoint method, α is a step length computed via a bracket line search, and P is a diagonal preconditioning matrix containing the discretized field P 1 1 , which is defined as where u i (x,z,t) is the synthetic wavefield excited by the ith STF.The main purpose of the preconditioner is to remove numerical artifacts caused by large amplitudes near the ground surface and to account for the geometric spreading of the wavefield.In Equation 5, H is the Hessian matrix; in practice, we approximate the Hessiangradient product, H∇ m χ, using a limited-memory Broyden-Fletcher-Goldfarb-Shanno (LBFGS) algorithm (D. C. Liu & Nocedal, 1989).Through exhaustive testing, we have heuristically estimated the number of LBFGS iterations needed for FWT to converge, after which model updates and misfit improvement become insignificant.
We note the number of LBFGS iterations used in the caption of Figure 4 (as well as in Figures 6 and 8).We used the same modeling strategy and optimization framework for the WD step of our workflow.
The FWT strategy used in this study was informed by our preliminary analysis of the data.We inverted surface waves and body waves separately, in different steps of the workflow, because the sensitivity with respect to surface waves is about two orders of magnitude higher than the sensitivity with respect to body waves (Figure 3).Hence it would require extreme scaling of the body waves to balance their contributions to model updates with those of the surface waves, which creates numerical artifacts.Separating the surface and body waves is also advantageous because of their different frequency contents.Using a multiscale approach (e.g., Bunks et al., 1995;G. X. Chen et al., 2019), we work through the frequency content of the surface waves gradually, focusing on their lower frequencies, while we step through frequencies of the body waves more aggressively to cover their broader bandwidth.
With these issues in mind, we chose to invert the surface waves first, using them to inform the upper portion of the Earth model.Then, in a quasi-layer-stripping approach, we updated the deeper part of the model using the body waves.During the surface wave step of the workflow, both Vp and Vs are updated, while only Vp is updated during the body wave step because the first arrivals and p-wave coda are primarily sensitive to Vp (Figure 3).We find that any Vs sensitivity shown in computed Fréchet derivatives for the first arrival or p-wave coda is likely a numerical artifact that degrades the fit of surface waves if incorporated into the model updates derived from the body waves (Figure 3).Throughout both the surface and body waves steps, we refrain from updating density, which we parameterize as a homogeneous field (e.g., Köhn et al., 2018;Luo & Schuster, 1991;Wang et al., 2023).We then go onto estimate source-time functions using the method described in subsection 3.3.The wave equation dispersion inversion step is discussed in detail in Section 3.5.FWT steps are explained in detail in Section 3.6.Please note that regularization during FWT is enforced by smoothing gradients with a Gaussian filter.During the surface wave FWT step, we allow model updates in Vp and Vs, although it should be noted that the model updates for Vp are relatively small.During the body wave FWT step, we only allow Vp to be updated and decrease the smoothing radius of the Gaussian filter every time higher frequency data is included.For the surface wave step, we use 15 LBFGS iterations for each frequency band.For the body wave step, we use 20 iterations for the first frequency band, and 40 iterations each for the last two bands.
Earth and Space Science This is because, with the exception of reflected phases, density primarily affects how amplitude changes with offset, but the NC norm is insensitive to such information (e.g., Borisov et al., 2020;Choi & Alkhalifah, 2012;X. Liu et al., 2022).
The preprocessing of waveforms in both steps included muting traces outside a particular offset range, applying a zero-phase Butterworth bandpass filter, normalizing all traces to have a maximum amplitude of 1, and muting various arrivals.In the surface wave inversion step, traces between 10-150 m offset were used, with 6-14, 6-18, and 6-22 Hz bandpass filters applied, while all phases arriving earlier than the higher mode were muted.In the body wave step, traces between 50-210 m offset were used, with 8-24, 8-40, and 8-56 Hz bandpass filters applied, and all phases arriving later than the p-wave coda were muted.To regularize the inversions, we smoothed the gradients by convolving them with 2D Gaussian kernels (e.g., Groos et al., 2017;Köhn et al., 2018;Modrak & Tromp, 2016;Smith et al., 2019).In the surface wave step, we used a smoothing radius of 10 m for all stages of the multiscale strategy, while for the body wave step, we used smoothing radii of 40, 20, and 10 m, decreasing the smoothing radius as we increased the frequency content during each stage of the multiscale strategy.

Workflow Validation With Synthetic Data
We benchmarked the FWT portion of our workflow (the last two boxes in Figure 4) by inverting synthetic data to illustrate the kinds of features that can be recovered.For this synthetic test, we used the same survey geometry and starting Vs and Vp models as in our real data case but added three velocity anomalies: a shallow high-velocity anomaly representing a corestone, a deeper high-velocity anomaly indicative of an area of bedrock with low fracture density, and a low-velocity zone characteristic of a fracture zone.Our FWT workflow recovers all three anomalies fairly well, although the shapes of the anomalies in the final models are not perfect (Figure 5).This is partly because the frequency content of the data was relatively low, as we used a 20 Hz Ricker for the STFs in the synthetic test.While using higher frequency data could enable better resolution, the relatively low frequency content in our synthetic test reflects what we observed in our field data (Figure 2) and estimated STFs (Figure S1 in Supporting Information S1).We also note that using a more accurate starting model could help better resolve the anomalies; however, we wanted to design a synthetic test that would reflect the effectiveness of our FWT workflow given the quality of the starting model we are generally capable of creating with FATT and WD.Nonetheless, this synthetic test bolsters confidence that we can trust relatively large-scale features (on the order of 10 m or larger) in our FWT results.

Surface Wave Step
After the FATT, source estimation, and WD, the low-frequency (6-14 Hz) surface wave data tends to fit within one wavelength but is not yet perfectly recovered (Figure S2 in Supporting Information S1).By the end of the first stage of the surface wave tomography, the phase information of Rayleigh waves is well represented by synthetics (Figure 6 and Figure S2 in Supporting Information S1).In the latter two stages, higher frequency data is progressively fit (6-18 Hz and 6-22 Hz).In these stages, only relatively small adjustments to the synthetic waveforms are needed to improve the model fits (Figure 6, Figures S3 and S4 in Supporting Information S1).Generally speaking, as the frequency content of the data being fit increases, diminishing returns in decreasing the misfit function are made (Figure 6).
In the resultant Vs model, we observe several features indicative of the increased resolution gained by performing FWT using surface waves (Figure 7).One such feature is a high-velocity zone around x = 50 m, where the 400-600 m/s velocity contours are bowed upward.The strongest vertical velocity gradients occur toward the far end of the line (x = 210-239 m).Note that the casing depths of the two boreholes correspond with shear-wave velocities of 400-500 m/s, suggesting that these velocities may be a good range to use for inferring the boundary between saprolite and fractured bedrock.It should also be noted that during the surface wave step, changes with respect to Vs are much larger than those made to Vp (Figure S5 in Supporting Information S1).

Body Wave Step
At the onset of the first stage of the multiscale strategy for the body waves, the low-frequency (8-24 Hz) data tend to fit reasonably well, implying that the FATT model and STF estimates provide a good initial parameterization Earth and Space Science 10.1029/2023EA003248 for performing FWT (Figure S6 in Supporting Information S1).In the ensuing stages of the multiscale strategy, we see that both the p-wave coda and first arrival are accurately fit by the synthetics, although the data fit degrades somewhat at further offsets where the signal to noise ratio is lower (Figure 8, Figures S6-S9 in Supporting Information S1).Convergence was slower for the body waves and required more iterations than during the surface wave step (Figure 8).
In the final Vp model, large updates can be observed, showing the impact of applying FWT to the body waves (Figure 9).Several novel features are observed in the final Vp model, including a high-velocity zone located at around 100 m, deep low-velocity zones located around x = 20 and 190 m, and various fine structures in the nearsurface.Generally speaking, vertical and lateral velocity gradients have increased substantially in several areas.
Interestingly, there appears to be more near-surface heterogeneity in the final Vp model than in the Vs model, and we discuss why this may be the case in Section 5.1.Several of the features we have noted were also observed by W. Wang, Chen, Lee, and Mu (2019), including deep low-velocity zones and more heterogeneity in Vp relative to Vs.

Comparison to Borehole Data
Two boreholes on our profile, BW1 and BW4, located at roughly x = 107 and 175 m, provide an opportunity to ground truth our FWT results.As summarized in Flinchum et al. (2022), the ∼10 m of the boreholes drilled through incompetent soil and saprolite were cased, and the deeper open holes were logged.This means that there is no ground truth in the upper 10 m where the FWT Vs model is well constrained (e.g., Figure 3), and at depths where Vs borehole logs exist, there is little to no difference between the initial and FWT Vs models (Figure S10 in Supporting Information S1).However, the borehole logs are an effective ground truth for the FWT Vp models, where the diving wave provides information on deep CZ structure (Figure 3).
The final Vp model shows much better agreement with the borehole logs than the initial model, demonstrating substantial gains from FWT (Figure 10).While the initial model created using FATT is too smooth and incorrectly estimates velocities at moderate depth (15-30 m), after applying FWT, this inconsistency is greatly reduced.At some depths, large discrepancies remain between the FWT model and the borehole log for BW1.We attribute this to the relatively low-frequency content of the data and the large discrepancy between the initial FATT model and the borehole log, which may have prevented convergence to the velocities measured downhole.Nonetheless, the borehole comparison suggests that both the existence of high-velocity and low-velocity zones in our final Vp model are true features rather than inversion artifacts.Although these low-velocity zones have not been directly observed in either borehole presented, the aforementioned synthetic tests support that we can recover such features using our FWT workflow.

Limitations, Uncertainties, and Outlook on Future Work
While our results are promising, some areas of improvement exist for our methodology, including potentially incorporating 3D modeling to more accurately recover geometric wavefield spreading.While more rigorous, incorporating 3D modeling would likely require a supercomputing cluster (e.g., Chow et al., 2020;W. Wang, Chen, Keifer, et al., 2019;W. Wang, Chen, Lee, & Mu, 2019), whereas limiting the technical overhead of FWT by implementing it on a workstation, as we have, makes the method accessible to more researchers.Furthermore, our workflow deals with unmodeled 3D effects by estimating STFs via stabilized deconvolution, which can help account for phase discrepancies caused by 2D approximations (e.g., Köhn et al., 2018).We also circumvent issues relating to amplitude mismatches by using the NC norm (e.g., Borisov et al., 2020;Choi & Alkhalifah, 2012;X. Liu et al., 2022).
Other areas of improvement for our workflow pertain to our parameterization of the Earth model.For example, including a reasonable estimate of anelasticity may help limit inversion artifacts (e.g., Borisov et al., 2020;Groos et al., 2014).In particular, we expect anelasticity to affect the surface wave tomography to a greater extent than it would the body wave tomography, as the surface waves travel significantly more cycles than the body waves.Nonetheless, given the complexities of parameterizing or inverting for a dynamic and heterogeneous near-surface anelasticity field (Askan et al., 2007), we leave this issue for future work.We also do not rigorously parameterize or invert for density in our workflow.In our modeling, density is parameterized as a homogenous scaler field (e.g., Köhn et al., 2018;Luo & Schuster, 1991;Y. Wang et al., 2023), however, changes in our parameterization of the density field could affect amplitude versus offset information (e.g., X. Liu et al., 2022).Since we exclusively use phase information in our tomography and normalize all the traces, our representation of density is not a major concern to us at this time.Given the inconsistent coupling of instruments in the data set we used, attempting to use amplitude information to constrain density would likely be ill-conceived, although future advances in instrumentation may someday make this a worthwhile pursuit (e.g., Yuan et al., 2015).Additionally, since CZ materials may exhibit significant seismic anisotropy (Eppinger et al., 2021;Novitsky et al., 2018), accounting for anisotropy may further improve FWT results, although this would likely require some prior information about the anisotropy of the study site or 3D, multi-component data coverage (e.g., Toyokuni & Zhao, 2021).
Another limitation in our FWT models relates to what extent the separately inverted Vp and Vs fields can be used to calculate Poisson's ratio in the CZ, given that the Vp FWT model is significantly more heterogeneous than the Vs model (Figures 7 and 9).One possible explanation is that the information contained in the surface waves varies from that of the body waves.Considering that anelasticity usually correlates with velocity (e.g., Askan et al., 2007;Borisov et al., 2020), the surface waves traveling primarily through lower-velocity material likely attenuate more than the body waves.This would result in the surface waves having a lower frequency content (e.g., Figure 2) and may cause models derived with them to lack high-wavenumber information.Alternatively, the Vp model may contain excess heterogeneity, namely inversion artifacts caused by the lower signal-to-noise ratio of the body waves.Furthermore, our sensitivity analysis and synthetic test show that the surface waves and body waves have different sensitivities and provide different scales of resolution.Figure 3 shows that the surface waves have sensitivity in the upper 10-20 m while the body waves have sensitivity at both shallower and deeper depths.This is reflected in the synthetic test where the deepest anomalies are recovered better for Vp than they are for Vs (e.g., Figure 5 at x ∼ 90 m and z ∼ 2,440 m).The synthetic test also shows that the surface waves provide a smoother model than the body waves, and this may be attributed to the different frequency content of the two wave types.
Future work could also make various theoretical advancements to our FWT workflow.For example, using source encoding could significantly reduce the computational cost of FWT in the CZ or reserve computational resources for incorporating 3D modeling and more data into workflows (e.g., Tromp & Bachmann, 2019).Additionally, revising the STF estimates repeatedly throughout the inversion as the Earth model is updated could more robustly contend with unmodeled anelastic and 3D effects (e.g., Borisov et al., 2020;Groos et al., 2017).More investigation into which misfit function is best for FWT in the CZ would also be beneficial.Looking into measurements that limit errors associated with source estimation and instrument response while simultaneously increasing resolution, such as the double difference measurement (e.g., Yuan et al., 2016) would be worthwhile.Trialing other misfit functions that capture traveltime differences of multiple events, such as the local traveltime inversion method proposed by Hu et al. (2020) could also be advantageous.Another promising branch of research is uncertainty quantification for FWT in the CZ, as these methods may help researchers to identify and avoid  interpreting inversion artifacts (e.g., Thurin et al., 2019).Additional improvements could be made via the introduction of horizontal component data and sources, as these may be particularly useful in better resolving Vs (e.g., J. Chen et al., 2017;Köhn et al., 2018).Incorporating these data acquisition methods could also lead to observing additional seismic phases such as refracted s-waves and converted waves.X. Liu et al. (2022) observed these phases using vertical sources and geophones, and subsequently employed FWT to build high-resolution models of Vp, Vs, and Vp-Vs ratios in the CZ.This strategy is advantageous because using body waves to constrain both Vp and Versus means that all the data have similar frequency content and sensitivity to shallow and deep CZ structure.Using horizontal sources and/or instruments could mean that the information-rich refracted swave and converted wave phases will be more consistently observed in the future.

Implications for Critical Zone Heterogeneity
One of the primary challenges in capturing and characterizing critical zone processes is the vast range in scales they span.At the smallest scales, chemical weathering occurs at the molecular and grain scale, driven by chemical reactions on individual mineral surfaces, often aided by symbiotic fungi at the micron scale (e.g., Brantley et al., 2017;Navarre-Sitchler et al., 2015;Sak et al., 2010).At larger scales, we might expect weathering to depend on climatic patterns that can vary at regional or watershed scales (e.g., Goodfellow et al., 2013).Other processes might be relevant at intermediate scales, including compositional heterogeneity, fracture zones, slopeaspect contrasts, or bedrock foliation (Callahan et al., 2022;Eppinger et al., 2021;Leone et al., 2020;Novitsky et al., 2018;West et al., 2019).This diverse set of processes acting across multiple scales creates heterogeneity in subsurface CZ structure, which is visible in outcrops (e.g., Dethier & Lazarus, 2006), corestones (Sak et al., 2010), and thin sections (e.g., Holbrook et al., 2019).Capturing such heterogeneity in the subsurface critical zone is a formidable challenge, for which improved geophysical methods like FWT are needed.
Our results show that critical zone structure is laterally heterogeneous at scales much smaller than can be attributed to large-scale forcing functions like climate or tectonic stress.For example, the depth to fast velocities associated with intact bedrock (Vp > 4,000) varies considerably throughout the FWT model.Throughout x = 80-120 m, the depth to intact bedrock is less than 30 m, and between x = 200-240, the depth to bedrock is even shallower, usually not more than 20 m.Contrastingly, between x = 0-80 m and x = 120-200 m, the depth to bedrock is consistently greater than 50 m.Contrasts at this horizontal scale cannot be the consequence of differing climate, and given the location of this profile along a ridgeline, it is similarly difficult to imagine other top-down processes (e.g., hydrology, vegetation) could produce such variability.Instead, we must seek bottom-up explanations for these changes, sourced in the local geology (e.g., composition or fractures).
Both the boreholes and the details of the full-waveform tomograms provide clues as to the causes of these strong lateral contrasts in critical zone structure.In particular, the drilling results at BW1 and BW4 combined with the FWT velocity model tell a story of two distinct weathering fronts at these locations.At BW1, we observe strong vertical velocity gradients in both the borehole log and FWT model and very few open fractures in the underlying bedrock (Figure 10).Meanwhile, at BW4, the vertical velocity gradients in the borehole log and FWT model are more diffuse, and more intensely fractured bedrock exists at depth.These results imply that the sharpness of the transition from weathered to unweathered materials depends on the fracture density of bedrock as it enters the CZ weathering engine.Indeed, the thickness of the fractured bedrock layer appears to be inversely correlated with the velocity of the underlying bedrock.In parts of the model with very fast (>4,500 m/s) bedrock velocities, there is a rapid transition to overlying saprolite, with little (or no?) weathered bedrock, while elsewhere slower deep bedrock underlies thick weathered bedrock layers-suggesting a bottom-up control on CZ architecture here (Figure 9).Such bottom-up controls could include lateral changes in composition (e.g., Bazilevskaya et al., 2013;Brantley et al., 2017), foliation (Leone et al., 2020), or fracture density (e.g., Novitsky et al., 2018).At our site, we suggest that changes in bedrock fracture density are most likely, given the observation of fracture zones in adjacent outcrops.
Given that the full waveform results show such heterogeneity, does this imply that the ray-based tomograms that have been the primary seismic tool for imaging the critical zone are wrong?To address this, we compared our FWT results with the FATT initial model.At a glance, the FWT model is much more detailed and heterogeneous than the FATT model (Figure 9).A comparison of the depth ranges of velocities associated with the saprolitebedrock transition (1.2 km/s), however, reveals that while the depth distributions are more variable in the FWT model, the average saprolite thicknesses are similar in the FWT and FATT results (Figure 11).The same can Earth and Space Science 10.1029/2023EA003248 EPPINGER ET AL. be said for the depth to intact bedrock (Figure 11).Thus, FATT accurately captures long-wavelength features in the CZ but misses smaller-scale heterogeneity.This point is further emphasized by the upper left panel of Figure 3, showing the banana-doughnut kernel for the first arrival.The large volume of the kernel implies that first arrival traveltimes are sensitive to the average velocity of a significant portion of the subsurface, and this detail is reflected in the smoothness of FATT models.These findings help contextualize previous conclusions based on FATT models, which have elucidated large-scale, first-order controls on CZ structure such as slope aspect (Befus et al., 2011), regional tectonic stresses (St. Clair et al., 2015), and foliation (Leone et al., 2020).Our findings show that FWT can build on this past research by unearthing the effects of smaller-scale processes.In other words, the average saprolite thickness at a site may reflect large-scale controls like climate or tectonic stress, while smallerscale lateral heterogeneity must have local causes, like variations in fracture density or composition.
Our results raise fundamental questions about the extent to which CZ architecture is controlled by large-scale forcing functions like climate, topography, and tectonic stress, versus local, smaller-scale characteristics of the bedrock.While past work has provided useful theories for the role of large-scale processes on CZ structure, our results suggest that smaller-scale factors also play an important role, as variability in bedrock characteristics over lateral scales of tens of meters imparts profound impacts on the overlying CZ architecture.We anticipate a concordance between the scale of forcing functions and their products.Seeking the signal of top-down processes like climate in CZ architecture will thus likely require comparing larger-scale averages across sites to filter out local variability (e.g., Callahan et al., 2022).We expect that future applications of the FWT workflow developed here will provide both new ideas and new hypothesis tests about the state and evolution of Earth's critical zone.

Conclusions
In this study, we present an FWT workflow specifically tailored to study weathering patterns in the CZ.Using existing and accessible open-source packages, we show how forward and adjoint modeling rooted in the spectral element method can be used to invert surface and body waves to constrain Vs and Vp.Our FWT results agree significantly better with borehole data than previously published FATT models.This, along with synthetic FWT experiments, bolsters confidence in our findings, which show remarkable heterogeneity in the CZ, previously undetectable using traveltime tomography.We hypothesize that local heterogeneity in Earth's weathering engine reflects local variations in bedrock composition and structure, including fracture density, foliation, and mineralogy.We suggest that FWT can be used to investigate a wide range of important CZ processes at smaller scales than previously possible.uploaded to a software repository (https://doi.org/10.5281/zenodo.10069477).Citations for these materials are included in the references section (Eppinger, 2023;Z. Liu, 2023 respectively).

Figure 1 .
Figure 1.Hill shade map of the BW study site modified from Flinchum et al. (2022), where various seismic refraction profiles collected by WyCEHG are demarcated with black lines.The seismic refraction profile used in this study is L29.The yellow dot indicates where x = 0 along the transect.The red stars show the locations of boreholes that were drilled and logged.In this work, we show borehole logs from BW1 and BW4, which are located directly on L29.

Figure 2 .
Figure 2. Left panel: a filtered (5-56 Hz) shot gather from a source located 0 m along the transect.Upper right panel: the geophone recording for the instrument located 120 m along the transect.The first arrival, p-wave coda, higher mode surface wave, and fundamental mode surface wave are highlighted.Each of these four phases is back-projected to construct two of the sensitivity kernels shown in Figure 3 (eight in total).One sensitivity kernel is with respect to p-wave velocity (Vp) while the other is with respect to shear-wave velocity (Vs).Bottom right panel: a spectrogram of the trace in the upper right panel, showing higher frequencies in the p-wave first arrival and lower frequencies in the Rayleigh wave.

Figure 3 .
Figure 3. Sensitivity kernels with respect to each model parameter (Vp and Vs) for each of the four highlighted phases in the upper right panel of Figure 2. The first column corresponds to sensitivity with respect to Vp, while the second column corresponds to sensitivity with respect to Vs.The first row shows the sensitivity of the first arrival, the second row shows the sensitivity of the p-wave coda, the third row shows the sensitivity of the higher mode surface wave, and the last row shows the sensitivity of the fundamental mode surface wave.

Figure 4 .
Figure 4.A flow chart of our FWT strategy showing both preliminary steps and FWT stages.The ray-based travel time tomography step was performed byFlinchum et al., 2022.We then go onto estimate source-time functions using the method described in subsection 3.3.The wave equation dispersion inversion step is discussed in detail in Section 3.5.FWT steps are explained in detail in Section 3.6.Please note that regularization during FWT is enforced by smoothing gradients with a Gaussian filter.During the surface wave FWT step, we allow model updates in Vp and Vs, although it should be noted that the model updates for Vp are relatively small.During the body wave FWT step, we only allow Vp to be updated and decrease the smoothing radius of the Gaussian filter every time higher frequency data is included.For the surface wave step, we use 15 LBFGS iterations for each frequency band.For the body wave step, we use 20 iterations for the first frequency band, and 40 iterations each for the last two bands.

Figure 5 .
Figure 5. Results from the synthetic FWT experiment.The left column has Vs models and the right column has Vp models.The first row shows the starting models, the second row shows the target models, and last row shows the inverted models.Velocity contours on the Vs and Vp models have intervals of 250 m/s and 500 m/s respectively.

Figure 6 .
Figure 6.Top panel: The evolution of the misfit function with each FWT iteration, segmented by each stage of the multiscale strategy.Middle panels: histograms of the coefficients, ∫ T û • û0 dt, for all traces before and after each stage of the multiscale strategy.Bottom panel: Preprocessed (body waves are muted and 6-22 Hz bandpass filtered) waveforms after surface wave FWT.

Figure 7 .
Figure 7. Top panel: the initial shear-wave velocity model derived with WD.Bottom panel: shear-wave velocity model after FWT.Velocity contours at 75 m/s intervals are also shown in white.The gray rectangles show the locations of borehole casings.Please note the limited elevation range of these plots.

Figure 8 .
Figure 8. Top panel: The evolution of the misfit function with each FWT iteration, segmented by each stage of the multiscale strategy.Middle panels: histograms of the correlation coefficients, ∫ T û • û0 dt, for all traces before and after each stage of the multiscale strategy.Bottom panel: Preprocessed (surface waves are muted and 8-56 Hz bandpass filtered) waveforms after body wave FWT.

Figure 9 .
Figure 9. Top panel: The initial p-wave velocity model, which is very close to that published by Flinchum et al. (2022), but with small changes due to surface wave FWT.Bottom panel: The p-wave velocity model after FWT.Velocity contours at 500 m/s intervals are also shown in white.The gray rectangles show the locations of borehole casings, while the gray lines show where the boreholes were logged.

Figure 10 .
Figure 10.Comparison of the FATT and FWT Vp models with the borehole logs from BW1 (left) and BW4 (right) and expanded views of the bedrock observed in optical logs over a 2 m depth range in each hole.Note higher fracture density visible in BW4, which corresponds to lower Vp in that hole.

Figure 11 .
Figure 11.Overlain histograms of the depth to saprolite (left) and intact bedrock (right) in the FATT (blue filled bars) and FWT (red hollow bars) models.The thick vertical lines indicate averages of the distributions displayed in the histograms.Note that in both cases, the averages of the FATT and FWT distributions are nearly identical.