Long Period Rayleigh Wave Focal Spot Imaging Applied to USArray Data

We demonstrate the effectiveness of seismic dense array surface wave focal spot imaging using USArray data from the western‐central United States. We study dispersion in the 60–310 s period range and assess the image quality of fundamental mode Rayleigh wave phase velocity maps. We apply isotropic spatial autocorrelation models to the time domain zero lag noise correlation wavefield data at distances of about one wavelength. Local estimates of the phase velocity, its uncertainty, and the regression quality imply overall better ZZ relative to ZR or RZ results. The extension of the depth resolution compared to passive surface wave tomography is demonstrated by the inversion of three clustered dispersion curves from different tectonic units. We observe anisotropic surface wave energy flux and the influence of body wave energy, but sensitivity tests at 60 s targeting the data range, correlation component, and processing choices show that the ZZ focal spots yield consistent high‐quality images compared to regional tomography results in the 60–150 s period range. In contrast, at 200–300 s the comparatively small scales of the imaged structures and the imperfect agreement with low‐resolution global tomography results highlight the persistent challenge to reconcile imaging results based on different data sources, theories, and techniques. Our study shows that surface wave focal spot imaging is an accurate, robust, local imaging approach. Better control over clean autocorrelation fields can further improve applications of this seismic imaging tool for increased resolution of the elastic structure below dense seismic arrays.


Introduction
The ambient seismic wavefield is predominantly generated by sources located on the surface of the Earth and is therefore mainly composed of surface waves.The first imaging approach that utilizes the ambient field data is the spatial autocorrelation or SPAC method (Aki, 1957), which is widely implemented for geotechnical engineering and near-surface exploration problems using an optimized size and geometry of the deployed SPAC arrays (Foti et al., 2018;Hayashi et al., 2022).The SPAC method is related to the cross-correlation approach (Prieto et al., 2009) that reconstructs the deterministic wavefield between two stations from the ambient field records (Shapiro & Campillo, 2004).Similar to earthquake tomography, this deterministic wavefield can be used for subsurface imaging on local, regional, and global scales (Lin et al., 2008;Lu et al., 2018;Retailleau et al., 2020;Sabra et al., 2005;Shapiro et al., 2005;Stehly et al., 2009;Yao et al., 2006).
The cross-correlation wavefield is equivalent to time-reversal experiment results (Derode, Larose, Campillo, & Fink, 2003;Derode, Larose, Tanter, et al., 2003).Considering a reference station in a dense seismic array such as the USArray used in this study, the time-distance representation of the surface wavefield at negative and positive times represents the converging and diverging waves, respectively (Figure 1), which form the basis for passive tomography applications.At zero lag time, the cross-correlation wavefield results from the interference pattern of the reconstructed converging and diverging waves.The equivalence of correlation and time reversal explains the emergence of the so-called time domain focal spot at the origin of a time-reversed surface wavefield, which differs from the propagating signals studied in seismic tomography.The zero lag time cross-correlation amplitude field is the focal spot, a term adopted from time-reversal acoustics (Fink, 2008).The focal spot is the time domain equivalent of the spatial autocorrelation or SPAC field in the frequency domain.The equivalence allows us to use the SPAC formulation or parametrization of the amplitude distribution (Aki, 1957;Haney et al., 2012) for our time domain analysis.
Applications in acoustics (Fink, 1999), passive elastography (Benech et al., 2009), and related medical imaging techniques (Zorgani et al., 2015) have long been using properties of the shear wave focal spot for imaging.The similarity between the multi-channel ultrasound transducers and modern dense seismic arrays (Giammarinaro et al., 2023) suggests to systematically explore the application of seismic focal spot imaging.The relevance is also implied by initial surface wave focal spot applications (Gallot et al., 2011;Roux et al., 2018), notably by the imaging results obtained from data recorded by 1,100 sensors deployed across a square 600 × 600 m 2 domain that show distributions of the local phase speed, azimuthal medium anisotropy, and a proxy for intrinsic attenuation (Hillers et al., 2016).These initial observations are further supported by results from numerical studies that demonstrate an overall robust estimation of the Rayleigh wave phase velocity from wavefields featuring biasing body wave energy or nonuniform background illumination (Giammarinaro et al., 2023), and good spatial resolution of lateral heterogeneities (Giammarinaro et al., 2024).
The deployment of dense arrays and access to the data sets have been stimulating the refinement and development of existing and new noise-based imaging methods, respectively, which have provided high-resolution images of tectonic features and structures in the crust and upper mantle.The USArray Transportable Array (TA) (Kerr, 2013) recorded an example influential data set.The ambient noise-based tomographic methods or variants applied to USArray data include ray theoretical approaches (Lin et al., 2008;Porter et al., 2016), eikonal tomography (Lin et al., 2009;Shen & Ritzwoller, 2016), and finite-frequency tomography (Lin & Ritzwoller, 2010;Zhao et al., 2020), whereas Ekström et al. (2009) and Ekström (2014) combined the SPAC method (Aki, 1957) with a tomographic approach to invert phase velocity estimates along ray paths.The interstation distance is scaled by the corresponding seismic wavelength.The energy is scaled by the maximum amplitude in each distance bin.The reconstructed wavefield is dominated by the familiar converging and diverging propagating surface waves at negative and positive lag times, respectively.In addition we highlight at zero lag time the azimuthally averaged spatial autocorrelation field, which is the time domain focal spot.The visible amplitude ridge between 1λ and 2λ corresponds to the second maximum of the J 0 Bessel function.The relatively wide period band illustrates dispersion of the propagating waves, but it degrades the focal spot shape (Giammarinaro et al., 2023).
For the western and central U.S. noise-based tomographic images using ray theory have been inverted using data from the period T range 8-40 s (Lin et al., 2008(Lin et al., , 2009;;Porter et al., 2016), while accounting for finite-frequency effects increases the long period limit up to T = 80 s (Zhao et al., 2020).Longer period surface waves can be considered by including larger interstation distances (Yang, 2014;Zhao et al., 2017) or by using earthquake data (Jin & Gaherty, 2015) in combination with ambient noise correlations (Lin & Ritzwoller, 2011).Very long period surface waves for T up to 180 s are studied using only earthquake data (Babikoff & Dalton, 2019).Most of these approaches are efficient for large sparse networks where the observation is obtained in the far field.Since data associated with interstation distances of less than two or three times the seismic wavelength λ are typically not included in ambient noise tomographic inversions (e.g., Chapter 5 in Nakata et al., 2019), surface waves with wavelengths equal to or longer than the dimension of the array are not considered, even if this propagation cutoffdistance can potentially be reduced to 1λ (Luo et al., 2015;Zhao et al., 2017).
Here we apply the focal spot imaging method to cross-correlation functions of continuous records from the western part of the contiguous U.S. to obtain fundamental mode Rayleigh wave velocity dispersion measurements between 60 and 310 s, thereby extending previous kilometer-scale focal spot applications (Hillers et al., 2016;Roux et al., 2018) to the continental scale.Our imaging technique estimates local phase velocities from narrow-band filtered focal spots using SPAC Bessel function parametrization (Haney et al., 2012), and the period dependent phase velocity maps are compiled by the iterative application of this approach using each station as reference location.That is, in contrast to tomography, we do not employ a 2D inversion for the phase velocity distributions on a regular grid.This inversion-free local imaging approach is computationally efficient and thus advantageous for dense array noise correlation data.The parameter estimation approach yields a measure of the local phase velocity uncertainty, which is important for the appraisal of the imaged features and for the evaluation of the downstream shear wave velocity model.At short periods around T = 60 s the comparison with tomographic phase velocity maps in the literature facilitates the assessment of the focal spot image quality and the resolution-uncertainty trade-off (Giammarinaro et al., 2023(Giammarinaro et al., , 2024;;Latallerie et al., 2022).At long periods up to T = 300 s our application explores the analysis of focal spots associated with wavelengths that approach the array dimension, which exceeds the current 150 s period limit using USArray ambient noise data (Zhao et al., 2017).We analyze the associated increase in vertical resolution using local autocorrelation field data from distances around one wavelength.
In this work we discuss focal spot-based Rayleigh wave phase velocity maps and their sensitivity to processing choices and parameters across a wide period range.The high similarity to tomography results at shorter periods demonstrates the relevance of the local imaging technique for the construction of 3D shear wave velocity models, that are, however, not computed in this work.In Section 2 we describe the application of the focal spot imaging method to USArray data.The implementation includes data selection, preprocessing, correlation, filtering, focal spot compilation, and the non-linear regression for phase velocity and uncertainty estimates.In Section 3 we present results on the nonuniform background illumination, dispersion curves in the 60-310 s period range, the period dependent lateral phase velocity distributions, and extensive sensitivity tests that investigate how these maps change in response to processing and parameter estimation choices.For the best phase velocity maps we evaluate the consistency of the imaged features to published results and the associated uncertainty estimates.We discuss the implications for the lateral and vertical resolution of the obtained results and the relation to tomography in Section 4, and conclude in Section 5 that focal spot imaging is a complementary seismic dense array surface wave imaging technique.

USArray Data Selection
We use USArray continuous records from the BK, CI, TA, and US network stations between 125°-90°W and 26°-50°N from 2004 to 2012.Most stations belong to the TA network with an average interstation distance of about 70 km.Starting from the West, the stations operated about two years at each location.In this rolling deployment the synchronicity of any reference station with other array stations scales with distance and azimuth (Figure 2).The sensors are three-component broadband seismographs and sample with a rate of 40 Hz.

Preprocessing and Noise Correlations
We download three-component continuous data recorded by the stations shown in Figure 2. We implement the standard signal processing operations that are typically applied in noise correlation imaging (Boué, Roux, et al., 2014).That is, the focal spot data processing does not require a specifically tuned processing chain, and correlation databases tailored for surface wave imaging can readily be used for focal spot imaging.Preprocessing aims to mitigate the effect of strong sources that are localized in space and persistent in time (Chapter 5 in Nakata et al., 2019).We decimate daily records to 1 Hz, remove the instrument response, and apply a bandpass filter.We cut the daily records into 4 hr segments and normalize the amplitude first in the frequency domain, with spectral whitening in the 2.85-340 s period band, and then in the time domain with amplitude clipping at three times the standard deviation of the distribution in each 4 hr segment.Tapering is performed to avoid spectral leakage.
We compute 4 hr time window cross-correlations in the vertical, North-South, and East-West ZNE system and stack them linearly.The relatively long maximum lag time of τ = 1,000 s is necessary for the dispersion analysis bandpass filtering up to T = 310 s.The cross-correlations are normalized by the square root of the energy in both traces.We exclude station pairs with less than 100 synchronous operation days.We rotate the cross-correlations from the ZNE to the vertical, radial, and transverse ZRT system along the great circle path.As shown in Section 2.3, the focal spot shapes are influenced by the anisotropic surface wave illumination, which can potentially be mitigated using data driven optimal rotation angles (Roux, 2009).
For the dispersion analysis we apply the same Gaussian narrow-band filter that was used in the numerical focal spot studies of Giammarinaro et al. (2023Giammarinaro et al. ( , 2024)).The comparatively wide bandpass filter used by Hillers et al. (2016) biases the focal spot reconstruction by causing an apparent attenuation effect (Giammarinaro et al., 2023).We zero pad the correlations before filtering to create more frequency samples and to better resolve the center frequencies f c .The actual f c values are those frequency values of the discrete Fourier transform of the signal that are closest to the nominal and reported analysis period in the 60-310 s range.

Focal Spots
Under the assumption of a diffuse wavefield, three-component cross-correlations of ambient noise approximate the 3 × 3 Green's tensor considering the combinations between the Z, R, and T components.The equivalence between the time domain focal spot and the spatial autocorrelation function in the frequency domain (Tsai & Moschetti, 2010) suggests that we can use the classical SPAC result (Aki, 1957) and its extension (Haney et al., 2012) for the focal spot parametrization in a lossless medium.In this work we do not consider the effect of anelastic attenuation on the spatial autocorrelation.At the here employed one wavelength distances the effect is small (Nakahara, 2012) and difficult to constrain considering the influence of directional incidence and body wave energy on the focal spot properties.
The ZZ, RR, and TT component autocorrelation fields are described by combinations of Bessel functions that allow wave speed estimates of the fundamental mode surface waves.The corresponding expressions are valid for a scalar wavefield, while the full Green's tensor including the non-diagonal terms extend these results to the study of a vector wavefield.Similar approaches facilitate the quantification of higher-mode surface wave dispersion curves (Wang et al., 2019).Here we use the SPAC formulas (Haney et al., 2012) to obtain fundamental mode Rayleigh wave velocity estimates from ZZ, ZR, and RZ focal spot data.This implementation estimates the phase velocity locally using short distances of about one wavelength around the origin, which is similar to passive elastography (Catheline et al., 2008;Gallot et al., 2011).The absence of Green's function near-field components at short distance is illustrated by the analogy to the time-reversal configuration, since the refocusing wavefield recorded by the time-reversal mirror does not include the evanescent near-field waves, and the autocorrelation is governed by the diffraction limit (Fink, 2006).
The here relevant correlation coefficients ϕ for the ZZ, ZR, and RZ components are (Haney et al., 2012) where P R (ω) denotes the Rayleigh wave power spectrum, R the horizontal-to-vertical displacement ratio of Rayleigh waves, ω the angular frequency, r the radial distance, c R the Rayleigh wave phase velocity, and J n denotes a Bessel function of the first kind of order n.The local Rayleigh wave properties are estimated from the ZZ, ZR, and RZ focal spots that follow the same parametrization as the SPAC expressions in Equation 1.At short distances, the here not considered RR and TT components contain both Rayleigh and Love wave energy (Haney et al., 2012), and an unbiased estimation of the two associated velocities requires a separation of the two components.
Figure 3b shows ZZ, ZR, and RZ component focal spots for the periods 100, 200, and 300 s together with the wavenumber k = ω/c R decomposition of the 100 s results in Figure 3a.The wavenumber domain images are obtained by the application of the 2D Fourier transform (2DFT) to the spatial autocorrelation amplitude distributions, and they can help analyze focal spot properties (Giammarinaro et al., 2023;Hillers et al., 2016).The focal spot shapes for all three components (Figure 3b) is not circular as the direction independent Bessel function model suggests.The distribution of the first one or two nodal lines delineate an ellipsoidal shape with a short axis oriented in the Northwest-Southeast direction.We understand this effect to result from anisotropic surface wave incidence (Giammarinaro et al., 2023) with a predominant direction oriented toward the North Pacific, which is a region of persistent noise generation in the considered frequency range (Ermert et al., 2017;Rhie & Romanowicz, 2004).A more subtle feature in the 100 s amplitude patterns (left column in Figure 3b) are the consistently reduced speed toward the Southwest across the Basin and Range province, where the ZR and RZ results also display a less coherent focal spot amplitudes.Such spatial variations of the autocorrelation field can potentially be used to increase the resolution.
The incidence direction is better seen in the wavenumber domain representation (Figure 3a) which shows increased energy in the Northwest-Southeast direction for all three components.Compared to the unidirectional and hence more extreme synthetic cases used to study biasing directivity effects (Giammarinaro et al., 2023), these spectral amplitude patterns are more complex, that is, they show energy at different directions with different strength.The directive yet overall more evenly distributed energy supports robust parameter estimation.The J 1 shape of the ZR and RZ data leads to the double-ring feature in the 2DFT results (Giammarinaro et al., 2023), and non-zero energy can be resolved at the center of the spectra.Different to beamforming, the peak energy does not indicate the actual wavenumber and hence speed of the surface waves.These are effects associated with the Fourier transform properties of a finite spatial distribution, and they becomes less significant when data from longer distances are included (Giammarinaro et al., 2023).The ZR and RZ focal spots in Figure 3b are more noisy compared to the ZZ data, which leads to overall higher uncertainties of the c R estimates using non-linear regression.

Phase Velocity and Uncertainty Estimates
We estimate the Rayleigh wave phase velocity by optimizing the fit between the focal spot data and the autocorrelation functions in Equation 1.We use the Python package (Giammarinaro & Hillers, 2022) developed for the numerical analysis by Giammarinaro et al. (2023), which uses the Levenberg-Marquardt algorithm to solve a non-linear least-squares regression problem (Aster et al., 2013).For a regular station distribution the number of focal spot samples scales with distance.This can put more weight on data away from the origin (Giammarinaro et al., 2023), however, the influence is small and we do not make use of the option to apply distance dependent weights.We reconstruct focal spots for each reference station in the array from the narrow-band filtered zero lag correlation amplitude fields.We perform the parameter estimation using the radial distance r dependent but direction independent amplitude distribution models to fit the focal spot data A(r, τ = 0) (Giammarinaro et al., 2023) (2) Note that while Equation 1 parametrizes the correlation coefficients in the frequency domain, Equation 2 provides the equivalent time domain focal spot Bessel function models.This equivalence results from the behavior of Fourier transform pairs at t = 0.In Equation 1, the theoretical amplitude factors are P R (ω) and ∓P R (ω)R for the ZZ and the ZR and RZ components, respectively.In the data regression, however, σ parametrizes the focal spot amplitude that depends on the Rayleigh wave power spectrum but also on processing effects that can influence the Rayleigh wave coherency on the vertical and radial component (Hillers et al., 2016).For a given frequency ω the two free parameters are the Bessel function argument wavenumber k-or the Rayleigh speed c R -and the amplitude factor σ. The data range r fit is the distance between the reference station and the most distant receiver from which data are included in the regression analysis.
Our focal spot approach has the analysis of the distance dependent autocorrelation amplitude in common with the frequency domain ESPAC method (Okada & Suto, 2003).This and other SPAC variants are typically implemented in geotechnical applications (Asten & Hayashi, 2018;Foti et al., 2011), where all wavelengths are analyzed over the same array aperture data range, which is usually on the 10-100 m scale.This full-array analysis style has also been adopted by Prieto et al. (2009) and Ekström (2014) at regional scales.Again, this differs from our frequency dependent data range that emphasizes a consistent local measurement across the USArray.
The optimal values for k or c R and σ are obtained in an iterative procedure.We perform the regression three times in which k and σ are simultaneously estimated.The actual autocorrelation value at r = 0 is not considered, it is typically on a different scale because the employed normalization affects cross-and autocorrelation differently.
The approach is illustrated in Figure 4 for ZZ, ZR, and RZ focal spot data at 100, 200, and 300 s, where the gray dots represent the data and the orange and black lines show the estimated models after the first and third iteration, respectively.In the first iteration we consider all amplitude data without distance restriction to obtain a data driven initial estimate of k 1 (orange model in Figure 4), which we need to compute r fit = 2πn/k in the following iterations.
The n value defines the short data range r fit .Here we explore values between 0.5λ and 1.5λ but present mostly results with r fit = 1.2λ.We do not need the σ 1 result associated with this first iteration.We perform the regression a second time now including data within the range r fit = 1.2λ.As said, σ is an amplitude factor linked to the energy The estimated models from the first and third iterations explained in Section 2.4 are plotted in orange and black, respectively.The scale of the two models is indicated by the correspondingly colored axis.The ZZ data (left column) show a better signal-to-noise ratio compared to the ZR (center column) and RZ (right column) data.The consistently higher ZZ amplitudes at short distances is attributed to the influence of coherent P wave energy.
of Rayleigh waves and to the processing-and normalization-affected coherency level of the autocorrelation field.
After the second iteration, we normalize the amplitudes by σ 2 to obtain consistent signal-to-noise ratio (SNR) and uncertainty estimates.The difference between k 1 and k 2 is typically a few percent, however, shorter r fit distances improve the lateral resolution (Giammarinaro et al., 2024).The third iteration (black model in Figure 4) is performed on the scaled amplitude data to access the normalized residual sum of squares (RSS) as SNR or goodness-of-fit proxy from the data misfit to the Bessel function models that is used to compute the uncertainty of the final k 3 estimate.Since r fit is not changed from the second to the third iteration, k 2 = k 3 , and σ 3 = 1 for all three components.
We apply the iterative parameter estimation to every array station and obtain thus estimates of the frequency dependent local Rayleigh wave speed (Hillers et al., 2016).We highlight that our approach has the advantage of yielding a straightforward estimate of the formal k uncertainty or standard error ɛ k at each reference station location and hence image pixel (Aster et al., 2013) where dof denotes the degrees of freedom, that is, the number of data points minus two, the number of fit parameters, and C k is the diagonal element of the parameter covariance matrix C associated with k.The error on the Rayleigh wave phase speed ε c R is then obtained by error propagation.For the example in Figure 4, the c R and ε c R uncertainties for the 100, 200, and 300 s ZZ data are 4.11 ± 0.03 km/s (0.8), 4.53 ± 0.03 km/s (0.6), and 5.25 ± 0.05 km/s (0.9), respectively, where the values in parenthesis are the errors in percent.These small errors in the 1% range imply overall well constrained c R estimates.For the ZR and RZ data the triplets are 3.99 ± 0.04 km/s (1.0), 4.26 ± 0.07 km/s (1.7), 5.36 ± 0.27 km/s (5.0), and 3.91 ± 0.05 km/s (1.3), 4.51 ± 0.09 km/s (2.0), 4.62 ± 0.21 km/s (4.5), respectively, all showing similarly good percentage values except for the longest period.For the 100-300 s ZZ component data in Figure 4 we estimate RSS values between 2 and 56 for the r fit = 1.2λ data range.It is the coherent elliptical shape of the ZZ spatial amplitude patterns (Figure 3b) that governs the data scatter around the azimuthal average model.In contrast, the obtained RSS values between 4 and 1,700 for the ZR and RZ data quantify the overall lower focal spot quality involving the radial component, in particular at 200 and 300 s.Although the first nodal line indicates a more circular shape of the amplitude patterns compared to the ZZ results (Figure 3b) that is compatible with synthetic observations on the impact of anisotropic illumination (Giammarinaro et al., 2023), the ZR and RZ data in Figure 4 show a significantly increased level of incoherent fluctuations around the best fitting model.
At short distances the ZZ amplitude samples in Figure 4 consistently exceed the estimated J 0 model associated with surface wave refocusing.The synthetic focal spot tests by Giammarinaro et al. (2023) suggest that this amplitude difference is caused by interfering body wave energy with near-vertical incidence.This phenomenon is observed in the focal spot analysis in the San Jacinto fault zone environment (Hillers et al., 2016).In this case the 600 × 600 m 2 array records body wave energy that is excited by the abundant microseismicity and then channeled along the low-velocity fault structure.Our USArray correlations contain teleseismic P wave energy excited by moderate and large earthquakes (Boué, Poli, et al., 2014;Lin et al., 2013).
As a result, a small-amplitude long wavelength sinc function associated with coherent P wave energy with high apparent velocity is superimposed on the J 0 surface wave model (Hillers et al., 2016).In the wavenumber domain the corresponding energy around the origin (Figure 3a) is superimposed on the small-k amplitude component associated with the 2DFT properties of the finite autocorrelation pattern (Giammarinaro et al., 2023).The error on the ZZ Rayleigh wave phase speed estimates depends on the energy ratio between surface and body waves (Giammarinaro et al., 2023).Interestingly, it also depends non-linearly on the data range r fit .Numerical experiments for azimuthally symmetric P wave incidence show that for r fit = 1.0λ and r fit = 1.5λ the c R value is first under-and then overestimated, respectively, by a few percent (Giammarinaro et al., 2023).This and the P wave energy minimization results discussed in Section 3.4.4imply that the potential bias associated with the r fit = 1.2λ parameter is not significantly affecting our main imaging results.In contrast, the vertical-radial component surface wave focal spots are insensitive to the P wave component.This can make them a preferred target for the analysis, however, as demonstrated, the radial component correlations are characterized by significantly lower data quality.
In the next Section 3, we first explore the anisotropic illumination effect (Section 3.1), discuss results of the dispersion analysis using the iterative parameter estimation (Section 3.2), and then assess the quality of the obtained phase velocity images (Section 3.3) and their sensitivity to processing choices and parameters (Section 3.4).

Background Illumination
Applications of correlation-based methods for passive imaging assume a diffuse ambient noise wavefield with an ideally isotropic distribution of seismic energy (Piña-Flores et al., 2021;Shapiro et al., 2000).Isotropic illumination is supported by a homogeneous distribution of the original noise sources and multiple scattering.However, low-frequency ambient noise is predominantly excited in the oceans or along coastlines, and the positive effect of multiple scattering (Derode, Larose, Campillo, & Fink, 2003) is frequency dependent and may be insufficient to mitigate strong, localized source effects.In the frequency range below 20 mHz or above 50 s period, which is of interest here, the noise wavefield, also referred to as seismic hum, is mainly generated by the interaction between ocean infragravity waves and the ocean floor (Ardhuin et al., 2015;Nishida, 2017).The North Pacific and Southern Oceans are two locations of persistent hum generation (Ermert et al., 2017;Nishida, 2017;Rhie & Romanowicz, 2004, 2006), and thus indicate a non-uniform noise source distribution at the relevant frequency range.
The direction to the noise sources can be estimated from the properties of the energy flux of the wavefield, which is typically determined by applying beamforming or frequency-wavenumber array analysis techniques to far-field observations (Nishida, 2017, and references within).However, even though information about the background illumination can help mitigate the bias on the phase velocity estimation (Igel et al., 2023), results of these methods are not routinely used for quality assessment in ambient noise tomography applications.The USArray focal spot shapes (Figure 3b) are influenced by anisotropic surface wave incidence (Section 2.3).We use the related wavenumber domain information (Figure 3a) to analyze the direction of the strong and weak Rayleigh wave energy incidence at each station for T = 100 and 200 s (Figure 5).
At any one station, the Rayleigh wave energy forms a heterogeneous ring-like feature in the wavenumber domain at radial wavenumbers smaller than k r = 0.02 rad/m (Figure 3a).For our analysis we associate the direction of the maximum and minimum value of the amplitude spectrum within the ring with the strong and weak energy incidence of the surface wavefield, respectively.We search for the maximum value along the ring and estimate the direction between that location and the center of the spectrum.Along a circle that is centered at the origin and that intersects with the maximum value we then identify the direction of the weakest incidence as the direction toward the minimum amplitude value.This procedure is implemented for each ZZ component focal spot on the array.In Figure 5 we plot the direction of strong and weak incidence at the location of each station as black and red lines, respectively.To enhance the display quality, the lengths are scaled in each panel individually which does not permit a comparison between strong and weak values.
We compute the average of all maximum directions and plot the resulting orange great circle direction using the average geographic location of the station coordinates as origin.As suggested by the orientation of the majority of the strong incidence indicating black lines, for both periods the average direction of maximum energy propagation points to the North Pacific, which is compatible with observations of hum excitation in this area.The strongest deviation from the average trend is observed in the West-Southwest and East-Northeast areas of USArray.The black lines in the California region are of equal length compared to the dominating regions West of the Rocky Mountains.In contrast, the energy flux in the Northeast is weaker.This indicates different mechanisms for the observed trend at the edges or corners of the study area.In the West, the direction toward the ocean between the coast and the Sierra Nevada implies a similarly strong flux albeit from a different direction.In contrast, in the Northeast the intra-continental location with the greatest distance to an oceanic source region leads to the smallest amplitudes in the wavenumber domain distributions.We consider these physical mechanisms related to the proximity to the source area and structural effects to be dominant.They dominate in comparison to 2DFT amplitude effects associated with the variable spatial sampling of the receiver station locations relative to the reference station locations in different parts of the array.
The dominant pattern in the distributions of the red indicated directions of weakest energy flow is the difference between the central region of the imaged area, where the lines are comparatively long, and the regions toward the image area boundaries to the West and East.This differs from the Southwest-Northeast edge area pattern of the strong incidence directions.We hypothesize that the consistent weak incidence line-length pattern is controlled by the array shape.In the center, the larger East-West array extension leads to an higher minimum spectral domain amplitude level compared to the reference stations along the edges that are associated with more narrow arrays.
These results indicate the anisotropic incidence of the wavefield, which implies that phase velocity estimates from these focal spots can be biased.However, if the phase velocity is estimated from a focal spot that is sufficiently reconstructed in all directions, the bias is not significant (Aki, 1957;Boschi et al., 2012;Giammarinaro et al., 2023;Nakahara, 2006;Tsai & Moschetti, 2010).This is compatible with the observed trend of increased error estimates along the edges of the array discussed in Section 4.3.Extensions to the isotropic SPAC model that consider the effect of anisotropic illumination on the autocorrelation fields (Nakahara, 2006) are suitable alternatives to mitigate potential biases (Tsarsitalidou et al., 2022) but are not applied in this work.

Dispersion Curves
Our goal is to demonstrate the effectiveness of the focal spot dispersion analysis for improving the depth resolution of surface wave array records.We apply the parameter estimation described in Section 2.4 to narrow-band filtered focal spots computed at each reference station in the period range 60-310 s, which results in dispersion curves at each station location (Figure 6).We use ZZ component data because of their better SNR compared to ZR and RZ data (Section 2.4, Figure 4).We use r fit = 1.2λ.Our 60-310 s range extends the T = 150 s long period limit obtained with a USArray noise tomography (Zhao et al., 2017) by a factor of two, which translates into a similarly extended depth resolution.
We use k-means clustering to collect the dispersion curves into k = 3 groups of equal variance by minimizing the sum of squares of the distances between each sample and the centroid (Pedregosa et al., 2011).This choice reflects a conservative estimate of the scale that we aim to resolve.The obtained dispersion patterns can be associated with large-scale tectonic provinces (Fenneman & Johnson, 1946), which supports our goal and implies that the results are not governed by imperfect wavefield properties or imaging configurations.Figure 6a shows the resulting three well separated average dispersion curves that are related to the algorithmically-defined groups of stations in the inset.The orange indicated stations to the West exhibit the lowest obtained phase velocities.They comprise the low-velocity regions of the Snake River Plain, and the Basin and Range and Colorado Plateau (Obrebski et al., 2011).The boundary between this group and the second green colored group follows the Rocky Mountains Front indicated by the manually inserted black line in the inset.This green colored group covers the Great Plains and the Mississippi Embayment.The third purple indicated group includes stations in the Northeast of the study area and is associated with the high velocities of the tectonic regions of the Mid-Continent Rift and Super Craton (Fenneman & Johnson, 1946).
We obtain a standard error estimate of the phase velocity value at each station and period (Equation 3).Through error propagation we compute period dependent error bars of the three average dispersion curves.For all three clusters, the error increases toward the low-and high-period limits of the dispersion curves, while it remains comparatively small at the central periods.At short periods the large error is related to the small number of samples per wavelength and an associated less well constrained estimate.At long periods the 1.2λ criterion translates to long r fit distances and an increased tendency to break the azimuthal average (Section 2.4).
We investigate the depth resolution of the dispersion curves for the period range explored by our focal spot data to inform future focal spot-based 3D shear wave velocity modeling.We invert the average dispersion curves using the Geopsy tool (Wathelet, 2008) that utilizes the Neighborhood Algorithm (Sambridge, 1999).Following Durand et al. (2017) we choose a five layer parametrization, where the topmost layer represents the crust and the four lower layers properties of the mantle.For the starting model, the interface between the fifth layer and the halfspace is at 900 km depth.We choose the v S values from Durand et al. (2017), and the not further relevant v P values are taken from the Preliminary Earth Model.For each interface depth and v S value we allow the Neighborhood Algorithm to explore a ±30% variability around the initial values.We evaluate 15,300 models in the inversion process (Figure S1 in Supporting Information S1).
For each of the three groups, the average of the 1,000 best models (solid line in Figure 6c) is used to compute the sensitivity kernels shown in Figure 6b (Haney & Tsai, 2017) using the same color convention as in Figures 6a and  6c. Figure 6c combines the average (solid lines) and the 1,000 best models (transparent lines) for each group.The choice of averaging the 1,000 best shear wave velocity models for the sensitivity kernels and the display in Figure 6c is arbitrary, yet it is sufficient for our discussion.The three synthetic dispersion curves shown by the thin lines in Figure 6a correspond to the single best inverted model.They show good agreement to the data for the full period range, except at the largest period T = 310 s.The three groups of sensitivity kernels shown for periods between 60 and 300 s indicate that the resolution of the focal spot dispersion measurements covers the depth range down to ∼800 km.The results highlight the complementarity of autocorrelation data from short, potentially subwavelength distances for passive dense array surface wave imaging.We discuss limits associated with the anisotropic incidence and the location dependent synchronicity pattern, and the features in the three average models in Figure 6c in the following sections.

The Reference Phase Velocity Map
In this section we discuss the imaging result using T = 60 s ZZ focal spot data (Figures 7a and 10a).This reference result is again used to demonstrate the practicality of the approach by comparing it to tomography results of the same period.In the next Section 3.4 we perform sensitivity tests to assess the robustness of our short distance imaging results by changing different processing parameters (Figures 7b-7g).We discuss the ZZ r fit = 1.2λ results for the full period range in Section 4, where we also systematically consider the obtained uncertainty and SNR estimates.
The focal spot dispersion measurements assign a c R , ε c R , and RSS estimate to each station location.The resulting phase velocity maps are compiled using a Voronoi tesselation governed by the USArray deployment geometry.
The Voronoi patterning reflects the data driven sampling of the focal spot image that is thus obtained without regularization or interpolation.This differs from ambient noise tomography imaging, where the size and shape of the grid cells is typically chosen to optimize the ray coverage in each cell and the solution to the inverse problem with its inherent smoothing or damping length scales (Schaefer et al., 2011).
The ZZ component 60 s reference image Figure 7a is obtained with the processing described in Section 2.2 and the parameter estimation from Section 2.4 using r fit = 1.2λ.As shown by numerical experiments (Giammarinaro et al., 2023(Giammarinaro et al., , 2024) ) and demonstrated in Section 3.4.1, a value around 1λ well trades off uncertainty and resolution.This, together with the test results discussed in Section 3.4, supports our choice to use the parametrization associated with the image in Figure 7a as reference throughout this study.For this combination of period, wave speed, data range, and average inter-station distance, the white circle with the 1.2λ radius indicates that the value in each pixel is constrained by about 50 amplitude samples.
The practicality of the technique consisting of a fairly standard noise imaging processing chain and parameter estimation is demonstrated, visually, by the convincing quality of the resulting 60 s image Figure 7a   Dalton (2019) (Section 4.3.1).The focal spot image convinces by its equivalent resolution, which includes the high similarity of the obtained velocity values in the 3.7-4.3km/s range, and the position and clarity of the imaged features, which are further supported by the small ε c R = 0.1 km/s scale (Section 4).On a large scale this includes the low velocities of the tectonically more active region in the West in contrast to the higher velocities indicative of the more stable area in the East.However, even smaller scale features such as the Yellowstone Snake River Plain, the Colorado Plateau, the Mississippi Embayment, and the velocity contrast along the Central Valley of California and the Sierra Nevada are well resolved in comparison to the tomography results.

Sensitivity Tests
In this section we investigate how the image quality is influenced by the choice of processing parameters for cross-correlation and parameter estimation (Figures 7b-7g and Figure S2 in Supporting Information S1).The order of the discussion follows the significance of the associated effect on the resulting image quality.We test how the r fit data range, the component, and processing choices associated with amplitude normalization and body wave minimization affect the focal spot images.We compare the obtained c R maps to the reference image Figure 7a, but do not evaluate systematic effects on the uncertainty or SNR.

Data Range r fit
In this first test we explore the imaging sensitivity to the r fit distance around the reference station that limits the data used to constrain the parameter estimation.Numerical simulations on lateral resolution show that spatial averaging or blurring effects scale approximately linearly with r fit in the range 0.5λ to 2λ (Giammarinaro et al., 2024).This implies short r fit values support high-resolution and high-contrast images.However, small r fit distances involve fewer data points per wavelength and because data are noisy, the estimated dispersion measurements and images tend to be noisy, too.We test the r fit distances 0.5λ, 1λ, and 1.5λ (Figures 7b-7d) and compare the results to the reference case with r fit = 1.2λ.We apply the reference parametrization ZZ component, T = 60 s, and standard deviation clipping.
Compatible with numerical experiments (Giammarinaro et al., 2024), we observe that the longer r fit distance results in smoother lateral velocity variations (Figure 7d).Larger r fit distances support a more stable estimation of the phase velocity, as more data points constrain the Bessel function models.However, while a long data range stabilizes the parameter estimation, it negatively influences the lateral resolution power of the here promoted dense array local imaging implementation by producing smoother distributions which is equivalent to a loss in amplitude contrast of small-scale velocity variations (Giammarinaro et al., 2024).The image quality is probably slightly enhanced for r fit = 1λ (Figure 7c) compared to the r fit = 1.2λ reference (Figure 7a), considering the increase in contrast along elongated features, for example, along the Snake River Plain and along the Central Valley and Sierra Nevada.At the shortest considered data range r fit = 0.5λ (Figure 7b), however, the image still features the large-scale structures but is dominated by incoherent fluctuations on the shortest length scale, in addition to an overall significant offset of the velocity values toward higher estimates.We explain this offset by the interference of P wave energy.The observed increased velocities are highly compatible with the numerical results of Giammarinaro et al. (2023) that show the most significant c R overestimation at r fit = 0.5λ across the range of the tested r fit values 0.25λ, 0.5λ, 1λ, and 1.5λ.
These results demonstrate the strong dependence of the image quality on the data range that is an important tuning parameter to trade off resolution and uncertainty.The r fit distance should be small enough to maintain a good resolution and at the same time large enough to mitigate systematic errors associated with body wave incidence or anisotropic surface wave illumination and random fluctuations (Giammarinaro et al., 2023(Giammarinaro et al., , 2024)).The parameter estimation at the end of the analysis chain allows for an efficient exploration of the r fit related effects compared to studying the effects of the more expensive seismogram processing discussed in Sections 3.4.3and 3.4.4.

Component
In the second test, we compile a 60 s phase velocity map from ZR component cross-correlations (Figure 7e).Similar to the reference case we include all available data in the database, we apply the same processing to the seismograms described in Section 2.2, and we use r fit = 1.2λ.Using the same color range as the reference Figure 7a shows that the ZR component-based velocity distribution in Figure 7e features consistently smaller values.The broad tectonic features are similarly resolved.However, the ZR image is characterized by more small-scale fluctuations, and details in the ZR results differ from the tomography supported ZZ images, in particular elements along the eastern limit of the imaged area.As mentioned, the ZR autocorrelation fields are insensitive to P wave energy (Giammarinaro et al., 2023), consistent with the notion that ZR correlations are generally more robust to interfering wavefield components (van Wijk et al., 2011).
We attribute the difference to the reference results to the overall reduced ZR component data quality.Combined with the smaller signal amplitude of the model J 1 function the overall higher ZR amplitude fluctuations (Figure 4) yield a significantly lower SNR and higher phase velocity uncertainties (Section 3.3).Hence despite the advantageous insensitivity to body wave interference, the noisy character of the radial component data requires improvement strategies for the estimates to be included in potentially combined ZZ, ZR, and RZ-based Rayleigh wave speed constraints.

Time Domain Normalization
Processing the continuous records prior to computing the correlation functions is an elementary part of ambient noise tomography applications.We analyze the influence of two alternative ambient noise preprocessing strategies on the focal spot image quality.We begin with testing an alternative time domain normalization scheme (Hillers et al., 2016).The choice of the time domain normalization affects the SNR of the reconstructed propwave, and different methods interfere to varying degree with the original amplitude values (Bensen et al., 2007).The reference image is reconstructed from data clipped at three times the standard deviation of the amplitude distribution in each 4 hr processing window (Section 2.2).For the test we apply the more aggressive one-bit clipping to the reference ZZ data case.
A comparison of the resulting image Figure 7f to the reference velocity distribution Figure 7a implies that the choice of the clipping threshold has little effect on the image quality.The visual inspection reveals no significant changes, the results are highly consistent down to each pixel.Hence although the coherency or signal quality of the correlation wavefield is sensitive to the amplitude processing, the reconstruction of the governing relative phase information is not affected and hence the resulting spatial velocity distribution does not change, compatible with the conclusions by Hillers et al. (2016).This inference is supported by additional tests using a subset of 100 stations deployed in a square area between 112.5°-105°W and 35°-41°N, using one year of data (Tsarsitalidou et al., 2021).Scaling the one-bit, one, and five times standard deviation clipping results by the reference three times standard deviation clipping values shows that the three corresponding c R value distributions are statistically equivalent to unity.For T = 80 s we obtain 1.0046 ± 0.0041, 1.0042 ± 0.0030, and 0.9983 ± 0.0014, respectively.However, there is a trend toward comparatively higher RSS estimates for the smaller amplitude clipping thresholds, and the ZR signal reconstruction is also significantly negatively affected by the one-bit clipping.In conclusion, manipulating time domain amplitudes to account for wavefield transients does affect the obtained data quality, but the ZZ focal spot images are insensitive to the details.These findings support focal spot imaging because good results do not critically depend on a critically tuned processing chain.This does not suggest that preprocessing is not important.Data processing can be optimized for a specific data application considering spatial and temporal scales of the acquisition and wavefield properties (Fichtner et al., 2020).Whereas anisotropic Rayleigh wave incidence can be accounted for using updated focal spot or autocorrelation models (Nakahara, 2006), optimized data processing can make important contributions to the separation of different wavefield components (Giammarinaro et al., 2023).Modern wavefield separation techniques (Mousavi & Beroza, 2022) can then also result in Love wave focal spots, which are otherwise unattainable because Love waves and Rayleigh waves contribute both to RR and TT component autocorrelation fields (Haney et al., 2012).

Removal of Earthquake Signals
The second alternative preprocessing strategy aims to minimize effects of body wave energy.The interference of P wave energy with the surface wave focal spot is suggested by the consistently elevated data points above the best fitting J 0 model (Section 2.4, Figure 4) that result in systematically higher c R values (Figures 7a and 7f) compared to the ZR imaging results (Figure 7e) discussed in the previous section.
In this last test, prior to cross-correlation, we exclude time segments including and after large earthquakes.This eliminates the steeply incident body wave energy associated with reverberations that continue to be part of the global wavefield long after the local surface wave coda energy has dissipated (Boué et al., 2013;Boué, Poli, et al., 2014).We use an empirical relation to calculate the length of the discarded time window t d that is based on the timing and moment magnitude M w of an event (Boué, Poli, et al., For example, for a M w 6.0 event t d = 25 hr, M w 7.0 leads to t d = 51 hr, and for M w 8.0 t d = 108 hr.We exclude windows after events with a minimum magnitude M w 6.0 that are listed in a global catalog (USGS, 2017).This M w 6.0 limit removes about 45% of the data available for cross-correlation.For alternative M w 5.5 or M w 6.5 thresholds the rates are 70% and 20%, respectively.We use all other reference case parameters, that is, we apply the preprocessing explained in Section 2.2, we analyze ZZ component focal spots, and we use r fit = 1.2λ.Inspection of the parameter estimation results associated with successively smaller M w thresholds shows a systematically reduced discrepancy between ZZ amplitude data and the J 0 Bessel function model (Figure S3 in Supporting Information S1), which suggests our approach does reduce interfering P wave energy.For the ZR and RZ results, however, the scatter increases significantly from the tested M w 6.5 to M w 5.5 limit, with corresponding changes in the estimated wavenumber and hence phase speed.This negative effect is exclusively controlled by the reduced temporal averaging considering the insensitivity of the radial component data to body wave energy (Giammarinaro et al., 2023).
Before we assess the impact of this strategy on the obtained image in Figure 7g we discuss the effects on the energy distribution in the wavenumber domain illustrated in Figure 8.For a reference station located at the center of the array, we show wavenumber domain narrow-band filtered focal spots at T = 60 s (Figures 8a and 8c) and T = 240 s (Figures 8b and 8d) obtained using the full data set (Figures 8a and 8b) and after removal of time windows associated with large earthquakes (Figures 8c and 8d).Comparing the test results in the lower row with the original upper row results shows two effects.First, we can discern a change in relative amplitude at the center part of the spectrum, and this change is stronger at the long period.It suggests that vertically incident P wave energy associated with small k values has been reduced.Second, in the case without the large event data, the azimuthal distribution of the Rayleigh wave energy has improved, it is now more isotropic.
The effect of the earthquake data removal on the focal spot image quality highlights the diverse effects of the different wavefield components and processing choices.Compared to the reference image Figure 7a the new T = 60 s results in Figure 7g have an equivalent quality in terms of the feature reconstruction and image coherence.There is perhaps a small systematic ∼0.1 km/s increase in the estimated c R values that is best seen by comparing the low-velocity regions in the western U.S. Here, then, the minimization of P wave energy using r fit = 1.2λ results in a change toward larger values, that is, the reference image has perhaps a weak tendency to underestimate c R , which is compatible with the discussion of the Giammarinaro et al. (2023) results in Section 3.3.
Together, the wavenumber domain focal spot properties and the image quality associated with earthquake data removal suggest a somewhat better performance.However, for the rolling TA deployment with an average twoyear stationing, data removal can significantly reduce temporal averaging, that is, the target improvement of removing P wave bias can be compensated or offset by a reduced surface wave focal spot quality.Overall, the choice of the threshold for the minimum magnitude should trade off the positive and negative effects associated with minimizing the biasing P wave energy and the data discard, respectively.For surface wavefields, the data range can be as short as the SNR or residual fluctuations allow.The presence of P wave energy requires r fit > 1λ for sufficiently good estimates (Giammarinaro et al., 2023).Alternative processing strategies that can result in similar improvements of the surface wave focal spot properties include iterative correlation approaches referred to as the C2 and C3 method, where the reconstructed ballistic surface wave arrival and the correlation coda waveforms are re-correlated, respectively (Stehly et al., 2008), with different implications for the long wavelength limit.

Sensitivity Test Summary
In Figure S2 in Supporting Information S1 we show the sensitivity test results Figure 7 that are scaled by the reference image Figure 7a.This alternative view illustrates some of the consistent c R effects discussed in Sections 3.4.1 to 3.4.4,for example, the consistently smaller ZRc R level by about 2%-4%, or the systematic effect of the earthquake removal in the 1%-2% range.In addition, the normalized images contain patterns that mirror imaged tectonic features (Figures S2b, S2d, and S2e in Supporting Information S1), which suggests systematically combined effects of processing parameters and structure related wavefield properties, in contrast to the spurious impression associated with large-scale linear features (Figures S2f and S2g in Supporting Information S1).
We conclude this section by comparing estimates from our seven different sensitivity tests to sample 60 s c R values at three locations in the Great Basin, the Columbia Basin, and the Williston Basin obtained with noise correlation dispersion analysis (Figure 9 in Lin et al., 2014).Rounded to the next 0.05 km/s interval considering the uncertainties, the c R values are 3.80, 3.90, and 4.05 km/s.We compare these to values from the single closest station in the maps in Figures 7a-7g using the same precision.The comparison shows that the agreement of our results to the tomography is within 0.5% for cases with longer r fit (Figures 7a, 7d, and 7f), and that the difference tends to increase for cases with short r fit and low data quality (Figures 7b, 7e, and 7g).The same trends are obtained using not rounded values, and averages from the five nearest stations to the sample locations.
Together with the scaled images (Figure S2 in Supporting Information S1) the sample location comparison suggests our ZZ reference data processing that we use for the dispersion analysis (Sections 3.2 and 4.3) and for the phase velocity maps (Section 4.3) produces quality results.Although our testing cannot conclusively reconcile the systematic ZZ and ZR differences, we conjecture that the variable data quality likely controls the observed different phase velocity levels.We establish that the ZZ results are based on considerably better data quality compared to the ZR and RZ autocorrelation fields, and that the influence of body wave energy on the image results can be neglected here.

Spatial Phase Velocity Variations in Relation to the Imaging Wavelength
In this section we analyze the spatial scale of the lateral phase velocity variations for the period range between 60 and 310 s and compare this size scale to the wavelength that is used to probe the lithosphere.For this we compute a period dependent ratio (Figure 9) that helps to interpret the obtained focal spot imaging results (Figure 10) in the next Section 4.
The range of wavelengths that can be analyzed with the focal spot technique depends on the minimum and maximum size of the reconstructed focal spot.As in seismic tomography, then, the period range of the seismic surface wave focal spot application is controlled by the array properties.Generally, the shortest and largest wavelengths are constrained by the interstation distance and the array aperture, respectively.In both cases, however, the data quality has a significant effect on the resolution.It is advantageous to achieve robust c R estimates with a small r fit data range as possible.For any given period, this length scale affects the parameter estimation and the lateral resolution immediately (Figures 7b-7d) (Giammarinaro et al., 2024), and it thus differs compared to the propagation distance limit employed in seismic surface wave tomography.For the USArray these limits facilitate the reconstruction of focal spots and associated phase velocity maps in the 60-310 s period range (Figure 10a) using a nominal data range value of r fit = 1.2λ.
In Figure 10a, the 60 s phase velocity results are corroborated by their similarity with established tomography results.Across the imaged domain the distribution is characterized by variations that are smooth compared to the wavelength that is indicated by the white line inside the circle with radius 1.2λ.Toward the long period limit, however, the results show variations on scales that are significantly smaller compared to the wavelength.To quantify this trend, we relate the scale of the spatial variations to a reference scale that we choose as the average wavelength for each period.This 1λ reference is a conservative interpretation of the diffraction limit and related to the typical resolution of surface wave tomography implementations (Barmin et al., 2001;Mordret et al., 2013).
For this we apply a Fourier decomposition to the obtained phase velocity maps.Before the 2DFT image analysis, we compute the positions of the station locations relative to a reference location at the center of the area and convert the Cartesian to Polar coordinates.We de-mean the phase velocity distributions at each period and compute the 2DFT.The resulting amplitude spectrum of the reference 60 s phase velocity distribution Figure 7a is shown in Figure 9b together with its azimuthal average in Figure 9c.The position of the large amplitude features in Figures 9b and 9c reflects the dominating pattern of the continental-scale velocity variations across the U.S. For periods between 60 and 200 s the spatial scale of the velocity variation is relatively constant with a maximum amplitude at k r ∼ 0.002 rad/km (Figure S4 in Supporting Information S1), but at longer periods the range extends to larger k r values and hence smaller wavelengths.The black circle in Figure 9b and the black line in Figure 9c indicate the average wavenumber computed from the average wavelength which is the considered lateral resolution reference.To estimate the power of the variations below and above that reference we compute the area below the average spectrum at radial wavenumbers k r that are smaller (blue area in Figure 9c) and larger (red area) than the mean wavenumber, respectively, and then take the ratio of the two integrals.
We estimate this ratio in the target period range between 60 and 310 s with a 10 s period interval.The ratio is influenced by the reference wavenumber and by the spectral amplitude distribution.The decreasing trend in Figure 9a hence reflects the decreasing k threshold, but also the redistribution of energy from smaller to larger k r .Comparatively large numbers as for the shortest period indicate relatively smooth velocity variations compared to the probing wavelength.At around T = 160 s the slope flattens significantly toward the small numbers at the long periods that reflect relatively rough distributions.We use these trends together with data quality markers, and parameter estimation and sensitivity tests for a comprehensive assessment of the obtained USArray focal spot imaging results.

General Observations
The analogy of modern dense seismic networks to ultrasonic transducers facilitates the application of elastography medical imaging concepts to seismic surface wave data.A high station density and clean correlation surface wavefields support the depth extension compared to ambient noise tomography and potentially sub-wavelength resolution of lateral medium variations.Numerical experiments targeting the lateral resolution power of the focal spot method demonstrate that the data range is the most important tuning parameter.Although the choice of r fit affects the contrast, the phase characteristics of the imaged velocity variations are robust (Giammarinaro et al., 2024).Overall, the application to dense local array data (Hillers et al., 2016) and the numerically investigated resolution properties (Giammarinaro et al., 2023(Giammarinaro et al., , 2024) ) demonstrate the effectiveness of the focal spot technique for a wide range of imaging applications including feature detection and accurate wave speed estimates, that is, the technique supports discovery mode and inference mode imaging (Tsai, 2023).Here we apply this method for the first time on a continental scale using USArray data to estimate dispersion and to image phase velocity distributions of the fundamental mode Rayleigh wave between 60 and 310 s period.
In this section we discuss the ZZ component results that are obtained with our reference processing which includes standard deviation amplitude clipping during the preprocessing stage.The ZZ focal spot data are analyzed using a model associated with isotropic propagation and a data range r fit = 1.2λ.Our data analysis implies that P wave energy interferes with the surface wave measurements.However, at T = 60 s, the sensitivity tests, a sample comparison to tomography phase speed estimates, and the high similarity of the images to equivalent tomography maps together demonstrate that our ZZ results are accurate and not systematically biased.We prefer to not use ZR results in our conclusive discussion due to the considerably lower SNR and the associated reduced image quality.Future focal spot applications to three-component data can potentially benefit from reconciling ZZ and ZR observations for yet better constrained estimates of the Rayleigh wave phase speed.
In the next section we continue the discussion from Section 3.2 about the dispersion results and the associated vertical resolution.Lateral resolution is then examined in a detailed discussion of the imaging results including uncertainty and quality estimates.We conclude with implications for resolving anelastic attenuation and medium anisotropy, and on the general significance of focal spot imaging.

Vertical Resolution
Our k-means cluster approach (Section 3.2) yields three average dispersion curves that show a general increase in the seismic velocities from West to East (Figure 6a).The uncertainty associated with each average dispersion curve suggests well separated clusters.The effectiveness of the focal spot method to increase the vertical resolution compared to array surface wave tomography (Section 3.2) is suggested by the good agreement between the data and the synthetic dispersion curves up to periods of 300 s.The sensitivity kernels (Figure 6b) associated with the averaged five-layer models imply a depth resolution down to the ∼800 km range.Averaged over the full period range, the decreasing error of 0.016, 0.015, and 0.014 km/s for the three clusters from West to East can be connected to the increasingly homogeneous conditions in the respective regions.By averaging the 1,000 best shear wave velocity models for the display in Figure 6c we consider a large enough number of sufficiently converged solutions.Details in the depth dependent profiles depend on this choice, and to a lesser degree also on the parameter sampling by the Neighborhood Algorithm.
From a suite of alternative average models associated with different random seeds we can identify three robust key features in Figure 6c.First, the crust and upper mantle in the top 300 km in the West shows slower velocities compared to the regions to the East.Second, the purple profile exhibits a low-velocity zone between 300 and 600 km, which is compatible with the relative low-velocity feature at the same depth along a North-South profile at about 92°West imaged by Porritt et al. (2014, their Figure 8).Third, the velocities in the mantle transition zone and lower mantle below 400-500 km depth are consistently higher in the East.These observations are averaged over relatively large areas, but the compatibility with results from other imaging techniques generally demonstrates the increased depth resolution of focal spot passive surface wave imaging for extended 3D shear wave velocity models.

Properties of Phase Velocity Maps
Figure 10a shows the focal spot obtained Rayleigh wave phase velocity maps at six periods.The phase velocities increase from around 4 km/s at T = 60 s to around 5.5 km/s at T = 300 s. Figure 10b illustrates one of the key advantages of the focal spot approach, that is, the access to the location dependent error estimate for the imaged area.As discussed for the dispersion results in Figure 6a the error is period dependent, with the largest values around 0.1 km/s at T = 60 s, smallest values around 0.035 km/s at T = 150 s, and then again larger values around 0.06 km/s toward the 300 s long period limit.Another measure of the data quality is the RSS value shown in Figure 10c.This shows a similar period dependence but here the values around 0.035 at T = 60 s are matched again at T = 200 s and then increase to 0.07 at T = 300 s.The best values of the RSS quality proxy are found at 100-150 s period.The RSS values are normalized to remove the effect of the number of samples that varies for a nominally constant range of r fit = 1.2λ for the considered wide period band.To facilitate the comparison with global tomography models we show in Figure 10d the low-pass filtered phase velocity distributions.We use a Gaussian kernel that has the same scale as the seismic wavelength at each period which is indicated by the white lines in the panels in Figure 6a. Figure 10e illustrates the original maps as perturbations in percent around the mean phase velocity estimate.
Supported by the information in Figures 10b-10e, the discussion of the phase velocity images in Figure 10a is divided into a part that focuses on the shorter period range where we can compare the results to surface wave tomography results, and a long period part for which no reference surface wave observations are available.For results between 60 and 180 s, we discuss a quantitative comparison between our focal spot images and the tomography results of Jin and Gaherty (2015), Zhao et al. (2017), and Babikoff and Dalton (2019) (Figure S5 in Supporting Information S1).

The Short Period Range 60-150 s
In the 60-150 s range we can compare our results to the images obtained using both earthquake (Babikoff & Dalton, 2019;Jin & Gaherty, 2015) and ambient noise data (Zhao et al., 2017).The obtained USArray phase velocity results present updated information about the crustal and mantle structure beneath the U.S. that is presented in numerous structural studies summarized by Shen and Ritzwoller (2016).The high similarity of our 60 s results in Figures 10a and 10e to the maps in Figure 9 in Jin and Gaherty (2015) and Figure 6 in Zhao et al. (2017) emphasizes, again, the equivalent focal spot imaging quality.Zhao et al. (2017) extend the long period limit of their ambient noise tomography by applying a comparatively short 1λ cutoff distance.This allows phase velocity estimates up to 150 s period.This 1λ threshold is the shortest propagation distance between two sensors considered for an unbiased speed estimate, whereas our 1.2λ data range is the maximum distance of the spatial autocorrelation field data used to constrain the local Bessel function model.The two approaches target different components of the surface wavefield (Figure 1), yet they are similarly sensitive to the physical medium properties.
At 100 s, we can consider the results of the two earthquake tomography studies by Jin and Gaherty (2015) and Babikoff and Dalton (2019).We find that the large-scale velocity structure of the generally slower velocities to the West of the Rocky Mountain Range, but also the comparatively higher velocities in South-East Texas, are equally resolved by our focal spot approach and by the ambient noise and earthquake tomography.Our 150 s results show a similar correspondence of the large-scale features in the images in Zhao et al. (2017).At smaller scale we highlight the consistently imaged high-velocity patch at the border of Oregon, California, and Nevada, and the North-South trending red indicated linear low-velocity feature to the North of the Texas Panhandle in Figure 10a.This feature is displayed in yellow color in the perturbation map Figure 10e, and, because of the domain size dependent mean value, in white color in the corresponding panel in Figure 6 in Zhao et al. (2017).The resolution of this one pixel wide feature in the 300 km depth range is noteworthy considering the indicated r fit data range and wavelength.It supports the notion that focal spot imaging has sub-wavelength resolution properties.at 150 s in the Washington and Oregon area differs between the noise-based solutions, which is likely influenced by the location close to the edge of the imaging domain.In addition to the Helmholtz tomography Babikoff and Dalton (2019) produce phase velocity maps between 25 and 180 s period that are computed from six different 3D shear wave velocity models.Differences between the obtained images are attributed to the variable sensitivity of the different methods and data sets.Comparing our focal spot to tomography results we can conclude that focal spot imaging is competitive in the 60-150 s range.
We compute the percentage difference (Figure S5 in Supporting Information S1) between the focal spot phase velocities obtained with the reference processing (Section 3.3) and results from three tomography studies (Babikoff & Dalton, 2019;Jin & Gaherty, 2015;Zhao et al., 2017).Spatial distributions of the percentage difference and the respective histograms show that the focal spot method tends to overestimate values compared to the other studies.The histograms exhibit an approximately Gaussian shape centered at +2% difference.As discussed in Section 3.4.4,a larger velocity is a likely effect of the P wave energy that interferes with the ZZ component surface wavefield.Stronger velocity differences are observed toward the edges, where our formal error estimates are also larger (Figure 10b).Apart these wavefield anatomy and configuration effects, we observe consistent small-scale features in the comparison maps that are associated with tectonic variations in the western part of the U.S. such as Yellowstone and the Snake River Plane.This is related to different contrasts across lateral variations obtained with our autocorrelation and the cross-correlation tomography implementations (Giammarinaro et al., 2024).
The uncertainty estimates (Figure 10b) are typically much smaller away from the domain edges.The size of the area that is affected by the edges increases with wavelength.However, we can extend our observation from Section 3.3 to the full image domain and conclude that the formal uncertainty is small compared to the imaged phase velocities.This implies sufficiently resolved speed variations, which supports the conclusion from the comparison to tomography results.
The normalized residual sum of squares RSS (Figure 10c) allows an alternative assessment of the regression results.Generally there is a tendency of the best RSS values to be located away from the edges.Compared to the uncertainty estimates, however, the RSS maps show more texture, which is governed by North-South oriented features in the here considered period band.The pattern of good RSS values-along the Rocky Mountains, along the northern and eastern image edge-and bad RSS values-to both sides of the Rocky Mountains, in southeast Texas-does not seem to be controlled by a single structural, wavefield, or deployment feature.Instead, the ZZ data scatter around the best fitting model appears to be influenced by a combination of various mechanisms that require further investigation to be more completely understood.
For the period band 60-150 s, the average uncertainty estimates (Figures 10b and 10c) decrease with increasing period.At the same time, the amplitude of the velocity variation decreases as well (Figure 10e), and together with the fact that in this period band the azimuthal average is almost complete for most stations, we interpret the large uncertainty at the smallest 60 s period to be controlled by stronger lateral heterogeneity.This is compatible with the effect of local scatterers on the focal spot reconstruction observed by Hillers et al. (2016).As summarized in the next Section 4.3.2, the error increase toward longer periods is mostly attributed to insufficient azimuthal averaging.
The low-pass filtered images (Figure 10d) are compatible with results from global tomography.The general West-East trend of increasing wave speed is similarly resolved by the hum-based images of Haned et al. (2015).In comparison to these low-resolution images the regional tomography-calibrated original maps (Figure 10a) show that the employed data range r fit = 1.2λ effectively increases the resolution.Figure 7d implies that a data range of r fit = 2λ or longer yield similarly low-pass filtered results.Longer data ranges degrade the amplitude contrast (Giammarinaro et al., 2024), but the resolution-uncertainty trade-off can help to improve estimates from noisy or biased data (Giammarinaro et al., 2023).We demonstrate this relation here, it is relevant for the discussion of the results at T ≥ 200 s for which only comparatively low-resolution global tomography results are available.
We iterate that the simple models (Equations 1 and 2) used for the c R estimation assume isotropic illumination or an equivalent complete azimuthal average.The general robust similarity to the surface wave tomography results even in regions close to the image boundary such as the velocity structure along the California Central Valley and the Sierra Nevada suggests that the images in the here considered 60-150 s range are not systematically biased by the directive surface wave energy flux (Figure 5).Incomplete averaging can influence the higher uncertainty estimates along the boundaries, but does not degrade the overall image quality here.
Journal of Geophysical Research: Solid Earth 10.1029/2023JB027417 autocorrelation data collected at distances on the order of one wavelength using models based on isotropic surface wave propagation.Giammarinaro et al. (2023) demonstrate that biasing body wave energy can be mitigated using longer data ranges for the surface wave assumption.Alternative strategies include the implementation of filters prior to or after correlation, or the consideration of the-potentially oblique-body wave incidence in the autocorrelation parametrization that can be further constrained by radial component data.As seismic energy redistribution through scattering at the longer periods T ≥ 200 s is less efficient compared to the crustal scales the associated focal spots are more influenced by the directive illumination.Together with the acquisition patterns and a smaller effective data range the obtained small-scale velocity variations require independent verification.However, the variations do not reflect a systematic connection to any of the considered control parameters including the relative sensor position, uncertainty, goodness-of-fit, or local flux pattern.
To further evaluate the consistency of the focal spot results we predict phase velocity maps (Babikoff & Dalton, 2019) from the global model SEISGLOB2 (Durand et al., 2017).We convert the S wave velocity perturbations relative to PREM to absolute v S values in km/s in each layer.The P wave speed v P in km/s and density ρ in units of g/m 3 are estimated using v P = 1.73v and ρ = 0.32v P + 0.77.We compute dispersion curves (Herrmann, 2013) on a 1°× 1°grid and compile period dependent phase velocity distributions.The predicted images (Figure S6 in Supporting Information S1) are low-pass filtered versions of the focal spot results, which is also evident from the similarity of the SEISGLOB2 predictions to the Gaussian-smoothed focal spot images in Figure 10d.That is, the broad patterning is consistent, with a correspondingly larger spatial scale of the phase velocity variations predicted from the global model.However, at short 60-150 s periods the small-scale features in the focal spot and the dense array tomography images (Section 4.3.1)are more similar to each other than to the SEISGLOB2 predictions.These observations highlight the need for controlling the imaging sensitivities associated with different acquisition patterns, wavefield components, theories, techniques, and algorithmic and parameter choices.
Our results suggest that the focal spot technique is a useful complementary imaging tool that can enhance the resolution of vertical and horizontal velocity variations of the mantle structure below the U.S. for the study of mantle dynamics and plate tectonics.Controlling cleaner surface wavefield properties the here applied r fit = 1.2λ data range can potentially decrease and in turn the lateral resolution can increase.The explicit consideration of the anisotropic illumination effect using updated SPAC expressions (Nakahara, 2006) can improve imaging results.We anticipate more significant enhancements at the longer periods compared to the shorter periods where sufficient azimuthal averaging has been shown to yield accurate estimates.Azimuth dependent surface wave speed variations are resolved in the 7-46 s range below the western U.S. (Lin & Ritzwoller, 2011;Lin et al., 2009;Zigone et al., 2015).A strong azimuthal Rayleigh wave speed dependence is inferred in the ∼1 km 2 -scale focal spot analysis in the San Jacinto fault zone environment (Hillers et al., 2016).This suggests that medium wave speed anisotropy at longer periods can be studied using USArray focal spot data, which requires a mitigation of the potentially interfering effects of anisotropic illumination (Giammarinaro et al., 2023;Nakahara, 2006;Roux et al., 2005).Prieto et al. (2009) apply SPAC to dense array data in California to estimate the anelastic structure.Nakahara (2012) and Magrini and Boschi (2021) confirm the application of the used ad-hoc exponential term that parametrizes the excess amplitude decay for weakly attenuating media.For this, however, it is essential to use long data ranges, hence constraining attenuation from noisy data at short one-wavelength focal spot distances appears challenging.The effectiveness of focal spot imaging to resolve first order wave speed variations compared to second order direction dependent and anelastic effects underlines universal principles of elastic and seismic wave propagation.

Conclusions
Focal spot imaging is an established medical imaging method.In contrast to tomography the local wave speed is estimated from properties of spatial autocorrelation fields.In seismology, the technique involves the application of SPAC models to dense array autocorrelation data at distances on the order of one wavelength, and this short scale allows an exploration of longer periods compared to tomography.In this work we apply focal spot imaging to noise correlation functions obtained from USArray data to estimate fundamental mode Rayleigh wave dispersion between periods of 60 and 310 s.The similarity to regional seismic tomography results in the 60-150 s range in this first large-scale survey demonstrates the effectiveness of passive seismic surface wave focal spot imaging.The local inversion-free imaging approach excels in dense array applications, which suggests its usefulness for the analysis of spatially dense data sampled by distributed acoustic sensing systems.
Our sensitivity tests show that the data range is an effective tuning parameter that trades off resolution and uncertainty.The test results at 60 s period reflect the better signal-to-noise ratio of the ZZ compared to ZR and RZ data, and that anisotropic illumination or interfering body wave energy do not critically influence the overall highquality, tomography-compatible ZZ images.Importantly, the local SPAC model regression yields estimates for the goodness-of-fit and the phase velocity error at each image pixel.Together with the depth sensitivity analysis results, this shows that 3D shear wave velocity models can be inverted using USArray focal spot phase velocity data and their uncertainties at periods up to 300 s.
The challenge to reconcile long period focal spot images in the 200-300 s range with results from comparatively low-resolution global tomography highlights the need to manage the data sensitivities, processing effects, and uncertainties.For focal spot imaging, better control over the combined effects of the array shape, relative position dependence, local flux pattern, wavefield component biases, and factors governing data quality and hence uncertainty can be rewarded with improved images of elastic material properties at previously unattainable depth ranges, and potentially also of anisotropy and anelastic properties.

Figure 1 .
Figure1.The ZZ cross-correlation wavefield associated with one USArray reference station at T = 65 s central period.The interstation distance is scaled by the corresponding seismic wavelength.The energy is scaled by the maximum amplitude in each distance bin.The reconstructed wavefield is dominated by the familiar converging and diverging propagating surface waves at negative and positive lag times, respectively.In addition we highlight at zero lag time the azimuthally averaged spatial autocorrelation field, which is the time domain focal spot.The visible amplitude ridge between 1λ and 2λ corresponds to the second maximum of the J 0 Bessel function.The relatively wide period band illustrates dispersion of the propagating waves, but it degrades the focal spot shape(Giammarinaro et al., 2023).

Figure 2 .
Figure 2. Illustration of the rolling Transportable Array deployment.Spatially dependent synchronicity patterns for three red indicated reference stations (a) M20A, (b) D06A, and (c) Z36A show that the deployment strategy results in station pairs that are predominantly oriented in the North-South compared to the East-West direction.

Figure 3 .
Figure 3. Representations of the spatial autocorrelation amplitude field or focal spot for the ZZ, ZR, and RZ components (top to bottom) at the reference station L15A.(a) Wavenumber spectra at T = 100 s.(b) Narrow-band filtered focal spots at T = 100 s (left column), T = 200 s (center column), and T = 300 s (right column).The location of L15A is indicated by the green cross.Amplitudes are scaled by the peak ZZ amplitude at each period.

Figure 4 .
Figure 4. Radial focal spot amplitude data associated with reference station R21A and the corresponding Bessel function models (Equations 1 and 2) illustrate the phase velocity estimation for (a) T = 100 s, (b) T = 200 s, and (c) T = 300 s.The estimated models from the first and third iterations explained in Section 2.4 are plotted in orange and black, respectively.The scale of the two models is indicated by the correspondingly colored axis.The ZZ data (left column) show a better signal-to-noise ratio compared to the ZR (center column) and RZ (right column) data.The consistently higher ZZ amplitudes at short distances is attributed to the influence of coherent P wave energy.

Figure 5 .
Figure 5. Spatial variation of the Rayleigh wave illumination pattern.Directions of strong (black lines) and weak (red lines) Rayleigh wave energy flux obtained from the ZZ component focal spot analysis in the wavenumber domain for (a) T = 100 s and (b) T = 200 s.For illustrative purposes the length of a line is scaled by the maximum length in each panel.The orange great circle represents the average maximum direction.

Figure 6 .
Figure 6.Results of the dispersion analysis.(a) Three average dispersion curves sampled with dT = 10 s interval for periods between T = 60 s and T = 310 s.The solid line is the observed cluster average dispersion curve, which is the target for the inversion of the 1D shear wave velocity model.The uncertainty estimate (Equation 3) is plotted at each period.The dashed lines indicate synthetic dispersion curves associated with the single best inversion result.The corresponding groups of stations are shown in the inset.U.S. physiographic divisions are colored in black.The Rocky Mountain Front is highlighted by the thick line.(b) Sensitivity kernels for six annotated periods associated with the average of the 1,000 best 1D shear wave velocity models.(c) The averages of the 1,000 best models obtained from the inversion of the three cluster average dispersion curves in panel (a) are plotted with solid lines.The corresponding 1,000 best models are plotted with transparent lines.
compared to tomography results, for example, to Figures 6 and 2f inLin and Ritzwoller (2011) andShen and Ritzwoller (2016), respectively, and by our formal comparison to results of Jin andGaherty (2015) andBabikoff and

Figure 7 .
Figure 7. Results from the sensitivity analysis (Section 3.4) using focal spot data at T = 60 s.All panels show fundamental mode Rayleigh wave phase velocity estimates.(a) Reference image compiled from ZZ component focal spot data using r fit = 1.2λ.We use the full data set and apply the preprocessing discussed in Section 2.2.Relevant information is summarized in each panel legend.Variations of the reference case using (b) r fit = 0.5λ, (c) r fit = 1.0λ,(d) r fit = 1.5λ,(e) ZR component data, and (f) an alternative preprocessing scheme.For (g) we exclude data containing global body wave reverberations excited by earthquakes with M w > 6.0.Most panels use the same color range to facilitate comparison.In each panel the white circle and the white line illustrate the applied data range r fit and the average wavelength, respectively.

Figure 8 .
Figure 8. Wavenumber spectra of ZZ component focal spots associated with the reference station R21A.Patterns for (a) T = 60 s and (b) T = 240 s are associated with focal spots that are reconstructed from all available data.For the corresponding patterns (c) and (d) we exclude data containing global body wave reverberations excited by earthquakes with M w > 6.0.Each spectrum is normalized by its maximum value.

Figure 9 .
Figure 9. Period dependent partitioning of phase velocity image variation scales.(a) The ratio quantifies the blue indicated area in relation to the red indicated area in panel (c), the relation of variations on scales that are smaller and larger, respectively, than the reference seismic wavenumber associated with the period dependent wavelength.(b) The 2DFT amplitude spectrum of the T = 60 s phase velocity map (Figure 7a) obtained from ZZ component focal spots.The black circle indicates the average seismic wavenumber at this period.(c) Azimuthal average of the spectrum in panel (b).The black line separating the blue from the red area corresponds to the circle in panel (b).

Figure 10 .
Figure 10.Focal spot imaging results between T = 60 s and T = 300 s obtained with ZZ component data and r fit = 1.2λ.(a) Phase velocity distributions.In each map the white circle and the white line illustrate the applied data range r fit and the average wavelength, respectively.The 60 s map is the same as Figure 7a.(b) The corresponding standard error estimate (Equation 3).(c) The normalized residual sum of squares (RSS) of the parameter estimation using the Bessel function model.(d) The low-pass filtered phase velocity distributions from panels (a) using a Gaussian kernel of size equal to the average seismic wavelength at each period.(e) The phase velocity distribution from panels (a) shown as perturbation around the average seismic velocity at each period.