3D Shear Wave Velocity Model of Salt Lake Valley via Rayleigh Wave Ellipticity across a Temporary Geophone Array

We construct a 3D shear velocity model of the Salt Lake Valley using Rayleigh waves excited by the 31 March 2020 M w 6.5 central Idaho earthquake recorded on a 168-station temporary nodal geophone network and the 49-station permanent regional network. The temporary array — deployed in response to the March 18 M w 5.7 Magna earth-quake — serendipitously recorded clear surface waves between 10 and 20 s period from the Idaho event at ∼ 500 km epicentral distance, from which we measure both Rayleigh wave phase velocity and ellipticity (H/V ratio). In addition, we employ multicomponent earthquake coda cross correlation to extend the measurements down to 5 s period. Because Rayleigh wave ellipticity features outstanding shallow sensitivity, we invert for a 3D upper crust V S model of the Salt Lake Valley. Our model shows basin structure in general agreement with and complements the current Community Velocity Model, which is mostly constrained by borehole and gravity measurements. Our model thus provides critical information for future earthquake hazard assessment studies, which require detailed shallow velocity structure.


Introduction
The Salt Lake basin is bounded by the Traverse Mountains to the south, Oquirrh Mountains to the west, and Wasatch Mountains to the east (Fig. 1a).There are three major normal fault zones in this area: the Salt Lake City segment (SLCS) of the range front west-dipping Wasatch fault zone (WFZ), the intrabasin east-dipping West Valley fault zone (WVFZ), and the range front west-dipping Oquirrh fault zone (OFZ; Hecker, 1993).The SLCS consists of three faults: the Warm Springs fault (WSF), East Bench fault (EBF), and Cottonwood fault (CF).The WSF and CF separate the Salt Lake basin from the Wasatch range to the east with prominent topographic scarps, whereas the EBF appears as a prominent intrabasin fault scarp with little known about fault geometry at depth and its association with other faults (Liberty, Clair, et al., 2021).The WVFZ and SLCS cut through the most densely populated area in Utah, and paleoseismic records reveal evidence of recurrent large magnitude earthquakes generated by both systems (DuRoss and Hylland, 2015).The ground motion from such a large earthquake would be further amplified by the low-velocity sediments in the basin and could cause major loss of lives and properties (Roten et al., 2011;Pankow et al., 2015).The 3D ground-motion simulations have been employed to assess probable seismic hazard in the area due to various earthquake scenarios (Olsen et al., 1995;Moschetti et al., 2017).
Because of the importance of the seismic velocity model on the evaluation of seismic hazard, various geophysical studies have been conducted in the Salt Lake Valley to determine the local basin structure (see references within Magistrale et al., 2008), including the determination of sediment depth from oil and gas wells, seismic exploration, gravity modeling, regional 3D crustal tomography, and analysis of regional earthquake travel time for P-wave structure near Moho.Magistrale et al. (2008) integrated the results from previous studies and constructed the initial version of the Wasatch front CVM.Since its establishment, the model has gone through several revisions, in which near-surface constraints, such as those derived from spatial autocorrelation (SPAC) microtremor arrays (Stephenson and Odum, 2010), geotechnical borehole logs, and other near-surface site response measurements were included in the model construction.A recent seismic study (Liberty, Clair, et al., 2021) using a land streamer has revealed shallow fluid and soil properties and identified faults and folds across Salt Lake City, though the results have not been incorporated into the CVM yet.Because of the limitations of both the seismic and borehole data, deeper CVM basin structures are mostly only constrained solely by earlier gravity measurements.
Traditional earthquake-based imaging is challenging in the Salt Lake area due to the overall low seismicity, dense urbanization, and high anthropogenic noise.However, the development of low-cost three-component temporary geophones opens up new opportunities for high-resolution surface-wave tomography either using ambient noise cross correlations (e.g., Gkogkas et al., 2021) or regional earthquakes.Traditional surface-wave tomography relies on dispersion measurements with peak depth sensitivity near one-third of the employed wavelengths (Lin et al., 2008).Rayleigh wave ellipticity, or H/ V ratio, has been shown to be much more sensitive to shallow crustal structure (Tanimoto and Rivera, 2008;Lin et al., 2014).Several recent studies have demonstrated the feasibility of highresolution 3D basin models constructed by jointly inverting Rayleigh wave dispersion and ellipticity measurements across densely distributed seismic arrays (Li et al., 2016;Berg et al., 2018;Liu et al., 2021).
In this study, we explore the passive data recorded by the temporary Magna aftershock array to image the Salt Lake basin structure.The M w 5.7 normal-faulting Magna earthquake occurred on 18 March 2020 near Magna, Utah, and was the largest earthquake within the Salt Lake Valley since the 1962 M L 5.2 Magna earthquake (Kleber et al., 2021).A 168-station nodal array was rapidly deployed across Salt Lake Valley following the mainshock to monitor and record the aftershock sequence (Pang et al., 2020).On 31 March 2020, an M w 6.5 earthquake occurred near Stanley, Idaho (Liberty, Lifton, et al., 2021;Pollitz et al., 2021), ∼500 km away from Salt Lake Valley, and the earthquake waveforms were recorded by the entire aftershock nodal array.In the Data and Rayleigh Wave H/V Ratio sections, we show that robust Rayleigh wave H/V ratios can be measured across Salt Lake Valley between 5 and 20 s period using the earthquake and its coda-wave cross correlations.Although not particularly sensitive to the shallow basin structure, homogeneous phase velocities between 5 and 20 s period are also obtained to constrain deeper structure and stabilize the inversion.We perform Markov chain Monte Carlo (MCMC) joint inversion of Rayleigh wave H/V ratios and phase velocities to construct the 3D basin shear wave velocity model (Shen et al., 2013).Though Rayleigh waves cannot resolve sharp discontinuities, our new 3D shear velocity model shares many common features with current CVM.This study demonstrates the potential of basin imaging based on earthquake coda interferometry.

Data
A total of 217 seismometers were used including 168 Zland three-component 5 Hz nodal geophones from the Magna aftershock array and 46 strong motion, two short period, and a broadband seismometer from the University of Utah Regional Seismic Network (Fig. 1).Clear Rayleigh wave moveout between 5 and 20 s period and 2.2 and 3.1 km/s group velocity can be observed following the M w 6.5 Stanley, Idaho, earthquake (Fig. 1b).We remove instrument responses from all sensors including flipping the vertical polarity of nodal seismometers due to the difference in the convention (nodes: positive if downward; others: positive if upward).
To extract shorter period signals, we also calculate the multicomponent cross correlations between all available station pairs using the earthquake coda signal following the direct arrival.Like ambient noise, cross correlation of earthquake coda can be used to retrieve surface-wave Green's functions (Campillo and Paul, 2003).The coda wave used is from 250 to 2250 s after the earthquake origin time; this is chosen empirically (Fig. S1, available in the supplemental material to this article) so that a longer window no longer improves the overall signal-to-noise ratio (SNR).We first divide the coda windows into 200-s-long segments with an overlap of 50 s and follow Wu et al. (2019) to temporally and spectrally normalize three-component coda waveforms in each segment to preserve the relative amplitude information between vertical and horizontal components, which is crucial for the H/V ratio measurement.After the preprocessing, we calculate the nine-component cross correlations among the north, east, and vertical components, stack the cross correlations from different segments, and rotate the horizontal motion into radial (R) and transverse (T) components (Lin et al., 2008).Example verticalvertical (ZZ) component cross correlations between station 401 and all other stations are shown as a record section in Figure 1c.Clear Rayleigh wave moveout between 1.5 and 4.5 km/s can be observed in both the positive and negative time lag.The stronger negative time lag signals indicate more energy is directed from the basin stations toward the mountain station 401 in the northeast.Similar Rayleigh wave signals can also be observed in the RR, RZ, and ZR components although with an overall lower SNR.Close examination of the cross correlations using beamforming reveals that the coda wavefield is semi-diffusive (Fig. S2) and hence satisfies the fundamental assumption behind associating the cross correlations with the empirical Green's functions (Campillo and Paul, 2003;Snieder, 2004).

Rayleigh wave H/V ratio
For each period and station, we measure Rayleigh wave ellipticity, or horizontal-to-vertical (H/V) amplitude ratios, using the vertical and radial component of the direct earthquake waveforms.Figure 2a,c show examples of earthquake waveforms filtered around 10 s period at two nodal stations, 243 and 033, which are located on the edge of the Oquirrh Mountains and near the center of Salt Lake Valley, respectively.At both locations, clear π=2 phase differences are observed between the vertical and radial components as expected for a retrograde Rayleigh wave particle motion.Figure 2e,g show the corresponding Rayleigh wave elliptical particle motions, which elongated vertically and horizontally reflecting no or thin sediment and thick sediment subsurface structures at the two locations, respectively (Berg et al., 2018).To determine the Rayleigh wave H/V ratio, we take the ratio of the vertical and radial envelope function maximums.Periods of direct H/V ratios range from 10 to 20 s, considering the number of stations that have robust measurements after selection criteria, including SNR of 5 and phase shift between vertical and radial components being within π=2 π=4.Similar H/V ratio measurements can be made by considering the coda cross correlations as the empirical Green's functions (Campillo and Paul, 2003).For each station pair, we can consider either station as the source and the other as the receiver.Two H/V ratio measurements can then be made at the receiver station using the amplitude ratio of the cross-correlation waveforms (i.e., ZR and ZZ or RR and RZ) excited either by a vertical or a horizontal force.
Figure 2b,d show the example cross-correlation waveforms around 10 s period for station pairs 401-243 and 401-033.Clear π=2 phase differences between the vertical and radial component again can be observed as with the direct arrival.The vertical component is ahead of the radial component due to the dominance and the time reversal of the negative lag cross-correlation signal (Snieder, 2004).The corresponding particle motion constructed using coda cross correlations at the two receivers (Fig. 2f,h) also reveal similar elliptical shapes compared to the direct arrival, demonstrating consistency between the two approaches.
For each period and each cross correlation, we determine the H/V ratio for the receiver station using the RR/RZ and ZR/ZZ amplitude ratios of both the causal and acausal parts of the correlograms.Several selection criteria are applied to retain only the most reliable measurements, including SNR greater than 5, interstation distance larger than half a wavelength, and the phase shift between vertical and radial components within π=2 π=4.The average and standard deviation of the mean of all measurements at the same receiver from different station pairs are then used to determine the final H/V ratio and its uncertainty.We iteratively remove spurious measurements outside three standard deviations of the mean.We have investigated taking the logarithm of H/V ratios before estimating the mean and its standard deviation in which no significant change is found.https://www.seismosoc.org/publications/the-seismic-record/For each period, we apply Gaussian lateral smoothing (Berg et al., 2018) to generate the final H/V ratio maps (Fig. 3).Periods of coda H/V ratios range from 5 to 10 s, considering the number of stations that have robust measurements, whereas the periods of direct H/V ratios range from 10 to 20 s.For the overlapping period of 10 s, highly correlated H/V ratio maps are observed although the direct H/V ratio is consistently smaller than the coda H/V ratio (note the color scale differences between Fig. 3b,c).This apparent discrepancy can be partially due to unaccounted off-great-circle and multipathing effects for direct H/V ratios and near field and not perfectly diffusive effects for coda H/V ratios.We note that Lin and Schmandt (2014) studied azimuthal anisotropy of H/V ratios across the USArray using noise cross correlations.They showed that 10 s H/V ratio measurements in northern Utah in north-northwest-south-southeast direction (same as the direction of the Idaho earthquake) are also systematically smaller than in other directions (∼8% peak to peak anisotropy) potentially due to the local stress and structure fabric orientation.To assess the uncertainty associated with the direct H/V ratios due to unaccounted for effects, we compare the direct and coda H/V ratio measurements at 10 s across the entire array.Here, we assume all direct H/V ratio measurements have the same uncertainty α.Based on the principle of uncertainty propagation, the difference between the direct (HV d ) and the mean coda (HV c ) H/V ratios at a receiver location (x i ) should have an uncertainty of E Q -T A R G E T ; t e m p : i n t r a l i n k -; d f 1 ; 3 1 4 ; 1 6 2 ; 1 in which σHV c x i is the uncertainty of the mean coda H/V ratio.Assuming the H/V ratio differences HV d x i − HV c x i are location independent and random, in this case, the normalized variance βx i : should have a normal distribution with the width equal to 1.We determine the uncertainty of our direct H/V ratio measurements (α) to be 0.211.We acknowledge that the assumption of random error might be incorrect as the variation can very well be systematic (e.g., due to azimuthal anisotropy).The uncertainty estimated here, nevertheless, allows us to combine the direct and coda H/V ratios in our revision described subsequently without imposing an arbitrary weighting factor.

Rayleigh wave phase velocity
Homogeneous phase velocity measurements between 5 and 20 s period for the entire study area are obtained to provide additional constraints on the deeper structure.We use frequencytime analysis (FTAN; Lin et al., 2008) to measure phase travel time for both the vertical direct arrival waveforms and ZZ coda cross correlations.For each period, we estimate the phase velocity and its uncertainty based on the linear regression of travel times and distances after 2π travel-time ambiguity corrections.
Zero travel time at zero distance is assumed for coda Rayleigh wave.For direct Rayleigh wave, we allow travel time at zero distance to vary to account for structure outside the area of interest.
The estimated homogeneous phase velocities increase with period as expected for typical surface-wave dispersion, while both direct and coda measurements are consistent within the overlapping periods of 9-12 s (Fig. 4c).Because of the large number of the coda cross correlations, the estimated phase velocity uncertainties based on coda cross correlations are rather small.Considering that we ignore the spatial variance of the phase velocity and to be conservative, we scale up all the uncertainty of coda phase velocity to 0.035 km/s, which matches the uncertainty of direct phase velocity at the overlapping periods (9-12 s).

Depth inversion
We use a nonlinear Bayesian MCMC method to jointly invert H/V ratios and phase velocities for shear wave velocity (Fig. 4; Shen et al., 2013).For each station location, we first extract a 1D V S reference model from the regional CVM (Magistrale et al., 2008), which we parameterize with four B-splines between 0 and 8 km depth.The empirical relationship from Brocher (2005) is used to determine V P and density from the V S model.Because Rayleigh wave phase velocities and H/V ratios have smoothly varying sensitivity with depth (Lin et al., 2014), we focus on resolving a smoothly variant V S model instead of sharp velocity interfaces such as the geotechnical layer and basement depth.We then explore the model space through 1000 iterations of random perturbations (up to 100% change in the starting model parameters) following the Metropolis algorithm, repeating the process from the starting model 6 times and characterizing model misfit by the χ 2 difference between the observations and model predictions (Herrmann and Ammon, 2004).Tests with more iterations do not show any significant change to the final result.We choose the relative weight between the H/V and phase velocity 2 to 1 because of the higher sensitivity of H/V ratios to basin structure and the fact that our phase velocity measurement is laterally invariant.As 5-20 s Rayleigh wave phase velocities are mostly sensitive to structure beneath ∼5 km depth, we effectively only use them to stabilize the inversion and obtain a reference V S structure beneath the basin.All models with misfit less than 1.5 of the minimum misfit are included in the posterior distribution.The final 1D model and the model uncertainty are calculated as the mean and standard deviation of the posterior.We then apply gaussian smoothing to all of the station-centered 1D models to create the final 3D model (Fig. 5).

Results and Discussion
The Rayleigh wave H/V ratio maps we observe (Fig. 3) show patterns consistent with known geological and basin structures (Kleber et al., 2021).Low H/V ratios are observed in the mountainous ranges including Wasatch Mountains, Oquirrh Mountains, and Traverse Mountains, and high H/V ratios are observed within the valley in particular in the area between the Warm Spring fault, East Bench fault, and West Valley fault zone.This is consistent with the shallow sensitivity of H/V ratio (Lin et al., 2014).The 10 s H/V ratios, for example, are most sensitive to V S contrast between 0-2 km and 2-10 km depth (Fig. S3).As the deeper velocity variation is likely small, the H/V ratio variation is expected to be mostly controlled by shallow basin V S structure.
Figure 4a shows the inverted 1D V S model for station 036 located between the East Bench fault and West Valley fault.Both the smoothed B-spline parameterized starting reference model and the original sharp contrast CVM model generate poor fits to the short-period H/V ratio and long-period phase velocity measurements.Compared to the starting reference model, the inverted model prefers a slower velocity in the shallowest depth consistent with the original CVM and a higher velocity between 2 and 5 km, which required to fit the elevated <8 s short-period H/V ratios (Fig. 4b).Although less sensitive to shallow basin structure, the presence of phase velocities (Fig. 4c) helps to constrain the deeper V S structure and stabilize the inversion.Despite the smoothed nature of our inversion methodology, the deepest portion of the basinbetween the West Valley and Wasatch faults-is sharply bounded and features very low velocities in the upper ∼0.5 km.This intrabasin graben structure agrees with previous gravity observation (Kleber et al., 2021).
At 0.5 km depth, the CVM features unrealistically sharp lateral boundaries likely due to the inconsistently spaced and relatively shallow borehole data that was used to create it.Our model captures the same general patterns of very low velocity in North Salt Lake and near the Great Salt Lake, but is laterally smooth.This is of critical importance to ground-motion simulations, wherein artificial boundaries can create unrealistic wavefields (e.g., Taborda et al., 2016).Despite the complete differences in data used, the two models share many common features; both models reveal (1) the major soft sediment within the intrabasin graben bounded by Warm Spring fault and the West Valley fault; (2) the thickening of the sediment toward the east, likely reflecting the surface tilting as the result of a higher slip rate along the Warm Spring fault; (3) a secondary graben (the Saltair graben discussed in Kleber et al., 2021)    (3) the overall faster sub-basin V S structure.Considering the complementary sensitivity of the data sets used in CVM and in our inversion, it is encouraging to see the overall consistency of the two models in which the discrepancy in model roughness is mostly expected.Although it is outside the scope of this study, the natural follow-up study would include combining the newly available gravity data set (Kleber et al., 2021), the previous data used in CVM construction, and our new seismic measurements and then performing a joint 3D inversion to update the CVM.

Conclusion
We constructed a 3D shear wave velocity model for the Salt Lake Valley from joint inversion of Rayleigh wave ellipticity and phase velocity, which provides new constraints to the 3D basin V S structure.Despite the differences in data type used in model construction, our new 3D model and CVM share many common features.The difference in vertical model roughness is expected considering how the models were parametrized.Our model is generally smoother and cannot resolve sharp discontinuities such as the shallow geotechnical layer (<50 m), which was constrained in CVM using borehole data.However, our model is laterally smooth, features higher basement velocities, and is self-consistent throughout the upper crust.
Our study demonstrates the potential of basin imaging based on earthquake coda interferometry.In an inland metropolitan basin, such as Salt Lake Valley, standard ambient noise tomography can be challenging due to the lack of short-period microseism and the presence of nondiffusive anthropogenic noise signals.With the deployment of a dense seismic array, however, a single significant earthquake in the region can produce sufficient coda energy that can be used to extract shortperiod surface waves.Moreover, the use of Rayleigh wave ellipticity improves the shallow sensitivity and allows basin structure to be imaged even when the sensitivity is below the basin for traditional dispersion-based surface-wave tomography.

Figure 1 .
Figure 1.(a) Station map of Magna aftershock nodal array.The gold star denotes the M w 5.7 Magna earthquake.The light blue circles are the three-component nodal geophones deployed across Salt Lake Valley between 18 March and 30 April in 2020 in response to the M w 5.7 Magna earthquake.The light red triangles are regional network seismometers from University of Utah Seismograph Stations (UUSS).Four example stations are shown by larger light blue circles and labeled; station 401 was used as the source station used for coda cross correlation shown in panel (c), stations 243 and 033 are two sample stations for H/V ratio example in Figure 2, and station 036 provides the example joint inversion result in Figure 4.The yellow star in the inset map shows the location of 31 March 2020 M w 6.5 Idaho earthquake.(b) Rayleigh wave on vertical component from M w 6.5 Idaho earthquake, filtered between 5 and 20 s.(c) Rayleigh wave from coda cross correlation on ZZ component, filtered between 5 and 20 s.

Figure 2 .
Figure 2. (a) Rayleigh wave from the Idaho earthquake at station 243 bandpassed around 10 s, normalized by the maximum absolute amplitude of both components.(b) Rayleigh wave from coda cross correlation at receiver station 243 on RR and RZ components, using station 401 as the source station, amplitude normalized.(c,d) Similar to panels (a) and (b) but for station 033.(e-h) Particle motions for panels (a-d).

Figure 3 .
Figure 3. (a) H/V map from coda Rayleigh wave at 6 s, H/V is measured at each station, the 2D background map is interpolated from station-wise measurements by gaussian lateral smoothing.(b) H/V map from coda Rayleigh wave at 10 s.(c) H/V map from direct Rayleigh wave at 10 s.(d) H/V map from direct Rayleigh wave at 14 s.

Figure 4 .
Figure 4. Example 1D Markov chain Monte Carlo (MCMC) joint inversion result at station 036.(a) Shear wave velocity versus depth showing original Community Velocity Model (CVM) as blue squares, the smoothed starting reference model as the red triangles (CVM parameterized using four B-splines), range of model space searched (green dashes), posterior models (cyan lines), and final mean model (white dots).(b) Rayleigh wave H/V dispersion curve including H/V measurements with uncertainties from coda cross correlations (black dots with error bars), H/V measurements with uncertainties from direct Rayleigh waves (gray dots with error bars), predicted H/V from CVM (blue squares), starting model (red triangles), posterior models (cyan lines), and final mean model (white dots).(c) Same as panel (b) but for Rayleigh wave phase velocity.

Figure 5 .
Figure 5. (a) The inverted shear velocity (V S ) model at 0.5 km depth.The red lines A-A′ and B-B′ show the location of two cross sections in panels (b) and (c).(b) The inverted V S model along cross section A-A′ for depth between 0 and 3 km.Relative elevation is shown on top of the cross section with same scale as depth.Approximate locations of two opposite dipping antithetic normal faults.West Valley fault zone (WVFZ) and Warm Springs fault (WSF) are identified by gray lines.(c) Same as panel (b) but for cross section B-B′.(d-f) Same as panels (a-c) but from the CVM.
west of the West Valley fault; and (4) the Warm Spring fault is