Cross-correlations of ambient seismic noise are widely used for seismic velocity imaging, monitoring, and ground motion analyses. A typical step in analyzing Noise Cross-correlation Functions (NCFs) is stacking short-term NCFs over longer time periods to increase the signal quality. Spurious NCFs could contaminate the stack, degrade its quality, and limit its use. Many methods have been developed to improve the stacking of coherent waveforms, including earthquake waveforms, receiver functions, and NCFs. This study systematically evaluates and compares the performance of eight stacking methods, including arithmetic mean or linear stacking, robust stacking, selective stacking, cluster stacking, phase-weighted stacking, time-frequency phase-weighted stacking, $N^{th}$-root stacking, and averaging after applying an adaptive covariance filter. Our results demonstrate that, in most cases, all methods can retrieve clear ballistic or first arrivals. However, they yield significant differences in preserving the phase and amplitude information. This study provides a practical guide for choosing the optimal stacking method for specific research applications in ambient noise seismology. We evaluate the performance using multiple onshore and offshore seismic arrays in the Pacific Northwest region. We compare these stacking methods for NCFs calculated from raw ambient noise (referred to as Raw NCFs) and from ambient noise normalized using a one-bit clipping time normalization method (referred to as One-bit NCFs). We evaluate six metrics, including signal-to-noise ratios, phase dispersion images, convergence rate, temporal changes in the ballistic and coda waves, relative amplitude decays with distance, and computational time. We show that robust stacking is the best choice for all applications (velocity tomography, monitoring, and attenuation studies) using Raw NCFs. For applications using One-bit NCFs, all methods but phase-weighted, time-frequency phase-weighted, and $N^{th}$-root stacking are good choices for seismic velocity tomography. Linear, robust, and selective stacking methods are all equally appropriate choices when using One-bit NCFs for monitoring applications. For applications relying on accurate relative amplitudes, both the robust and cluster stacking methods perform well with One-bit NCFs. The evaluations in this study can be generalized to a broad range of time-series analysis that utilizes data coherence to perform ensemble stacking. Another contribution of this study is the accompanying open-source software, which can be used for general purposes in time-series stacking.

Volcanic arcs consist of many distinct vents that are ultimately fueled by the common process of melting in the subduction zone mantle wedge. Seismic imaging of crustal scale magmatic systems can provide insight into how melt is organized in the deep crust and eventually focused beneath distinct vents as it ascends and evolves. Here we investigate the crustal-scale structure beneath a section of the Cascades arc spanning four major stratovolcanoes: Mt. Hood, Mt. St. Helens, Mt. Adams, and Mt. Rainier, based on ambient noise interferometry measurements from 234 seismographs. Simultaneous inversion of Rayleigh and Love wave dispersion better constrain the isotropic shear velocity (Vs) and identify the unusual occurrence of radially anisotropic structures. Isotropic Vs shows two sub-parallel low-Vs zones at ~15-30 km depth with one connecting Mt. Rainier to Mt. Adams, and another connecting Mt. St. Helens to Mt. Hood, which are interpreted as deep crustal magma reservoirs containing up to ~2.5-6% melt, assuming near-equilibrium melt geometry. Negative radial anisotropy is prevalent in this part of the Cascadia margin, but is interrupted by positive radial anisotropy extending vertically beneath Mt. Adams and Mt. Rainier at ~10-30 km depth and weaker positive anisotropy beneath Mt. St. Helens with a west dipping. The positive anisotropy regions are adjacent to rather than co-located with the isotropic low-Vs anomalies. Ascending melt that stalled and mostly crystallized in sills with possible compositional difference from the country rock may explain the near-average Vs and positive radial anisotropy adjacent to the active deep crustal magma reservoirs.

Adjoint tomography (i.e., full-waveform inversion) has been recently applied to ambient seismic noise and teleseismic P waves separately to unveil fine-scale lithospheric structures beyond the resolving ability of traditional ray-based traveltime tomography. In this study, we propose a joint inversion scheme that alternates between frequency-dependent traveltime inversions of ambient noise surface waves and waveform inversions of teleseismic P waves to take advantage of their complementary sensitivities to the Earth’s structure. We apply our method to ambient noise empirical Green’s functions from 60 virtual sources, direct P and scattered waves from 11 teleseismic events recorded by a dense linear array (～ 7 km station spacing) and other regional stations (～ 40 km average station spacing) in central California. To evaluate the performance of the method, we compare tomographic results from ambient noise adjoint tomography, full-waveform inversion of teleseismic P waves, and the joint inversion of the two data sets. Both applications to practical field data sets and synthetic checkerboard tests demonstrate the advantage of the joint inversion over individual inversions as it combines the complementary sensitivities of the two independent data sets towards a more unified model. The 3D model from our joint inversion not only shows major features of velocity anomalies and discontinuities in agreement with previous studies. but also reveals small-scale heterogeneities which provide new constraints on the geometry of the Isabella Anomaly and mantle dynamic processes in central California. The proposed joint inversion scheme can be applied to other regions with similar array deployments for high-resolution lithospheric imaging.